Generalized Additive Mixed Models
We start by showing how galamm can be used to estimated
generalized additive mixed models.
Gaussian Responses
The cognition dataset contains simulated data with
measurements of abilities in three cognitive domains.
head(cognition)
#> id domain x timepoint item trials y
#> 1 1 1 0.06475113 1 11 1 0.16788973
#> 2 1 1 0.06475113 1 12 1 0.08897838
#> 3 1 1 0.06475113 1 13 1 0.03162123
#> 4 1 1 0.15766278 2 11 1 0.46598362
#> 5 1 1 0.15766278 2 12 1 0.84564656
#> 6 1 1 0.15766278 2 13 1 0.20549872For this first example, we focus only on the first item measured for the first domain.
Each subject in this dataset has been measured eight times, and we can plot the measurements as follows:
plot(dat$x, dat$y, type = "n", xlab = "x", ylab = "y")
for(i in unique(dat$id)) {
dd <- dat[dat$id == i, ]
lines(dd$x, dd$y, col = "gray")
}
points(dat$x, dat$y, pch = 20, lwd = .05)
We use a generalized additive mixed model with random intercepts per subject to estimate the function relating \(x\) to \(y\). In terms of the model framework outlined in the introductory vignette, we model the \(i\)th response from the \(j\)th subject with
\[ y_{ij} = f(x_{ij}) + \eta_{j} + \epsilon_{ij} \]
where \(f(x_{ij})\) is a smooth function to be estimated, \(\eta_{j} \sim N(0, \psi)\) is a random intercept, and \(\epsilon_{ij} \sim N(0, \phi)\) is a residual term.
This model can be estimated using gamm4 as follows:
The package gamm4 uses lme4 to fit the
underlying model, and the resulting model has two components.
mod_gamm4$mer contains the mixed model representation,
whereas in mod_gamm4$gam the fixed and random effects
corresponding to spline coefficients have been converted into single
smooth terms. We can look at the model summary for each:
summary(mod_gamm4$mer)
#> Linear mixed model fit by maximum likelihood ['lmerMod']
#>
#> AIC BIC logLik deviance df.resid
#> 3025.2 3052.1 -1507.6 3015.2 1595
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.93755 -0.65215 0.00612 0.62654 3.14289
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.8551 0.9247
#> Xr s(x) 2.0341 1.4262
#> Residual 0.2501 0.5001
#> Number of obs: 1600, groups: id, 200; Xr, 8
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> X(Intercept) 1.26938 0.06657 19.067
#> Xs(x)Fx1 -0.15814 0.20156 -0.785
#>
#> Correlation of Fixed Effects:
#> X(Int)
#> Xs(x)Fx1 0.000
summary(mod_gamm4$gam)
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> y ~ s(x)
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.26938 0.06657 19.07 <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(x) 6.681 6.681 324.9 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.253
#> lmer.REML = 3015.2 Scale est. = 0.25012 n = 1600We can also plot the estimated smooth term:
In contrast, invoking the plot function on the mixed
model part gives us a diagnostic plot.
With galamm we use similar argument, but the
random specification is now part of the model formula.
As opposed to gamm4, galamm gives a single
summary. As can be seen, smooth terms are both reported as random
effects, and in a separate line under the header “Approximate
significance of smooth terms:”. Reassuringly, the results from fitting
the model with gamm4 and with galamm are
essentially equally, even though they use somewhat different
computational algorithms.
summary(mod)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ s(x) + (1 | id)
#> Data: dat
#>
#> AIC BIC logLik deviance df.resid
#> 3025.2 3052.1 -1507.6 3015.2 1595
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.93755 -0.65215 0.00612 0.62654 3.14290
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.8551 0.9247
#> Xr s(x) 2.0346 1.4264
#> Residual 0.2501 0.5001
#> Number of obs: 1600, groups: id, 200; Xr, 8
#>
#> Fixed effects:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.2694 0.06657 19.0672 4.731e-81
#> s(x)Fx1 -0.1582 0.20236 -0.7818 4.343e-01
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(x) 6.681 6.681 324.9 <2e-16The plot function now gives us a diagnostic plot, which
by inspection can be seen to be almost identical to the plot produced
from the mixed model part of the gamm4 model.
In order to plot the smooth term, we use
plot_smooth.
The plot_smooth function is a thin wrapper around the
plot.gam function provided by the mgcv package
(Wood 2017).
This means that the arguments used by plot.gam can be used
also here, as see with the examples below:
Binomial Responses
In the cognition dataset, the responses relating to domain 2 are binomially distributed. We will use the first trial to illustrate how such data can be modeled.
Again we can fit this model using gamm4.
We can look at the summary output as before.
summary(mod_gamm4$mer)
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#>
#> AIC BIC logLik deviance df.resid
#> 983.7 1005.2 -487.8 975.7 1596
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -8.2547 0.0786 0.1946 0.4248 0.8213
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.2306 0.4802
#> Xr s(x) 0.9694 0.9846
#> Number of obs: 1600, groups: id, 200; Xr, 8
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> X(Intercept) 2.8115 0.2034 13.824 < 2e-16 ***
#> Xs(x)Fx1 -1.4110 0.4242 -3.326 0.00088 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Correlation of Fixed Effects:
#> X(Int)
#> Xs(x)Fx1 -0.440
summary(mod_gamm4$gam)
#>
#> Family: binomial
#> Link function: logit
#>
#> Formula:
#> y ~ s(x)
#>
#> Parametric coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.8115 0.1662 16.91 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(x) 2.14 2.14 113.9 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.116
#> glmer.ML = 908.96 Scale est. = 1 n = 1600And we can plot the smooth term. The diagnostic plot is not very useful in the binomial case, so we omit it.
Again the galamm syntax is similar, but it puts the
random effect specification into the model formula.
The estimates are very similar, although not identical. The
difference in deviance is due to differences in the way deviance is
defined. The call deviance(mod_gamm4$mer) gives the same
value as in the summary for the model fitted with galamm.
summary(mod)
#> GALAMM fit by maximum marginal likelihood.
#> Formula: y ~ s(x) + (1 | id)
#> Data: dat
#>
#> AIC BIC logLik deviance df.resid
#> 983.7 1005.2 -487.8 908.8 1596
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -8.2237 0.0792 0.1947 0.4245 0.8226
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> id (Intercept) 0.2316 0.4813
#> Xr s(x) 0.9387 0.9688
#> Number of obs: 1600, groups: id, 200; Xr, 8
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.808 0.1956 14.354 9.989e-47
#> s(x)Fx1 -1.406 0.4134 -3.402 6.697e-04
#>
#> Approximate significance of smooth terms:
#> edf Ref.df Chi.sq p-value
#> s(x) 2.124 2.124 114 <2e-16
