How to add random intercept and random bias to a GAMM model in R
I am trying to specify both random intercept and random bias in a single fixed effect GAMM model.
I have successfully simulated a random intercept model using the below code in the library mgcv
, but now I can't figure out what the syntax is for random skew inside a function gamm()
:
M1 = gamm(dur ~ s(dep, bs="ts", k = 4), random= list(fInd = ~1), data= df)
If I were to use random intercept and slope in a linear mixed effects model, I would write it like this:
M2 = lme(dur ~ dep, random=~1 + dep|fInd, data=df)
The accompanying documentation gamm()
states that random members should be specified in the form list
as in lme()
, but I cannot find a single interpreted example that includes both slope and hooking points. Any advice / solution would be much appreciated.
source to share
The function gamm4
in the package gamm4
contains a way to do this. You define random intercept and slope just like in style lmer
. In your case:
M1 = gamm4(dur~s(dep,bs="ts",k=4), random = ~(1+dep|fInd), data=df)
Here is the gamm4 documentation: https://cran.r-project.org/web/packages/gamm4/gamm4.pdf
source to share
Here is the syntax gamm()
for injecting correlated random intercept and slope effects using a dataset sleepstudy
.
library(nlme)
library(mgcv)
data(sleepstudy,package='lme4')
# Model via lme()
fm1 <- lme(Reaction ~ Days, random= ~1+Days|Subject, data=sleepstudy, method='REML')
# Model via gamm()
fm1.gamm <- gamm(Reaction ~ Days, random= list(Subject=~1+Days), data=sleepstudy, method='REML')
VarCorr(fm1)
VarCorr(fm1.gamm$lme)
# Both are identical
# Subject = pdLogChol(1 + Days)
# Variance StdDev Corr
# (Intercept) 612.0795 24.740241 (Intr)
# Days 35.0713 5.922103 0.066
# Residual 654.9424 25.591843
The syntax for introducing uncorrelated random intercept and slope effects is the same for lme()
and gamm()
.
# Model via lme()
fm2 <- lme(Reaction ~ Days, random= list(Subject=~1, Subject=~0+Days), data=sleepstudy, method='REML')
# Model via gamm()
fm2.gamm <- gamm(Reaction ~ Days, random= list(Subject=~1, Subject=~0+Days), data=sleepstudy, method='REML')
VarCorr(fm2)
VarCorr(fm2.gamm$lme)
# Both are identical
# Variance StdDev
# Subject = pdLogChol(1)
# (Intercept) 627.5690 25.051328
# Subject = pdLogChol(0 + Days)
# Days 35.8582 5.988172
# Residual 653.5838 25.565285
This answer also shows you how to inject multiple random effects into lme()
.
source to share