Example 3.3 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-77)
Source:R/Exam3.3.R
Exam3.3.RdExam3.3 use RCBD data with fixed location effect and different forms of estimable functions are shown in this example.
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 for Gaussian data
#-----------------------------------------------------------------------------------
data(DataSet3.2)
DataSet3.2$trt <- factor(x = DataSet3.2$trt, level = c(3,0,1,2))
DataSet3.2$loc <- factor(x = DataSet3.2$loc, level = c(8, 1, 2, 3, 4, 5, 6, 7))
levels(DataSet3.2$trt)
#> [1] "3" "0" "1" "2"
levels(DataSet3.2$loc)
#> [1] "8" "1" "2" "3" "4" "5" "6" "7"
Exam3.3.lm1 <- lm(formula = Y~ trt + loc, data = DataSet3.2)
summary( Exam3.3.lm1 )
#>
#> Call:
#> lm(formula = Y ~ trt + loc, data = DataSet3.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.7750 -0.7875 0.3625 1.0813 2.3000
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 28.2500 0.9945 28.407 < 2e-16 ***
#> trt0 -2.1375 0.8481 -2.520 0.01988 *
#> trt1 -2.9500 0.8481 -3.478 0.00224 **
#> trt2 0.2875 0.8481 0.339 0.73798
#> loc1 -0.9750 1.1994 -0.813 0.42538
#> loc2 -3.2500 1.1994 -2.710 0.01312 *
#> loc3 -1.6750 1.1994 -1.397 0.17713
#> loc4 -0.3500 1.1994 -0.292 0.77329
#> loc5 0.8250 1.1994 0.688 0.49907
#> loc6 -0.3000 1.1994 -0.250 0.80492
#> loc7 -3.5750 1.1994 -2.981 0.00713 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.696 on 21 degrees of freedom
#> Multiple R-squared: 0.6818, Adjusted R-squared: 0.5303
#> F-statistic: 4.5 on 10 and 21 DF, p-value: 0.001795
#>
#-------------------------------------------------------------
## Individual least squares treatment means
#-------------------------------------------------------------
library(emmeans)
(Lsm3.3 <- emmeans(object = Exam3.3.lm1, specs = ~trt))
#> trt emmean SE df lower.CL upper.CL
#> 3 27.1 0.6 21 25.8 28.3
#> 0 24.9 0.6 21 23.7 26.2
#> 1 24.1 0.6 21 22.9 25.4
#> 2 27.4 0.6 21 26.1 28.6
#>
#> Results are averaged over the levels of: loc
#> Confidence level used: 0.95
#---------------------------------------------------
## Pairwise treatment means estimate
#---------------------------------------------------
contrast(object = Lsm3.3 , method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> trt3 - trt0 2.138 0.848 21 2.520 0.0856
#> trt3 - trt1 2.950 0.848 21 3.478 0.0111
#> trt3 - trt2 -0.287 0.848 21 -0.339 0.9862
#> trt0 - trt1 0.812 0.848 21 0.958 0.7742
#> trt0 - trt2 -2.425 0.848 21 -2.859 0.0430
#> trt1 - trt2 -3.237 0.848 21 -3.817 0.0051
#>
#> Results are averaged over the levels of: loc
#> P value adjustment: tukey method for comparing a family of 4 estimates
#---------------------------------------------------
## Revpairwise treatment means estimate
#---------------------------------------------------
contrast(object = Lsm3.3, method = "revpairwise")
#> contrast estimate SE df t.ratio p.value
#> trt0 - trt3 -2.138 0.848 21 -2.520 0.0856
#> trt1 - trt3 -2.950 0.848 21 -3.478 0.0111
#> trt1 - trt0 -0.812 0.848 21 -0.958 0.7742
#> trt2 - trt3 0.287 0.848 21 0.339 0.9862
#> trt2 - trt0 2.425 0.848 21 2.859 0.0430
#> trt2 - trt1 3.237 0.848 21 3.817 0.0051
#>
#> Results are averaged over the levels of: loc
#> P value adjustment: tukey method for comparing a family of 4 estimates
#-------------------------------------------------------
## LSM Trt0 (This term is used in Walter Stroups' book)
#-------------------------------------------------------
contrast(
object = emmeans(object = Exam3.3.lm1, specs = ~ trt)
, list(trt = c(0, 1, 0, 0))
)
#> contrast estimate SE df t.ratio p.value
#> trt 24.9 0.6 21 41.605 <.0001
#>
#> Results are averaged over the levels of: loc
library(phia)
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1)))
#>
#> 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.95
#>
#> Std. Error:
#>
#> 0.5996899
#>
#>
#> Linear hypothesis test:
#> (Intercept) + trt0 + 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 ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 5040.437
#> 2 21 60.417 1 4980.02 1730.962 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## LSM Trt0 alt(This term is used in Walter Stroups' book)
#-------------------------------------------------------
# contrast(
# object = emmeans(object = Exam3.3.lm1, specs = ~ trt + loc)
# , list(
# trt = c(0, 1, 0, 0)
# , loc = c(1, 0, 0, 0, 0, 0, 0, 0)
# )
# )
#
#
# list3.3.2 <-
# list(
# trt = c("0" = 1 )
# , loc = c("1" = 0, "2" = 0,"3" = 0,"4" = 0,"5" = 0,"6" = 0,"7" = 0)
# )
# testFactors(model = Exam3.3.lm1, levels = list3.3.2)
#-------------------------------------------------------
## Trt0 Vs Trt1
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(0, 1, -1, 0))
)
#> contrast estimate SE df t.ratio p.value
#> trt 0.812 0.848 21 0.958 0.3489
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1, "1" = -1)))
#>
#> 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.8125
#>
#> Std. Error:
#>
#> 0.8480896
#>
#>
#> Linear hypothesis test:
#> trt0 - trt1 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 63.05812
#> 2 21 60.41750 1 2.640625 0.91783 0.34895
#> ------
#-------------------------------------------------------
## average Trt0 + Trt1
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(0, 1/2, 1/2, 0))
)
#> contrast estimate SE df t.ratio p.value
#> trt 24.5 0.424 21 57.880 <.0001
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5 , "1" = 0.5)))
#>
#> 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.4240448
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.5 trt0 + 0.5 trt1 + 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 ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 9698.748
#> 2 21 60.417 1 9638.331 3350.105 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## average Trt0+2+3
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(1/3, 1/3, 0, 1/3))
)
#> contrast estimate SE df t.ratio p.value
#> trt 26.5 0.346 21 76.454 <.0001
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"2" = 1/3,"3" = 1/3)))
#>
#> 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:
#>
#> 26.47083
#>
#> Std. Error:
#>
#> 0.3462311
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.333333333333333 trt0 + 0.333333333333333 trt2 + 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 ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 16877.338
#> 2 21 60.417 1 16816.92 5845.249 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## Trt 2 Vs 3 difference
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(-1, 0, 0, 1))
)
#> contrast estimate SE df t.ratio p.value
#> trt 0.287 0.848 21 0.339 0.7380
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("2" = 1,"3" = -1)))
#>
#> 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.2875
#>
#> Std. Error:
#>
#> 0.8480896
#>
#>
#> Linear hypothesis test:
#> trt2 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 60.74812
#> 2 21 60.41750 1 0.330625 0.11492 0.73798
#> ------
#-------------------------------------------------------
## Trt 1 Vs 2 difference
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(0, 0, 1, -1))
)
#> contrast estimate SE df t.ratio p.value
#> trt -3.24 0.848 21 -3.817 0.0010
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"2" = -1)))
#>
#> 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:
#>
#> -3.2375
#>
#> Std. Error:
#>
#> 0.8480896
#>
#>
#> Linear hypothesis test:
#> trt1 - trt2 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 102.3431
#> 2 21 60.4175 1 41.92562 14.57257 0.0010045 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## Trt 1 Vs 3 difference
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(-1, 0, 1, 0))
)
#> contrast estimate SE df t.ratio p.value
#> trt -2.95 0.848 21 -3.478 0.0022
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("1" = 1,"3" = -1)))
#>
#> 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:
#>
#> -2.95
#>
#> Std. Error:
#>
#> 0.8480896
#>
#>
#> Linear hypothesis test:
#> trt1 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 95.2275
#> 2 21 60.4175 1 34.81 12.09931 0.0022437 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## Average trt0+1 vs Average Trt2+3
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(-1/2, 1/2, 1/2, -1/2))
)
#> contrast estimate SE df t.ratio p.value
#> trt -2.69 0.6 21 -4.481 0.0002
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 0.5,"1" = 0.5,"2" = -0.5,"3" = -0.5)))
#>
#> 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:
#>
#> -2.6875
#>
#> Std. Error:
#>
#> 0.5996899
#>
#>
#> Linear hypothesis test:
#> 0.5 trt0 + 0.5 trt1 - 0.5 trt2 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 118.1988
#> 2 21 60.4175 1 57.78125 20.08369 0.00020551 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#-------------------------------------------------------
## Trt1 vs Average Trt0+1+2
#-------------------------------------------------------
contrast(
emmeans(object = Exam3.3.lm1, specs = ~trt)
, list(trt = c(1/3, 1/3, -1, 1/3))
)
#> contrast estimate SE df t.ratio p.value
#> trt 2.33 0.692 21 3.370 0.0029
#>
#> Results are averaged over the levels of: loc
testFactors(model = Exam3.3.lm1, levels = list(trt = c("0" = 1/3,"1" = -1,"2" = 1/3,"3" = 1/3)))
#>
#> 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:
#>
#> 2.333333
#>
#> Std. Error:
#>
#> 0.6924623
#>
#>
#> Linear hypothesis test:
#> - 2.08166817117217e - 17 ((Intercept) + 0.333333333333333 trt0 - trt1 + 0.333333333333333 trt2 - 6.93889390390723e - 18 loc1 - 6.93889390390723e - 18 loc2 - 6.93889390390723e - 18 loc3 - 6.93889390390723e - 18 loc4 - 6.93889390390723e - 18 loc5 - 6.93889390390723e - 18 loc6 - 6.93889390390723e - 18 loc7 = 0
#>
#> Model 1: restricted model
#> Model 2: Y ~ trt + loc
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 22 93.08417
#> 2 21 60.41750 1 32.66667 11.35433 0.0028974 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------