Extracting standard errors from random effects of class GAMM in r

7
votes

I have built a generalized additive mixed effects model with a fixed effect and a random intercept effect (that is a categorical variable). After running the model I am able to extract the random intercepts per each of my categories using ranef(m1$lme)$x[[1]]. However, when I try to extract the standard errors of the random effects using se.ranef(m1$lme), the function does not work. Other attempts using se.ranef(m1) and se.ranef(m1$gam) dont work either. I don't know if this is because these function only apply to models of the class lmer?

Can any one help me so that I can pull out my standard errors of my random intercept from a class "gamm"? I would like to use the random intercepts and standard errors to plot the Best Linear Unbiased Predictors of my gamm model.

My initial model is of the form: gamm(y ~ s(z), random = list(x = ~1), data = dat).

library(mgcv)
library(arm)

example <- gamm(mag ~ s(depth), random = list(stations = ~1), data = quakes)
summary(example$gam)

#Family: gaussian 
#Link function: identity 

#Formula:
 # mag ~ s(depth)

#Parametric coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  5.02300    0.04608     109   <2e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Approximate significance of smooth terms:
#  edf Ref.df     F p-value    
#s(depth) 3.691  3.691 43.12  <2e-16 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#R-sq.(adj) =  0.0725   
#Scale est. = 0.036163  n = 1000

ranef(example$lme)$stations[[1]] # extract random intercepts
#se.ranef(example$lme) # extract se of random intercepts - Problem line - doesn't work?
1
rrandom

1 Answers

9
votes

I'm not an expect on the inner workings of nlme::lme but I don't think it is easy to get what you want from that model --- the ranef() method doesn't allow for the posterior or conditional variance of the random effects to be returned, unlike the method for models fitted by lmer() and co.

Two options spring to mind

  1. fit the model using gamm4:gamm4(), where the mixed model face of the object should work with se.ranef(), or
  2. fit the model with gam() using the random effect basis.

Setup

library("mgcv")
library("gamm4")
library("arm")

library("ggplot2")
library("cowplot")
theme_set(theme_bw())

Option 1: gamm4::gamm4()

This is a straight translation of your model to the syntax required for gamm4::gamm4()

quakes <- transform(quakes, fstations = factor(stations))

m1 <- gamm4::gamm4(mag ~ s(depth), random = ~ (1 | fstations),
                   data = quakes)
re1 <- ranef(m1$mer)[["fstations"]][,1]
se1 <- se.ranef(m1$mer)[["fstations"]][,1]

Note I convert stations to a factor as mgcv::gam needs a factor to fit a random intercept.

Option 2: mgcv::gam()

For this we use the random effects basis. The theory of penalised spline models shows that if we write the math down in a particular way, the model has the same form as a mixed model, with the wiggly parts of the basis acting as random effects and the infinitely smooth parts of the basis used as fixed effects. The same theory allows the reverse process; we can formulate a spline basis that is fully penalised, which is the equivalent of a random effect.

m2 <- gam(mag ~ s(depth) + s(fstations, bs = "re"),
          data = quakes, method = "REML")

We also need to do a little more work to get the "estimated" random effects and standard errors. We need to predict from the model at the levels of fstations. We also need to pass in values for the other terms in the models but as the model is additive we can ignore their effect and just pull out the random effect.

newd <- with(quakes, data.frame(depth = mean(depth),
                                fstations = levels(fstations)))
p <- predict(m2, newd, type = "terms", se.fit = TRUE)
re2 <- p[["fit"]][ , "s(fstations)"]
se2 <- p[["se.fit"]][ , "s(fstations)"]

How do these options compare? 

re <- data.frame(GAMM = re1, GAM = re2)
se <- data.frame(GAMM = se1, GAM = se2)

p1 <- ggplot(re, aes(x = GAMM, y = GAM)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  coord_equal() +
  labs(title = "Random effects")

p2 <- ggplot(se, aes(x = GAMM, y = GAM)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  coord_equal() +
  labs(title = "Standard errors")

plot_grid(p1, p2, nrow = 1, align = "hv")

enter image description here

The "estimates" are the equivalent, but their standard errors are somewhat larger in the GAM version.