CODEWITHSUNDEEP
rlme4mgcvrandom-effects
olyashevska
olyashevska

Reputation: 447

Variance-covariance matrix for random effect from GAM with mgcv package

Random effect and variance-covariance matrix of random effect with lme4 package are extracted as following:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.rr <- ranef(fm1,condVar=TRUE)
fm1.pv <- attr(rr[[1]],"postVar")

I wonder how I can do this with mgcv? 'gam.vcomp' function does extract the estimated variance components, but not for each level of random effect.

library(mgcv)
fm2 <- gam(Reaction ~ Days + s (Subject, bs="re"), data = sleepstudy, method = "REML")
gam.vcomp(fm2)

Upvotes: 1

Views: 1411

Answers (2)

JWilliman
JWilliman

Reputation: 3883

You can also use the package gamm4, which is based on the gamm package but using lme4 underneath. The model would be fitted as:

fm3 <- gamm4(Reaction ~ Days, random = ~ (1|Subject), data = sleepstudy)

Random effects and variance-covariance matrix of random effects can be obtained following the normal lme4 procedure.

fm3.rr <- ranef(fm3$mer,condVar=TRUE)
fm3.pv <- attr(fm3.rr[[1]],"postVar")[1,1,]

However gamm4 can be much slower than gam so read the help file to see when it best suits your need.

Upvotes: 1

olyashevska
olyashevska

Reputation: 447

library(lme4)
data(sleepstudy)

fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.rr <- ranef(fm1,condVar=TRUE)$Subject[,1]
fm1.pv <- sqrt(attr(ranef(fm1,condVar=TRUE) [['Subject']],"postVar")[1,1,])

library(mgcv)
fm2 <- gam(Reaction ~ Days + s (Subject, bs="re"), 
data = sleepstudy,   method = "REML")

To extract random effect for each Subject 

idx <-grep("Subject", names(coef(fm2)))
fm2.rr<-coef(fm2)[idx]
attributes(fm2.rr)<-NULL

We can see that random effects in both models are identical as expected.

To extract variance-covariance matrix for random effect and calculate an error we use parameter Vp which is a Bayesian posterior covariance matrix:

fm2.pv <-sqrt(diag(fm2$Vp))[idx]

Or frequentist estimated covariance matrix Ve 

fm2.pv <-sqrt(diag(fm2$Ve))[idx]

We can see that random effect errors estimated with mgcv slightly differ that those estimated with lme4 model. Errors based on a Bayesian posterior covariance matrix are larger, whereas based on a frequentist matrix are smaller.

Upvotes: 3

Related Questions