<jupyter_start><jupyter_text>>>> Work in process <<<  Analyzing ordinal data with metric models: What could possibly go wrong?<NAME>., & <NAME>. (2018, April 5). Ordinal Data Analysis. Retrieved from http://osf.io/53ce9PyMC3 implementation of the ordinal probit model.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pymc3 as pm
import theano.tensor as tt

from IPython.display import Image
from theano.compile.ops import as_op
from scipy.stats import norm

%matplotlib inline

color = '#87ceeb'
f_dict = {'size':14}<jupyter_output><empty_output><jupyter_text>Function `gammaShRaFromModeSD` ported from https://osf.io/eujgd/<jupyter_code># Calculate Gamma shape and rate from mode and sd.
def gammaShRaFromModeSD(mode, sd):
    rate = (mode + np.sqrt( mode**2 + 4 * sd**2 ) ) / ( 2 * sd**2 )
    shape = 1 + mode * rate
    return(shape, rate)
# Read data file directly from https://osf.io.
# Selection of movie ratings from Amazon.com.
data = pd.read_csv('https://osf.io/38ue2/download', skiprows=1, usecols=np.arange(1,8))
data.info()
data.head()
# Convert the precentages to counts with help of N column. Order the columns.
y = (data.iloc[:,2:7]
     .divide(100)
     .multiply(data.iloc[:,1], axis=0)
     .round()
     .astype(int)
     .iloc[:,::-1])
# Strip "Prcnt" from the column names
y.columns = y.columns.map(lambda x: x[:5])
y.head()
z = y.sum(axis=1)
z.head()
Ncases = data.index.size
Ncases
# Number of outcomes
nYlevels = y.columns.size
nYlevels
gammaShRa = gammaShRaFromModeSD(3,3)
gammaShRa<jupyter_output><empty_output><jupyter_text>Lower panel of figure 1<jupyter_code>Image('Ordered_Probit_Model.png')<jupyter_output><empty_output><jupyter_text>*A latent scale on the horizontal axis is divided into subintervals with thresholds marked bydashed lines. The cumulative normal probability in the subintervals is the probability of theordinal values. The cumulative normal probability mass within each interval is indicatedby the height of the corresponding bar, with numerical scale indicated on the right verticalaxis.* (<NAME>, 2018, April 5)  Model - JAGS (Kruschke & Liddell, 2018, April 5)  https://osf.io/mvsc3/ Note: model contains a hierarchical structure on the standard deviations across movies, but not on the means of the movies. ``` THE *ORDERED PROBIT* MODEL:modelString = "  model {    for ( i in 1:Ncases ) {      y[i, ] ~ dmulti( pr[i,1:nYlevels] , z[i] )      pr[i,1] <- pnorm( thresh[1] , mu[x[i]] , 1/sigma[x[i]]^2 )      for ( k in 2:(nYlevels-1) ) {        pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[x[i]] , 1/sigma[x[i]]^2 )                           - pnorm( thresh[k-1] , mu[x[i]] , 1/sigma[x[i]]^2 ) )      }      pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[x[i]] , 1/sigma[x[i]]^2 )    }    for ( j in 1:Ncases ) {       mu[j] ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )      sigma[j] ~ dgamma( sigmaSh , sigmaRa )    }    sigmaSh <- 1 + sigmaMode * sigmaRa    sigmaRa <- ( ( sigmaMode + sqrt( sigmaMode^2 + 4*sigmaSD^2 ) )                 / ( 2*sigmaSD^2 ) )    sigmaMode ~ dgamma( gammaShRa[1] , gammaShRa[2] )     sigmaSD ~ dgamma( gammaShRa[1] , gammaShRa[2] )     for ( k in 2:(nYlevels-2) ) {   1 and nYlevels-1 are fixed, not stochastic      thresh[k] ~ dnorm( k+0.5 , 1/2^2 )    }  }"  close quote for modelString```  Model - PyMC3<jupyter_code># Thresholds, masking the the inner two values. 
thresh = [k + .5 for k in range(1, nYlevels)]
thresh_obs = np.ma.asarray(thresh)
thresh_obs[1:-1] = np.ma.masked

print('thresh:\t\t{}'.format(thresh))
print('thresh_obs:\t{}'.format(thresh_obs))
@as_op(itypes=[tt.dvector, tt.dvector, tt.dvector], otypes=[tt.dmatrix])
def outcome_probabilities(theta, mu, sigma):
    out = np.empty((nYlevels, Ncases))
    n = norm(loc=mu, scale=sigma)
    lbound = np.repeat(0, Ncases)
    
    # Thresholded cumulative normal probabilities.
    # Four thresholds (theta values) define the 5 outcome probabilities.
    out[0,:] = n.cdf(theta[0])        
    out[1,:] = np.max([lbound, n.cdf(theta[1]) - n.cdf(theta[0])], axis=0)
    out[2,:] = np.max([lbound, n.cdf(theta[2]) - n.cdf(theta[1])], axis=0)
    out[3,:] = np.max([lbound, n.cdf(theta[3]) - n.cdf(theta[2])], axis=0)
    out[4,:] = 1 - n.cdf(theta[3])

    return out
with pm.Model() as ordinal_model_multi_groups:    
    # Latent means (rating) of the movies           
    mu = pm.Normal('mu', mu=(1+nYlevels)/2.0, tau=1.0/(nYlevels)**2, shape=Ncases)
    
    # Latent standard deviations of the ratings. Hierarchical structure.
    sigmaSD = pm.Gamma('sigmaSD', gammaShRa[0], gammaShRa[1])
    sigmaMode = pm.Gamma('sigmaMode', gammaShRa[0], gammaShRa[1])
    sigmaRa = pm.Deterministic('sigmaRa', ((sigmaMode + pm.math.sqrt(sigmaMode**2 + 4*sigmaSD**2)) / (2*sigmaSD**2)))
    sigmaSh = pm.Deterministic('sigmaSh', 1 + sigmaMode*sigmaRa)
    sigma = pm.Gamma('sigma', sigmaSh, sigmaRa, shape=Ncases)
    
      
    # Latent thresholds between the ratings (ordinal values)
    theta = pm.Normal('theta', mu=thresh, tau=1/np.repeat(2**2, len(thresh)),
                      shape=len(thresh), observed=thresh_obs)#, testval=thresh)
    
    # Cumulative normal probabilities for ratings (ordinal values)
    pr = outcome_probabilities(theta, mu, sigma)
    
    # Likelihood
    out = pm.Multinomial('out', n=z.values.reshape(-1,1), p=pr.T, observed=y.values)  
with ordinal_model_multi_groups:
    step1 = pm.Slice([mu, theta, sigma])
    step2 = pm.NUTS([sigmaMode, sigmaSD])
    trace = pm.sample(10000, step=[step1, step2], cores=1, progressbar=True)    
pm.traceplot(trace, ['mu', 'theta_missing', 'sigmaSD', 'sigmaMode', 'sigma']);
pm.plot_posterior(trace['theta_missing'][:,0], point_estimate='mode');
pm.plot_posterior(trace['theta_missing'][:,1], point_estimate='mode');
pm.plot_posterior(trace['sigmaMode'], point_estimate='mode');
pm.plot_posterior(trace['sigmaSD'], point_estimate='mode');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1555443208628244, 'fraction_numerical': 0.02510953825412875, 'mean_word_length': 3.9623745819397995, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 23, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '192495', 'license': 'MIT', 'n_tokens_mistral': 2315, 'n_tokens_neox': 2138, 'n_words': 522, 'path': 'Notebooks/Ordinal Model_Kruschke_Liddell.ipynb', 'repo_name': 'steinszzh/DAV-5400', 'size': '5958'}