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

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

#-------------------------------------------------------------
## 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)
summary(Exam3.2.glm)
#> 
#> Call:
#> glm(formula = F/N ~ trt, family = quasibinomial(link = "logit"), 
#>     data = DataSet3.1)
#> 
#> 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
#> 
library(parameters)
model_parameters(Exam3.2.glm)
#> Parameter   | Log-Odds |   SE |         95% CI |  t(18) |      p
#> ----------------------------------------------------------------
#> (Intercept) |    -2.25 | 0.22 | [-2.70, -1.84] | -10.37 | < .001
#> trt [1]     |     1.36 | 0.26 | [ 0.87,  1.89] |   5.29 | < .001
#> 
#> Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
#>   computed using a Wald t-distribution approximation.

#-------------------------------------------------------------
## Individula least squares treatment means
#-------------------------------------------------------------
library(emmeans)
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'
emmeans(object  = Exam3.2.glm, specs   = "trt")
#>  trt emmean    SE  df asymp.LCL asymp.UCL
#>  0   -2.247 0.217 Inf     -2.67    -1.822
#>  1   -0.883 0.140 Inf     -1.16    -0.608
#> 
#> Results are given on the logit (not the response) scale. 
#> Confidence level used: 0.95 
emmeans(object  = Exam3.2.glm, specs   = "trt", type = "response")
#>  trt   prob     SE  df asymp.LCL asymp.UCL
#>  0   0.0956 0.0187 Inf    0.0646     0.139
#>  1   0.2926 0.0290 Inf    0.2392     0.352
#> 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the logit scale 

#---------------------------------------------------
## Over all mean
#---------------------------------------------------
library(phia)
#> Loading required package: car
#> Loading required package: carData
list3.2 <-   list(trt = c("0" = 0.5, "1" = 0.5 ))
testFactors(model  =  Exam3.2.glm, levels =  list3.2 )
#> 
#> Call: base::tryCatch(model = base::withCallingHandlers({     NULL     base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp/RtmpST5oiP/callr-fun-96ab599ed3c1"),          base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv,          quote = TRUE), file = "/tmp/RtmpST5oiP/callr-res-96ab4ebdc98f",          compress = FALSE)     base::flush(base::stdout())     base::flush(base::stderr())     NULL     base::invisible() }, error = function(e) {     {         callr_data <- base::as.environment("tools:callr")$`__callr_data__`         err <- callr_data$err         if (FALSE) {             base::assign(".Traceback", base::.traceback(4), envir = callr_data)             utils::dump.frames("__callr_dump__")             base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`,                  envir = callr_data)             base::rm("__callr_dump__", envir = .GlobalEnv)         }         e <- err$process_call(e)         e2 <- err$new_error("error in callr subprocess")         class <- base::class         class(e2) <- base::c("callr_remote_error", class(e2))         e2 <- err$add_trace_back(e2)         cut <- base::which(e2$trace$scope == "global")[1]         if (!base::is.na(cut)) {             e2$trace <- e2$trace[-(1:cut), ]         }         base::saveRDS(base::list("error", e2, e), file = base::paste0("/tmp/RtmpST5oiP/callr-res-96ab4ebdc98f",              ".error"))     } }, interrupt = function(e) {     {         callr_data <- base::as.environment("tools:callr")$`__callr_data__`         err <- callr_data$err         if (FALSE) {             base::assign(".Traceback", base::.traceback(4), envir = callr_data)             utils::dump.frames("__callr_dump__")             base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`,                  envir = callr_data)             base::rm("__callr_dump__", envir = .GlobalEnv)         }         e <- err$process_call(e)         e2 <- err$new_error("error in callr subprocess")         class <- base::class         class(e2) <- base::c("callr_remote_error", class(e2))         e2 <- err$add_trace_back(e2)         cut <- base::which(e2$trace$scope == "global")[1]         if (!base::is.na(cut)) {             e2$trace <- e2$trace[-(1:cut), ]         }         base::saveRDS(base::list("error", e2, e), file = base::paste0("/tmp/RtmpST5oiP/callr-res-96ab4ebdc98f",              ".error"))     } }, callr_message = function(e) {     base::try(base::signalCondition(e)) }), error = function(e) {     NULL     if (FALSE) {         base::try(base::stop(e))     }     else {         base::invisible()     } }, interrupt = function(e) {     NULL     if (FALSE) {         e     }     else {         base::invisible()     } }) 
#> 
#> Term (Intercept) 
#> 
#> Adjusted mean:
#>           
#> 0.1729334 
#> 
#> Std. Error of link function:
#>           
#> 0.1290583 
#> 
#> 
#> Linear hypothesis test:
#> (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
#> ------

#---------------------------------------------------
## Repairwise treatment means estimate
#---------------------------------------------------
contrast(emmeans(object  = Exam3.2.glm, specs   = "trt"))
#>  contrast    estimate    SE  df z.ratio p.value
#>  trt0 effect   -0.682 0.129 Inf  -5.287  <.0001
#>  trt1 effect    0.682 0.129 Inf   5.287  <.0001
#> 
#> Results are given on the log odds ratio (not the response) scale. 
#> P value adjustment: fdr method for 2 tests 
contrast(emmeans(object  = Exam3.2.glm, specs   = "trt", type = "response"))
#>  contrast    odds.ratio     SE  df null z.ratio p.value
#>  trt0 effect      0.505 0.0652 Inf    1  -5.287  <.0001
#>  trt1 effect      1.978 0.2553 Inf    1   5.287  <.0001
#> 
#> P value adjustment: fdr method for 2 tests 
#> Tests are performed on the log odds ratio scale
Example 3.2 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-73)

References

See also

Author

Examples
