commit 3b3a2ecbd59c5758d1e9afcf3f6ef749bf22a2aa
parent dc432e2ab8444be6eca9052e232bd19ebc025c03
Author: eamoncaddigan <eamon.caddigan@gmail.com>
Date: Sat, 22 Aug 2015 22:46:56 -0400
Model changes to prepare for the addition of a subject parameter and eventual intervention parameter.
Diffstat:
2 files changed, 32 insertions(+), 18 deletions(-)
diff --git a/Jags-Yord-Xnom1grp-Mnormal.R b/Jags-Yord-Xnom1grp-Mnormal.R
@@ -36,10 +36,10 @@ genMCMC = function( datFrm, yName , qName,
thresh[, nYlevels-1] = nYlevels-1 + 0.5
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
- q = q,
- nQlevels = nQlevels,
+ x1 = q,
+ Nx1Lvl = nQlevels,
y = y,
- nYlevels = nYlevels,
+ NyLvl = nYlevels,
thresh = thresh,
Ntotal = Ntotal
)
@@ -48,25 +48,38 @@ genMCMC = function( datFrm, yName , qName,
modelString = "
model {
for (i in 1:Ntotal) {
- y[i] ~ dcat(pr[i,1:nYlevels])
- pr[i,1] <- pnorm(thresh[q[i], 1] , mu[q[i]] , 1/sigma[q[i]]^2)
- for (k in 2:(nYlevels-1)) {
- pr[i,k] <- max( 0 , pnorm(thresh[q[i], k] , mu[q[i]] , 1/sigma[q[i]]^2 )
- - pnorm(thresh[q[i], k-1] , mu[q[i]] , 1/sigma[q[i]]^2 ))
+ y[i] ~ dcat(pr[i,1:NyLvl])
+ pr[i,1] <- pnorm(thresh[x1[i], 1] , mu[i] , 1/sigma[x1[i]]^2)
+ for (k in 2:(NyLvl-1)) {
+ pr[i,k] <- max( 0 , pnorm(thresh[x1[i], k] , mu[i] , 1/sigma[x1[i]]^2 )
+ - pnorm(thresh[x1[i], k-1] , mu[i] , 1/sigma[x1[i]]^2 ))
}
- pr[i,nYlevels] <- 1 - pnorm(thresh[q[i], nYlevels-1] , mu[q[i]] , 1/sigma[q[i]]^2)
+ pr[i,NyLvl] <- 1 - pnorm(thresh[x1[i], NyLvl-1] , mu[i] , 1/sigma[x1[i]]^2)
+ mu[i] <- a0 + a1[x1[i]]
}
- # Unique mu, sigma, and thresh vector for each question
- for (j in 1:nQlevels) {
- mu[j] ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
- sigma[j] ~ dunif( nYlevels/1000 , nYlevels*10 )
+ a0 ~ dnorm((1+NyLvl)/2, 1/(NyLvl)^2)
+ for (j in 1:Nx1Lvl) {
+ # Constant sigma for beta1, we're treating all Qs as independent
+ a1[j] ~ dnorm(0.0 , 1/(NyLvl)^2)
- # Threshold distributions. 1 and nYlevels-1 are fixed, not stochastic
- for (k in 2:(nYlevels-2)) {
+ # Sigma for normal CDF, unique for each Q.
+ sigma[j] ~ dunif(NyLvl/1000 , NyLvl*10)
+
+ # Threshold distributions. 1 and NyLvl-1 are fixed, not stochastic
+ for (k in 2:(NyLvl-2)) {
thresh[j, k] ~ dnorm( k+0.5 , 1/2^2 )
}
}
+
+ # Convert a0,a[] to sum-to-zero b0,b[] :
+ for (j in 1:Nx1Lvl) {
+ m[j] <- a0 + a1[j]
+ }
+ b0 <- mean(m[1:Nx1Lvl])
+ for (j in 1:Nx1Lvl) {
+ b1[j] <- m[j] - b0
+ }
}
" # close quote for modelString
# Write out modelString to a text file
@@ -76,7 +89,7 @@ genMCMC = function( datFrm, yName , qName,
initsList = NULL
#-----------------------------------------------------------------------------
# RUN THE CHAINS
- parameters = c( "mu" , "sigma" , "thresh" )
+ parameters = c("b0", "b1", "sigma", "thresh")
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
runJagsOut <- run.jags( method=runjagsMethod ,
diff --git a/antivax-attitudes.R b/antivax-attitudes.R
@@ -78,9 +78,9 @@ mcmcCoda <- genMCMC(datFrm = modelData,
saveName = fileNameRoot)
# Display diagnostics of chain, for specified parameters:
-# (everything)
+# betas, sigmas
parameterNames <- varnames(mcmcCoda)
-parameterNames <- parameterNames[grepl("^mu", parameterNames) | grepl("^sigma", parameterNames)]
+parameterNames <- parameterNames[grepl("^b", parameterNames) | grepl("^sigma", parameterNames)]
for (parName in parameterNames) {
diagMCMC(codaObject = mcmcCoda,
parName = parName,
@@ -89,6 +89,7 @@ for (parName in parameterNames) {
}
# Display posterior information:
+# XXX - broken by model changes
plotMCMC(mcmcCoda,
datFrm = modelData,
yName = "response",
