Example with interactions
This file demonstrates how to run a model with multiple predictors taht interact. A good understanding of centering and effect coding from regression will help. My focus will be on the nlme syntax
We begin by loading the necessary R libraries.
The Basics - Loading data, preparing for analysis
library(nlme) # loads the library used for nonlinear multilevel modeling commands
library(lattice) # loads the library that includes xyplotLoad data and set up for analysis using effect coding (contr.sum) and centered predictors (by subtracting the mean).
myd<-read.csv("Young E1 DelayDiscounting.dat")
myd$Cluster<-as.factor(myd$Cluster)
contrasts(myd$Sex)=contr.sum(2)
contrasts(myd$Cluster)=contr.sum(3)
myd$c.VGExp<-myd$VGExp-mean(myd$VGExperience, na.rm=T) # create centered experience variableRun the analysis
Run a model with the full interaction between Sex and (centered) VG Experience. For starting values we need 4: male-intercept, female-intercept, male-VG slope, female-VG slope. Note that by using effect coding, the interpretation is different, but thinking of it this way allows you to count up the necessary number of starting values.
myr<-nlme(Indifference.point~1000/(1+exp(logk)*Delay), fixed=logk~Sex*c.VGExp, random=logk~1|Subject, data=subset(myd, !is.na(c.VGExp)), start=c(logk=-6, 0, 0, 0))
summary(myr)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Indifference.point ~ 1000/(1 + exp(logk) * Delay)
## Data: subset(myd, !is.na(c.VGExp))
## AIC BIC logLik
## 2233.021 2251.693 -1110.511
##
## Random effects:
## Formula: logk ~ 1 | Subject
## logk.(Intercept) Residual
## StdDev: 2.505964 143.1833
##
## Fixed effects: logk ~ Sex * c.VGExp
## Value Std.Error DF t-value p-value
## logk.(Intercept) -6.746734 0.8316339 128 -8.112626 0.0000
## logk.Sex1 1.385526 0.8316339 128 1.666029 0.0982
## logk.c.VGExp 0.274983 0.1754461 128 1.567338 0.1195
## logk.Sex1:c.VGExp -0.068744 0.1754461 128 -0.391822 0.6958
## Correlation:
## lg.(I) lgk.S1 l..VGE
## logk.Sex1 -0.772
## logk.c.VGExp -0.751 0.803
## logk.Sex1:c.VGExp 0.803 -0.751 -0.844
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.869146311 -0.438504679 -0.002750696 0.256652846 2.771505684
##
## Number of Observations: 166
## Number of Groups: 35As in regression based on centered and effect-coded variables, the intercept is the grand mean k. The regression equation is -6.75 + 1.39*Sex + .27*c.VGExp - .069*Sex*c.VGExp. To derive predictions, Sex is -1 for males and +1 for females (see the contrasts).
Each p-value evaluates each term: Intercept > 0, Sex1 dummy > 0 (main effect of sex), VG Exp > 0 (main effect of VG Exp), Sex1:c.VGExp > 0 (interaction).
You can also use the emmeans library to derive for group means and (supposedly) group slopes.
library(emmeans)
emmeans(myr, "Sex", param="logk")## NOTE: Results may be misleading due to involvement in interactions## Sex emmean SE df lower.CL upper.CL
## f -5.361208 0.5553082 128 -6.45998 -4.262436
## m -8.132260 1.5464251 128 -11.19213 -5.072394
##
## Confidence level used: 0.95To compare the effects of video game experience separately for each sex, you should be able to use the emtrends command, but I can’t get it to work for this model. So, I’ll use the more laborious multcomp library. Although the non-significant Sex*VGExp interaction indicates these slopes don’t differ, I’ll demonstrate how to estimate the two slopes.
The four elements in the matrix correspond to the four fixed effects in the model summary. Since we are testing slopes, the 1 in the third column is the main effect for slope and the +1 or -1 in the fourth column indicates the sex adjustments.
# The following isn't working
#emtrends(myr, "Sex", var="c.VGExp", param="logk") # Alternatively you can use the following syntax
library(multcomp)
contrast.matrix <- rbind(
'Male VG Exp slope' = c(0, 0, 1, 1),
'Female VG Exp slope' = c(0, 0, 1, -1)
)
comps <- glht(myr, contrast.matrix)
summary(comps)##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: nlme.formula(model = Indifference.point ~ 1000/(1 + exp(logk) *
## Delay), data = subset(myd, !is.na(c.VGExp)), fixed = logk ~
## Sex * c.VGExp, random = logk ~ 1 | Subject, start = c(logk = -6,
## 0, 0, 0))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Male VG Exp slope == 0 0.20624 0.09679 2.131 0.0651 .
## Female VG Exp slope == 0 0.34373 0.33285 1.033 0.5125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Note that the the effect of videogame experience on the discounting rate is a bit weaker for the males, but there’s a lot of overlap.
par(mar=c(4,10,4,2)) # necessary to change margins to avoid truncation of labels
plot(comps)
Interacting categorical variables
Now things get complicated for my 2 x 3 design. Six parameters are estimated. One of the challenges is counting the number of parameters for the starting values. When they’re categorical, it’s easy: 2(Sex) x 3(Cluster) = 6.
myr<-nlme(Indifference.point~1000/(1+exp(logk)*Delay), fixed=logk~Sex*Cluster, random=logk~1|Subject, data=subset(myd, !is.na(c.VGExp)&!is.na(Cluster)), start=c(logk=-6, 0, 0, 0, 0, 0))
summary(myr)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Indifference.point ~ 1000/(1 + exp(logk) * Delay)
## Data: subset(myd, !is.na(c.VGExp) & !is.na(Cluster))
## AIC BIC logLik
## 2177.139 2201.79 -1080.569
##
## Random effects:
## Formula: logk ~ 1 | Subject
## logk.(Intercept) Residual
## StdDev: 2.517503 145.5232
##
## Fixed effects: logk ~ Sex * Cluster
## Value Std.Error DF t-value p-value
## logk.(Intercept) -6.100571 0.7346577 122 -8.303964 0.0000
## logk.Sex1 0.217885 0.7346577 122 0.296581 0.7673
## logk.Cluster1 0.969568 1.0648788 122 0.910496 0.3644
## logk.Cluster2 -1.226881 0.9231293 122 -1.329046 0.1863
## logk.Sex1:Cluster1 -0.335137 1.0648788 122 -0.314718 0.7535
## logk.Sex1:Cluster2 -0.674706 0.9231293 122 -0.730890 0.4662
## Correlation:
## lg.(I) lgk.S1 lgk.C1 lgk.C2 l.S1:C1
## logk.Sex1 -0.004
## logk.Cluster1 0.070 -0.678
## logk.Cluster2 -0.335 0.169 -0.373
## logk.Sex1:Cluster1 -0.678 0.070 -0.472 0.425
## logk.Sex1:Cluster2 0.169 -0.335 0.425 0.129 -0.373
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.752183568 -0.513404742 -0.002461481 0.270338455 2.727017331
##
## Number of Observations: 161
## Number of Groups: 34To derive the group estimates, you can again use emmeans. The following also shows the brainless post hoc approach (Tukey’s) testing for all (15!) possible pairwise comparisons.
emmeans(myr, ~Sex*Cluster, param="logk")## Sex Cluster emmean SE df lower.CL upper.CL
## f 1 -5.248255 0.5966088 122 -6.429302 -4.06720833
## m 1 -5.013752 2.5513227 122 -10.064350 0.03684652
## f 2 -7.784273 1.5660770 122 -10.884479 -4.68406746
## m 2 -6.870631 1.0756504 122 -8.999989 -4.74127416
## f 3 -4.615529 2.5500953 122 -9.663698 0.43263893
## m 3 -7.070986 1.3145705 122 -9.673310 -4.46866235
##
## Confidence level used: 0.95pairs(emmeans(myr, ~Sex*Cluster, param="logk"))## contrast estimate SE df t.ratio p.value
## f,1 - m,1 -0.2345035 2.620151 122 -0.089 1.0000
## f,1 - f,2 2.5360182 1.675870 122 1.513 0.6565
## f,1 - m,2 1.6223764 1.230027 122 1.319 0.7740
## f,1 - f,3 -0.6327258 2.618955 122 -0.242 0.9999
## f,1 - m,3 1.8227308 1.443620 122 1.263 0.8046
## m,1 - f,2 2.7705216 2.993634 122 0.925 0.9392
## m,1 - m,2 1.8568798 2.768803 122 0.671 0.9848
## m,1 - f,3 -0.3982223 3.607247 122 -0.110 1.0000
## m,1 - m,3 2.0572343 2.870077 122 0.717 0.9796
## f,2 - m,2 -0.9136418 1.899900 122 -0.481 0.9967
## f,2 - f,3 -3.1687440 2.992588 122 -1.059 0.8965
## f,2 - m,3 -0.7132874 2.044674 122 -0.349 0.9993
## m,2 - f,3 -2.2551022 2.767672 122 -0.815 0.9643
## m,2 - m,3 0.2003545 1.698564 122 0.118 1.0000
## f,3 - m,3 2.4554566 2.868986 122 0.856 0.9561
##
## P value adjustment: tukey method for comparing a family of 6 estimatesIt’s better to do targeted comparisons. In this situation, I might want to know male/female differences for each cluster (3 comparisons instead of 15).
pairs(emmeans(myr, ~Sex|Cluster, param="logk"))## Cluster = 1:
## contrast estimate SE df t.ratio p.value
## f - m -0.2345035 2.620151 122 -0.089 0.9288
##
## Cluster = 2:
## contrast estimate SE df t.ratio p.value
## f - m -0.9136418 1.899900 122 -0.481 0.6315
##
## Cluster = 3:
## contrast estimate SE df t.ratio p.value
## f - m 2.4554566 2.868986 122 0.856 0.3938If you compare these results to the prior 15-comparisons output, you’ll notice that the p-values drop because you aren’t penalized for those 12 other comparisons. They’re still not significant, but this is the proper approach for this data set.
ANOVA?
It would be great to get an ANOVA table for this situation due to the three-level cluster variable, but the anova command uses the wrong sum-of-squares (sequential). Easy solution:
anova(myr, type="marginal")## numDF denDF F-value p-value
## logk.(Intercept) 1 122 68.95581 <.0001
## logk.Sex 1 122 0.08796 0.7673
## logk.Cluster 2 122 0.98293 0.3772
## logk.Sex:Cluster 2 122 0.46767 0.6276