Example 3.9 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-118)

Exam3.9 used to differentiate conditional and marginal binomial models with and without interaction for S2 variable.
References

Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
Author

Muhammad Yaseen (myaseen208@gmail.com)
Adeela Munawar (adeela.uaf@gmail.com)
Examples

#-----------------------------------------------------------------------------------
## Binomial conditional GLMM without interaction, logit link
#-----------------------------------------------------------------------------------
library(MASS)
DataSet3.2$trt <- factor( x  =  DataSet3.2$trt )
DataSet3.2$loc <- factor( x  =  DataSet3.2$loc )

Exam3.9.fm1   <-
  glmmPQL(
      fixed    =  S2/Nbin~trt
    , random   = ~1|loc
    , family   =  quasibinomial(link = "logit")
    , data     =  DataSet3.2
    , niter    = 10
    , verbose  = TRUE
  )
#> iteration 1
#> iteration 2
#> iteration 3
#> iteration 4
summary(Exam3.9.fm1)
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: DataSet3.2 
#>   AIC BIC logLik
#>    NA  NA     NA
#> 
#> Random effects:
#>  Formula: ~1 | loc
#>         (Intercept) Residual
#> StdDev:   0.7710488  0.33184
#> 
#> Variance function:
#>  Structure: fixed weights
#>  Formula: ~invwt 
#> Fixed effects:  S2/Nbin ~ trt 
#>                  Value Std.Error DF   t-value p-value
#> (Intercept) -1.7997839 0.4591969 21 -3.919416  0.0008
#> trt1         0.2823215 0.4752384 21  0.594063  0.5588
#> trt2         0.8788371 0.4502862 21  1.951730  0.0644
#> trt3         1.0548226 0.4461060 21  2.364511  0.0278
#>  Correlation: 
#>      (Intr) trt1   trt2  
#> trt1 -0.561              
#> trt2 -0.596  0.571       
#> trt3 -0.603  0.576  0.610
#> 
#> Standardized Within-Group Residuals:
#>        Min         Q1        Med         Q3        Max 
#> -2.0036140 -0.5699977 -0.2243692  0.4913423  2.4790311 
#> 
#> Number of Observations: 32
#> Number of Groups: 8 
library(parameters)
model_parameters(Exam3.9.fm1)
#> # Fixed Effects
#> 
#> Parameter   | Log-Odds |   SE |         95% CI | t(21) |      p
#> ---------------------------------------------------------------
#> (Intercept) |    -1.80 | 0.46 | [-2.69, -0.91] | -3.92 | < .001
#> trt [1]     |     0.28 | 0.48 | [-0.64,  1.21] |  0.59 | 0.559 
#> trt [2]     |     0.88 | 0.45 | [ 0.00,  1.75] |  1.95 | 0.064 
#> trt [3]     |     1.05 | 0.45 | [ 0.19,  1.92] |  2.36 | 0.028 
#> 
#> # Random Effects
#> 
#> Parameter           | Coefficient
#> ---------------------------------
#> SD (Intercept: loc) |        0.77
#> SD (Residual)       |        0.33
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.

#-------------------------------------------------------------
##  treatment means
#-------------------------------------------------------------
library(emmeans)
emmeans(object = Exam3.9.fm1, specs = ~trt, type = "response")
#>  trt  prob     SE df lower.CL upper.CL
#>  0   0.142 0.0559  7   0.0529    0.329
#>  1   0.180 0.0646  7   0.0722    0.382
#>  2   0.285 0.0833  7   0.1315    0.511
#>  3   0.322 0.0881  7   0.1546    0.552
#> 
#> Degrees-of-freedom method: containment 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the logit scale 
emmeans(object = Exam3.9.fm1, specs = ~trt, type = "link")
#>  trt emmean    SE df lower.CL upper.CL
#>  0   -1.800 0.459  7    -2.89  -0.7140
#>  1   -1.517 0.438  7    -2.55  -0.4813
#>  2   -0.921 0.409  7    -1.89   0.0458
#>  3   -0.745 0.403  7    -1.70   0.2089
#> 
#> Degrees-of-freedom method: containment 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95 
emmeans(object = Exam3.9.fm1, specs = ~trt, type = "logit")
#>  trt emmean    SE df lower.CL upper.CL
#>  0   -1.800 0.459  7    -2.89  -0.7140
#>  1   -1.517 0.438  7    -2.55  -0.4813
#>  2   -0.921 0.409  7    -1.89   0.0458
#>  3   -0.745 0.403  7    -1.70   0.2089
#> 
#> Degrees-of-freedom method: containment 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95 

##--- Normal Approximation
library(nlme)
Exam3.9fm2 <-
  lme(
      fixed       = S2/Nbin~trt
    , data        = DataSet3.2
    , random      = ~1|loc
    , method      = c("REML", "ML")[1]
  )

Exam3.9fm2
#> Linear mixed-effects model fit by REML
#>   Data: DataSet3.2 
#>   Log-restricted-likelihood: 2.381826
#>   Fixed: S2/Nbin ~ trt 
#> (Intercept)        trt1        trt2        trt3 
#>  0.16041667  0.03958333  0.14375000  0.17916667 
#> 
#> Random effects:
#>  Formula: ~1 | loc
#>         (Intercept)  Residual
#> StdDev:   0.1218281 0.1659259
#> 
#> Number of Observations: 32
#> Number of Groups: 8 
model_parameters(Exam3.9fm2)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |        95% CI | t(21) |     p
#> ----------------------------------------------------------------
#> (Intercept) |        0.16 | 0.07 | [ 0.01, 0.31] |  2.20 | 0.039
#> trt [1]     |        0.04 | 0.08 | [-0.13, 0.21] |  0.48 | 0.638
#> trt [2]     |        0.14 | 0.08 | [-0.03, 0.32] |  1.73 | 0.098
#> trt [3]     |        0.18 | 0.08 | [ 0.01, 0.35] |  2.16 | 0.043
#> 
#> # Random Effects
#> 
#> Parameter           | Coefficient
#> ---------------------------------
#> SD (Intercept: loc) |        0.12
#> SD (Residual)       |        0.17
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.

emmeans(object  = Exam3.9fm2, specs = ~trt)
#>  trt emmean     SE df lower.CL upper.CL
#>  0    0.160 0.0728  7  -0.0117    0.333
#>  1    0.200 0.0728  7   0.0279    0.372
#>  2    0.304 0.0728  7   0.1321    0.476
#>  3    0.340 0.0728  7   0.1675    0.512
#> 
#> Degrees-of-freedom method: containment 
#> Confidence level used: 0.95 


##---Binomial GLMM with interaction
Exam3.9fm3   <-
  glmmPQL(
      fixed       =  S2/Nbin~trt
    , random      = ~1|trt/loc
    , family      =  quasibinomial(link = "logit")
    , data        =  DataSet3.2
    , niter = 10
    , verbose = TRUE
  )
#> iteration 1
#> iteration 2
#> iteration 3
#> iteration 4
#> iteration 5
#> iteration 6
#> iteration 7
#> iteration 8
#> iteration 9
#> iteration 10
summary(Exam3.9fm3)
#> Warning: NaNs produced
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: DataSet3.2 
#>   AIC BIC logLik
#>    NA  NA     NA
#> 
#> Random effects:
#>  Formula: ~1 | trt
#>         (Intercept)
#> StdDev: 1.76876e-05
#> 
#>  Formula: ~1 | loc %in% trt
#>          (Intercept)  Residual
#> StdDev: 2.946513e-06 0.4265118
#> 
#> Variance function:
#>  Structure: fixed weights
#>  Formula: ~invwt 
#> Fixed effects:  S2/Nbin ~ trt 
#>                  Value Std.Error DF   t-value p-value
#> (Intercept) -1.4485522 0.4379601 28 -3.307498  0.0026
#> trt1         0.2349559 0.5998313  0  0.391703     NaN
#> trt2         0.6971934 0.5791864  0  1.203746     NaN
#> trt3         0.7330440 0.5794759  0  1.265012     NaN
#>  Correlation: 
#>      (Intr) trt1   trt2  
#> trt1 -0.730              
#> trt2 -0.756  0.552       
#> trt3 -0.756  0.552  0.571
#> 
#> Standardized Within-Group Residuals:
#>        Min         Q1        Med         Q3        Max 
#> -1.8229824 -0.8793369 -0.4149872  0.5574551  2.2264946 
#> 
#> Number of Observations: 32
#> Number of Groups: 
#>          trt loc %in% trt 
#>            4           32 
model_parameters(Exam3.9fm3)
#> Warning: NaNs produced
#> Warning: NaNs produced
#> Warning: NaNs produced
#> # Fixed Effects
#> 
#> Parameter   | Log-Odds |   SE |         95% CI |     t | df |     p
#> -------------------------------------------------------------------
#> (Intercept) |    -1.45 | 0.44 | [-2.29, -0.61] | -3.31 | 28 | 0.003
#> trt [1]     |     0.23 | 0.60 |                |  0.39 |  0 |      
#> trt [2]     |     0.70 | 0.58 |                |  1.20 |  0 |      
#> trt [3]     |     0.73 | 0.58 |                |  1.27 |  0 |      
#> 
#> # Random Effects
#> 
#> Parameter           | Coefficient
#> ---------------------------------
#> SD (Intercept: trt) |    1.77e-05
#> SD (Intercept: loc) |    2.95e-06
#> SD (Residual)       |        0.43
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.
emmeans(object = Exam3.9fm3, specs = ~trt)
#> Warning: NaNs produced
#>  trt emmean    SE df lower.CL upper.CL
#>  0   -1.449 0.438  3    -2.84  -0.0548
#>  1   -1.214 0.410  0      NaN      NaN
#>  2   -0.751 0.379  0      NaN      NaN
#>  3   -0.716 0.379  0      NaN      NaN
#> 
#> Degrees-of-freedom method: containment 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95 


##---Binomial Marginal GLMM(assuming compound symmetry)
Exam3.9fm4   <-
  glmmPQL(
      fixed       =  S2/Nbin~trt
    , random      = ~1|loc
    , family      =  quasibinomial(link = "logit")
    , data        =  DataSet3.2
    , correlation =  corCompSymm(form = ~1|loc)
    , niter       = 10
    , verbose     = TRUE
  )
#> iteration 1
summary(Exam3.9fm4)
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: DataSet3.2 
#>   AIC BIC logLik
#>    NA  NA     NA
#> 
#> Random effects:
#>  Formula: ~1 | loc
#>          (Intercept)  Residual
#> StdDev: 0.0007465953 0.4412173
#> 
#> Correlation Structure: Compound symmetry
#>  Formula: ~1 | loc 
#>  Parameter estimate(s):
#>       Rho 
#> 0.3817414 
#> Variance function:
#>  Structure: fixed weights
#>  Formula: ~invwt 
#> Fixed effects:  S2/Nbin ~ trt 
#>                  Value Std.Error DF   t-value p-value
#> (Intercept) -1.6551311 0.4544041 21 -3.642421  0.0015
#> trt1         0.2688368 0.4854471 21  0.553792  0.5856
#> trt2         0.8275968 0.4605683 21  1.796903  0.0867
#> trt3         0.9899796 0.4564203 21  2.169008  0.0417
#>  Correlation: 
#>      (Intr) trt1   trt2  
#> trt1 -0.608              
#> trt2 -0.686  0.577       
#> trt3 -0.701  0.583  0.624
#> 
#> Standardized Within-Group Residuals:
#>        Min         Q1        Med         Q3        Max 
#> -1.6252165 -0.6846589 -0.3489664  0.6677812  2.6023420 
#> 
#> Number of Observations: 32
#> Number of Groups: 8 
model_parameters(Exam3.9fm4)
#> # Fixed Effects
#> 
#> Parameter   | Log-Odds |   SE |         95% CI | t(21) |     p
#> --------------------------------------------------------------
#> (Intercept) |    -1.66 | 0.45 | [-2.54, -0.77] | -3.64 | 0.002
#> trt [1]     |     0.27 | 0.49 | [-0.68,  1.21] |  0.55 | 0.586
#> trt [2]     |     0.83 | 0.46 | [-0.07,  1.72] |  1.80 | 0.087
#> trt [3]     |     0.99 | 0.46 | [ 0.10,  1.88] |  2.17 | 0.042
#> 
#> # Random Effects
#> 
#> Parameter           | Coefficient
#> ---------------------------------
#> SD (Intercept: loc) |    7.47e-04
#> SD (Residual)       |        0.44
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.
emmeans(object  = Exam3.9fm4, specs  = ~trt)
#>  trt emmean    SE df lower.CL upper.CL
#>  0   -1.655 0.454  7    -2.73  -0.5806
#>  1   -1.386 0.417  7    -2.37  -0.4005
#>  2   -0.828 0.362  7    -1.68   0.0296
#>  3   -0.665 0.352  7    -1.50   0.1675
#> 
#> Degrees-of-freedom method: containment 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95
Example 3.9 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-118)

References

See also

Author

Examples
