Linear and Generalized Mixed Models
Matt Chesney
11/15/2022
Introduction
This is a tutorial to walk through of how to complete Mixed regression in R. The tutorial is in two parts, the first is mixed regression utilizing linear methods and part 2 is generalized linear regression.
Mixed models
What are mixed models? introduce fixed and random effects and how they differ.
In an experiment there is a unit that the experiment is being conducted on, this is the experimental unit. For the Experimental unit there is some measured attribute of interest, this is the independent variable or our Y.
A fixed effect is a treatment that is applied to our experimental unit. The treatment’s effect is ascertained by a shift in the average value of Y, or the independent variable being monitored.
The NULL hypothesis is that the treatment has no effect on the mean. \(\beta\) is our Fixed effect.
\[H_0: \beta=0\] \[H_a: \beta \neq 0\]
Examples:
1. Comparing tensile strengths of alloy blends.
2. Comparing weight gain in pigs from several food types.
3. Comparing maximum speeds of cars with different fuel additives.
4. Comparing the blood pressure change in a patient over several days of medication.
A random effect is also applied to the experimental unit. The treatment’s effect is determined by a difference in the variation in Y due to the application of the treatment. A random effect is assigned when the population of treatments is assumed to be homogeneous.
The NULL hypothesis is that the treatment has no effect on the variance. \(\sigma_{R}\) is the variance attributed to our treatment.
\[H_0: \sigma_{R}=0\] \[H_a: \sigma_{R}\neq 0\] Examples:
1. A factory’s various production lines effect on variation tensile strength.
2. A pig farm’s pig pens’ effect on variation in weight gain.
3. A professional driver’s effect on variation in top speeds of cars.
4. A medical subject’s effect on blood pressure medication treatment.
The mixed model contains both fixed effects and random effects.There are several ways this relationship can manifest. The fixed effect can have a consistent effect within the random effect, intercept changes but slope remains the same. The fixed effect random effect can only vary in the presence of the fixed effect, intercept constant and slope varies. The fixed effect can vary with respect to the random effect, slope and intercept change. Each of these options will be examined in our sleep data example.
Examples:
1. Comparing tensile strengths of alloy blends from several production lines in a factory.
2. Comparing weight gain in pigs from several food types living in one of several randomly assigned pens
3. Comparing maximum speeds of cars with different fuel additives driven by randomly assigned professional drivers.
4. Comparing the average blood pressure reduction after a treatment over several days applied to randomly selected subjects.
We will now look at a dataset that has one random variable and one fixed variable.
Introduction to the Data
The sleep study dataset is available as part of the lme4 package. It is from a study conducted by Belenkyn et all (2003).
The response variable is the reaction time of the subjects to stimuli as recorded by researchers Reaction. The subjects are the random effect, Subject. The number of days the subject underwent sleep deprivation is the fixed effect Days.
Below is the structure of the sleep study data:
str(sleepstudy)## 'data.frame': 180 obs. of 3 variables:
## $ Reaction: num 250 259 251 321 357 ...
## $ Days : num 0 1 2 3 4 5 6 7 8 9 ...
## $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...The first few lines of data. These are all observations for subject 308.
head(sleepstudy)## Reaction Days Subject
## 1 249.5600 0 308
## 2 258.7047 1 308
## 3 250.8006 2 308
## 4 321.4398 3 308
## 5 356.8519 4 308
## 6 414.6901 5 308Plot of sleep study data, color coded by patient numbers:
ggplot(sleepstudy,aes(x=Days, y=Reaction, group=Subject, color=Subject)) +
geom_point(size=1) + theme_classic()
The Model
This example will follow the convention of using Greek symbols to denote fixed effects and capital letters to denote random effects.
The variables in our model are: \(y_{ij}\) is each day’s reaction time (j) of our subject (i), \(c_i\) our random intercept term caused by the subject, \(\beta\) is the Days coefficient, and \(D\) is our days value. .
Model 1, Only intercept varies [Fig 1]
\[y_{ij} = \alpha + c_{i} + \beta \times D_{j} + e_{ij}\] Model 1 assumptions: \(C_i \underset{i.i.d.}{\thicksim}(0,\sigma^2_c)\)
\(e_ij\underset{i.i.d.}{\thicksim}(0,\sigma^2)\)
Model 2 has all the same terms as model 1 along with \(d_i\) as our random slope term due to subject and a correlation term for the intercept and slope effect due to subject : \(\rho_{c_j,d_i} \ne 0\) [Fig 2]
\[y_{ij} = \alpha + c_{i} + (\beta_+d_{i})\times D_j + e_{ij}\]
Model 2 assumptions: \([\begin{smallmatrix}c \\ d\end{smallmatrix}] \underset{i.i.d.}{\thicksim} N(0,\sum)\) where: \(\sum=\big(\begin{smallmatrix}\sigma^2_c & \rho_{c,d}\\\rho_{c,d} & \sigma^2_d\end{smallmatrix}\big)\)
\(e\) \((0,\sigma^2)\)
c and d are distributed jointly as indicated and not correlated with e, e is not correlated with any other factor.
Model 3, intercept and slope vary and random terms are not correlated: \(\rho_{(c,d)} = 0\) [Fig 3] ()
\[y_{ij} = \alpha + c_{i} + (\beta+d_{i})\times D_j+ e_{ij}\]
Model 3 assumptions: \(C_d\underset{i.i.d.}{\thicksim}(0,\sigma^2_d)\)
\(C_i\underset{i.i.d.}{\thicksim}(0,\sigma^2_c)\)
\(e_ij\underset{i.i.d.}{\thicksim}(0,\sigma^2)\)
c,d, and e are all independent and not correlated with each other
LME4 Package
To analyze the mixed effect model we will utilize the R package LME4. The relevant syntax for the model is as follows:
| Model Type | Syntax | Abbreviated Syntax | |
|---|---|---|---|
| Random intercept | (1|RandomEffect) |
||
| Random Slope effect | (0+FixedEffect|RandomEffect) |
||
| Random intercept and slope with random effect correlated | 1+FixedEffect+(1+FixedEffect|RandomEffect) |
FixedEffect+(FixedEffect|RandomEffect |
|
| Random intercept and slope with random effect uncorrelated | 1+FixedEffect+(1|RandomEffect)+(0+FixedEffect|RandomEffect) |
FixedEffect+(FixedEffect||RandomEffect) |
We will utilize the abbreviated syntax in this demonstration.
REML
Restricted Maximum likelihood is recommended for the final model, however REML does not work with model selection techniques such as AIC and BIC. If it is necessary to compare models REML should not be used. Once the best model is selected REML can be used with that model for more accurate confidence intervals.This can be achieved in the LME4 package by adding REML=False as an option. The default in LME4 is true.
For model comparison we will use the option REML=F, once the model is selected we will use REML
Model 1 [Fig 1]
[MtrxEx]
model1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy,REML=F)
summary(model1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
##
## AIC BIC logLik deviance df.resid
## 1802.1 1814.9 -897.0 1794.1 176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2347 -0.5544 0.0155 0.5257 4.2648
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 1296.9 36.01
## Residual 954.5 30.90
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.4051 9.5062 26.45
## Days 10.4673 0.8017 13.06
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.380Output for Model 1
The summary of the model gives us AIC, BIC, and log likelihood scores along with deviance and degrees of freedom. After the scaled residuals we get parameter estimates. Listed first are estimates for the random effects. The \(c_i\) standard deviation is estimated at 36.01 and the \(e_{ij}\) is estimated at 30.9. Next are the estimates for the fixed effects of intercept and the Days term of 251.405 and 10.467 respectively.
Confidence intervals for Model 1
The lme4 package contains a confint.mermod() function that computes confidence intervals for fixed and random effects. The function is shown below with all of its parameters.
# S3 method for merMod credit: https://www.rdocumentation.org/packages/lme4/versions/1.1-29/topics/confint.merMod
confint(object, parm, level = 0.95,
method = c("profile", "Wald", "boot"), zeta,
nsim = 500,
boot.type = c("perc","basic","norm"),
FUN = NULL, quiet = FALSE,
oldNames = TRUE, ...)By default all parameters will be output at a 95% interval utilizing the profile likelihood method. The function will return confidence intervals for all random effects followed by fixed effects in the order that they appear in the summary output. In order to extract a specific confidence interval simply select its numerical order followed by a comma to get the entire row of output. Below we show the entire output and then only the fixed effects extracted.
NOTE: oldNames=F will return more informative names rather than .sig01, .sig02, etc.
confint.merMod(model1, oldNames = F)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|Subject 26.007120 52.93598
## sigma 27.813847 34.59105
## (Intercept) 231.992326 270.81788
## Days 8.886551 12.04802confint.merMod(model1)[3:4,]## Computing profile confidence intervals ...## 2.5 % 97.5 %
## (Intercept) 231.992326 270.81788
## Days 8.886551 12.04802Visulization of Model 1
The lines in this graph are the model for each of the Subjects.
\[ \hat{y_{ij}}=\hat{\alpha}+\hat{c_{i}}+\hat{\beta}*Days_{j}\]
ggplot(NULL) +
geom_point(data=sleepstudy,aes(x=Days, y=Reaction, group=Subject, color=Subject))+geom_line(data=sleepstudy%>%mutate(pred_dist=fitted(model1)),aes(x=Days, y=pred_dist, group=Subject, color=Subject))+ theme_classic()
Model 2 [Fig 2]
model2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,REML=F)
summary(model2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
## Data: sleepstudy
##
## AIC BIC logLik deviance df.resid
## 1763.9 1783.1 -876.0 1751.9 174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9416 -0.4656 0.0289 0.4636 5.1793
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 565.48 23.780
## Days 32.68 5.717 0.08
## Residual 654.95 25.592
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.405 6.632 37.907
## Days 10.467 1.502 6.968
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.138Output for Model 2
As in model 1 we get estimates of our \(c_i\) term for the intercept of 23.780 for st dev. In addition to the output received in model 1, we also get an estimate for the \(d_i\) term in the random effects listed in the Days row with std dev of 5.717. We also get the \(\rho\) estimate for \((c_i , d_i)\) which is .08. Finally our fixed effects estimates follow with values of 251.405 and 10.467 for intercept of Days respectively.
Confidence intervals for Model 2
For out fixed effects we have:
confint.merMod(model2)[5:6,]## Computing profile confidence intervals ...## 2.5 % 97.5 %
## (Intercept) 237.680695 265.12951
## Days 7.358653 13.57592In this model we also have our correlation term for our random effects with a point estimate of .08, the 95% confidence interval is shown below:
confint.merMod(model2)[2,]## Computing profile confidence intervals ...## 2.5 % 97.5 %
## -0.4815004 0.6849861NOTE: When there is a correlation term for the random effects the confint.merMod() function places the correlation term in between the two random effects it applies to.
Visualization of Model 2
The lines on this plot represent each of the subjects with the correlation term following the equation: \[\hat{y_{ij}} = \hat{\alpha_{intercept}} + \hat{c_{i}} + (\hat{\beta} + \hat{d_i})*Days_{j} \]
ggplot(NULL) +
geom_point(data=sleepstudy,aes(x=Days, y=Reaction, group=Subject, color=Subject))+geom_line(data=sleepstudy%>%mutate(pred_dist=fitted(model2)),aes(x=Days, y=pred_dist, group=Subject, color=Subject))+ theme_classic()
Model 3 [Fig 3]
model3 <- lmer(Reaction ~ Days + (Days||Subject), sleepstudy,REML=F)
summary(model3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
## Data: sleepstudy
##
## AIC BIC logLik deviance df.resid
## 1762 1778 -876 1752 175
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9535 -0.4673 0.0239 0.4625 5.1883
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 584.27 24.172
## Subject.1 Days 33.63 5.799
## Residual 653.12 25.556
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.405 6.708 37.48
## Days 10.467 1.519 6.89
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.194Output for Model 3
The output for model 3 is in the same format as model 2 but without the correlation term.
Confidence intervals fixed effects for Model 3
confint.merMod(model3)[4:5,]## Computing profile confidence intervals ...## 2.5 % 97.5 %
## (Intercept) 237.572148 265.23806
## Days 7.334067 13.60051Data visulization for Model 3
Each line represents the subject without a random effects correlation term according to the equation: \[ \hat{y_{ij}}= \hat{\alpha}+\hat{c_{i}} +(\hat{\beta}+\hat{d_i})*Days_{j}\]
ggplot(NULL) +
geom_point(data=sleepstudy,aes(x=Days, y=Reaction, group=Subject, color=Subject))+geom_line(data=sleepstudy%>%mutate(pred_dist=fitted(model3)),aes(x=Days, y=pred_dist, group=Subject, color=Subject))+ theme_classic()
Model Comparison
Here we compare the AIC, BIC, and log likelihood scores for each of the models.
## Data: sleepstudy
## Models:
## model1: Reaction ~ Days + (1 | Subject)
## model3: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
## model2: Reaction ~ Days + (Days | Subject)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model1 4 1802.1 1814.8 -897.04 1794.1
## model3 5 1762.0 1778.0 -876.00 1752.0 42.0754 1 8.782e-11 ***
## model2 6 1763.9 1783.1 -875.97 1751.9 0.0639 1 0.8004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC and BIC select the uncorrelated model ( model 3), log likelihood selects the correlated model (model 2) . We will go with model 3 since the confidence interval for the correlation term in model 2 indicated that it was not likely significant.
Refit with REML
First, as mentioned previously we will refit the selected model with REML=T.
model3 <- lmer(Reaction ~ Days + (Days||Subject), sleepstudy,REML=T)REML Confidence intervals
We have our refit confidence intervals for our selected model below.
confint.merMod(model3,oldNames = F)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|Subject 15.258647 37.786471
## sd_Days|Subject 3.964074 8.769159
## sigma 22.880555 28.787598
## (Intercept) 237.572148 265.238062
## Days 7.334067 13.600505Bootstrap Confidence intervals
Next we use the same function to calculate the bootstrap confidence intervals for all parameters. The default number of iterations is 500. The bootstrap method utilized is a non parametric type leveraging the percentile method ( running the requisite number of iterations then selecting the percentile requested for the specific interval of interest). We will use 1000 iterations for all of our confidence intervals in this demonstration.
Further information on the specifics of this bootstrap method is described in the boot.ci and bootMer packages’ documentation.
Bootstrap method 1: percent
In this method, the number requisite number of sims is run, placed in numerical order, then the quantiles and the appropriate percentage values are selected.This method is appropriate when the test statistic is unbiased and homoscedastic and the standard error of the bootstrap statistic and sample statistics are the same.
confint.merMod(model3, method = 'boot',boot.type = 'perc', oldNames = F, nsim = 1000)## Computing bootstrap confidence intervals ...##
## 1 message(s): boundary (singular) fit: see ?isSingular
## 1 warning(s): Model failed to converge with max|grad| = 0.00274894 (tol = 0.002, component 1)## 2.5 % 97.5 %
## sd_(Intercept)|Subject 13.496320 35.599046
## sd_Days|Subject 3.539577 8.388968
## sigma 22.594612 28.533273
## (Intercept) 237.133900 265.286650
## Days 7.451359 13.507604Bootstrap method 2: basic
In this method, the number requisite number of sims is run, placed in numerical order, then the percentage requested is selected. This method is appropriate when the test statistic is unbiased and homoscedastic and the bootstrap statistic can be transformed into a normal distribution.
confint.merMod(model3, method = 'boot',boot.type = 'basic', oldNames = F, nsim = 1000)## Computing bootstrap confidence intervals ...##
## 1 message(s): boundary (singular) fit: see ?isSingular
## 1 warning(s): Model failed to converge with max|grad| = 0.00241023 (tol = 0.002, component 1)## 2.5 % 97.5 %
## sd_(Intercept)|Subject 14.907629 36.46922
## sd_Days|Subject 3.701596 8.31672
## sigma 22.481408 28.66672
## (Intercept) 238.283358 264.71251
## Days 7.548358 13.26260Bootstrap method 3: norm
In this method, the number requisite number of sims is run, then the standard error of the sample of samples is taken to calculate the confidence interval. This is good if the statistic in question is assumed to be normally distributed and unbiased.
confint.merMod(model3, method = 'boot',boot.type = 'norm', oldNames = F, nsim = 1000)## Computing bootstrap confidence intervals ...##
## 1 message(s): boundary (singular) fit: see ?isSingular## 2.5 % 97.5 %
## sd_(Intercept)|Subject 14.258734 37.287065
## sd_Days|Subject 3.787987 8.394926
## sigma 22.556438 28.545110
## (Intercept) 237.814143 264.680483
## Days 7.325546 13.539793Extracting Design matrix
The matrix instantiating for mixed models can be generalized as: \[ Y= X \beta+Zu+e\] Where Y is the response variable a [nx1] matrix, X is a [n x p] design matrix where n = number of obs and p = number of fixed parameters. \(\beta\) is a [px1] matrix with the values for the fixed parameters. Z is the random design matrix [gxn] matrix, where g is the number of grouping variables. \(e\) is a matrix of residuals from the model of [nx1]. u is a [gx1] matrix containing the estimates for the random effects for each group member. There will be a u for each random effect, in our case we have a u for intercept and a u for the slope effect.
Each element of the design matrix for model 3 can be extracted as seen below:
First we have the X matrix, we extract this using the getME(model,"X") function:
# Fixed design matrix
Xi=getME(model3,"X")
head(Xi)## (Intercept) Days
## 1 1 0
## 2 1 1
## 3 1 2
## 4 1 3
## 5 1 4
## 6 1 5Shown below are the estimates for the random effects in out model or U1 and U2.
#Function to retrieve random effects from the model
ranef(model3)## $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"In order to extract the u elements we can utilize the getME(model,"b"). This will give us all the random effects in a single vector, we will need to split these into each element for intercept and the slope effect.
#Ui is a composite list of all random effects
Ui=getME(model3,"b")Next we will extract the Random design matrices, one for the intercept and one for the slope.
#Ztlist is two design matrices for the random effects, intercept and days
melist=getME(model3,'Ztlist')
dim(melist$`Subject.(Intercept)`)## [1] 18 180dim(melist$Subject.Days)## [1] 18 180We will place these into two separate variables, Zint for intercept and ZDays for the slope effect. We will do the same for the U matrix. Uint and UDays are shown.
Zint=melist$`Subject.(Intercept)`
ZDays=melist$Subject.Days
Uint=Ui[1:18]
UDays=Ui[19:36]
Uint## [1] 1.5126648 -40.3738728 -39.1810279 24.5189244 22.9144470 9.2219759
## [7] 17.1561243 -7.4517382 0.5787623 34.7679030 -25.7543312 -13.8650598
## [13] 4.9159912 20.9290332 3.2586448 -26.4758468 0.9056510 12.4217547UDays## [1] 9.3234970 -8.5991757 -5.3877944 -4.9686503 -3.1939378 -0.3084939
## [7] -0.2872078 1.1159911 -10.9059754 8.6276228 1.2806892 6.7564064
## [13] -3.0751356 3.5122123 0.8730514 4.9837910 -1.0052938 1.2584037Here we extract the beta matrix using getME(model,"beta")
#beta are the fixed effect values
beta=getME(model3,"beta")We will also extract the response variable vector with getME(model,"y")
# y is the response variable matrix
y=getME(model3,"y")Referencing our original model in Fig[2] we will use Zint and Uint to create \(C_i\) and ZDays and UDays to create \(D_i\). We are multiplying \(Zu\) for each random effect. We will have to add \(\beta\) to Udays prior to multiplying by ZDays in order to get the correct value for the combined slope term.
Ci=Uint%*%Zint
Di=(UDays+beta[2])%*%ZDaysBelow we add the intercept estimate to the Ci matrix to get the subject specific intercept. We then add both int and \(D_i\) together to get the model estimate yest. Finally we subtract the model estimates from the observations and take out the residuals to demonstrate we have correctly extracted the matrix. Notice that all values are effectively zero.
#intercept+Ci
int= Ci+beta[1]
yest=int+Di
head(y-yest-residuals(model3))[1:180]## [1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [6] 0.000000e+00 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00
## [11] 0.000000e+00 0.000000e+00 2.842171e-14 0.000000e+00 0.000000e+00
## [16] 2.842171e-14 -2.842171e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [21] 0.000000e+00 0.000000e+00 2.842171e-14 -2.842171e-14 0.000000e+00
## [26] 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [31] 0.000000e+00 5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [36] 5.684342e-14 -5.684342e-14 0.000000e+00 0.000000e+00 5.684342e-14
## [41] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [46] 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [51] 0.000000e+00 0.000000e+00 5.684342e-14 0.000000e+00 0.000000e+00
## [56] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [61] 0.000000e+00 -5.684342e-14 0.000000e+00 -5.684342e-14 -5.684342e-14
## [66] -5.684342e-14 -5.684342e-14 -5.684342e-14 -5.684342e-14 0.000000e+00
## [71] 0.000000e+00 2.842171e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [76] 5.684342e-14 0.000000e+00 -5.684342e-14 0.000000e+00 5.684342e-14
## [81] 0.000000e+00 0.000000e+00 2.842171e-14 0.000000e+00 0.000000e+00
## [86] 2.842171e-14 -2.842171e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [91] 0.000000e+00 0.000000e+00 5.684342e-14 0.000000e+00 0.000000e+00
## [96] 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [101] 0.000000e+00 2.842171e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [106] 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [111] 0.000000e+00 2.842171e-14 0.000000e+00 0.000000e+00 0.000000e+00
## [116] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [121] 0.000000e+00 0.000000e+00 5.684342e-14 0.000000e+00 5.684342e-14
## [126] 0.000000e+00 -5.684342e-14 0.000000e+00 0.000000e+00 5.684342e-14
## [131] 0.000000e+00 0.000000e+00 0.000000e+00 -5.684342e-14 -5.684342e-14
## [136] 0.000000e+00 -1.136868e-13 -5.684342e-14 0.000000e+00 0.000000e+00
## [141] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [146] 0.000000e+00 -5.684342e-14 0.000000e+00 5.684342e-14 5.684342e-14
## [151] 0.000000e+00 2.842171e-14 2.842171e-14 0.000000e+00 5.684342e-14
## [156] 0.000000e+00 0.000000e+00 5.684342e-14 0.000000e+00 0.000000e+00
## [161] 0.000000e+00 0.000000e+00 0.000000e+00 -5.684342e-14 -5.684342e-14
## [166] 0.000000e+00 -5.684342e-14 -5.684342e-14 -5.684342e-14 0.000000e+00
## [171] 0.000000e+00 0.000000e+00 5.684342e-14 0.000000e+00 0.000000e+00
## [176] 5.684342e-14 0.000000e+00 0.000000e+00 0.000000e+00 5.684342e-14Nested Example
We will be using the dataset from John Starkweather at the Univeristy of North Texas. The experimental unit are 50 students within 4 types of classes within 6 schools for a total of 1200 observations. The response variable is the extroversion of students extro as predicted by the fixed effects of agreeableness agree, openness to new experiences open, and social engagement social along with the fixed categorical effects of class type. The random effects are the schools and the classes within schools.
nest <- read.table("http://bayes.acs.unt.edu:8083/BayesContent/class/Jon/R_SC/Module9/lmm.data.txt",
header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE)
head(nest)## id extro open agree social class school
## 1 1 63.69356 43.43306 38.02668 75.05811 d IV
## 2 2 69.48244 46.86979 31.48957 98.12560 a VI
## 3 3 79.74006 32.27013 40.20866 116.33897 d VI
## 4 4 62.96674 44.40790 30.50866 90.46888 c IV
## 5 5 64.24582 36.86337 37.43949 98.51873 d IV
## 6 6 50.97107 46.25627 38.83196 75.21992 d IBelow we show the data looking at one of the fixed effects and coloring by first class type then school. We see a much more clear separation by school in the data with respect to the open variable.
ggplot(NULL) +
geom_point(data=nest,aes(x=class, y=extro, group=school, color=school))
The Model
In this example the model options are simplified. The fixed effects open (\(O_k\)), agree (\(A_k\)), and social (\(S_k\)) are assumed not to have random interactions with the school and class within school. The class variable is categorical and therefore only impacts the slope and can have no further interaction with random effects. So although we have more variables our model options are limited to random slope effects. Therefore the model we will be analyzing is listed below: \[y_{ijk} = \alpha +b_i+ c_{j}(b_i) + \beta*O_{k}+\gamma*A_{k} +\delta*S_{k}+ \tau_j\times X_j+ e_{ijk}\] \(y_{ijk}\) is the extroversion response variable i being school and j being class. \(\alpha\) is the intercept, \(b\) is the random effect for school, \(c(b)\) is the random effect for class nested within school, \(\beta\) is the coefficient for openness, \(\gamma\) is the coefficient for agreeableness, \(\delta\) is the coefficient for social engagement, and \(\tau_j\times X_j\) is the class specific intercept shift mutiplied by a fixed effect design matrix, and \(e_{ijk}\) represents the remaining residual.
Model Assumptions:
\(b_i \underset{i.i.d.}{\thicksim}N(0, \sigma^2_{b_i})\)
\(c_j(b) \underset{i.i.d.}{\thicksim}N(0, \sigma^2_{c(b)})\)
\(O_k,A_k,S_k,e_{ijk} \underset{i.i.d.}{\thicksim}N(0, \sigma^2)\)
A note on nesting and correlation
When a variable is nested within a variable the correlation effect is pooled with the higher effect. This is because it is impossible to discern correlation between a nested effect and the higher order effect since they would be necessarily confounded. In our case, the class within school and school interaction cannot be ascertained since class b within school a is not the same as class b within school b.
LME4 with nested models
In LME4 we can denote nesting with the / symbol. So school/classindicates school is nested within class. With LME4 when a nested variable is included as a random effect, the parent variable (in this case school), will also be included as a random effect in the model.
nestModel=lmer(extro~open+agree+social+class+(1|school/class), REML = T)Next we will use a function called step from the lmerTest library to evaluate each variable. This function utilizes backwards elimination and eliminates variables according to comparative analysis, dropping one variable at a time and comparing AIC scores of the models. We will need the REML=F option on in order to run this function effectively. The step output is listed below.
## Warning: package 'lmerTest' was built under R version 4.1.3##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepNOTE: There is a known error in ‘LmerTest’ that causes an error during the Rstudio Knit process. The resulting model is shown below.
The function indicated that the variables of social, open, and agree do not contribute to the model and are thus eliminated leaving us with a reduced model of:
\[y_{ijk} = \alpha +b_i+ c_{j}(b_i) + \tau_j+ e_{ijk}\]
nestReduced=lmer(extro ~ class + (1 | school/class))Below are the confidence intervals for each of the parameters in the original model:
confint.merMod(nestModel,oldNames = F)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|class:school 1.150044059 2.239429608
## sd_(Intercept)|school 5.496421925 17.908789836
## sigma 0.944381622 1.023898354
## (Intercept) 48.848703467 65.919025409
## open -0.003586661 0.015860294
## agree -0.018896611 0.003430591
## social -0.003100553 0.004156909
## classb 0.196539758 3.913060839
## classc 1.846703019 5.563158376
## classd 3.807484711 7.523983935All of the factors that were eliminated by step traverse 0 at the 95% confidence interval.
We will compare the reduced model to the orginal model using anova function.
anova(nestModel,nestReduced)## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## nestReduced: extro ~ class + (1 | school/class)
## nestModel: extro ~ open + agree + social + class + (1 | school/class)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nestReduced 7 3526.5 3562.1 -1756.2 3512.5
## nestModel 10 3529.0 3579.9 -1754.5 3509.0 3.4142 3 0.3321According to AIC/BIC scoring, the reduced model is the best model.
Below is the summary for the reduced final model. ### Output Nested Model
summary(nestReduced)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: extro ~ class + (1 | school/class)
##
## REML criterion at convergence: 3503.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.9751 -0.3342 0.0078 0.3483 10.6423
##
## Random effects:
## Groups Name Variance Std.Dev.
## class:school (Intercept) 2.8844 1.6983
## school (Intercept) 95.1804 9.7560
## Residual 0.9687 0.9842
## Number of obs: 1200, groups: class:school, 24; school, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 57.4166 4.0432 5.2288 14.201 2.24e-05 ***
## classb 2.0486 0.9838 15.0000 2.082 0.05485 .
## classc 3.6985 0.9838 15.0000 3.759 0.00189 **
## classd 5.6563 0.9838 15.0000 5.749 3.84e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) classb classc
## classb -0.122
## classc -0.122 0.500
## classd -0.122 0.500 0.500Random Effects
As shown in the summary we first have the variation for the random effects listed with class within school first at 1.6983 Std. Dev, school next at 9.7560, and residual effect last at .9842. In this example the majority of the variation is accounted for by the school effect.
Fixed Effects
The fixed effects are shown next starting out with the intercept at 57.4166 and the class intercept estimates of 2.049, 3.699, 5.656 for A,B,C respectively. Notice class B is borderline insignificant at the 95% level. Depending on the goals of the researcher it may make sense to keep or eliminate class B. In this case since class b in one of 3 hierarchical categorical variables, two of which are significant, It may make sense to keep them all for increased model interpretability.
Confidence Intervals
confint.merMod(nestReduced,oldNames = F)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|class:school 1.1501768 2.239738
## sd_(Intercept)|school 5.4964683 17.909994
## sigma 0.9457512 1.025383
## (Intercept) 48.8990897 65.934092
## classb 0.1900934 3.907041
## classc 1.8400749 5.557022
## classd 3.7978019 7.514749As can be seen from the confidence intervals the fixed effects for class b is close to but does not traverse zero.
Bootstrap Confidence Intervals
We will now look at bootstrap confidence intervals first showing our profile likelihood confidence intervals. The intervals are very similar as can be seen.
confint.merMod(nestReduced,oldNames = F)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## sd_(Intercept)|class:school 1.1501768 2.239738
## sd_(Intercept)|school 5.4964683 17.909994
## sigma 0.9457512 1.025383
## (Intercept) 48.8990897 65.934092
## classb 0.1900934 3.907041
## classc 1.8400749 5.557022
## classd 3.7978019 7.514749## Warning: package 'glmm' was built under R version 4.1.3## Loading required package: trust## Loading required package: mvtnorm## Loading required package: parallel## Loading required package: doParallel## Loading required package: foreach## Loading required package: iterators## Warning: package 'glmmTMB' was built under R version 4.1.3## Warning: package 'LaplacesDemon' was built under R version 4.1.3##
## Attaching package: 'LaplacesDemon'## The following objects are masked from 'package:mvtnorm':
##
## dmvt, rmvtGenarlized Mixed models
Sometimes our mixed model data may not be from a normal distribution. Mixed models are already very complex and challenging to implement as intended in R, extra caution is advised when adding non-normality parameters.
There are many different packages in R for working with generalized mixed models. Each of them have their advantages and disadvantages. We will be using two packages; the glmmTMB package and the glmer package.
Data
We will be using the salamander dataset from McCullagh et all 1986. The objective of this experiment was to breed salamander pairs. Each experimental unit is a pair of salamanders. There are two types of salamanders, White side and rough butt. The variable cross indicates the mating pair composition. W/R white side female with a rough butt male.
##Model We are interested in the probability of a successful mating outcome as modeled by the log odds ratio. Our model uses this log odds ratio as the response and our predictors include and intercept term, the fixed categorical variable Cross, and the two random effects of Male and Female. The individual males and females will have random effects on the mating outcome. The Cross composition will be a fixed categorical effect consisting of 3 possibilities. The variable \(\beta\) is a 1x3 vector of the model estimated fixed value estimates for the mating pairs. The variable X is a 3 value dummy variable that has one indicator value indicating the type of mating pair corresponding to cross (one pair is calculated into the intercept term). F is the random effect for a specific female. M is the random effect for the specific male.
\[log(P/(1-P))=\alpha + (\beta \times X) +F+M\] ## Model assumptions
\[F\underset{i.i.d.}{\thicksim}N(0, \sigma^2_{F})\] \[M\underset{i.i.d.}{\thicksim}N(0, \sigma^2_{M})\] \(F and M\) are independent of one another
head(salamander)## Mate Cross Female Male
## 1 1 R/R 10 10
## 2 1 R/R 11 14
## 3 1 R/R 12 11
## 4 1 R/R 13 13
## 5 1 R/R 14 12
## 6 1 R/W 15 28inputing the model with glmmTMB and glmer
The syntax for both packages is the same as previously disused in part 1. In our case we have one categorical fixed effect Y~Fixed and two random effects +(1|randomeEffect).This is the same syntax previously used and is shown below.
We will compare the effectiveness of each package by splitting up our data set. We will train the data set on 70% of our salamander obs and test it on the remaining 30%.
newsal=salamander[,2:4]
trainSplit = createDataPartition(y = salamander[,1], p = 0.7, list = FALSE)
salamanderTrain=salamander[trainSplit,]
salTestX=salamander[-trainSplit,][2:4]
salTestY=salamander[-trainSplit,][1]
sal=glmmTMB(Mate ~ Cross+( 1 | Female) + (1| Male), data = salamanderTrain,
family= binomial(link = 'logit'))
sal2=glmer(Mate ~ Cross+( 1 | Female) + (1| Male), data = salamanderTrain,
family= binomial(link = 'logit'))Output of packages
Below we will look at the outputs for each package.
Summary output glmmTMB
summary(sal)## Family: binomial ( logit )
## Formula: Mate ~ Cross + (1 | Female) + (1 | Male)
## Data: salamanderTrain
##
## AIC BIC logLik deviance df.resid
## 318.5 339.7 -153.2 306.5 246
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Female (Intercept) 0.4658 0.6825
## Male (Intercept) 0.3929 0.6268
## Number of obs: 252, groups: Female, 60; Male, 60
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0188 0.3570 2.854 0.00432 **
## CrossR/W -0.9487 0.4520 -2.099 0.03582 *
## CrossW/R -2.5271 0.5514 -4.583 4.59e-06 ***
## CrossW/W -0.1856 0.4830 -0.384 0.70079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fixed Effects
We have the categorical fixed effect indicating the type of salamander pair. W/W is included in the intercept estimate 1.019. The estimates for R/W, W/R, and W/W are -0.949, -2.527, -0.186 ``
Random Effects
The two random effects stadard deviations, F and M, are estimated to be 0.6824607 for Female and 0.6261523 for Male.
Confidence Intervals
Below are the confidence intervals for each of the factors in the model:
confint(sal)## 2.5 % 97.5 % Estimate
## (Intercept) 0.3191037 1.71843964 1.0187717
## CrossR/W -1.8344825 -0.06282041 -0.9486515
## CrossW/R -3.6078274 -1.44630910 -2.5270683
## CrossW/W -1.1323570 0.76113579 -0.1856106
## Std.Dev.(Intercept)|Female 0.2872755 1.62127509 0.6824607
## Std.Dev.(Intercept)|Male 0.2375406 1.65382918 0.6267787Summary output for GLMER
summary(sal2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Mate ~ Cross + (1 | Female) + (1 | Male)
## Data: salamanderTrain
##
## AIC BIC logLik deviance df.resid
## 318.5 339.7 -153.2 306.5 246
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8517 -0.7651 0.4224 0.6339 2.0273
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female (Intercept) 0.4657 0.6825
## Male (Intercept) 0.3928 0.6268
## Number of obs: 252, groups: Female, 60; Male, 60
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0188 0.3570 2.854 0.00432 **
## CrossR/W -0.9486 0.4520 -2.099 0.03582 *
## CrossW/R -2.5271 0.5514 -4.583 4.59e-06 ***
## CrossW/W -0.1856 0.4830 -0.384 0.70080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) CrsR/W CrsW/R
## CrossR/W -0.691
## CrossW/R -0.664 0.458
## CrossW/W -0.684 0.529 0.460Fixed Effects
We have the categorical fixed effect indicating the type of salamander pair. W/W is included in the intercept estimate 1.019. The estimates for R/W, W/R, and W/W are -0.949, -2.527, -0.186 ``
Random Effects
The two random effects Standard deviations, F and M, are estimated to be 0.682457 for Female and 0.6267659 for Male.
Confidence Intervals
Below are the confidence intervals for each of the factors in the model:
confint.merMod(sal2,method="Wald")## 2.5 % 97.5 %
## .sig01 NA NA
## .sig02 NA NA
## (Intercept) 0.3191113 1.71843090
## CrossR/W -1.8344717 -0.06282553
## CrossW/R -3.6078106 -1.44632333
## CrossW/W -1.1323473 0.76113416NOTE: There is no built in method for estimating confidence intervals for GLMER. We will compare with bootstrap later.
Comparing the two packages prediction capabilty
Below we utilize the package LaplacesDemon to leverage the invlogit function. We use each of the tow models to compare their prediction capability. As can be seen they perform precisely the same. While there are subtle differences between the packages they will often yield very similar results as we have seen here.
salpred=round((invlogit(predict(sal,newdata=salTestX))))
boolglmm=salpred==salTestY
sum(boolglmm)/length(boolglmm)## [1] 0.6388889sal2pred=round((invlogit(as.numeric(predict(sal2,newdata=salTestX)))))
boolglmer=sal2pred==salTestY
sum(boolglmer)/length(boolglmer)## [1] 0.6388889Bootstrap intervals
Before conducting our bootstrap intervals we will refit the models with all available data.
#refit with all data
sal2=glmer(Mate ~ Cross+( 1 | Female) + (1| Male), data = salamander,
family= binomial(link = 'logit'))
sal=glmmTMB(Mate ~ Cross+( 1 | Female) + (1| Male), data = salamander,
family= binomial(link = 'logit'))##
## Attaching package: 'boot'## The following object is masked from 'package:LaplacesDemon':
##
## logit## The following object is masked from 'package:lattice':
##
## melanomaWe will now go through bootstraps for each of the two packages.
Bootstrap for glmmTMB
The glmmTMB package does not include a built in boot strap function. It is possible to write our own as we have below. The FUN function samples from the original data set with replacement from the original number of observations and refits the model. Next the values of interest from our model are returned.
Next we use the lme4::bootmer for mixed model boot straps and put the sims into the variable b1. Finally we output the percentile bootstrap confidence intervals.
#NOTE: Takes considerble time to run
FUN <- function(model=sal){
fm1R <- refit(sal, simulate(sal)[[1]])
return( c(fm1R$obj$report()$sd[[1]],fm1R$obj$report()$sd[[2]],as.numeric(fm1R$fit$par[1:4])))
}
b1 <- lme4::bootMer(sal,FUN, nsim=100, .progress="txt")## ================================================================================ if (requireNamespace("boot")) {
cat("95% Bootstrap confidence intervals for glmmTMB\n")
for ( i in 1: length(namesBoot)){
cat(c(namesBoot[i], ": "))
print(boot.ci(b1,type="perc",index=i)[4]$percent[4:5])
}
}## 95% Bootstrap confidence intervals for glmmTMB
## Female : [1] 0.3332123 1.3936561
## Male : [1] 0.2983972 1.5842521
## Intercept : [1] 0.2686178 1.9475829
## CrossR/W : [1] -1.6585578 0.1240604
## CrossW/R : [1] -4.205659 -1.860600
## CrossW/W : [1] -1.145193 1.083512