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.
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
  # , weights
  # , control
    , 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 

#-------------------------------------------------------------
##  treatment means
#-------------------------------------------------------------
library(lsmeans)
(Lsm3.9fm1   <-
  lsmeans::lsmeans(
      object  = Exam3.9.fm1
    , specs   = "trt"
    , link=TRUE
    # , ...
  )
)
#>  trt     lsmean        SE df  lower.CL    upper.CL
#>  0   -1.7997839 0.4591969  7 -2.885612 -0.71395572
#>  1   -1.5174624 0.4381757  7 -2.553583 -0.48134141
#>  2   -0.9209468 0.4088300  7 -1.887676  0.04578253
#>  3   -0.7449613 0.4033992  7 -1.698849  0.20892630
#> 
#> 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
    , weights     = NULL
  # , subset
    , method      = "REML" #c("REML", "ML")
    , na.action   = na.fail
  # , control     = list()
    , contrasts   = NULL
    , keep.data   = TRUE
  )
(Lsm3.9fm2    <-
  lsmeans::lsmeans(
      object  = Exam3.9fm2
    , specs   = "trt"
    # , ...
  )
)
#>  trt    lsmean         SE df    lower.CL  upper.CL
#>  0   0.1604167 0.07277836  7 -0.01167680 0.3325101
#>  1   0.2000000 0.07277836  7  0.02790654 0.3720935
#>  2   0.3041667 0.07277836  7  0.13207320 0.4762601
#>  3   0.3395833 0.07277836  7  0.16748987 0.5116768
#> 
#> 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
  # , weights
  # , control
    , 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.768771e-05
#> 
#>  Formula: ~1 | loc %in% trt
#>          (Intercept)  Residual
#> StdDev: 2.946519e-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 
(Lsm3.9fm3    <-
  lsmeans::lsmeans(
      object  = Exam3.9fm3
    , specs   = "trt"
  # , ...
  )
)
#>  trt     lsmean        SE df  lower.CL    upper.CL
#>  0   -1.4485522 0.4379601  3 -2.842337 -0.05476769
#>  1   -1.2135963 0.4098640  0 -2.517966  0.09077395
#>  2   -0.7513588 0.3790090  0 -1.957535  0.45481710
#>  3   -0.7155082 0.3794512  0 -1.923091  0.49207499
#> 
#> 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)
  # , weights
  # , control
    , 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.00074658 0.4412173
#> 
#> Correlation Structure: Compound symmetry
#>  Formula: ~1 | loc 
#>  Parameter estimate(s):
#>       Rho 
#> 0.3817415 
#> 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 
(Lsm3.9fm4  <-
  lsmeans::lsmeans(
      object  = Exam3.9fm4
    , specs   = "trt"
    # , ...
  )
)
#>  trt     lsmean        SE df  lower.CL    upper.CL
#>  0   -1.6551311 0.4544041  7 -2.729626 -0.58063610
#>  1   -1.3862944 0.4169109  7 -2.372132 -0.40045667
#>  2   -0.8275344 0.3624894  7 -1.684685  0.02961675
#>  3   -0.6651516 0.3521450  7 -1.497842  0.16753898
#> 
#> 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

Examples

Contents

Author
