This package was specifically for mixed models using mgcv. At present it has been almost entirely superseded by the mixedup package, which provides all the same functionality and more, while the visualizations were originally from the visibly package anyway. That only leaves the prediction functionality as unique to this package, so I leave it here as I may still use it for that, and I’m contemplating moving the GAM specific visuals from visibly at some point.

# Package Description

## Introduction

The goal of gammit is to provide a set of functions to aid using mgcv (possibly solely) for mixed models. Lately I’ve been using it in lieu of lme4, especially the bam function, for GLMM with millions of observations and multiple random effects. It’s turning out very useful in this sense (see this post for details), but I’d like some more/different functionality with the results. Furthermore, mgcv just has some nice things going on for such models anyway, like the ability to add other smooth terms, alternative distributions for the target variable, etc., so I’m looking to make it easier for me to get some things I want when I use it.

At present there are four functions: extract_vc, extract_ranef, extract_fixed, summary_gamm, and predict_gamm.

## Installation

You can install the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("m-clark/gammit")


## Example

This example demonstrates the summary_gamm function with comparison to the corresponding lme4 model.

library(mgcv)
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(lme4)

Attaching package: 'lme4'
The following object is masked from 'package:nlme':

lmList
library(gammit)

lmer_model = lmer(Reaction ~  Days + (Days || Subject), data=sleepstudy)

ga_model = gam(Reaction ~  Days + s(Subject, bs='re') + s(Days, Subject, bs='re'),
data=sleepstudy,
method = 'REML')

summary(lmer_model)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
Data: sleepstudy

REML criterion at convergence: 1743.7

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.9626 -0.4625  0.0204  0.4653  5.1860

Random effects:
Groups    Name        Variance Std.Dev.
Subject   (Intercept) 627.57   25.051
Subject.1 Days         35.86    5.988
Residual              653.58   25.565
Number of obs: 180, groups:  Subject, 18

Fixed effects:
Estimate Std. Error t value
(Intercept)  251.405      6.885  36.513
Days          10.467      1.560   6.712

Correlation of Fixed Effects:
(Intr)
Days -0.184

summary_gamm(ga_model)

Variance components:
group    effect variance     sd sd_2.5 sd_97.5 var_prop
Subject Intercept  627.571 25.051 16.085  39.015    0.477
Subject      Days   35.858  5.988  4.025   8.908    0.027
Residual            653.582 25.565 22.792  28.676    0.496

Fixed Effects:
Term Estimate Std. Error t value Pr(>|t|)
(Intercept)  251.405      6.885  36.513    0.000
Days   10.467      1.560   6.712    0.000

Extract the variance components with extract_vc.

data.frame(VarCorr(lmer_model))
grp        var1 var2      vcov     sdcor
1   Subject (Intercept) <NA> 627.56905 25.051328
2 Subject.1        Days <NA>  35.85838  5.988187
3  Residual        <NA> <NA> 653.58350 25.565279

extract_vc(ga_model)
group    effect variance     sd sd_2.5 sd_97.5 var_prop
1  Subject Intercept  627.571 25.051 16.085  39.015    0.477
2  Subject      Days   35.858  5.988  4.025   8.908    0.027
3 Residual            653.582 25.565 22.792  28.676    0.496

Extract the random effects with extract_ranef.

ranef(lmer_model)
$Subject (Intercept) Days 308 1.5126648 9.3234970 309 -40.3738728 -8.5991757 310 -39.1810279 -5.3877944 330 24.5189244 -4.9686503 331 22.9144470 -3.1939378 332 9.2219759 -0.3084939 333 17.1561243 -0.2872078 334 -7.4517382 1.1159911 335 0.5787623 -10.9059754 337 34.7679030 8.6276228 349 -25.7543312 1.2806892 350 -13.8650598 6.7564064 351 4.9159912 -3.0751356 352 20.9290332 3.5122123 369 3.2586448 0.8730514 370 -26.4758468 4.9837910 371 0.9056510 -1.0052938 372 12.4217547 1.2584037 with conditional variances for "Subject" extract_ranef(ga_model) # A tibble: 36 x 7 group_var effect group value se lower_2.5 upper_97.5 <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> 1 Subject Intercept 308 1.51 13.3 -24.5 27.5 2 Subject Intercept 309 -40.4 13.3 -66.4 -14.3 3 Subject Intercept 310 -39.2 13.3 -65.2 -13.2 4 Subject Intercept 330 24.5 13.3 -1.51 50.5 5 Subject Intercept 331 22.9 13.3 -3.11 48.9 6 Subject Intercept 332 9.22 13.3 -16.8 35.2 7 Subject Intercept 333 17.2 13.3 -8.87 43.2 8 Subject Intercept 334 -7.45 13.3 -33.5 18.6 9 Subject Intercept 335 0.579 13.3 -25.4 26.6 10 Subject Intercept 337 34.8 13.3 8.74 60.8 # … with 26 more rows Extract the fixed effects extract_fixef. fixef(lmer_model) (Intercept) Days 251.40510 10.46729 extract_fixed(ga_model) # A tibble: 2 x 7 term value se t p lower_2.5 upper_97.5 <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 Intercept 251. 6.88 36.5 0 238. 265. 2 Days 10.5 1.56 6.71 0 7.39 13.5 ## Prediction There are a couple of ways to do prediction, and the main goal for gammit was to make it easy to use the lme4 style to include random effects or not. mgcv already has this functionality as well, so the functionality of predict_gamm is mostly cosmetic. One benefit here is to provide standard errors for the prediction also. head(predict_gamm(ga_model)) prediction 1 252.9178 2 272.7086 3 292.4994 4 312.2901 5 332.0809 6 351.8717 Add standard errors. head(data.frame(predict_gamm(ga_model, se=T))) prediction se 1 252.9178 12.410220 2 272.7086 10.660891 3 292.4994 9.191224 4 312.2901 8.153871 5 332.0809 7.724998 6 351.8717 8.003034 compare = data.frame( gam_original = predict_gamm(ga_model)$prediction,
gam_fe_only   = predict_gamm(ga_model, re_form = NA)$prediction, gam_fe_only2 = predict_gamm(ga_model, exclude = c('s(Subject)', "s(Days,Subject)"))$prediction,
lme4_fe_only  = predict(lmer_model, re.form = NA))

gam_original gam_fe_only gam_fe_only2 lme4_fe_only
1     252.9178    251.4051     251.4051     251.4051
2     272.7086    261.8724     261.8724     261.8724
3     292.4994    272.3397     272.3397     272.3397
4     312.2901    282.8070     282.8070     282.8070
5     332.0809    293.2742     293.2742     293.2742
6     351.8717    303.7415     303.7415     303.7415

Along with that, one can still use include/exclude for other smooth terms as above. Unfortunately, some options do not yet work with bam objects, but this is to due to the functionality in predict.gam from mgcv and should change in the near future.