Example 9.2 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup (p-276)
Source:R/Exam9.2.R
Exam9.2.RdExam9.2 Two way random effects nested model
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(DataSet9.2)
DataSet9.2$a <- factor(x = DataSet9.2$a)
DataSet9.2$b <- factor(x = DataSet9.2$b)
library(lmerTest)
Exam9.2lmer <- lmer(y ~ (1|b/a), data = DataSet9.2)
summary(Exam9.2lmer)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: y ~ (1 | b/a)
#> Data: DataSet9.2
#>
#> REML criterion at convergence: 109.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -1.2000 -0.4615 -0.0562 0.4228 2.4853
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> a:b (Intercept) 2.6008 1.6127
#> b (Intercept) 0.3577 0.5981
#> Residual 1.3614 1.1668
#> Number of obs: 28, groups: a:b, 14; b, 2
#>
#> Fixed effects:
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 19.2571 0.6429 1.0000 29.96 0.0212 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Exam9.2lmer2 <- lm(y ~ a + b %in% a, data = DataSet9.2)
summary(Exam9.2lmer2)
#>
#> Call:
#> lm(formula = y ~ a + b %in% a, data = DataSet9.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.1500 -0.4625 0.0000 0.4625 2.1500
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.705e+01 8.251e-01 20.665 6.9e-12 ***
#> a2 1.650e+00 1.167e+00 1.414 0.17918
#> a3 1.500e+00 1.167e+00 1.286 0.21945
#> a4 3.000e-01 1.167e+00 0.257 0.80083
#> a5 3.000e+00 1.167e+00 2.571 0.02219 *
#> a6 4.400e+00 1.167e+00 3.771 0.00207 **
#> a7 1.000e-01 1.167e+00 0.086 0.93292
#> a1:b2 3.000e-01 1.167e+00 0.257 0.80083
#> a2:b2 1.350e+00 1.167e+00 1.157 0.26663
#> a3:b2 1.500e+00 1.167e+00 1.286 0.21945
#> a4:b2 3.600e+00 1.167e+00 3.085 0.00806 **
#> a5:b2 -1.721e-15 1.167e+00 0.000 1.00000
#> a6:b2 1.700e+00 1.167e+00 1.457 0.16719
#> a7:b2 5.500e-01 1.167e+00 0.471 0.64464
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.167 on 14 degrees of freedom
#> Multiple R-squared: 0.8258, Adjusted R-squared: 0.664
#> F-statistic: 5.104 on 13 and 14 DF, p-value: 0.00233
#>
##--- Over all mean
library(phia)
list9.2 <- list(a = c("1" = 1/7,"2" = 1/7
, "3" = 1/7,"4" = 1/7
, "5" = 1/7,"6" = 1/7
, "7" = 1/7
))
phia::testFactors(model = Exam9.2lmer2, levels = list9.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:
#>
#> 19.25714
#>
#> Std. Error:
#>
#> 0.220505
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.142857142857143 a2 + 0.142857142857143 a3 + 0.142857142857143 a4 + 0.142857142857143 a5 + 0.142857142857143 a6 + 0.142857142857143 a7 + 0.0714285714285714 a1:b2 + 0.0714285714285714 a2:b2 + 0.0714285714285714 a3:b2 + 0.0714285714285714 a4:b2 + 0.0714285714285714 a5:b2 + 0.0714285714285714 a6:b2 + 0.0714285714285714 a7:b2 = 0
#>
#> Model 1: restricted model
#> Model 2: y ~ a + b %in% a
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 15 10402.51
#> 2 14 19.06 1 10383.45 7626.879 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#---BLUP Estimates
coef <- unlist(ranef(Exam9.2lmer)$a)
BLUPa <- NULL
for(i in 1:length(coef)){
BLUPa[i] <- (mean(DataSet9.2$y) + coef[i])
}
print(BLUPa)
#> [1] 17.72837 17.52509 19.03610 19.66501 18.91721 19.66501 17.96614 20.37832
#> [9] 20.10606 19.66501 21.21564 22.12195 17.80762 17.80249
#---BLUP Estimates Narrow
BLUPaNar <- NULL
for( i in 1:length(coef)) {
BLUPaNar[i] <- (mean(DataSet9.2$y) + coef[i])
}
BLUPaNar
#> [1] 17.72837 17.52509 19.03610 19.66501 18.91721 19.66501 17.96614 20.37832
#> [9] 20.10606 19.66501 21.21564 22.12195 17.80762 17.80249