######################################################
#Lecture 6 Script
#Missing data
#Anastasia Ejova 24/11/2014
#R101/PSY532
######################################################

#CREATING AN SS DATA SET WITH MISSING VALUES 
#(5% of values not measured pre-experimentally and not denoting SeqCond will be missing)
#The data set is already in your workspace under the name SS3. This code just shows 
#you how it was created.

#Reduce SS to all variables except those that should have some missing values. Also
#remove the variables that make up the total illusion of control measures. We will
#re-introduce them later. For now, we will assume that the illusion of control measures 
#(PostNaturalIoC and PostSupIoC) are missing because of some missing values in their components 
#(the PostList variables). Note that, if your data set includes both a total score (e.g., PostSupIoC)
#and the items constituting it (e.g., PostListLuckyPlay, PostListDeserve, etc.), you will
#have a "singularity" problem when using the norm and mi packages:
#https://lists.gking.harvard.edu/pipermail/amelia/2013-July/001030.html
drops <- c("PNo","PostStrategyDescr",
           "PostOutcomeSeqOther", "ChanceCat",
           "PostHowManySingle", "PostListSkill",
           "PostListLuckyMoments", "PostListPlayers", "PostListStrategy", 
           "PostListLuckyPlay", "PostListLuck", "PostListLuckyPerson",
           "PostListSoccerKnowl", "PostListPractice", "PostListComputerGames",
           "PostListDeserve", "PostListGoalieMoves", "PostListKickDirn", "PreDBC_Total") 

SS3 <- SS[,!(names(SS) %in% drops)]

alpha <- .05 #probability of missing data
N <- 334

rbinom(n = N, size = 1, prob=alpha) #As an illustration of what we are doing in lines 34-40

r.PostNaturalIoC.mcar<-rbinom(n = N, size = 1, prob=alpha)
r.PostSupIoC.mcar<-rbinom(N,1,prob=alpha)
r.LChance.mcar<-rbinom(N,1,prob=alpha)
r.Strat.mcar<-rbinom(N,1,prob=alpha)
r.Goal.mcar<-rbinom(N,1,prob=alpha)
r.Single.mcar<-rbinom(N,1,prob=alpha)
r.Seq.mcar<-rbinom(N,1,prob=alpha)

SS3$PostStrategyPresent<-as.numeric(SS3$PostStrategyPresent)
SS3$PostOutcomeSeq<-as.numeric(SS3$PostOutcomeSeq)

#x.mcar=x1 if not missing, 99999 if missing:
SS3$PostNaturalIoC<-SS3$PostNaturalIoC*(1-r.PostNaturalIoC.mcar)+r.PostNaturalIoC.mcar*99999  
SS3$PostSupIoC<-SS3$PostSupIoC*(1-r.PostSupIoC.mcar)+r.PostSupIoC.mcar*99999  
SS3$PostListCHANCE<-SS3$PostListCHANCE*(1-r.LChance.mcar)+r.LChance.mcar*99999
SS3$PostStrategyPresent<-SS3$PostStrategyPresent*(1-r.Strat.mcar)+r.Strat.mcar*99999
SS3$Post100RoundsGoal<-SS3$Post100RoundsGoal*(1-r.Goal.mcar)+r.Goal.mcar*99999
SS3$PostHowManySingleWins<-SS3$PostHowManySingleWins*(1-r.Single.mcar)+r.Single.mcar*99999
SS3$PostOutcomeSeq<-SS3$PostOutcomeSeq*(1-r.Seq.mcar)+r.Seq.mcar*99999

SS3[SS3 == 99999] <- NA

SS3$PostStrategyPresent<-as.factor(SS3$PostStrategyPresent)
SS3$PostOutcomeSeq<-as.factor(SS3$PostOutcomeSeq)

levels(SS3$PostStrategyPresent) <- levels(SS$PostStrategyPresent)
levels(SS3$PostOutcomeSeq) <- levels(SS$PostOutcomeSeq)

#Add interaction terms
#Since we know from our own analyses and past research that they predict illusion of control
SS3$SeqCapt <- interaction(SS3$SeqCond,SS3$PostHowManySingleCaptionType, drop = TRUE) 
SS3$SeqBels <- interaction(SS3$SeqCond,SS3$PreDBC_IOC, drop = TRUE) 
SS3$SeqBelsCapt <- interaction(SS3$SeqCond,SS3$PreDBC_IOC,SS3$PostHowManySingleCaptionType, 
                               drop = TRUE) 

######################################################

#EXPLORING THE MISSING DATA PATTERNS
#install.packages("mice")
library(mice)
md.pattern(SS3)
#Each row represents a missing data pattern (1 = observed, 0 = missing).
#There are 239 (out of 334) rows that are complete, and
#15 rows in which only PostListCHANCE is missing (etc., etc.)
#The hightest number of missing values (N = 21) is in PostListCHANCE.

#To see the table better
miss <- md.pattern(SS3)
write.table(miss, file = "MissTable.txt", sep = ",", dec = ".") 

######################################################

#DETERMINING WHETHER DATA IS MISSING COMPLETELY AT RANDOM (MCAR)
#It should be in this case, because we've created a data set with MCAR from a 
#data set with no missing values.

#Little's MCAR test
#install.packages("BaylorEdPsych")
#install.packages("mvnmle") #necessary for LittleMCAR function
library(BaylorEdPsych)
library("mvnmle")
littleMCAR <- LittleMCAR(SS3) #be prepared, the calculations take some time

littleMCAR$chi.square
  #342.6976

littleMCAR$df
  #343

littleMCAR$p.value
  #0.494451 Non-significant result implies that the data is MCAR

#--

#Visual method
#install.packages("VIM")
library(VIM)
marginplot(SS3[, c("PreDBC_Sup", "PostSupIoC")]) 
#Red points: Missing on X or Y, but not both
#Blue points: Not missing X or Y
#Under MCAR, the red and blue boxplots should mostly overlap
#The plot can also be useful to determine the variables in the data set responsible
#for the missing data, as assumed under the MAR assumption. Here, prior beliefs
#in supernatural phenomena do not predict missingness in scores on the illusion of 
#supernatural control.
#http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm 
#van Buuren & Karin Groothuis-Oudshoornmice (2011) "mice step-by-step" p10

######################################################

#MULTIPLE IMPUTATION USING the norm package
#install.packages("norm")
library(norm)

#Step 1 (on Slide 8)
SS4 <- data.matrix(SS3) #norm accepts only data matrices

s <- prelim.norm(SS4) #do preliminary manipulations

thetahat <- em.norm(s) #find the maximum likelihood estimate for a starting value

rngseed(1234567) #set random number generator seed

theta <- da.norm(s, start = thetahat, steps=20, showits=TRUE) 
#Take 20 steps to create a theta matrix to insert into imputation function at line 126
#What is this theta matrix? Let's take a look:

getparam.norm(s,theta)
#mu is the mean parameter necessary for defnining a regression equation
#sigma is a set of covariance parameters

#Create 5 imputed data sets
imputed1 <- imp.norm(s, theta, SS4)
DBCTot <- imputed1[,"PreDBC_Sup"] + imputed1[,"PreDBC_IOC"] #Can calculate scale totals
                                                  #(e.g., PreDBC_Total) by adding items
                                                  #within the imputed dataset. Line 135
                                                  #adds the calculated total to the imputed
                                                  #dataset.
imputed1 <- cbind(imputed1, DBCTot)

theta <- da.norm(s, start = thetahat, steps=20, showits=TRUE)
imputed2 <- imp.norm(s, theta, SS4)
DBCTot <- imputed2[,"PreDBC_Sup"] + imputed2[,"PreDBC_IOC"]
imputed2 <- cbind(imputed2, DBCTot)

theta <- da.norm(s, start = thetahat, steps=20, showits=TRUE)
imputed3 <- imp.norm(s, theta, SS4)
DBCTot <- imputed3[,"PreDBC_Sup"] + imputed3[,"PreDBC_IOC"]
imputed3 <- cbind(imputed3, DBCTot)

theta <- da.norm(s, start = thetahat, steps=20, showits=TRUE)
imputed4 <- imp.norm(s, theta, SS4)
DBCTot <- imputed4[,"PreDBC_Sup"] + imputed4[,"PreDBC_IOC"]
imputed4 <- cbind(imputed4, DBCTot)

theta <- da.norm(s, start = thetahat, steps=20, showits=TRUE)
imputed5 <- imp.norm(s, theta, SS4)
DBCTot <- imputed5[,"PreDBC_Sup"] + imputed5[,"PreDBC_IOC"]
imputed5 <- cbind(imputed5, DBCTot)

imputed1 <- data.frame(imputed1)
imputed2 <- data.frame(imputed2)
imputed3 <- data.frame(imputed3)
imputed4 <- data.frame(imputed4)
imputed5 <- data.frame(imputed5)

#Steps 2 and 3 on Slide 8
#install.packages("mitools")
library(mitools)
imputedDA <- imputationList(list(imputed1, imputed2, imputed3, imputed4, imputed5))
lmDA <- with(imputedDA, lm(PostNaturalIoC ~ SeqCond*PreDBC_IOC)) #run ANOVA
summary(MIcombine(lmDA)) #view pooled ANOVA results
  #lower and upper confidence intervals suggest significant effect of background beliefs
  #but not success-slope

#Compare to lm result for full SS data set:
lmFull <- aov(PostNaturalIoC ~ SeqCond*PreDBC_IOC, data = SS)
summary(lmFull)

#Run generalised linear model... and this is where we encounter a problem
#See Slide 12. Because we have some negative values for Post100RoundsGoal in one of the 
#imputed data sets, the negative binomial and poisson distributions cannot be used to
#describe the distribution of the data.
library(MASS)
glmnbDA <- with(imputedDA, glm.nb(Post100RoundsGoal ~ 
                                 SeqCond*PostHowManySingleCaptionType*DBCTot))

######################################################

#MULTIPLE IMPUTATION THROUGH CHAINED EQUATIONS
#library(mice); already loaded in the data exploration section (line 73)

#Imputations are generated according to the default method, which, for numerical data, 
#is predictive mean matching (pmm). The entries imp$VisitSequence and imp$PredictorMatrix
#are algorithmic options that will be discusseD later. The default number of 
#multiple imputations is equal to m = 5.

?mice #for more on "method", see van Buuren & Karin Groothuis-Oudshoornmice (2011) p16
#Below we will just use the default methods.

SS5 <- SS3[,-(17:18)] #algorigthm cannot cope with the complex interaction terms, 
  #created on lines 67 and 68. These were interactions between a numeric and categorical
  #variable, resulting in many levels (200 or more).
  #Need to probably divide background beliefs into levels and create simpler interactions

imputedCE <- mice(data = SS5, 
                  m = 5, #no of imputed datasets created: default
                  method = c("", "", "", "", "", "", "", #first 7 vars have no missing values
                             "pmm", "pmm", "pmm", #default method for next three vars (numeric)
                             "logreg", #default logistic regr for StrategyPresent (2 factors)
                             "pmm", #Post100RoundsGoal
                             "", #no missing for CaptionType
                             "pmm", #PostHowManySingleWins
                             "polyreg", #PostOutcomeSeq: unordered, >2 factors
                              ""), #last variable (interaction term) has no missing values
                  print=FALSE, #to avoid lenghty print-out on console
                  seed = 23109) #set seed to some large number

print(imputedCE) #confirms methods

#Diagnostic checking
#(1) Visual scan of all imputed values for each variable in all 5 generated data sets
imputedCE$imp$PostListCHANCE #nothing >10, <0
imputedCE$imp$PostNaturalIoC
imputedCE$imp$PostSupIoC
imputedCE$imp$PostStrategyPresent #factor names respected
imputedCE$imp$Post100RoundsGoal
imputedCE$imp$PostHowManySingleWins #nothing >100, <0
imputedCE$imp$PostOutcomeSeq #factor names respected

#(2) Graph of imputed vs. existing values
stripplot(imputedCE) #can't draw just for variables with missing values

#(3) Variable pairs of interest
xyplot(imputedCE, PostSupIoC ~ PreDBC_Sup|.imp)

#(4) Distributions of imputed and original values - check of MCAR/MAR
#van Buuren & Karin Groothuis-Oudshoornmice (2011) Section 4.5
densityplot(imputedCE)

#Data analysis with imputed values
lmCE <- with(imputedCE, lm(PostNaturalIoC ~ SeqCond*PreDBC_IOC)) 
#with is a mice rather than mitools function

#print(pool(lmCE)) #Pooled based on Rubin's rules

summary(pool(lmCE)) #Significant difference betweeen ascending and descending condition
#which is being used as the base level in dummy coding

#Could check result with 50 datasets:
imputedCEf <- mice(data = SS5, 
                  m = 50,
                  method = c("", "", "", "", "", "", "", #first 7 vars have no missing values
                             "pmm", "pmm", "pmm", #default method for next three vars (numeric)
                             "logreg", #default logistic regr for StrategyPresent (2 factors)
                             "pmm", #Post100RoundsGoal
                             "", #no missing for CaptionType
                             "pmm", #PostHowManySingleWins
                             "polyreg", #PostOutcomeSeq: unordered, >2 factors
                             ""), #last variable (interaction term) has no missing values
                  print=FALSE,
                  seed = 23109)

lmCEf <- with(imputedCEf, lm(PostNaturalIoC ~ SeqCond*PreDBC_IOC))
#print(pool(lmCEf))
summary(pool(lmCEf))

#What about glm with negative binomial distribution?
#using 5 datasets
library(MASS)
glmnbCE <- with(imputedCE, glm.nb(Post100RoundsGoal ~ SeqCond *
                                    PostHowManySingleCaptionType * PreDBC_IOC))
summary(pool(glmnbCE))

#The calculation is possible, but we can't calculate PreDBC_Total variable...

#Introducing a realistic complication - scale totals
#Let's create a data set with missing PostList items, which results in us
#not being able to calculate the scale totals

#This is already saved as SS3a in your Workspace, but the code for creating it 
#is below
drops3a <- c("PNo","PostStrategyDescr",
             "PostOutcomeSeqOther", "ChanceCat", 
             "PostHowManySingle", "PostStrategyType", "PostNaturalIoC",
             "PostSupIoC", "PreDBC_Total") 
SS3a <- SS[,!(names(SS) %in% drops3a)]

#We will assume that only the post-experimental data is missing, with everyone having
#filled in the pre-experimental info
#Note that we did not include List totals as a result of many missing List subscale responses
#http://www.niehs.nih.gov/about/visiting/events/pastmtg/assets/docs_k_m/missingdatar.r
alpha <- .05 #probability of missing data
N <- 334

r.LSkill.mcar<-rbinom(N,1,prob=alpha) #yes/no for whether PostListskill is missing
r.LMoments.mcar<-rbinom(N,1,prob=alpha)
r.LPlayers.mcar<-rbinom(N,1,prob=alpha)
r.LStrategy.mcar<-rbinom(N,1,prob=alpha)
r.LPlay.mcar<-rbinom(N,1,prob=alpha)
r.LLuck.mcar<-rbinom(N,1,prob=alpha)
r.LPers.mcar<-rbinom(N,1,prob=alpha)
r.LKnowl.mcar<-rbinom(N,1,prob=alpha)
r.LPrac.mcar<-rbinom(N,1,prob=alpha)
r.LComp.mcar<-rbinom(N,1,prob=alpha)
r.LDeserve.mcar<-rbinom(N,1,prob=alpha)
r.LGoalie.mcar<-rbinom(N,1,prob=alpha)
r.LKick.mcar<-rbinom(N,1,prob=alpha)
r.LChance.mcar<-rbinom(N,1,prob=alpha)
r.Strat.mcar<-rbinom(N,1,prob=alpha)
r.Goal.mcar<-rbinom(N,1,prob=alpha)
r.Single.mcar<-rbinom(N,1,prob=alpha)
r.Seq.mcar<-rbinom(N,1,prob=alpha)

SS3a$PostStrategyPresent<-as.numeric(SS3a$PostStrategyPresent)
SS3a$PostOutcomeSeq<-as.numeric(SS3a$PostOutcomeSeq)

SS3a$PostListSkill<-SS3a$PostListSkill*(1-r.LSkill.mcar)+r.LSkill.mcar*99999  #x.mcar=x1 if not missing, 99999 if missing
SS3a$PostListLuckyMoments<-SS3a$PostListLuckyMoments*(1-r.LMoments.mcar)+r.LMoments.mcar*99999
SS3a$PostListPlayers<-SS3a$PostListPlayers*(1-r.LPlayers.mcar)+r.LPlayers.mcar*99999
SS3a$PostListStrategy<-SS3a$PostListStrategy*(1-r.LStrategy.mcar)+r.LStrategy.mcar*99999
SS3a$PostListLuckyPlay<-SS3a$PostListLuckyPlay*(1-r.LPlay.mcar)+r.LPlay.mcar*99999
SS3a$PostListLuck<-SS3a$PostListLuck*(1-r.LLuck.mcar)+r.LLuck.mcar*99999
SS3a$PostListLuckyPerson<-SS3a$PostListLuckyPerson*(1-r.LPers.mcar)+r.LPers.mcar*99999
SS3a$PostListSoccerKnowl<-SS3a$PostListSoccerKnowl*(1-r.LKnowl.mcar)+r.LKnowl.mcar*99999
SS3a$PostListPractice<-SS3a$PostListPractice*(1-r.LPrac.mcar)+r.LPrac.mcar*99999
SS3a$PostListComputerGames<-SS3a$PostListComputerGames*(1-r.LComp.mcar)+r.LComp.mcar*99999
SS3a$PostListDeserve<-SS3a$PostListDeserve*(1-r.LDeserve.mcar)+r.LDeserve.mcar*99999
SS3a$PostListGoalieMoves<-SS3a$PostListGoalieMoves*(1-r.LGoalie.mcar)+r.LGoalie.mcar*99999
SS3a$PostListKickDirn<-SS3a$PostListKickDirn*(1-r.LKick.mcar)+r.LKick.mcar*99999
SS3a$PostListCHANCE<-SS3a$PostListCHANCE*(1-r.LChance.mcar)+r.LChance.mcar*99999
SS3a$PostStrategyPresent<-SS3a$PostStrategyPresent*(1-r.Strat.mcar)+r.Strat.mcar*99999
SS3a$Post100RoundsGoal<-SS3a$Post100RoundsGoal*(1-r.Goal.mcar)+r.Goal.mcar*99999
SS3a$PostHowManySingleWins<-SS3a$PostHowManySingleWins*(1-r.Single.mcar)+r.Single.mcar*99999
SS3a$PostOutcomeSeq<-SS3a$PostOutcomeSeq*(1-r.Seq.mcar)+r.Seq.mcar*99999

SS3a[SS3a == 99999] <- NA

SS3a$PostStrategyPresent<-as.factor(SS3a$PostStrategyPresent)
SS3a$PostOutcomeSeq<-as.factor(SS3a$PostOutcomeSeq)

levels(SS3a$PostStrategyPresent) <- levels(SS$PostStrategyPresent)
levels(SS3a$PostOutcomeSeq) <- levels(SS$PostOutcomeSeq)

#Add interaction term 
#Since we intend to include it in subsequent analyses
SS3a$SeqCapt <- interaction(SS3a$SeqCond,SS3a$PostHowManySingleCaptionType,
                            drop = TRUE) 

#Add scale totals
SS3a$PreDBC_Total <- SS3a$PreDBC_Sup + SS3a$PreDBC_IOC 

SS3a$PostNaturalIoC <- (SS3a$PostListSkill + SS3a$PostListPlayers + 
                          SS3a$PostListStrategy + 
                          SS3a$PostListSoccerKnowl + SS3a$PostListPractice + 
                          SS3a$PostListComputerGames + SS3a$PostListGoalieMoves + 
                          SS3a$PostListKickDirn)/8
SS3a$PostSupIoC <- (SS3a$PostListLuckyMoments + SS3a$PostListLuckyPlay +
                      SS3a$PostListLuck + SS3a$PostListLuckyPerson + 
                      SS3a$PostListDeserve)/5

#Using predictor matrix
ini <- mice(SS3a, maxit = 0, print = FALSE) #dry run of mice
ini$method
meth <- ini$method #create methods vector
meth["PreDBC_Total"] <-"~I((PreDBC_Sup + PreDBC_IOC))"
meth["PostNaturalIoC"] <- "~I((PostListSkill + PostListPlayers + \n
  + PostListStrategy + PostListSoccerKnowl + PostListPractice + \n
  + PostListComputerGames + PostListGoalieMoves + PostListKickDirn)/8)" #I operator for passive imputation
meth["PostSupIoC"] <-"~I((PostListLuckyMoments + PostListLuckyPlay + \n
    + PostListLuck + PostListLuckyPerson + PostListDeserve)/5)"

pred <- ini$predictorMatrix
pred["PreDBC_Total",] <- 0 #Do not use totals as predictors
pred["PostNaturalIoC",] <- 0
pred["PostSupIoC",] <- 0

imputedCEa <- mice(SS3a, predictorMatrix = pred, method = meth, maxit = 20,
                   seed = 9212, print = FALSE)

head(complete(imputedCEa), action = 1, 10) 
#Can check that values have been computed properly; action = 1 (default) shows
#first completed dataset

glmnbCEa <- with(imputedCEa, 
                 glm.nb(Post100RoundsGoal ~ SeqCond*PostHowManySingleCaptionType * PreDBC_Total))
summary(pool(glmnbCEa))
