antivax-attitudes

Reanalyses of data from Horne, Powell, Hummel & Holyoak (2015)
git clone https://git.eamoncaddigan.net/antivax-attitudes.git
Log | Files | Refs | README | LICENSE

antivax-attitudes.R (4182B)

      1 # Note: this code has been pulled over to antivax-attitudes.Rmd and subsequently
      2 # modified. I'm keeping the file around for archival purposes.
      3 
      4 # Yeah, the standard Hadley stack. /sigh
      5 library(readxl)
      6 library(tidyr)
      7 library(dplyr)
      8 library(ggplot2)
      9 
     10 
     11 # Gather and clean the data -----------------------------------------------
     12 
     13 # Generates warnings, but only for the Ps who didn't do day 2
     14 expData <- read_excel("Vacc_HPHH_publicDataset.xlsx", sheet = 2)
     15 
     16 # Just a bit of data cleaning
     17 expData.clean <- expData %>%
     18   # Add a subject number
     19   mutate(subject_number = 1:nrow(.)) %>%
     20   # Just exclude Ps who didn't do day 2 and failed the attention checks
     21   filter(Returned == 1,
     22          `AttentionCheck_PostTest (if = 4 then include)` == 4,
     23          `AttentionChecks_Sum(include if = 4)` == 4,
     24          Paid_Attention == 1)
     25 
     26 # Get all the dependent measures into a DF
     27 questionnaireData <- expData.clean %>%
     28   # pull out the columns and use BETTER NAMES (jeez Zach)
     29   select(subject_number,
     30          intervention = Condition,
     31          pretest.healthy = Healthy_VaxscalePretest,
     32          posttest.healthy = Healthy_VaxscalePosttest,
     33          pretest.diseases = Diseases_VaxScalePretest,
     34          posttest.diseases = Diseases_VaxScalePosttest,
     35          pretest.doctors = Doctors_VaxScalePreTest,
     36          posttest.doctors = Doctors_VaxScalePostTest,
     37          pretest.side_effects = Sideeffects_VaxScalePreTest,
     38          posttest.side_effects = Sideeffects_VaxScalePostTest,
     39          pretest.plan_to = Planto_VaxScalePreTest,
     40          posttest.plan_to = Planto_VaxScalePostTest) %>%
     41   # reverse-code the approrpiate columns
     42   mutate(pretest.diseases = 7 - pretest.diseases,
     43          posttest.diseases = 7 - posttest.diseases,
     44          pretest.side_effects = 7 - pretest.side_effects,
     45          posttest.side_effects = 7 - posttest.side_effects) %>%
     46   # "tidy" the data
     47   gather("question", "response", -subject_number, -intervention) %>%
     48   separate(question, c("interval", "question"), sep = "\\.") %>% 
     49   mutate(interval = factor(interval, c("pretest", "posttest"), ordered = TRUE),
     50          question = factor(question))
     51 
     52 
     53 # Some plots --------------------------------------------------------------
     54 
     55 # Check out the distribution of responses before and after the intervention
     56 p1 <- ggplot(questionnaireData, aes(x = question, y = response, fill = interval)) +
     57   geom_violin() + 
     58   facet_grid(intervention ~ .)
     59 print(p1)
     60 
     61 # Look at each subject's change for each question
     62 p2 <- ggplot(questionnaireData, aes(x = interval, y = response, group = subject_number)) + 
     63   geom_line(alpha = 0.2, position = position_jitter(w = 0.15, h = 0.15)) + 
     64   facet_grid(intervention ~ question)
     65 print(p2)
     66 
     67 
     68 # Bayesian analysis of survey data ----------------------------------------
     69 
     70 # Fit a model to each question using pre-intervention data. 
     71 modelData <- questionnaireData
     72 
     73 source("Jags-Yord-Xnom1grp-Mnormal.R")
     74 fileNameRoot = "antivax-mcmc"
     75 
     76 mcmcCoda <- genMCMC(datFrm = modelData,
     77                     yName = "response",
     78                     x1Name = "question",
     79                     x2Name = "intervention",
     80                     x3Name = "interval",
     81                     numSavedSteps = 15000,
     82                     thinSteps = 10,
     83                     saveName = fileNameRoot)
     84 
     85 # Display diagnostics of chain, for specified parameters:
     86 # betas, sigmas
     87 parameterNames <- varnames(mcmcCoda)
     88 parameterNames <- parameterNames[grepl("^b", parameterNames) | grepl("^sigma", parameterNames)]
     89 for (parName in parameterNames) {
     90   diagMCMC(codaObject = mcmcCoda,
     91            parName = parName, 
     92            saveName = fileNameRoot,
     93            saveType = "png")
     94 }
     95 
     96 # Display posterior information:
     97 plotMCMC(mcmcCoda,
     98          datFrm = modelData,
     99          yName = "response",
    100          qName = "question",
    101          compVal = 3.5, 
    102          saveName = fileNameRoot,
    103          saveType = "png")
    104 
    105 # This right here is the good stuff:
    106 mcmcMat <- as.matrix(mcmcCoda)
    107 plotPost((mcmcMat[, "b2b3[3,2]"] - mcmcMat[, "b2b3[3,1]"]) - 
    108            (mcmcMat[, "b2b3[1,2]"] - mcmcMat[, "b2b3[1,1]"]))
    109 # It shows that post-pre for the "disease risk" is greater than post-pre for
    110 # "autism correction", supporting the authors' findings.