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

Exam3.2 used binomial data, two treatment samples

References

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

See also

DataSet3.1

Examples

#-------------------------------------------------------------
## Linear Model and results discussed in Article 1.2.1 after Table1.1
#-------------------------------------------------------------
data(DataSet3.1)
DataSet3.1$trt <- factor(x =  DataSet3.1$trt)
Exam3.2.glm <-
  glm(
         formula    =  F/N~trt
       , family     =  quasibinomial(link = "logit")
       , data       =  DataSet3.1
       , weights    =  NULL
    #  , subset
    #  , na.action
       , start      =  NULL
    #  , etastart
    #  , mustart
    #  , offset
    #  , control    =  list(...)
    #  , model      =  TRUE
       , method     =  "glm.fit"
    #  , x          =  FALSE
    #  , y          =  TRUE
       , contrasts  =  NULL
    #  , ...
  )
summary(Exam3.2.glm)
#> 
#> Call:
#> glm(formula = F/N ~ trt, family = quasibinomial(link = "logit"), 
#>     data = DataSet3.1, weights = NULL, start = NULL, method = "glm.fit", 
#>     contrasts = NULL)
#> 
#> Deviance Residuals: 
#>      Min        1Q    Median        3Q       Max  
#> -0.44825  -0.16321   0.07221   0.09450   0.25513  
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  -2.2473     0.2168 -10.366 5.12e-09 ***
#> trt1          1.3646     0.2581   5.287 5.01e-05 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for quasibinomial family taken to be 0.04062774)
#> 
#>     Null deviance: 2.21652  on 19  degrees of freedom
#> Residual deviance: 0.92684  on 18  degrees of freedom
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 5
#> 

#-------------------------------------------------------------
## Individula least squares treatment means
#-------------------------------------------------------------
library(lsmeans)
(Lsm3.2    <-
  lsmeans::lsmeans(
    object  = Exam3.2.glm
    , specs   = "trt"
    # , ...
  )
)
#>  trt     lsmean        SE df asymp.LCL  asymp.UCL
#>  0   -2.2472830 0.2167875 NA -2.672179 -1.8223873
#>  1   -0.8826749 0.1400978 NA -1.157262 -0.6080883
#> 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95 
OddsRatioMean3.2 <-  1/(1 + exp(-summary(Lsm3.2)[c("lsmean")] ) )
#---------------------------------------------------
## Over all mean
#---------------------------------------------------
library(phia)
list3.2<-   list(trt=c("0" = 0.5,"1" = 0.5 ))
(Test3.2 <-
  testFactors(
      model  =  Exam3.2.glm
    , levels =  list3.2 )
)
#> 
#> Call: pkgdown::build_site(pkg = ".", examples = TRUE, run_dont_run = FALSE,      seed = 12345, mathjax = TRUE, lazy = FALSE, override = list(),      preview = NA) 
#> 
#> Term (Intercept) 
#> 
#> Adjusted mean:
#>           
#> 0.1729334 
#> 
#> Std. Error of link function:
#>           
#> 0.1290583 
#> 
#> Linear hypothesis test
#> 
#> Hypothesis:
#> (Intercept)  + 0.5 trt1 = 0
#> 
#> Model 1: restricted model
#> Model 2: F/N ~ trt
#> 
#>   Res.Df Df    Chisq Pr(>Chisq)    
#> 1     19                           
#> 2     18  1 147.0432 < 2.22e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#---------------------------------------------------
## Pairwise treatment means estimate
#---------------------------------------------------
contrast(object = Lsm3.2 , method = "pairwise")
#>  contrast  estimate        SE df z.ratio p.value
#>  0 - 1    -1.364608 0.2581167 NA  -5.287  <.0001
#> 
#> Results are given on the log odds ratio (not the response) scale. 
#---------------------------------------------------
## Repairwise treatment means estimate
#---------------------------------------------------
## contrast( object = Lsm3.2 , method = "repairwise")