joshd<-read.csv("Beckmann timing data.dat") summary(joshd) # center and effect code - avoids non-essential multicollinearity joshd$c.trial<-joshd$Trial - mean(joshd$Trial) joshd$c.duration<-joshd$Duration - mean(joshd$Duration) contrasts(joshd$Audio.Visual) = contr.sum(2) # All three predictors are within-subject, so all three considered in RE structure # Given complexity, will start by trying to figure out the RE structure myr <- glmer(Long.~c.trial*c.duration*Audio.Visual+(c.trial+c.duration+Audio.Visual|Subject), data=joshd, family=binomial) summary(myr) # Now the real work starts! Simplifying to ensure convergence. Convergence issues often due to too much RE flexibility. # Indeed, the REs are highly correlated in this analysis. So, which ones to eliminate? # I explored simpler models until the following converged and had the lowest AIC. # If you have difficulties, this is where the art comes in. BUT, always indicate in paper what you did. myr2 <- glmer(Long.~c.trial*c.duration*Audio.Visual+(c.duration-1|Subject), data=joshd, family=binomial) summary(myr2) # BTW, if you compare the parameter estimates and their SEs from 'myr' and 'myr2', they are very similar which # gives us greater confidence that the simpler RE structure is acceptable. # Now, what does all of this output mean? # Much more likely to choose 'long' as duration increases (duh!) # More likely to choose 'long' if Audio.Visual = 1 than if it's -1 (Wait, which is which?) contrasts(joshd$Audio.Visual) # Ah, so more likely to choose 'long' if the stimulus also has a sound change. # Finally, the duration slope is steeper as the trials increase (not surprising) # We could use lsmeans and lstrends to obtain simple effects library(lsmeans) lsmeans(myr2, "Audio.Visual") # To turn into P(long), use inv.logit in the boot library, or use the following but don't trust the SE because # it should be asymmetrical around the mean lsmeans(myr2, "Audio.Visual", transform="response") # This simply confirms the very similar duration slopes for the two conditions lstrends(myr2, ~ Audio.Visual, var="c.duration") # To examine the slope for particular values of c.trial - about the 25th percentile and 75th percentile lstrends(myr2, ~ 1|c.trial, var="c.duration", at=list(c.trial=c(-10, 10))) # much greater sensitivity as the experiment proceeds # okay, now let's do some plots ggplot(joshd, aes(x = c.duration+mean(Duration), y = fitted(myr2), linetype=Audio.Visual)) + geom_smooth(span=1)+xlab("Duration")+ylab("P(Long)") # Note, this is not a publishable graph because the error region is not correct and the smooth distorts the logistic fit # MORE WORK to do the best graph - note all variables in the model must be in newdat newdat <- expand.grid( c.duration = seq(0.5, 1.7, .005) - .99375 , Audio.Visual = c("Sound", "Static") , c.trial = c(12, 24, 36, 48) - 24.5 , Long. = 0 ) contrasts(newdat$Audio.Visual)=contr.sum(2) # ABSOLUTELY NECESSARY to match analysis modelToPlot<-myr2 mm <- model.matrix(terms(modelToPlot),newdat) newdat$Long. <- mm %*% fixef(modelToPlot) pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm)) newdat <- data.frame( newdat , PLong = inv.logit(newdat$Long.) , plo = inv.logit(newdat$Long.-sqrt(pvar1)) , phi = inv.logit(newdat$Long.+sqrt(pvar1)) ) # Plotting behavior as a function of Delay and Level pd <- position_dodge(width = 0) # need to uncenter c.trial in new column for plotting - must be done before ggplot command! newdat$Trial<-newdat$c.trial+24.5 # I'm changing the order of things here to ensure that the line is plotted on top of my error "bars" g1 <- ggplot(newdat, aes(x=c.duration+.99375, y= PLong, linetype=Audio.Visual)) g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=0, position = pd, col="gray")+ xlab("Duration")+ylab("P(Long)\n")+theme_bw()+theme(text = element_text(size=16))+facet_grid(.~Trial)+geom_line(position = pd, size=1, colour="black")+ ylim(0,1) + xlim(.5, 1.8) # IF TIME, an RT analysis using gamma ggplot(joshd, aes(x = c.duration+mean(Duration), y = RT, linetype=Audio.Visual)) + geom_smooth(span=1)+xlab("Duration")+ylab("Mean(RT)") myr <- glmer(RT~c.trial*c.duration*Audio.Visual+c.trial*I(c.duration^2)*Audio.Visual+(c.trial+c.duration+Audio.Visual|Subject), data=joshd, family=Gamma(link="log")) summary(myr) ggplot(joshd, aes(x = c.duration+mean(Duration), y = fitted(myr), linetype=Audio.Visual)) + geom_smooth(span=1)+xlab("Duration")+ylab("Predicted(RT)")