-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathOneOddGroupModelComp2E_Stan.R
executable file
·177 lines (154 loc) · 6.69 KB
/
OneOddGroupModelComp2E_Stan.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# OneOddGroupModelComp2E.R
# Accompanies the book:
# Kruschke, J. K. (2014). Doing Bayesian Data Analysis:
# A Tutorial with R, JAGS, and Stan. 2nd Edition. Academic Press / Elsevier.
# Adapted for Stan by Joe Houpt
graphics.off()
rm(list=ls(all=TRUE))
source("DBDA2E-utilities.R")
#require(rjags)
#require(runjags)
fileNameRoot="OneOddGroupModelComp2E-"
#------------------------------------------------------------------------------
# THE DATA.
# Randomly generated fictitious data.
# For each subject, specify the condition s/he was in,
# the number of trials s/he experienced, and the number correct.
npg = 20 # number of subjects per group
ntrl = 20 # number of trials per subject
cond_of_subj = c( rep(1,npg) , rep(2,npg) , rep(3,npg) , rep(4,npg) )
n_trl_of_subj = rep( ntrl , 4*npg )
set.seed(47405)
condMeans = c(.40,.50,.51,.52)
nCorrOfSubj = c( rbinom(npg,ntrl,condMeans[1]) , rbinom(npg,ntrl,condMeans[2]) ,
rbinom(npg,ntrl,condMeans[3]) , rbinom(npg,ntrl,condMeans[4]) )
n_cond = length(unique(cond_of_subj))
n_subj = length(cond_of_subj)
# jitter the data to be as close as possible to desired condition means:
for ( cIdx in 1:n_cond ) {
nToAdd = round(condMeans[cIdx]*npg*ntrl)-sum(nCorrOfSubj[cond_of_subj==cIdx])
if ( nToAdd > 0 ) {
for ( i in 1:nToAdd ) {
thisNcorr = ntrl
while ( thisNcorr == ntrl ) {
randSubjIdx = sample(which(cond_of_subj==cIdx),size=1)
thisNcorr = nCorrOfSubj[randSubjIdx]
}
nCorrOfSubj[randSubjIdx] = nCorrOfSubj[randSubjIdx]+1
}
}
if ( nToAdd < 0 ) {
for ( i in 1:abs(nToAdd) ) {
thisNcorr = 0
while ( thisNcorr == 0 ) {
randSubjIdx = sample(which(cond_of_subj==cIdx),size=1)
thisNcorr = nCorrOfSubj[randSubjIdx]
}
nCorrOfSubj[randSubjIdx] = nCorrOfSubj[randSubjIdx]-1
}
}
}
show( aggregate( nCorrOfSubj , by=list(cond_of_subj) , FUN=mean ) / ntrl )
# Package the data:
dataList = list(
n_cond = n_cond ,
n_subj = n_subj ,
cond_of_subj = cond_of_subj ,
n_trl_of_subj = n_trl_of_subj ,
n_corr_of_subj = nCorrOfSubj
)
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Let Stan do it...
#------------------------------------------------------------------------------
# RUN THE CHAINS.
parameters = c("omega","kappa","omega0","theta","model_prob_1")
adaptSteps = 1000 # Number of steps to "tune" the samplers.
burnInSteps = 5000 # Number of steps to "burn-in" the samplers.
nChains = 3 # Number of chains to run.
numSavedSteps=12000 # Total number of steps in chains to save.
thinSteps=10 # Number of steps to "thin" (1=keep every step).
# Translate to C++ and compile to DSO:
stanDso <- stan_model( file="OneOddGroupModelComp.stan" )
stanFit <- sampling( object=stanDso ,
data = dataList ,
pars = parameters , # optional
chains = nChains ,
iter = ( ceiling(numSavedSteps/nChains)*thinSteps
+burnInSteps ) ,
warmup = burnInSteps ,
thin = thinSteps ,
init = "random" ) # optional
# Or, accomplish above in one "stan" command; note stanDso is not separate.
# For consistency with JAGS-oriented functions in DBDA2E collection,
# convert stan format to coda format:
codaSamples = mcmc.list( lapply( 1:ncol(stanFit) ,
function(x) { mcmc(as.array(stanFit)[,x,]) } ) )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
save( codaSamples , file=paste(fileNameRoot,"Mcmc.Rdata",sep="") )
save( stanFit , file=paste(fileNameRoot,"StanFit.Rdata",sep="") )
save( stanDso , file=paste(fileNameRoot,"StanDso.Rdata",sep="") )
#-------------------------------------------------------------------------------
# Display diagnostics of chain:
parameterNames = varnames(codaSamples) # get all parameter names
show(parameterNames)
for ( parName in c("model_prob_1","omega[1]","omega0","kappa[1]","theta[1]") ) {
diagMCMC( codaSamples , parName=parName ,
saveName=fileNameRoot , saveType="eps" )
}
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.
mcmcMat = as.matrix(codaSamples,chains=TRUE)
xLim=c(0.35,0.75)
# Display the model index
pM1 = mcmcMat[, "model_prob_1" ]
pM2 = 1 - pM1
string1 =paste("p( Diff Omega M1 | D )=",round(mean(pM1),3),sep="")
string2 =paste("p( Same Omega M2 | D )=",round(mean(pM2),3),sep="")
openGraph(10,4)
nStepsToPlot = 1000
plot( 1:nStepsToPlot , pM1[1:nStepsToPlot] , type="l" , lwd=2 , ylim=c(0,1),
xlab="Step in Markov chain" , ylab="Model Index (1, 2)" ,
main=paste(string1,", ",string2,sep="") , col="skyblue" )
saveGraph(file=paste0(fileNameRoot,"model_prob_1"),type="eps")
# Display the omega0 posterior
omega0sample = mcmcMat[, "omega0" ]
openGraph()
#layout( matrix(1:2,nrow=2) )
plotPost( omega0sample , main="Posterior for Same Omega" ,
xlab=expression(omega[0]) , xlim=xLim )
saveGraph(file=paste0(fileNameRoot,"Omega0"),type="eps")
# Display the omega[j] posterior
omega1sample = mcmcMat[, "omega[1]" ]
omega2sample = mcmcMat[, "omega[2]" ]
omega3sample = mcmcMat[, "omega[3]" ]
omega4sample = mcmcMat[, "omega[4]" ]
openGraph(10,5)
layout( matrix(1:4,nrow=1,byrow=T) )
plotPost( omega1sample , main="Posterior for Diff Omega" ,
xlab=expression(omega[1]) , xlim=xLim )
plotPost( omega2sample , main="Posterior for Diff Omega" ,
xlab=expression(omega[2]) , xlim=xLim )
plotPost( omega3sample , main="Posterior for Diff Omega" ,
xlab=expression(omega[3]) , xlim=xLim )
plotPost( omega4sample , main="Posterior for Diff Omega" ,
xlab=expression(omega[4]) , xlim=xLim )
saveGraph(file=paste0(fileNameRoot,"OmegaCond"),type="eps")
# Display the differences of omega[j]'s
omegaSample = rbind( omega1sample , omega2sample , omega3sample , omega4sample )
openGraph(10,5)
layout( matrix(1:6,nrow=2,ncol=3,byrow=T) )
xmin = -0.25
xmax = 0.25
for ( i in 1:3 ) {
for ( j in (i+1):4 ) {
plotPost( omegaSample[i,]-omegaSample[j,] , compVal=0.0 ,
xlab=bquote(omega[.(i)]-omega[.(j)]) ,
#breaks=unique( c( min(c(xmin,omegaSample[i,]-omegaSample[j,])),
# seq(xmin,xmax,len=20),
# max(c(xmax,omegaSample[i,]-omegaSample[j,])) )) ,
main="" , xlim=c(xmin,xmax) )
}
}
saveGraph(file=paste0(fileNameRoot,"OmegaDiff"),type="eps")