library(lme4) library(lsmeans) library(lattice) library(ggplot2) myd <- read.csv("S-D Mixed and Num icons entropy.dat") summary(myd) # Note renaming of columns when illegal characters appear # I'm going to drop the Mixed type to avoid trying to discuss a nested design! myd <- subset(myd, Type!="Mixed") myd$Type <- factor(myd$Type) # necessary to tell it there are now only two legal levels contrasts(myd$Type)=contr.sum(2) # to do effect coding and thus match JMP's output contrasts(myd$Type) myd$c.numOfIcons <- myd$X..of.icons - mean(myd$X..of.icons) # to avoid non-essential multicollinearity myr <- lmer(X..judged.different~Type*c.numOfIcons + (Type*c.numOfIcons|Bird), data=myd) summary(myr) xyplot(fitted(myr)+X..judged.different~X..of.icons|Type, group=Bird, data=myd, type=c("a", "p")) # To obtain Tukey's https://cran.r-project.org/web/packages/lsmeans/vignettes/using-lsmeans.pdf lsmeans(myr, "Type") # to examine condition differences pairs(lsmeans(myr, "Type")) # Tukey's (uninteresting here) lstrends(myr, ~ Type, var="c.numOfIcons") # to examine slope differences # for nesting problem lsmeans(myr, c("Type", "X..of.icons"), at=list(Type="Mixed", X..of.icons=16)) pairs(lsmeans(myr, c("Type", "X..of.icons"), at=list(X..of.icons=16))) # only compares for common value of 16 icons ############### For Repeated measures logistic regression example ###################### myd <- read.csv("S-D Mixed disaggregated.dat") summary(myd) myr<-glmer(ChoseDifferent~entropy+(entropy|Bird), family=binomial, data=myd, weights=Freq) summary(myr) xyplot(fitted(myr)~entropy, group=Bird, data=myd, type="a",auto.key=T, ylab="P(Different)") # Note, can't just plot ChoseDifferent here as a DV because I was forced to use 'weights' because I don't # have trial-by-trial data. But, I can do this with ggplot. ggplot(myd, aes(x = entropy, y = ChoseDifferent)) + geom_smooth(aes(weight = myd$Freq, linetype=Bird), se = FALSE, size = 1.5) ggplot(myd, aes(x = entropy, y = fitted(myr))) + geom_smooth(aes(weight = myd$Freq, linetype=Bird), se = FALSE, size = 1.5) ranef(myr) coef(myr) # In general, no reason to test residuals for a logistic regression, either their distribution or homogeneity hist(residuals(myr)) # Deviance residuals ############## PLOTTING with error bars newdat <- expand.grid( entropy = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4) , ChoseDifferent = 0 ) modelToPlot<-myr mm <- model.matrix(terms(modelToPlot),newdat) newdat$ChoseDifferent <- mm %*% fixef(modelToPlot) pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm)) newdat <- data.frame( newdat , PDiff = inv.logit(newdat$ChoseDifferent) , plo = inv.logit(newdat$ChoseDifferent-sqrt(pvar1)) , phi = inv.logit(newdat$ChoseDifferent+sqrt(pvar1)) ) # Plotting behavior as a function of Delay and Level pd <- position_dodge(width = 0) g1 <- ggplot(newdat, aes(x=entropy, y= PDiff))+geom_line(position = pd)+geom_point(position = pd)+ ylim(0,1) g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.2, position = pd)+ xlab("Entropy")+ylab("P(Different)\n")+theme_bw()+theme(text = element_text(size=18)) ### Error region newdat <- expand.grid( entropy = seq(0,4,.01) , ChoseDifferent = 0 ) modelToPlot<-myr mm <- model.matrix(terms(modelToPlot),newdat) newdat$ChoseDifferent <- mm %*% fixef(modelToPlot) pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm)) newdat <- data.frame( newdat , PDiff = inv.logit(newdat$ChoseDifferent) , plo = inv.logit(newdat$ChoseDifferent-sqrt(pvar1)) , phi = inv.logit(newdat$ChoseDifferent+sqrt(pvar1)) ) # Plotting behavior as a function of Delay and Level pd <- position_dodge(width = 0) g1 <- ggplot(newdat, aes(x=entropy, y= PDiff))+geom_line(position = pd)+ ylim(0,1) g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=0, position = pd, col="gray", alpha=.5)+ xlab("Entropy")+ylab("P(Different)\n")+theme_bw()+theme(text = element_text(size=18)) # Model comparison myr1<-glmer(ChoseDifferent~entropy+(entropy|Bird), family=binomial, data=myd, weights=Freq) myr2<-glmer(ChoseDifferent~entropy+(1|Bird), family=binomial, data=myd, weights=Freq) AIC(myr1, myr2) anova(myr1, myr2)