# Load the necessary packages
library(MASS)
library(lme4)
## Warning: package 'lme4' was built under R version 4.2.3
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.2.3
# Load the bacteria dataset
data("bacteria")
str(bacteria)
## 'data.frame': 220 obs. of 6 variables:
## $ y : Factor w/ 2 levels "n","y": 2 2 2 2 2 2 1 2 2 2 ...
## $ ap : Factor w/ 2 levels "a","p": 2 2 2 2 1 1 1 1 1 1 ...
## $ hilo: Factor w/ 2 levels "hi","lo": 1 1 1 1 1 1 1 1 2 2 ...
## $ week: int 0 2 4 11 0 2 6 11 0 2 ...
## $ ID : Factor w/ 50 levels "X01","X02","X03",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ trt : Factor w/ 3 levels "placebo","drug",..: 1 1 1 1 3 3 3 3 2 2 ...
summary(bacteria)
## y ap hilo week ID trt
## n: 43 a:124 hi:122 Min. : 0.000 X03 : 5 placebo:96
## y:177 p: 96 lo: 98 1st Qu.: 2.000 X04 : 5 drug :62
## Median : 4.000 X05 : 5 drug+ :62
## Mean : 4.455 X07 : 5
## 3rd Qu.: 6.000 X08 : 5
## Max. :11.000 X09 : 5
## (Other):190
# Fit a GLM using logistic regression
glm_model <- glm(y ~ trt + week, data = bacteria, family = binomial)
summary(glm_model)
##
## Call:
## glm(formula = y ~ trt + week, family = binomial, data = bacteria)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2899 0.3885 0.5400 0.7027 1.1077
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.54629 0.40555 6.279 3.42e-10 ***
## trtdrug -1.10667 0.42519 -2.603 0.00925 **
## trtdrug+ -0.65166 0.44615 -1.461 0.14412
## week -0.11577 0.04414 -2.623 0.00872 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 217.38 on 219 degrees of freedom
## Residual deviance: 203.81 on 216 degrees of freedom
## AIC: 211.81
##
## Number of Fisher Scoring iterations: 4
# Fit a GLMM using the glmer function from the lme4 package
glmm_model <- glmer(y ~ trt + week + (1 | ID), data = bacteria, family = binomial)
summary(glmm_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.8 224.7 -98.9 197.8 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8175 0.1755 0.2958 0.4171 1.2930
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.314 1.146
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.14392 0.62249 5.051 4.41e-07 ***
## trtdrug -1.32014 0.64240 -2.055 0.03988 *
## trtdrug+ -0.79544 0.65198 -1.220 0.22245
## week -0.14369 0.05099 -2.818 0.00484 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.618
## trtdrug+ -0.580 0.492
## week -0.591 0.118 0.085
# Default method in glmer is adaptive Gauss-Hermite quadrature
glmm_quad_model <- glmer(y ~ trt + week + (1 | ID), data = bacteria, family = binomial, nAGQ = 10)
summary(glmm_quad_model)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.4 224.4 -98.7 197.4 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8008 0.1717 0.2884 0.4129 1.3316
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.446 1.202
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.16561 0.62870 5.035 4.78e-07 ***
## trtdrug -1.32455 0.65735 -2.015 0.0439 *
## trtdrug+ -0.80488 0.66745 -1.206 0.2279
## week -0.14553 0.05136 -2.834 0.0046 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.615
## trtdrug+ -0.582 0.491
## week -0.587 0.112 0.081
# Laplace approximation method in glmer
glmm_laplace_model <- glmer(y ~ trt + week + (1 | ID), data = bacteria, family = binomial, nAGQ = 1)
summary(glmm_laplace_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.8 224.7 -98.9 197.8 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8175 0.1755 0.2958 0.4171 1.2930
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.314 1.146
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.14392 0.62249 5.051 4.41e-07 ***
## trtdrug -1.32014 0.64240 -2.055 0.03988 *
## trtdrug+ -0.79544 0.65198 -1.220 0.22245
## week -0.14369 0.05099 -2.818 0.00484 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.618
## trtdrug+ -0.580 0.492
## week -0.591 0.118 0.085
# Fit a GLMM using the glmmPQL function from the MASS package
glmmPQL_model <- glmmPQL(y ~ trt + week, random = ~ 1 | ID, family = binomial, data = bacteria)
## iteration 1
## iteration 2
## iteration 3
## iteration 4
## iteration 5
## iteration 6
summary(glmmPQL_model)
## Linear mixed-effects model fit by maximum likelihood
## Data: bacteria
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 1.325243 0.7903088
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: y ~ trt + week
## Value Std.Error DF t-value p-value
## (Intercept) 3.0302276 0.4791396 169 6.324310 0.0000
## trtdrug -1.2176812 0.6160113 47 -1.976719 0.0540
## trtdrug+ -0.7886376 0.6193895 47 -1.273250 0.2092
## week -0.1446463 0.0392343 169 -3.686730 0.0003
## Correlation:
## (Intr) trtdrg trtdr+
## trtdrug -0.622
## trtdrug+ -0.609 0.464
## week -0.481 0.050 0.030
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.2868074 0.2039043 0.3140333 0.5440835 1.9754065
##
## Number of Observations: 220
## Number of Groups: 50
# Print summaries to compare
summary(glm_model)
##
## Call:
## glm(formula = y ~ trt + week, family = binomial, data = bacteria)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2899 0.3885 0.5400 0.7027 1.1077
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.54629 0.40555 6.279 3.42e-10 ***
## trtdrug -1.10667 0.42519 -2.603 0.00925 **
## trtdrug+ -0.65166 0.44615 -1.461 0.14412
## week -0.11577 0.04414 -2.623 0.00872 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 217.38 on 219 degrees of freedom
## Residual deviance: 203.81 on 216 degrees of freedom
## AIC: 211.81
##
## Number of Fisher Scoring iterations: 4
summary(glmm_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.8 224.7 -98.9 197.8 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8175 0.1755 0.2958 0.4171 1.2930
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.314 1.146
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.14392 0.62249 5.051 4.41e-07 ***
## trtdrug -1.32014 0.64240 -2.055 0.03988 *
## trtdrug+ -0.79544 0.65198 -1.220 0.22245
## week -0.14369 0.05099 -2.818 0.00484 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.618
## trtdrug+ -0.580 0.492
## week -0.591 0.118 0.085
summary(glmm_quad_model)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.4 224.4 -98.7 197.4 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8008 0.1717 0.2884 0.4129 1.3316
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.446 1.202
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.16561 0.62870 5.035 4.78e-07 ***
## trtdrug -1.32455 0.65735 -2.015 0.0439 *
## trtdrug+ -0.80488 0.66745 -1.206 0.2279
## week -0.14553 0.05136 -2.834 0.0046 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.615
## trtdrug+ -0.582 0.491
## week -0.587 0.112 0.081
summary(glmm_laplace_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ trt + week + (1 | ID)
## Data: bacteria
##
## AIC BIC logLik deviance df.resid
## 207.8 224.7 -98.9 197.8 215
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8175 0.1755 0.2958 0.4171 1.2930
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.314 1.146
## Number of obs: 220, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.14392 0.62249 5.051 4.41e-07 ***
## trtdrug -1.32014 0.64240 -2.055 0.03988 *
## trtdrug+ -0.79544 0.65198 -1.220 0.22245
## week -0.14369 0.05099 -2.818 0.00484 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtdrg trtdr+
## trtdrug -0.618
## trtdrug+ -0.580 0.492
## week -0.591 0.118 0.085
summary(glmmPQL_model)
## Linear mixed-effects model fit by maximum likelihood
## Data: bacteria
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 1.325243 0.7903088
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: y ~ trt + week
## Value Std.Error DF t-value p-value
## (Intercept) 3.0302276 0.4791396 169 6.324310 0.0000
## trtdrug -1.2176812 0.6160113 47 -1.976719 0.0540
## trtdrug+ -0.7886376 0.6193895 47 -1.273250 0.2092
## week -0.1446463 0.0392343 169 -3.686730 0.0003
## Correlation:
## (Intr) trtdrg trtdr+
## trtdrug -0.622
## trtdrug+ -0.609 0.464
## week -0.481 0.050 0.030
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.2868074 0.2039043 0.3140333 0.5440835 1.9754065
##
## Number of Observations: 220
## Number of Groups: 50