Example 3.5 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-85)
Source:R/Exam3.5.R
Exam3.5.RdExam3.5 fixed location, factorial treatment structure, Gaussian response
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
data(DataSet3.2)
DataSet3.2$A <- factor(x = DataSet3.2$A)
DataSet3.2$B <- factor(x = DataSet3.2$B)
DataSet3.2$loc <- factor(x = DataSet3.2$loc, level = c(8, 1, 2, 3, 4, 5, 6, 7))
Exam3.5.lm <- lm(formula = Y~ A + B +loc, data = DataSet3.2)
Exam3.5.lm
#>
#> Call:
#> lm(formula = Y ~ A + B + loc, data = DataSet3.2)
#>
#> Coefficients:
#> (Intercept) A1 B1 loc1 loc2 loc3
#> 25.981 -0.550 2.688 -0.975 -3.250 -1.675
#> loc4 loc5 loc6 loc7
#> -0.350 0.825 -0.300 -3.575
#>
##---a0 marginal mean
library(emmeans)
contrast(
object = emmeans(object = Exam3.5.lm, specs = ~ B)
, list(trt = c(1, 0))
)
#> contrast estimate SE df t.ratio p.value
#> trt 24.5 0.416 22 58.974 <.0001
#>
#> Results are averaged over the levels of: A, loc
library(phia)
testFactors(model = Exam3.5.lm, levels = list(B = c("0" = 1,"1" = 0) ))
#>
#> 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:
#>
#> 24.54375
#>
#> Std. Error:
#>
#> 0.4161811
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.5 A1 + 0.125 loc1 + 0.125 loc2 + 0.125 loc3 + 0.125 loc4 + 0.125 loc5 + 0.125 loc6 + 0.125 loc7 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ A + B + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 23 9699.299
#> 2 22 60.969 1 9638.331 3477.901 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
##---b0 marginal mean
testFactors(model = Exam3.5.lm, levels=list(B = c("0" = 1, "1" = 0)))
#>
#> 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:
#>
#> 24.54375
#>
#> Std. Error:
#>
#> 0.4161811
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.5 A1 + 0.125 loc1 + 0.125 loc2 + 0.125 loc3 + 0.125 loc4 + 0.125 loc5 + 0.125 loc6 + 0.125 loc7 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ A + B + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 23 9699.299
#> 2 22 60.969 1 9638.331 3477.901 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
##---Simple effect of A on B0
testInteractions(model = Exam3.5.lm, custom = list(B = c("0" = 1,"1" = 0)), across = "B")
#> Warning: number of columns of result, 6, is not a multiple of vector length 5 of arg 2
#> F Test:
#> P-value adjustment method: holm
#> Value SE Df Sum of Sq F Pr(>F)
#> B1 -2.6875 0.5886 1.000 57.781 20.85 0.0001513 ***
#> Residuals 22.0000 60.969
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##---Simple effect of B on A0
testInteractions(model = Exam3.5.lm, custom = list(A = c("0" = 1, "1" = 0)), across = "A")
#> Warning: number of columns of result, 6, is not a multiple of vector length 5 of arg 2
#> F Test:
#> P-value adjustment method: holm
#> Value SE Df Sum of Sq F Pr(>F)
#> A1 0.55 0.5886 1.000 2.42 0.8732 0.3602
#> Residuals 22.0000 60.969
##---Simple Effect of A over B
testInteractions(model = Exam3.5.lm, fixed = "A", across = "B")
#> Warning: number of columns of result, 6, is not a multiple of vector length 5 of arg 2
#> F Test:
#> P-value adjustment method: holm
#> Value SE Df Sum of Sq F Pr(>F)
#> 0 -2.6875 0.5886 1.000 57.781 20.85 0.0003027 ***
#> 1 -2.6875 0.5886 1.000 57.781 20.85 0.0003027 ***
#> Residuals 22.0000 60.969
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##---Simple Effect of B over A
testInteractions(model = Exam3.5.lm, fixed = "B", across = "A")
#> Warning: number of columns of result, 6, is not a multiple of vector length 5 of arg 2
#> F Test:
#> P-value adjustment method: holm
#> Value SE Df Sum of Sq F Pr(>F)
#> 0 0.55 0.5886 1.000 2.42 0.8732 0.7204
#> 1 0.55 0.5886 1.000 2.42 0.8732 0.7204
#> Residuals 22.0000 60.969
#-------------------------------------------------------------
## Individula least squares treatment means
#-------------------------------------------------------------
emmeans(object = Exam3.5.lm, specs = ~A*B)
#> A B emmean SE df lower.CL upper.CL
#> 0 0 24.8 0.51 22 23.8 25.9
#> 1 0 24.3 0.51 22 23.2 25.3
#> 0 1 27.5 0.51 22 26.4 28.6
#> 1 1 27.0 0.51 22 25.9 28.0
#>
#> Results are averaged over the levels of: loc
#> Confidence level used: 0.95