library(languageR)
################ Reads in Data #############################
d<-read.table("fakedata.txt", header=TRUE)
#renames factors
d$IV1<-ifelse(d$IV1==1, "silent", "noise")
d$IV2<-ifelse(d$IV2==2, "word", ifelse(d$IV2==3, "legal", "illegal")) 

d$NoiseCond<-as.factor(d$IV1)
d$WordCond<-as.factor(d$IV2)
d$WordCond<-as.factor(d$WordCond)
d$NoiseCond<-as.factor(d$NoiseCond)
d$Freq<-as.numeric(d$Freq)

##below lets you look at condition means for fake data set
attach(d)
RTmeans<-aggregate(RT, list(WordCond), FUN=mean)
RTmeansWN<-aggregate(RT, list(WordCond, NoiseCond),FUN= mean)


########### Sets up Treatment Coding ############################

d$WordCond.Treatment<-d$WordCond


# R automatically assigns levels alphabetically, this isn't always what you'll want, so you can reassign the order of the levels as shown below...


d$WordCond.Treatment<-factor(d$WordCond.Treatment, levels=c("word","legal","illegal")) #reorders levels to put "word" in baseline position (1st in list)

# R's default is to set coding scheme to Treatment, so here you don't need to do anything else now that the levels are ordered appropriately.

#More generally, if you just want to specify which level is the baseline you can do the following: 

 #contrasts(d$WordCond.Treatment)<-contr.treatment(3, base=3)
 
#This says set the contrasts to treatment coding with 3 levels, with the 3rd level being the base condition


lin.Treatment<-lmer(RT ~ WordCond.Treatment + (1|Subject) + (1|Item), data=d) #linear model

log.Treatment<-lmer(Response~WordCond.Treatment +(1|Subject) + (1|Item), data=d, family="binomial" )

###another way to change the base condition that lets you see what levels are being tested

lin.Treatment.relevel<-lmer(RT ~ relevel(WordCond, ref="legal") + (1|Subject) + (1|Item), data=d)

##this example sets the base level to legal nonwords.  Notice, each coefficient is different here because the comparisions being tested are different (i.e., "word vs. legal", "illegal vs. legal")



########### Sets up Effects Coding ############################

d$WordCond.Effects<-d$WordCond

# If you want to make your outputs more readable (easier to figure out what each Ci is doing), then you can manually create the contrast matrix

contrasts(d$WordCond.Effects)<-cbind("illegal.v.GM"= c(1, 0, -1),  "legal.v.GM"= c(0, 1, -1))  #renames Cis to give indication of what is being tested... C1 = illegal vs. grandmean, C2= legal vs.grandmean

#The basic command to create effects coding variables is contr.sum(n), where n = the number of levels of your factor.

lin.Effects<-lmer(RT ~ WordCond.Effects + (1|Subject) + (1|Item), data=d) #output is exactly the same as before, but now the levels of output give description of contrasts

log.Effects<-lmer(Response ~ WordCond.Effects + (1|Subject) + (1|Item), data=d, family="binomial") 

########### Sets up Helmert Coding ############################

# R has a default command for creating helmert contrast weights contr.helmert(); however, these weights don't result in directly interpretable coefficients. So, we'll need to set them by hand.


#Let's make each Ci equal the difference between the means we're testing.
d$WordCond.Helm.Reg<-d$WordCond

contrasts(d$WordCond.Helm.Reg)<-cbind("leg.vs.ill"= c(-.5, .5, 0),"word.vs.nons"=c (-(1/3), -(1/3), (2/3)) )  #renames Cis to give indication of what is being tested... C1 = illegal vs. legal, C2= word vs.nonwords(mean of other two levels) 


lin.Helm.Reg<-lmer(RT ~ WordCond.Helm.Reg + (1|Subject) + (1|Item), data=d) # This model allows for directly interpretable Bs - difference between conditions

log.Helm.Reg<-lmer(Response ~ WordCond.Helm.Reg + (1|Subject) + (1|Item), data=d, family="binomial") # This model allows for directly interpretable Bs - difference between conditions


########### Sets up Polynomial Coding ############################

d$WordCond.Poly<-d$WordCond

contrasts(d$WordCond.Poly)<-contr.poly(3) #This specifies that you want to do polynomial coding (tests linear and quadratic components for this example with 3-levels)

#Orthogonal coding scheme, but interpretations of Cis are different because you're not testing for differences among group means

lin.Poly<-lmer(RT ~ WordCond.Poly + (1|Subject) + (1|Item), data=d)


# Output shows .L for linear component and .Q for quadratic.  If the .L coeff is significant, then the regression requires a linear component, if the .Q coeff is significant, the effect has a quadratic component.

#########Rsquare comparisons#######################

rsq.Treatment<-cor(d$RT, fitted(lin.Treatment))^2
rsq.Effects<-cor(d$RT, fitted(lin.Effects))^2
rsq.Helm<-cor(d$RT, fitted(lin.Helm.Reg))^2
rsq.Poly<-cor(d$RT, fitted(lin.Poly))^2



###############CATEGORICAL INTERACTION##########################
#These are linear models....

#d$NoiseCond (default is to treatment code )
l.Treatment.int<-lmer(RT ~ NoiseCond*WordCond.Treatment + (1|Subject) + (1|Item), data=d)

###############CATEGORICAL X CONTINUOUS INTERACTION##############


l.Treatment.cont<-lmer(RT ~ Freq*WordCond.Treatment + (1|Subject) + (1|Item), data=d)


####centered frequency variable

d$c.Freq<-d$Freq-mean(d$Freq)

l.Treatment.cont.c<-lmer(RT ~ c.Freq*WordCond.Treatment + (1|Subject) + (1|Item), data=d)


p.Treatment.cont<-pvals.fnc(l.Treatment.cont)
p.Treatment.cont.c<-pvals.fnc(l.Treatment.cont.c)


rsq.Treatment<-cor(d$RT, fitted(l.Treatment.cont))^2
rsq.Treatment.c<-cor(d$RT, fitted(l.Treatment.cont.c))^2




#########Non-Orthogonal Coding and Collinearity ################

##set up coding schemes by hand...creates new coding variables that have the appropriate coding scheme... here we're setting up the equivalent of treatment coding with base of "word"

d$Legal.v.Word<- ifelse(d$WordCond == "legal", 1, 0)
d$Illegal.v.Word<- ifelse(d$WordCond == "illegal", 1, 0)

nonOrth<-lmer(RT~Legal.v.Word+Illegal.v.Word + (1|Subject)+ (1|Item), data=d)

#compare nonOrth to the previous model we ran with treatment coding using contr.treatment() and "word" as our baseline.
summary(lin.Treatment)

##run models excluding one of the coding predictors
#use model comparison to determine if the predictor representing the comparison significantly improves model fit


nonOrth1<-lmer(RT~Legal.v.Word + (1|Subject)+ (1|Item), data=d)
nonOrth2<-lmer(RT~Illegal.v.Word + (1|Subject)+ (1|Item), data=d)

anova(nonOrth, nonOrth1)
anova(nonOrth, nonOrth2)

###so what should we determine here?


##more complex example of setting up a model by hand with interaction of two categorical variables

d$Silent.v.Noise<-ifelse(d$NoiseCond == "silent", 1, 0)


###full model (compare to l.Treatment.int)
nonOrthInt<-lmer(RT~ Silent.v.Noise + Legal.v.Word + Illegal.v.Word + Silent.v.Noise:Legal.v.Word + Silent.v.Noise:Illegal.v.Word + (1|Subject) + (1|Item), data=d)

###reduced models removing each of the interaction terms
nonOrthInt2<-lmer(RT~ Silent.v.Noise + Legal.v.Word + Illegal.v.Word + Silent.v.Noise:Legal.v.Word  + (1|Subject) + (1|Item), data=d)

nonOrthInt3<-lmer(RT~ Silent.v.Noise + Legal.v.Word + Illegal.v.Word  + Silent.v.Noise:Illegal.v.Word + (1|Subject) + (1|Item), data=d)

###compares the fully crossed model to the reduced models
anova(nonOrthInt, nonOrthInt2)
anova(nonOrthInt, nonOrthInt3)