## MIXED-EFFECTS REGRESSION LAB SESSION, MARCH 31, 2011, 9:15 - 11:00
## Solutions to exercises
## Martijn Wieling, http://www.martijnwieling.nl


# EXERCISE 1: add the random intercept for Transcriber and store the model in
#             mixedModel3. Print the results and test if the inclusion of this 
#             intercept is necessary. 
mixedModel3 = lmer(LD ~ GeoGAM + PopCnt + PopAge + PopIncome + IsNoun + WordFreq + VCratio + (1|Location) + (1|Word) + (1|Transcriber), data=dialectdata)
print(mixedModel3,cor=F)
anova(mixedModel2,mixedModel3)


# EXERCISE 2: Check if a random slope for PopAge is necessary and visualize
#             the random slopes. 
mixedModel7 = lmer(LD ~ GeoGAM + PopCnt + PopAge + PopIncome + IsNoun + WordFreq + VCratio + (1|Location) + (1|Word) + (0+PopCnt+PopIncome|Word) + (0+PopAge|Word) + (1|Transcriber), data=dialectdata)
print(mixedModel7,cor=F)
anova(mixedModel6,mixedModel7)
myXYPlot(dialectdata,"LD","PopAge","Word",mixedModel6,mixedModel7)


# EXERCISE 3: Check if the random slope for PopAge should be correlated with
#             the other random slopes (this will take a few minutes to run)
mixedModel8 = lmer(LD ~ GeoGAM + PopCnt + PopAge + PopIncome + IsNoun + WordFreq + VCratio + (1|Location) + (1|Word) + (0+PopCnt+PopIncome+PopAge|Word) + (1|Transcriber), data=dialectdata)
print(mixedModel8,cor=F)
anova(mixedModel7,mixedModel8)


# EXERCISE 4: Use your last model to see if a random slope
# for IsNoun is necessary to include. Also check if it has to be correlated
# with the random intercept. First think about if the slope of IsNoun varies
# per Word or Location?
mixedModel9 = lmer(LD ~ GeoGAM + PopCnt + PopAge + PopIncome + IsNoun + WordFreq + VCratio + (1|Location) + (0+IsNoun|Location) + (1|Word) + (0+PopCnt+PopIncome+PopAge|Word) + (1|Transcriber), data=dialectdata)
print(mixedModel9,cor=F)
anova(mixedModel8,mixedModel9)

timestart = proc.time();
mixedModel10 = lmer(LD ~ GeoGAM + PopCnt + PopAge + PopIncome + IsNoun + WordFreq + VCratio + (1+IsNoun|Location) + (1|Word) + (0+PopCnt+PopIncome+PopAge|Word) + (1|Transcriber), data=dialectdata)
duration = proc.time() - timestart; duration; 

print(mixedModel10,cor=F)
anova(mixedModel9,mixedModel10)


# EXERCISE 5: visualize the by-word random slopes for PopCnt and PopIncome
myCoefPlot(mixedModel10,"Word","PopCnt")
myCoefPlot(mixedModel10,"Word","PopIncome")


# EXERCISE 6: visualize the all remaining pairs of by-word random slopes
myCoefPlot2(mixedModel10,"Word","PopCnt","PopAge")
myCoefPlot2(mixedModel10,"Word","PopCnt","PopIncome")


# EXERCISE 7: If interested, take a look at the files
#             datacreation.R and myFunctions.R
# No answer necessary.