1

Hi I've started fitting simple gams with random intercepts using mgcv and bs='re', but I can't seem to find how to extract the conditional modes/random effects/BLUPs for each level of my factor. ie the equivalent of ranef() in lmer. My googling has drawn a blank as well. Is this possible?

Mike Dunbar
  • 11

1 Answers1

0

The easiest way is to either

  1. identify which model coefficients are associated with a specific random effect smooth, or
  2. evaluate the random effect smooth at the levels of the grouping factor that you want.

Neither of these is especially difficult with mgcv, but my gratia package does make doing them somewhat easier. Using a simple example:

load_mgcv()
data(sleepstudy, package = "lme4")
m <- gam(Reaction ~ Days + s(Subject, bs = "re") +
           s(Days, Subject, bs = "re"),
         data = sleepstudy, method = "REML")

We can extract the coefficients for a selected smooth using smooth_coefs():

> library("gratia")
> smooth_coefs(m, "s(Subject)")                                                                             
 s(Subject).1  s(Subject).2  s(Subject).3  s(Subject).4  s(Subject).5  s(Subject).6  s(Subject).7 
    1.5127209   -40.3739742   -39.1811090    24.5189267    22.9144576     9.2219858    17.1561444 
 s(Subject).8  s(Subject).9 s(Subject).10 s(Subject).11 s(Subject).12 s(Subject).13 s(Subject).14 
   -7.4517412     0.5786996    34.7679974   -25.7543565   -13.8650381     4.9159796    20.9290802 
s(Subject).15 s(Subject).16 s(Subject).17 s(Subject).18 
    3.2586540   -26.4758514     0.9056463    12.4217779

(if you want their indices among the vector of coefficients, use smooth_coef_indices().)

A tidier output is provided by smooth_estimates(), which implements the second idea:

> smooth_estimates(m)
# A tibble: 1,818 × 7
   smooth     type          by        est    se Subject  Days
   <chr>      <chr>         <chr>   <dbl> <dbl> <fct>   <dbl>
 1 s(Subject) Random effect NA      1.51   13.3 308        NA
 2 s(Subject) Random effect NA    -40.4    13.3 309        NA
 3 s(Subject) Random effect NA    -39.2    13.3 310        NA
 4 s(Subject) Random effect NA     24.5    13.3 330        NA
 5 s(Subject) Random effect NA     22.9    13.3 331        NA
 6 s(Subject) Random effect NA      9.22   13.3 332        NA
 7 s(Subject) Random effect NA     17.2    13.3 333        NA
 8 s(Subject) Random effect NA     -7.45   13.3 334        NA
 9 s(Subject) Random effect NA      0.579  13.3 335        NA
10 s(Subject) Random effect NA     34.8    13.3 337        NA
# ℹ 1,808 more rows
# ℹ Use `print(n = ...)` to see more rows

I plan to add a ranef() method but haven't got round to it just yet.

Gavin Simpson
  • 47,626