commit ddc983158d7b79042ea81d6ea68070e56304050a
parent 0d6df9df4add1f4b4502cb942e5b61820c8a1aae
Author: eamoncaddigan <eamon.caddigan@gmail.com>
Date: Sun, 30 Aug 2015 14:31:08 -0400
Tweaks
Diffstat:
| M | antivax-attitudes.Rmd | | | 186 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1 file changed, 186 insertions(+), 0 deletions(-)
diff --git a/antivax-attitudes.Rmd b/antivax-attitudes.Rmd
@@ -87,3 +87,189 @@ print(p2)
```
The above figure shows pre- and post-intervention responses to each question. Each line represents represents a single participant's responses before and after the intervention to a single question. Lines are colored by the magnitude of the change in response; blue lines indicate more agreement (toward a more pro-vaccine stance) and red lines indicate less agreement (a more anti-vaccine stance).
+
+```{r, echo=FALSE}
+# Get the data ready for JAGS
+x1 <- as.numeric(as.factor(questionnaireData[["question"]]))
+Nx1Lvl <- max(x1)
+x2 <- as.numeric(as.factor(questionnaireData[["intervention"]]))
+Nx2Lvl <- max(x2)
+x3 <- as.numeric(as.factor(questionnaireData[["interval"]]))
+Nx3Lvl <- max(x3)
+y <- as.numeric(questionnaireData[["response"]])
+Ntotal <- length(y)
+nYlevels <- max(y)
+
+# Threshold 1 and nYlevels-1 are fixed; other thresholds are estimated.
+# This allows all parameters to be interpretable on the response scale.
+thresh <- matrix(data = NA, nrow = Nx1Lvl, ncol = nYlevels-1)
+thresh[, 1] <- 1 + 0.5
+thresh[, nYlevels-1] <- nYlevels-1 + 0.5
+# Specify the data in a list, for later shipment to JAGS:
+dataList <- list(
+ x1 = x1,
+ Nx1Lvl = Nx1Lvl,
+ x2 = x2,
+ Nx2Lvl = Nx2Lvl,
+ x3 = x3,
+ Nx3Lvl = Nx3Lvl,
+ y = y,
+ NyLvl = nYlevels,
+ thresh = thresh,
+ Ntotal = Ntotal
+)
+
+# Prepare the model for JAGS
+modelString <- "
+ model {
+ for (i in 1:Ntotal) {
+ # Thresholded cummulative normal distribution
+ 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,NyLvl] <- 1 - pnorm(thresh[x1[i], NyLvl-1] , mu[i] , 1/sigma[x1[i]]^2)
+
+ # mu ~ x1*x2*x3
+ mu[i] <- a0 + a1[x1[i]] + a2[x2[i]] + a3[x3[i]] +
+ a1a2[x1[i], x2[i]] + a1a3[x1[i], x3[i]] + a2a3[x2[i], x3[i]] +
+ a1a2a3[x1[i], x2[i], x3[i]]
+ }
+
+ a0 ~ dnorm((1+NyLvl)/2, 1/(NyLvl)^2)
+
+ for (j1 in 1:Nx1Lvl) {
+ # Constant sigma for beta1, we're treating all Qs as independent
+ a1[j1] ~ dnorm(0.0, 1/(NyLvl)^2)
+
+ # Sigma for normal CDF, unique for each x1.
+ sigma[j1] ~ dunif(NyLvl/1000, NyLvl*10)
+
+ # Threshold distributions. 1 and NyLvl-1 are fixed, not stochastic
+ for (k in 2:(NyLvl-2)) {
+ thresh[j1, k] ~ dnorm(k+0.5, 1/2^2)
+ }
+ }
+
+ # Constant sigma for beta2, the interventions are independent
+ for (j2 in 1:Nx2Lvl) {
+ a2[j2] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+
+ # Constant sigma for beta3 (for now)
+ for (j3 in 1:Nx3Lvl) {
+ a3[j3] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+
+ # Interaction terms also has homogenous variance (for now)
+ for (j1 in 1:Nx1Lvl) {
+ for (j2 in 1:Nx2Lvl) {
+ a1a2[j1, j2] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+ }
+ for (j1 in 1:Nx1Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ a1a3[j1, j3] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+ }
+ for (j2 in 1:Nx2Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ a2a3[j2, j3] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+ }
+ for (j1 in 1:Nx1Lvl) {
+ for (j2 in 1:Nx2Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ a1a2a3[j1, j2, j3] ~ dnorm(0.0, 1/(NyLvl)^2)
+ }
+ }
+ }
+
+ # Compute cell means
+ for (j1 in 1:Nx1Lvl) {
+ for (j2 in 1:Nx2Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ m[j1, j2, j3] <- a0 + a1[j1] + a2[j2] + a3[j3] +
+ a1a2[j1, j2] + a1a3[j1, j3] + a2a3[j2, j3] +
+ a1a2a3[j1, j2, j3]
+ }
+ }
+ }
+
+ # Convert a0, a1[], a2[], &c. to sum-to-zero b0, b1[], b2[], &c.
+ b0 <- mean(m[1:Nx1Lvl, 1:Nx2Lvl, 1:Nx3Lvl])
+ for (j1 in 1:Nx1Lvl) {
+ b1[j1] <- mean(m[j1, 1:Nx2Lvl, 1:Nx3Lvl]) - b0
+ }
+ for (j2 in 1:Nx2Lvl) {
+ b2[j2] <- mean(m[1:Nx1Lvl, j2, 1:Nx3Lvl]) - b0
+ }
+ for (j3 in 1:Nx3Lvl) {
+ b3[j3] <- mean(m[1:Nx1Lvl, 1:Nx2Lvl, j3]) - b0
+ }
+ for (j1 in 1:Nx1Lvl) {
+ for (j2 in 1:Nx2Lvl) {
+ b1b2[j1, j2] <- mean(m[j1, j2, 1:Nx3Lvl]) - (b0 + b1[j1] + b2[j2])
+ }
+ }
+ for (j1 in 1:Nx1Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ b1b3[j1, j3] <- mean(m[j1, 1:Nx2Lvl, j3]) - (b0 + b1[j1] + b3[j3])
+ }
+ }
+ for (j2 in 1:Nx2Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ b2b3[j2, j3] <- mean(m[1:Nx1Lvl, j2, j3]) - (b0 + b2[j2] + b3[j3])
+ }
+ }
+ for (j1 in 1:Nx1Lvl) {
+ for (j2 in 1:Nx2Lvl) {
+ for (j3 in 1:Nx3Lvl) {
+ b1b2b3[j1, j2, j3] <- m[j1, j2, j3] - (b0 + b1[j1] + b2[j2] + b3[j3] +
+ b1b2[j1, j2] + b1b3[j1, j3] + b2b3[j2, j3])
+ }
+ }
+ }
+ }
+" # close quote for modelString
+print(modelString)
+# Write out modelString to a text file
+writeLines(modelString , con="TEMPmodel.txt")
+
+# Tell JAGS which parameters to return
+parameters <- c("b0", "b1", "b2", "b3", "b1b2", "b1b3", "b2b3", "b1b2b3",
+ "sigma", "thresh")
+
+# JAGS parameters. We'll let it iniaialize itself
+initsList <- NULL
+adaptSteps <- 500 # Number of steps to "tune" the samplers
+burnInSteps <- 1000
+numSavedSteps <- 15000
+thinSteps <- 10
+nChains <- 4
+fileNameRoot <- "antivax-mcmc"
+
+# Since running JAGS takes forever, we'll skip redoing it every time we knit.
+saveName <- paste0(fileNameRoot, "-coda.Rdata")
+if (file.exists(saveName)) {
+ load(saveName)
+} else {
+ runJagsOut <- run.jags(method="parallel",
+ model="TEMPmodel.txt",
+ monitor=parameters,
+ data=dataList,
+ #inits=initsList,
+ n.chains=nChains,
+ adapt=adaptSteps,
+ burnin=burnInSteps,
+ sample=ceiling(numSavedSteps/nChains),
+ thin=thinSteps,
+ summarise=FALSE,
+ plots=FALSE)
+ codaSamples <- as.mcmc.list(runJagsOut)
+ save(codaSamples, file=saveName)
+}
+```
+
