library(lme4) # loads the library containing the lmer and glmer commands for multilevel modeling library(lattice) # use for basic graphing ############################################################################################################ # Load a data set. The data is in CSV format, so feel free to open it to examine the structure. # R will look for the data in the working directory. If your file isn't there, you either need to # move it there (use "getwd()" to find out what it is) or change the working directory (it's in your menu) ############################################################################################################ ############################################################################################################ # I'm assuming that your data is in the long format, in comma-delimited form, and has no missing values. # If you have missing values, you may need to use na.strings in the read.csv call. # This data set always has SS delays that are zero. The data were simulated to include people with different # k values from a Kirby questionnaire - not ideal for the approach being advocated. ############################################################################################################ myd<-read.csv("Example one group data glmer.dat") # The following is the basic multilevel model without an intercept myr1<-glmer(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial) # This version assumes a linear relationship for both predictors myr2<-glmer(Waited~I(LLmagn/SSmagn)+I(LLdelay+1)-1+(I(LLmagn/SSmagn)+I(LLdelay+1)-1|Subject), data=myd, family=binomial) # Look at the results summary(myr1) summary(myr2) AIC(myr1, myr2) # To compare the models # Note that a likelihood ratio test could normally be performed using "anova(myr1, myr2)," but an LR test is only valid # when the two models are "nested" meaning that one contains a subset of the terms used in the other model. This is not # the case for these two models. # Look at the results. Note that plotting binomial data isn't very interesting, but I'll try it anyway in the second graph. # The second graph will jitter the 0/1 values to make them easier to see. # These graphs aren't super pretty. Better control can be obtained using ggplot which will be demonstrated in a separate file. xyplot(fitted(myr1)~log(LLmagn/SSmagn)|as.factor(Subject), data=myd, type="a", ylab="P(Waited)") xyplot(fitted(myr1)+jitter(Waited)~log(LLmagn/SSmagn)|as.factor(Subject), data=myd, type=c("a", "p"), distribute.type=T, ylab="P(Waited)") ############################################################################################################ # Let's examine individual differences. I'll focus on the myr1 model. ############################################################################################################ # An aside - for this simulated data, the Kirby survey was used assuming a hyperbolic discounting function. This means that # the delay and magnitude predictors are highly correlated (r = -.97). Thus, the individual subject beta weights are also # highly correlated (r = 1.0). This confirms the problems with the Kirby for identifying the independent contributions of # delay and magnitude. coef(myr1) # shows the individual regression weights # Plots the combinations; each dot is an individual. Not very interesting for this simulated data. plot(coef(myr1)$Subject[,1] ~ coef(myr1)$Subject[,2], ylab="log(LLmagn/SSmagn)", xlab="log(LLdelay+1)") # A histogram of each of the regression weights hist(coef(myr1)$Subject[,1], xlab="log(LLmagn/SSmagn) slope") hist(coef(myr1)$Subject[,2], xlab="log(LLdelay+1) slope") ############################################################################################################ # Now, let's use data with non-zero SS delays. ############################################################################################################ myd<-read.csv("Simulated NonZero data.dat") # The following model uses the log ratios discussed in the manuscript. myrNZ<-glmer(Waited~log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1+(log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1|Subject), data=myd, family=binomial) # Look at the results summary(myrNZ) # Look at the results. xyplot(fitted(myrNZ)~log(LLmagn/SSmagn)|as.factor(Subject), group=cut(log((LLdelay+1)/(SSdelay+1)), 3), data=myd, type="a", ylab="P(Waited)", auto.key=T) # The above plot mixes lots of combinations of continous values of delay and magnitude, so it's particularly ugly. # This analysis really would be better served by a ggplot of the model fit. However, I'm going to # demonstrate another alternative using a bubble plot where the color of each dot indicates the likelihood # of waiting for each combination of the two ratios. In the example below, green = wait, red = don't. myd$fit<-fitted(myrNZ) # grabs the P(Waiting) predicted by the model xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd, col=rgb(1-sqrt(myd$fit), sqrt(myd$fit), 0), # This sets the color in RGB notation. panel=function(x,y,..., col, subscripts) { panel.xyplot(x,y,col=col[subscripts]) } ) # Here I'm plotting the same data but where the size of the dot corresponds to the P(Waiting). # The problem here is that the small dots are too small for most purposes. Dot sizes are # designated using the "cex" parameter. Here I'm changing the range using a power function. xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd, cex=myd$fit^.25, panel=function(x,y,..., cex, subscripts) { panel.xyplot(x,y,cex=cex[subscripts], col="black") } ) # So, here I'm adding a constant to the point size (designed with "cex") to increase the minimum # size and changing the power to .10 to shrink the range of circle sizes. xyplot(log(LLmagn/SSmagn)~log((LLdelay+1)/(SSdelay+1))|as.factor(Subject), data=myd, cex=myd$fit^.1+.2, panel=function(x,y,..., cex, subscripts) { panel.xyplot(x,y,cex=cex[subscripts], col="black") } )