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) Loading required package: nlme This is mgcv 1.8-31. For overview type 'help("mgcv-package")'. library(lme4) Loading required package: Matrix 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)) head(compare) 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.