######################################################
#Lecture 2 Script (with code for Seminar starting at line 121 - ANOVA)
#Anastasia Ejova 06/10/2014
#R101/PSY532
######################################################

#-----------------------------------------------------
#T-TESTS - quick demonstration
#-----------------------------------------------------
library(lsr)
independentSamplesTTest(
  formula = PostSupIoC ~ ThinkAloud, # formula specifying outcome and group variables
  data = SF,
  var.equal = FALSE, one.sided = FALSE, conf.level = .95) #These last three are defaults

#Other examples
oneSampleTTest(SS$Age, mu=18)

pairedSamplesTTest(
  formula = ~ PlayerChange_First24 + PlayerChange_Last24, # one-sided formula listing the two variables
  data = SS2)

#-----------------------------------------------------
#REGRESSION
#-----------------------------------------------------
#Descriptive statistics
#Graph
  #install.packages(ggplot2)
  library(ggplot2)
  qplot(data = SF, x = WinFreqPerc, y = PostSupIoC, geom = c("point", "smooth"), 
        method = "lm", se = FALSE)
  
  log10(50)

  qplot(data = SF, x = log10(WinFreqPerc), y = PostSupIoC, geom = c("point", "smooth"), 
      method = "lm", se = FALSE)

  SF$LogWinFreqPerc <- log10(SF$WinFreqPerc)

#Correlation table
  SF$PreSoccerInterest <- as.integer(SF$PreSoccerInterest)
  
  SFsub <- subset(SF, select = c(LogWinFreqPerc,
                                 PreSoccerInterest,
                                 PreDBC_Sup:PreDBC_Total,
                                 PostSupIoC,
                                 EndCredits))

  correlate(SFsub, 
          corr.method = "spearman", #default is pearson
          test = TRUE, #default is FALSE, but we wish to test significance
          p.adjust.method = "none") #lsr package
#--

#Running the analysis - hierarchical regression
#Ensure that all steps of the hierarchical regression will be conducted on the 
#same participants
#i.e., remove those who are missing a value on any of the variables in the model
  SF$SumPreds <- SF$PreSoccerInterest + SF$PreDBC_Sup + SF$LogWinFreqPerc
  SFhr<- SF[!is.na(SF$SumPreds),]

#Fit the original model
  model0 <- lm(PostSupIoC ~ PreDBC_Sup + PreSoccerInterest, data = SFhr)
  summary(model0)

#Fit a model with one or more additional predictors with a different name
  model1 <- lm(PostSupIoC ~ PreDBC_Sup + PreSoccerInterest + LogWinFreqPerc, data = SFhr)
  summary(model1)

#Conduct ANOVA to determine whether the second model is significantly better than the first
  anova (model0, model1)
anova(model1, model0)
  
#Can then add another predictor and include three arguments in the anova function

#Additional useful functions
  standardCoefs(model1) #standardCoefs is a function in the lsr package
  AIC (model0, model1)
#--

#Checking assumptions
#Normality of residuals
  hist(residuals(model1), breaks = 20) #histogram with multiple breaks
  plot(model1, which = 2) #Q-Q plot - residuals are normal if points lie on the diagonal
  shapiro.test(residuals(model1))

#Constant variance of residuals
  #Lack of influential points
  plot(model1, which = 4)
  plot(model1, which = 5)

    #Running regression without influential points:
    summary(lm(PostSupIoC ~ PreDBC_Sup + PreSoccerInterest + LogWinFreqPerc, data = SFhr,
            subset = -c(5, 38, 74))) # in this case, the findings remain the same, which is good
  
  #Same variance/standard deviation across all parts of the fitted line (i.e., for all fitted values)
  plot(model1, which = 3)
  #install.packages(car)
  library(car) #for each value of the predictor, 
  #with tests of whether the visible curvature is significant
  ncvTest(model1)

    #If there is significant heteroscedasticity (i.e., if ncvTest result is signficant)
    #install.packages("lmtest")
    library(lmtest) #also, car package must be loaded
    coeftest(model1, vcov= hccm)
  
#Linearity of the relationship
  #Plot of fitted values (predicted by model) against observed values
  fitted <- fitted.values(model1)
  plot( x = fitted,
        y = SFhr$PostSupIoC,
        xlab = "Fitted Values",
        ylab = "Observed Values")
  
  #Plots of residuals against fitted values and values of the predictors
  plot(model1, which = 1) #for whole model
  residualPlots(model1) #function from car package
  
#-----------------------------------------------------
#ANOVA: INDEPENDENT SAMPLES
#-----------------------------------------------------
#Example 1: Success-frequency as a single categorical predictor
  library(psych)
  describeBy(SF$PostSupIoC, group = SF$WinFreq) #note that 1/16 mean is 2.11
  
  anova_SFSingle <- aov(PostSupIoC ~ WinFreq, data = SFhr) 
    #Note the use of WinFreq - a categorical variable rather than the numeric
    #WinFreqPerc or LogWinFreqPerc.
  
  summary(anova_SFSingle) 
    #Tells us the result of F-test.
  
  summary.lm(anova_SFSingle)
    #Tells us the results of the regression component (without F-test).
    #The coding used in the regression is the default in R:
        options("contrasts") #Treatment contrasts
    #For WinFreq this means:
        contrasts(SFhr$WinFreq) #Exactly what we have on the lecture slide (slide 18)
  
    #Based on what we know from our describeBy table, we can see that:
      #b0 = Mean of 1/16 win frequency group
      #b1 = Mean of 1/8 group minus mean of 1/16 group
      #b2 = Mean of 1/4 group minus mean of 1/16 group
      #etc.
    
    #The results suggest that b4 (difference between 1/16 group and 1/2 group)
    #is significantly different from 0.
  
  #The F-test as a model comparison
    anova_null <- aov(PostSupIoC ~ 1, data = SFhr) #intercept only model
    summary(anova_null)
    anova_full1 <- aov(PostSupIoC ~ WinFreq, data = SFhr)
    anova(anova_null, anova_full1) #F-value corresponds to F-value for anova_SFSingle
      #Back to lecture - other possible contrasts
  
  #Other options for contrasts - we will see the use for these more in Example 3
    #(1) Use existing contrast matrices: contr.helmert and contr.sum
    options(contrasts = c("contr.helmert", "contr.poly")) #change default contrasts
                                                          #for all analyses, then re-run
                                                          #ANOVA
      ?contr.helmert 
      #According to help file: Helmert contrasts contrast the second level with the 
      #first, the third with the average of the first two, and so on
      contrasts(SFhr$WinFreq) 
      anova_SFSingle <- aov(PostSupIoC ~ WinFreq, data = SFhr) 
      summary.lm(anova_SFSingle)
    
    options(contrasts = c("contr.sum", "contr.poly"))
      contrasts(SFhr$WinFreq) 
      anova_SFSingle <- aov(PostSupIoC ~ WinFreq, data = SFhr) 
      summary.lm(anova_SFSingle)
        #Can see that intercept is close to the grand mean - in fact, it is the weighted
        #grand mean: sum of (group sample size * group mean) for all groups
        #divided by total sample size
          mean(SFhr$PostSupIoC) #2.82

          #Based on what we know from our describeBy table, we can see that:
          #b0 = (weighted) grand mean
          #b1 = Mean of 1/16 group minus weighted grand mean
          #b2 = Mean of 1/8 group minus weighted grand mean
          #...
          #b4 = Mean of 1/3 group minus weighted grand mean
    
    #(2) Change the reference group for whatever contrast is being used
    options(contrasts = c("contr.treatment", "contr.poly"))
    SFhr$WinFreq <- relevel(SFhr$WinFreq, ref = "1/2" )
    
    #(3) Manually set orthogonal contrasts - see rules for doing this in lecture,
    #and an example in Example 3
#--

#Example 2: Success-frequency as a predictor in an ANCOVA

anova_SFHyp <- aov(PostSupIoC ~ PreDBC_Sup + PreSoccerInterest + WinFreq, data = SFhr)
  
summary(anova_SFHyp)
contrasts(SFhr$WinFreq)  
summary.lm(anova_SFHyp)
library(lsr)
etaSquared(anova_SFHyp)
#--

#Example 3 (main example, with all steps): Hypothesis 2 from SS data
#Descriptive statistics
  interaction.plot(x.factor = SS$SeqCond, trace.factor = SS$PostHowManySingleCaptionType, 
                   response = SS$PostHowManySingleWins)
  describeBy(SS$PostHowManySingleWins, group = list(SS$SeqCond, 
                                                    SS$PostHowManySingleCaptionType))
  
  cor.test(SS$PreDBC_Sup, SS$PostHowManySingleWins, method = "spearman",
           exact = TRUE) #asking for test of correlation's significance
      #Can see moderate significant correlation
      #Note: this function is part of the base package

#Running the analysis
  #F-test - Type II Sums of Squares
  anova_SSHyp2 <- lm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                       PreDBC_Total, data = SS)
  library(car)
  Anova(anova_SSHyp2, type = 2) #Run as Type II sum of squares
  etaSquared(anova_SSHyp2)
  
  #Planned comparisons/contrasts
  options(contrasts = c("contr.treatment", "contr.poly")) #reset all contrasts to defaults
  contrasts(SS$SeqCond)
  #SS$SeqCond <- relevel(SS$SeqCond, ref = "Descending" ) #Since we're interested in 
                                          #comparing the Descending condition to 
                                          #all the others,
                                          #I would have needed to do this, had the
                                          #Descending condition not been the reference
                                          #category already.
  
  #anova_SSHyp2 <- lm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                       #PreDBC_Total, data = SS) 
  #Anova(anova_SSHyp2, type = 2) #Would have also needed to re-run the ANOVA.
  summary.lm(anova_SSHyp2) #See lecture "Reporting the analysis" for the interpretation.
  
    #What if I were interested in the two conditions with many early wins relative to
    #all others? An illustration of manually set orthogonal contrasts:
    levels(SS$SeqCond)
    contrasts(SS$SeqCond) <- cbind(c(-1, -1, 1, 1)) #Wish to compare the Desc and
                                                    #U-shaped conditions to the other
                                                    #two.
    contrasts(SS$SeqCond) #R has created an orthogonal contrasts matrix with 
                          #the above matrix as the first column
    anova_SSHyp2 <- lm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                     PreDBC_Total, data = SS)
    summary.lm(anova_SSHyp2) #Given the contrasts matrix, we are interested just in 
                              #whether the first comparison (SeqCond1) is
                              #significantly different from 0. And it is not! For the
                              #interaction, we are also interested only in the first
                              #coefficient, where the difference we are interested in
                              #does not vary significantly across CaptionType.
       
  #Estimated marginal means - useful for graphs etc. after ANCOVA
  #install.packages("effects")
  library(effects)
  eff <- effect(term = "SeqCond*PostHowManySingleCaptionType", mod = anova_SSHyp2)
  eff
  summary(eff) #provides confidence intervals around estimates

#Checking assumptions
  #Normality of residuals
  hist(residuals(anova_SSHyp2))
  
  shapiro.test(residuals(anova_SSHyp2))
               
  #Homogeneity of variance
  leveneTest(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType, data = SS) #in car package
    #violated in this case
    #but cannot run kruskal.test() and oneway.test() with this complex design

  #Homogeniety of regression (ANCOVA)
  HRS <- aov(PreDBC_Total ~ SeqCond*PostHowManySingleCaptionType*PreDBC_Total, data = SS)
  Anova(HRS, type = 2)
  
  #Independence between covariate and predictor(s)
  Ind <- aov(PreDBC_Total ~ SeqCond*PostHowManySingleCaptionType, data = SS)
  Anova(Ind, type = 2)

#-----------------------------------------------------
#REPEATED MEASURES ANOVA
#-----------------------------------------------------
#Example 1: Hypothesis 1a from Lecture 1 (player profile changes and dir entropy over
#time and across success-slopes)
  #Restructuring data into long form
  #install.packages(reshape)
  library(reshape)
  SS2subset <- SS2[,c("PNo", "SeqCond", "PreDBC_Total", "PlayerChange_First24",
                   "PlayerChange_Last24", "DirEntropy_First24",
                   "DirEntropy_Last24")]
  
  SS2long <- reshape(SS2subset, idvar = c("SeqCond", "PNo", "PreDBC_Total"),
                        varying = rbind(c("PlayerChange_First24",  
                                          "PlayerChange_Last24"),
                                        c("DirEntropy_First24", 
                                          "DirEntropy_Last24")),
                        v.names = c("PlayerChanges", "DirEntropy"),
                        timevar = "Time",
                        times = c(1,2),
                        direction = "long")
  
  #Descriptive statistics
  #install.packages(sciplot)
  library(sciplot)
  lineplot.CI(x.factor = SS2long$Time, response = SS2long$PlayerChanges, group = SS2long$SeqCond)
  
  library(psych)
  describeBy(SS2long$PlayerChanges, group = list(SS2long$SeqCond, SS2long$Time))
  
  cor.test(SS2$PlayerChange_First24, SS2$PreDBC_Total, method = "spearman", exact = TRUE)
  cor.test(SS2$PlayerChange_Last24, SS2$PreDBC_Total, method = "spearman", exact = TRUE)
  cor.test(SS2$DirEntropy_First24, SS2$PreDBC_Total, method = "spearman", exact = TRUE)
  cor.test(SS2$DirEntropy_Last24, SS2$PreDBC_Total, method = "spearman", exact = TRUE)
    #Some significant correlations here, so PreDBC_Total should be included as a covariate
  
  #Running the analysis
    class(SS2long$Time)
    SS2long$Time <- as.factor(SS2long$Time)
    SS2long$PNo <- as.factor(SS2long$PNo)
    
    #Centre the covariate, as recommended by Baguley p651
    SS2long$PreDBC_TotalCentred <- SS2long$PreDBC_Total -  mean(SS2$PreDBC_Total)
    
    #Analysis of player changes
    #install.packages(nlme)
    library(nlme)
    lme_SSHyp1aPlayers <- lme(fixed = PlayerChanges ~ PreDBC_TotalCentred + SeqCond*Time, 
              random = ~1|PNo, data=SS2long) #see Baguley p660
    anova(lme_SSHyp1aPlayers, type = "marginal") #Type 2 sums of squares again

    #Interpretation: The interaction effect is significant, which is what we expected to see.

    #Analysis of kick direction entropy
    lme_SSHyp1aDir <- lme(fixed = DirEntropy ~ PreDBC_TotalCentred + SeqCond*Time, 
                              random = ~1|PNo, data=SS2long) #see Baguley p660
    anova(lme_SSHyp1aDir, type = "marginal")
    #All effects significant, including the interaction that was significant in the 
    #player changes analysis.
  
  #Checking assumptions
  #In this case, we will not examine sphercity corrections because there are only two levels of
  #the repeated-measures predictor variable. To check the normality of residuals assumption:
  
  residuals(lme_SSHyp1aPlayers)
  hist(as.numeric(residuals(lme_SSHyp1aPlayers)), breaks = 20)
  shapiro.test(residuals(lme_SSHyp1aPlayers)) #Shapiro test significant, and the residuals are somewhat skewed
    #but they are acceptable
  
  residuals(lme_SSHyp1aDir)
  hist(as.numeric(residuals(lme_SSHyp1aDir)), breaks = 20)
  shapiro.test(residuals(lme_SSHyp1aDir))
  
#--

#Example 2 
#Uses SF2 database, which contains three additional variables (and no LogWinFreqPerc):
  #SF2$PlayerChangesPerc_First30 - % of first 30 rounds on which there were player profile changes
  #SF2$PlayerChangesPerc_Mid40 - % of the next 40 rounds on which there were player profile changes
  #SF2$PlayerChangesPerc_Last30 - % of last 30 rounds on which there were player profile changes

  #Restructuring data into long form
  SF2subset <- SF2[,c("PNo", "WinFreq", "PreDBC_Total", "PlayerChangesPerc_First30",
                    "PlayerChangesPerc_Mid40",
                    "PlayerChangesPerc_Last30")]

  SF2long <- reshape(SF2subset, idvar = c("WinFreq", "PNo", "PreDBC_Total"),
                   varying = rbind(c("PlayerChangesPerc_First30",  
                                     "PlayerChangesPerc_Mid40",
                                     "PlayerChangesPerc_Last30")),
                   v.names = c("PlayerChangesPerc"),
                   timevar = "Time",
                   times = c(1,2,3),
                   direction = "long")

  #Descriptive statistics
  lineplot.CI(x.factor = SF2long$Time, response = SF2long$PlayerChangesPerc, 
              group = SF2long$WinFreq)

  describeBy(SF2long$PlayerChangesPerc, 
             group = list(SF2long$WinFreq, SF2long$Time))
             
  cor.test(SF2$PlayerChangesPerc_First30, SF2$PreDBC_Total, method = "spearman", exact = TRUE)
  cor.test(SF2$PlayerChangesPerc_Mid40, SF2$PreDBC_Total, method = "spearman", exact = TRUE)
  cor.test(SF2$PlayerChangesPerc_Last30, SF2$PreDBC_Total, method = "spearman", exact = TRUE)
    #No significant correlations here, so PreDBC_Total will not be included as a covariate
  
  #Running the analysis and checking the sphercity assumption
  #install.packages("ez")
  #Can use this package if there are no covariates. Its covariance analyses are still
  #under development. It is particularly useful for providing sphercity corrections.
  library(ez)
  SF2long$Time <- as.factor(SF2long$Time)
  SF2long$PNo <- as.factor(SF2long$PNo)
  
  anova_SFPlayerChanges <- ezANOVA(SF2long, dv = .(PlayerChangesPerc), wid = .(PNo), 
                                    within = .(Time), 
                                    between = .(WinFreq), 
                                    type = 2, #Type II sums of squares
                                    detailed = TRUE) #print out df and sums of squares
  #.() = is equivalent to list(), as we have been doing for describeBy. 
  #Each .() can contain numerous comma-separated elements.
  
  anova_SFPlayerChanges
    #Given the epsilon values, would report p[GG]. The ges is an effect size
    #that can be interpreted in the same way as partial eta squared (see LSR p.502).
    #Interpretation: Significant effect of time observed, 
    #but not the expected significant interaction.
  
  #Checking the normality of residuals assumption and further examining the effects of time
    lmeCheck <- lme(fixed = PlayerChangesPerc ~ WinFreq*Time, 
                random = ~1|PNo, data=SF2long)
    anova(lmeCheck, type = "marginal")
    hist(as.numeric(residuals(lmeCheck)),breaks = 20)
    shapiro.test(residuals(lmeCheck))
  
    #To further investigate the effects of time
    summary(lmeCheck)  #No significant differences between Time 1 vs 2, or Time 1 vs 3
    
    #What about Time 2 vs 3
    SF2long$Time <- relevel(SF2long$Time, ref = "3") 
    lmeCheck2 <- lme(fixed = PlayerChangesPerc ~ WinFreq*Time, 
              random = ~1|PNo/WinFreq, data=SF2long)
    summary(lmeCheck2) #No significant differences between Time 2 and 3 either
   
    #Interpretation: There is some decrease in the number of player changes made over time, but 
    #it is not visible in pairwise comparisons.  
  
#-----------------------------------------------------
#BOOTSTRAPPING - quick demonstration
#-----------------------------------------------------
#install.packages("simpleboot")
library(simpleboot)

library(car)
#http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Bootstrapping.pdf
  
#install.packages(Hmisc)
library(Hmisc)

#Independent samples t-test
  v1 <- subset(SS, PostStrategyPresent == "No", select = c("PostNaturalIoC"), drop = T)
  v2 <- subset(SS, PostStrategyPresent == "Yes", select = c("PostNaturalIoC"), drop = T)
  tb <- two.boot(sample1 = v1, sample2 = v2, #simpleboot function; v1 and v2 must be vectors
                 FUN = mean, 
                 R = 1000) #number of bootstrap replicates
  boot.ci(tb) #BCa usually recommended as the statistic to report

#--

#Mulitple regression
  regr1 <- lm(PostSupIoC ~ PreDBC_Sup + PreSoccerInterest + LogWinFreqPerc, data = SFhr)
  #This is model1 from line 67
  summary(regr1)

  regr1b <- Boot(object = regr1, #car functiion
               f=coef, labels=names(coef(regr1)),
               R=999, method=c("case"))
  summary(regr1b)
  confint(regr1b)

#--

#ANOVA: factorial or single-predictor
  anova1 <- lm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                 PreDBC_Total, data = SS)
  #anova_SSHyp2 from line 216
  summary(anova1)

  anova1b <- Boot(object = anova1, 
                f=coef, labels=names(coef(anova1)),
                R=999, method=c("case"))
  summary(anova1b)
  confint(anova1b)

#--

#Repeated-measures ANOVA with single predictor termed 'time'
  SF2long$Time <- as.integer(SF2long$Time)
  SF2long$PNo <- as.integer(SF2long$PNo)
  anovaRMBoot <- rm.boot(time = SF2long$Time, y = SF2long$PlayerChangesPerc,
                    subset = SF2long$SeqCond == "Ascending",
                    id=SF2long$PNo, 
                    bootstrap.type='x fixed',
                    nk=0, B=500)
  plot(anovaRMBoot)
  quantile(anovaRMBoot$Coef[,2], c(.05, .95))



  
  