######################################################
#Lecture 4 Script (including seminar)
#Skewed, categorical and clustered data (Seminar material starts line,
#but pay attention also to new line 25 and 62)
#Anastasia Ejova 27/10/2014
#R101/PSY532
######################################################

#The workspace is in the Lecture 4 folder.

#-----------------------------------------------------
#Skewed data: GENERALIZED LINEAR MODELLING
#-----------------------------------------------------
#Descriptive statistics (as for ANOVA/ANCOVA in Lecture 2)
#--

#Running the analysis
?glm #canonical link functions specified if you follow the family link

#Original Hypothesis 2 test - with normal (i.e., gaussian) random component
original <- glm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType +
                  PreDBC_Total, family = gaussian(link = "identity"), data = SS)
summary(original)
library(car)
Anova(original, type =2)

#Attempt with Poisson random component
glmp <- glm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
              PreDBC_Total, family = poisson(link = "log"), data = SS)
summary(glmp) #Compared to original: reduction in residual deviance but increase in AIC.

#Check for overdispersion of residuals in standard Poisson model
#install.packages("boot")
library(boot)
library(ggplot2)
stdres <- glm.diag(glmp)$rp #vector of standardized residuals from boot package
fitted <- glmp$fitted.values
df <- data.frame(stdres, fitted)

diagplot <- ggplot(aes(x = fitted, y = stdres), data = df) + 
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = 2, linetype = "dashed") +
  geom_hline(yintercept = -2, linetype = "dashed")
diagplot

#Since there is overdispersion (i.e., since many residuals fall outside the lines),
#an overdispersed Poisson model should be better
glmpo <- glm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
              PreDBC_Total, family = quasipoisson (link = "log"), data = SS) 
summary(glmpo) #Indeed, we see a further reductiion in residual deviance, but not possible
                #to calculate AIC.

#Would a negative binomial random component work better?
library(MASS) #for glm.nb function
glmnb <- glm.nb(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                  PreDBC_Total, data = SS)
summary(glmnb) #Yes, a better option than original: lower residual deviance and lower AIC

#For reporting results, need an ANOVA table:
Anova(glmnb, type = 2)
#--

#Diagnostics
#Influential points
plot(cooks.distance(glmnb), type = "h")

  #None are very large, but if they were, you'd need to run the model without the problematic
  #cases, which you can identify using the following command:
  which(cooks.distance(glmnb)>3*mean(cooks.distance(glmnb)))

#Normality of residuals and homogeneity of variance
#Four plots for this purpose are included in the boot package
library(boot)
glmnb.diag <- glm.diag(glmnb) #to compute the information needed for the plot
glm.diag.plots(glmnb, glmnb.diag, iden = T) #draw plots, incorporating information computed at the previous
                                  #line and by the original glm function
  
  ?glm.diag.plots
  #Left top corner plot provides a test of the homogeneity of variance - plot of residuals against
  #predicted values
  #Right top corner plot is a q-q plot of residuals: residuals are normal if they lie on the
  #diagonal line
  #Right bottom corner plot is the same as plot(cooks.distance(glmnb), type = "h") except in that
  #the dashed line is added to suggest which points might be influential (i.e., those above it)

#You can also plot different kinds of information about residuals and their variance using the
#plot function in the base package
plot(glmnb,1)
plot(glmnb,2)
plot(glmnb,3)
plot(glmnb,4)
plot(glmnb,5)
plot(glmnb,6)

#--

#-----------------------------------------------------
#Categorical data: CHI-SQUARE TESTS
#-----------------------------------------------------
#Descriptive statistics - but also provided as part of the output
freqtable <- xtabs(~SeqCond + PostStrategyPresent, data = SS)
freqtable

#Running the test
library(lsr)
associationTest(~SeqCond + PostStrategyPresent, data = SS)

#-----------------------------------------------------
#Categorical data: LOGISTIC REGRESSION
#-----------------------------------------------------
#Select only complete cases for variables we are interested in
CG1$SumPreds <- as.numeric(CG1$PathCat) +
  as.numeric(CG1$PlayFreq) + as.numeric(CG1$GameNoCat) + as.numeric(CG1$Online)
CG1nm<- CG1[!is.na(CG1$SumPreds),]

#Descriptive statistics
summary(CG1nm$PlayFreq) #Categories not even, so will reduce number of categories.
  #Will make a big category "Weekly or more" containing 53 people
  library(plyr)
  ?mapvalues
  CG1nm$PlayFreq <- mapvalues(CG1nm$PlayFreq, from = c("1-2 times per week", "3-4 times per week",
                                                   "Almost every day"), 
                            to = c("Weekly", "Weekly", "Weekly"))
  summary(CG1nm$PlayFreq)
   
summary(CG1nm$GameNoCat)

summary(CG1nm$Online)

summary(CG1nm$PathCat)
#--

#Analysis
#Step 1
model1 <- glm(PathCat ~ PlayFreq, family=binomial(link="logit"), data = CG1nm)
summary(model1)

#To assess improvement in model fit - contribution of PlayFreq as a whole:
anova(model1, test = "Chisq") #"Chisq", "LRT", "Rao", "F" or "Cp"; see recommendations made in
                              #the help file

#Step 2
model2 <- glm(PathCat ~ PlayFreq + GameNoCat, family=binomial(link="logit"), data = CG1nm)
summary(model2)

anova(model1, model2, test = "Chisq") #Marginally significant additional contribution of
                                    #variety of games played

#Step 3
model3 <- glm(PathCat ~ PlayFreq + GameNoCat + Online, family=binomial(link="logit"), 
              data = CG1nm)
summary(model3) 

anova(model1, model2, model3, test = "Chisq") #No significant contribution of whether one of
                                          #the games was played online

#Interpretation of coefficients in final model (model3)
summary(model3) #Logit scale coefficients with Wald z-tests
  #The p-values associated with the Wald tests tell you whether the predictor is significant.
  #The "Estimates" (called "log odds") can be interpreted in the following way:
    #For the first coefficient and for all of them: The probability (note, different to "odds") 
    #of being a pathological gambler if playing once per month is .25 * the Estimate more (or less
    #depending on the sign) than the probability of being a pathological gambler if playing once
    #total. So, for the first coefficient (-0.76936), the probability of being a pathological
    #gambler if playing less than once a month is surprisingly 19% LESS (.25*(-.76) = -.19) than
    #at the reference point for play frequency (i.e., when the person has played only once in total
    #in the last 12 months). This difference is not significant, however, and is often an
    #overestimate (Baguley p678).

exp(model3$coefficients) #Coefficients expressed as odds ratios

  #0.1274872: The odds of pathology when all the coefficients are at reference level.

  #0.4633090: The odds of gambling pathology when play frequency is less than once a month 
  #as a proportion of the odds of gambling pathology when play frequency is once total.  
  #So, the probability of being a pathological 
  #gambler when playing less than once a month is surprisingly LESS than the probability of being a 
  #pathological gambler when playing only once a year, but the diference is not significant
  #as the Wald z tests indicate. 
  #Note, we are assuming that all the other predictors are fixed at some level, which means
  #that we are assuming no relationship between predictors - more on this later.

  #1.1166979: The odds of gambling pathology when play frequency is once a month 
  #as a proportion of the odds of gambling pathology when play frequency is once total.
  #So, the probability of being a pathological 
  #gambler when playing once a month is higher than the probability of being a 
  #pathological gambler when playing only once a year, but the diference is again not significant. 
  #Note, we are assuming that all the other predictors are fixed at some level.
  
  #... Moving to the coefficients significantly different from 0 (on the logit scale)

  # 2.9338080 (marginally significant: p = .07): The odds of gambling pathology when play frequency 
  #is weekly
  #as a proportion of the odds of gambling pathology when play frequency is once total.
  #So, the probability of being a pathological 
  #gambler when playing weekly or more often is (significantly) higher than the probability of being a 
  #pathological gambler when playing only once a year. 
  #Note, we are assuming that all the other predictors are fixed at some level.
  
  #2.0642397 (marginally significant: p = .07): Gambling pathology has more than twice the
  #odds of existing when two or more games are played, compared to when only one game style
  #is played. Again, all other predictors are assumed to be unchanging.

  #For continuous predictors: The intercept is the odds ratio of pathology when all coeffs
  #are at 0, and the coefficient is the amount of change in the odds of the outcome variable being
  #equal to one value (e.g., gambling pathology) as opposed to the alternative (e.g., no gambling
  #pathology) for each unit change in the outcome variable. See Baguley text p679 for suggestion on
  #how to compare the coefficients of continuous and categorical predictors that are in the same
  #model. For continuous predictors, the
  #"Estimates" (log odds) obtained through summary(model3) also need to be interpreted
  #"per unit change" rather than "with change from reference point to the listed category".

  #Some helpful plots
  plot(CG1nm$PlayFreq, model3$fitted.values,
     ylab = "Probability of gambling pathology", xlab = "Play frequency")
  
  plot(CG1nm$GameNoCat, model3$fitted.values,
     ylab = "Probability of gambling pathology", xlab = "Game variety")
#--

#Diagnostics
#Detecting collinearity (i.e., relationships between predictors)
library(lsr)
associationTest(~PlayFreq + GameNoCat, data = CG1nm)
associationTest(~PlayFreq + Online, data = CG1nm)
associationTest(~GameNoCat + Online, data = CG1nm)
#It is evident that all the predictors are correlated - which makes sense! Afterall,
#if you have played only once, you are less likely to have played more than one game!

#An option would be to include just two predictors (Online being one of them), and 
#test for the effect of their interaction as well.

model4 <- glm(PathCat ~ GameNoCat * Online, family=binomial(link="logit"), data = CG1nm)
summary(model4) #Still no significant effect of Online, and the interaction effect is not significant.

library(car)
Anova(model4, type = 2)

#Detecting outliers/influential points
boxplot(cooks.distance(model3)) #The Cook's distances are all very small in this case
#However, imagine that they were bigger and I wanted to run the model without the cases with the
#two highest Cook's distances.

row.names(CG1nm[cooks.distance(model3) >=.06,]) #the high Cook's distances belong to participants
                                                #15 and 42
model3a <- glm(PathCat ~ PlayFreq + GameNoCat + Online, family=binomial(link="logit"), 
               data = CG1nm, subset = -c(15,42))
summary(model3a) #Same overall pattern, with GameNoCat actually becoming significant
anova(model3a, test = "Chi")

#Checking for sparse data - how much are you predicting over the null model?
threshold <- mean(model3$fitted.values) #Threshold for deciding whether a fitted value
                                        #should go to 1 or 0. This is a good default option.
                                        #However, there is a whole literature on this - e.g.,
                                        #http://www.talkstats.com/showthread.php/52814-ROC-cutoff-value?p=147952&viewfull=1#post147952

basic <- table(model3$y, model3$fitted.values > threshold)
basic
prop.table(basic) #approx 60% of cases classified correctly: 43.6 + 15.7

#Is this a lot? In other words, what is the null model's probability of being correct?
model0 <- glm(PathCat ~ 1, family=binomial(link="logit"), data = CG1nm)
head(model0$fitted.values,20) #See the first 20 fitted values for null model
#All equal to 0.213198. Let's assume that we can classify everyone then as a "no risk".
#In this situation, we would be correct in 155 out of 197 (79%!) of cases.

#-----------------------------------------------------
#Clustered data: ZERO-INFLATED MODELLING
#-----------------------------------------------------
#This section and the next follow up on Assignment 2, where we found that the assumptions for the ANCOVA we were 
#conducting were not met. What we noticed is that prior beliefs (PreDBC_IoC) interacted with the main predictor 
#(SeqCond) in affecting the
#the illusion of control (PostNaturalIoC). We decided to split prior beliefs into two categories and treat them as a 
#categorical variable in an ANOVA with two predictors: SeqCond AND categories of PreDBC_IOC. The assumptions for that
#ANOVA weren't met either, so, with zero-inflated modelling and multilevel modelling we trial alternative approaches.

#Example 1: SS Hypothesis 1
#Descriptive statistics
#Plot to decide whether a zero-inflated model is appropriate
qplot(SS$PostNaturalIoC, geom = "histogram", binwidth = 1)
stem(SS$PostNaturalIoC) #Stem and leaf plot could also be useful

#Preparation
#Outcome variable must be in integer form, with zeroes evident
SS$PostNaturalIoCInt <- as.integer(SS$PostNaturalIoC)

#Running the analysis
library(pscl)
?zeroinfl
modelziNat1 <- zeroinfl(PostNaturalIoCInt ~ SeqCond*PreDBC_Total, data = SS) #zero-infl poisson regr.
# y ~ x1 + x2 equivalent to y ~ x1 + x2 | x1 + x2 where the second set of predictors
# relates to the point mass (zeros) component of the model (note, as per help file,
# there are two sources of zeros under this modelling approach)

summary(modelziNat1)
#We see a significant effect only of prior beliefs in influencing the illuision of control in the "non-zero" category.

modelziNat2 <- zeroinfl(PostNaturalIoCInt ~ SeqCond*PreDBC_Total, dist = c("negbin"), data = SS)
#zero-infl negative binomial regression
summary(modelziNat2)

#As can be seen from the summaries, the models have the same Log Likelihood. So I will
#report the slightly simpler model (Poisson with 16 Df). Otherwise, I would have chosen
#the model with the Log Likelihood closer to zero.

#--

#Example 2: SS Hypothesis 1 - "It was all chance"
#Descriptive statistics/preparation
#Plot to decide whether a zero-inflated model is appropriate
qplot(SS$PostListCHANCE, geom = "histogram", binwidth = 1)
stem(SS$PostListCHANCE) #Stem and leaf plot could also be useful

#Break the variable into categories, defining the zero point (in this case, the zero is the response
#of value 10, indicating complete agreement that "It was all chance")
SS$ChanceCat[SS$PostListCHANCE == 10] <- 0
SS$ChanceCat[SS$PostListCHANCE >= 8 & SS$PostListCHANCE <= 9] <- 1
SS$ChanceCat[SS$PostListCHANCE >= 5 & SS$PostListCHANCE <= 7] <- 2
SS$ChanceCat[SS$PostListCHANCE >= 0 & SS$PostListCHANCE <= 4] <- 3

#install.packages("pscl")
library(pscl)
?zeroinfl
modelziChance1 <- zeroinfl(ChanceCat ~ SeqCond + PreDBC_Total, data = SS) #zero-infl poisson regr.
# y ~ x1 + x2 equivalent to y ~ x1 + x2 | x1 + x2 where the second set of predictors
# relates to the point mass (zeros) component of the model (note, as per help file,
# there are two sources of zeros under this modelling approach)
summary(modelziChance1)

modelziChance2 <- zeroinfl(ChanceCat ~ SeqCond + PreDBC_Total, dist = c("negbin"), data = SS)
#zero-infl negative binomial regression
summary(modelziChance2)

#-----------------------------------------------------
#Clustered data: MULTILEVEL MODELLING
#-----------------------------------------------------
SS$CatPreDBC_IOC <- cut(x = SS$PreDBC_IOC, breaks = 2, labels = c("Low", "High"))
summary(SS$CatPreDBC_IOC)

#Descriptive statistics
library(psych)
describeBy(x = SS$PostNaturalIoC, group = list(SS$SeqCond, SS$CatPreDBC_IOC))

describeBy(x = SS$PreDBC_IOC, group = SS$CatPreDBC_IOC)

#--

#Running the analysis
library(lme4)

#Since the PostNaturalIoC means in the baseline Descending group could be different
#across belief categories (1.74 in low, 2.74 in high), we fit a varying slope, varying 
#intercept model
Hyp1mm <- lmer(PostNaturalIoC ~ SeqCond +  (1 + SeqCond | CatPreDBC_IOC), data = SS)
print(Hyp1mm)
#For guidance on model notation, see: http://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet and
#Gelman text Ch13

#Interpretation:
#estimated standard deviation of the prior belief intercepts: .70
#estimated standard deviation of the prior-belief slopes: 1.83
#estimated correlation between intercepts and slopes is... very high!
  #So we should be fitting a simpler (varying intercept) model.
  #Look here for discussion: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q3/016780.html
  #Fitting the simpler model...

Hyp1mm <- lmer(PostNaturalIoC ~ SeqCond +  (1 | CatPreDBC_IOC), data = SS)
print(Hyp1mm) #Interpretation still as in lines 295-6 above

coef(Hyp1mm) #Estimated differences between named condition and Desc condition for each prior belief level
#Do not differ across belief levels in this case, since no slope variability as a function of the 
#random variable is assumed 

ranef(Hyp1mm) #Estimated amount to be added to the differences betw named condition and Desc condition for each
#belief level if we are to try to predict the group differences from this model. So, the predicted difference
#between the Asc and U-shaped conditiions among people with low prior illusion of control beliefs, based on this
#model is:
#0.4859199 (from coef(Hyp1mm)) PLUS
#-0.8936404  (from ranef(Hyp1mm)) = -0.4077205

#Hypothesis-testing approaches
#(1) Confidence intervals around group differences on the fixed factor. The differences are 
#significant if the confidence interval does not include 0. Analogous to summary.lm(modelname).
#A significant difference between the Ascending and Descending conditions is evident here.
confint(Hyp1mm, method="Wald")

#(2)Wald test of a fixed factor's significance when there is more than one fixed factor. For the detailed maths, see: 
#http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/6749
Anova(Hyp1mm, type = 2)

#(3)Model comparison
Hyp1mm0 <- lmer(PostNaturalIoC ~  1| CatPreDBC_IOC, data = SS) #Model according to which only prior beliefs matter
print(Hyp1mm0)
anova(Hyp1mm0, Hyp1mm, test = "Chi") #Comparison of Hyp1mm0 with the full model.

#--

#Diagnostics - We'll check these for Hyp1mm, but, really, these
#should be checked for the best model of the data you decide on (e.g., Hyp1mm0 - the simpler
#model that doesn't perform worse)
#Normality of residuals
hist(residuals(Hyp1mm), breaks = 20) #Quite okay
shapiro.test(residuals(Hyp1mm)) #Significant, but this is often the case

#Homogeniety of variance
plot(Hyp1mm) #quite okay - even spread around 0 line
