library(nlme)
library(lattice)

# This file demonstrates how to include an interaction.  A good understanding of centering and effect
# coding from regression will help.  My focus will be on the nlme syntax. 

################################################################################
# Load data and set up for analysis using effect coding (contr.sum) and centered
# predictors (by subtracting the mean).
################################################################################

myd<-read.csv("Young E1 DelayDiscounting.dat")
myd$Cluster<-as.factor(myd$Cluster)
contrasts(myd$Sex)=contr.sum(2)
contrasts(myd$Cluster)=contr.sum(3)
myd$c.VGExp<-myd$VGExp-mean(myd$VGExperience, na.rm=T)  # create centered experience variable

# Run a model with the full interaction between Sex and (centered) VG Experience. 
# For starting values we need 4: male-intercept, female-intercept, male-VG slope, female-VG slope.
# Note that by using effect coding, the interpretation is different, but thinking of it this way 
# allows you to count up the necessary number of starting values. 
myr<-nlme(Indifference.point~1000/(1+exp(logk)*Delay), fixed=logk~Sex*c.VGExp, random=logk~1|Subject, data=subset(myd, !is.na(c.VGExp)), start=c(logk=-6, 0, 0, 0))
summary(myr)

# As in regression based on centered and effect-coded variables, the intercept is the grand mean k.  The regression
# equation is -6.75 + 1.39*Sex + .27*c.VGExp - .069*Sex*c.VGExp.  To derive predictions, Sex is -1 for males 
# and +1 for females (see the contrasts).  

# Each p-value evaluates each term: Intercept > 0, Sex dummy > 0 (main effect of sex), VG Exp > 0 
# (main effect of VG Exp), Sex1:c.VGExp > 0 (interaction)

################################################################################
# Multiple comparisons (post hocs and planned comparisons).
################################################################################
# The multcomp library can be used for specific tests or to obtain predicted k values with SE for graphing purposes. 
# The contrast.matrix allows you to specify specific combinations of predictors for specific tests. 
# Each element in the "c" list (see below) corresponds to one of the coefficients in the model "summary" table. 
# For example, c(1,-1,0,0) means that you are using the equation above (the -6.75... one) with the 
# intercept (always present here) plus -1 for Sex, 0 for c.VGExp, and 0 for sex*c.VGExp.  This produces the predicted value when
# Sex is -1 (male) and the VG Experience is average (c.VGExp is 0).  If you need to compute comparisons,
# you subtract the two vectors.  For example, to compare male VERSUS female when VG Experience is 5:
# c(1,-1, 5-11.59, -1*(5-11.59)) - c(1, 1, 5-11.59, +1*(5-11.59)) => c(0, -2, 0, 13.18). 
# *** Note the need to 'center' the VGExp value of 5 by subtracting the 11.59 mean from 5.
# See statistics texts on contrast coding. 
library(multcomp)
contrast.matrix <- rbind(
  'Male at ave VG Exp' = c(1, -1, 0, 0),          # 11.59 is the mean VG Experience score
  'Female at ave VG Exp' = c(1, 1, 0, 0), 
  'Male at VG Exp of 5' = c(1, -1, 5-11.59, -1*(5-11.59)), 
  'Female at VG Exp of 5' = c(1, 1, 5-11.59, +1*(5-11.59)),
  'Female vs. Male at VG Exp = 5' = c(0, -2, 0, 13.18)
)
comps <- glht(myr, contrast.matrix)
summary(comps)



# Now things get complicated: interactions between categorical variables. Six parameters are estimated.
# One of the challenges is counting the number of parameters for the starting values.  When they're
# categorical, it's easy: 2(Sex) x 3(Cluster) = 6. 
myr<-nlme(Indifference.point~1000/(1+exp(logk)*Delay), fixed=logk~Sex*Cluster, random=logk~1|Subject, data=subset(myd, !is.na(c.VGExp)&!is.na(Cluster)), start=c(logk=-6, 0, 0, 0, 0, 0))
summary(myr)

# It would be great to get an ANOVA table for this situation due to the three-level cluster variable,
# but the anova command uses the wrong sum-of-squares (sequential).  Easy solution:
anova(myr, type="marginal")