Example 9.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup (p-273)
Source:R/Exam9.1.R
Exam9.1.RdExam9.1 One-way random effects only 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.1)
DataSet9.1$a <- factor(x = DataSet9.1$a)
##---Random effects model
library(lmerTest)
Exam9.1lmer <- lmer( y ~ 1 + (1|a), data = DataSet9.1)
summary(Exam9.1lmer)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: y ~ 1 + (1 | a)
#> Data: DataSet9.1
#>
#> REML criterion at convergence: 101.4
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -1.37988 -0.56281 0.05619 0.52242 1.44495
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> a (Intercept) 5.511 2.348
#> Residual 1.531 1.237
#> Number of obs: 24, groups: a, 12
#>
#> Fixed effects:
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 15.9000 0.7232 11.0000 21.98 1.93e-10 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##---fixed effects model
Exam9.1lmer2 <- lm(y ~ a, data = DataSet9.1)
summary(Exam9.1lmer2)
#>
#> Call:
#> lm(formula = y ~ a, data = DataSet9.1)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.55 -0.85 0.00 0.85 1.55
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 16.1500 0.8749 18.460 3.55e-10 ***
#> a2 1.2000 1.2373 0.970 0.35125
#> a3 -4.0000 1.2373 -3.233 0.00718 **
#> a4 0.2000 1.2373 0.162 0.87427
#> a5 -1.7500 1.2373 -1.414 0.18266
#> a6 1.7000 1.2373 1.374 0.19456
#> a7 3.9000 1.2373 3.152 0.00834 **
#> a8 2.2500 1.2373 1.819 0.09401 .
#> a9 1.3500 1.2373 1.091 0.29665
#> a10 -3.5500 1.2373 -2.869 0.01411 *
#> a11 -3.2500 1.2373 -2.627 0.02211 *
#> a12 -1.0500 1.2373 -0.849 0.41269
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.237 on 12 degrees of freedom
#> Multiple R-squared: 0.8826, Adjusted R-squared: 0.775
#> F-statistic: 8.201 on 11 and 12 DF, p-value: 0.0005143
#>
#---------------------------------------------------
## Over all mean narrow( page # 274)
#---------------------------------------------------
library(emmeans)
library(phia)
list9.1 <- list(a = c( "1" = 1/12,"2" = 1/12
, "3" = 1/12,"4" = 1/12
, "5" = 1/12,"6" = 1/12
, "7" = 1/12,"8" = 1/12
, "9" = 1/12,"10" = 1/12
, "11" = 1/12,"12" = 1/12
))
phia::testFactors(model = Exam9.1lmer2, levels = list9.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:
#>
#> 15.9
#>
#> Std. Error:
#>
#> 0.2525564
#>
#>
#> Linear hypothesis test:
#> (Intercept) + 0.0833333333333333 a2 + 0.0833333333333333 a3 + 0.0833333333333333 a4 + 0.0833333333333333 a5 + 0.0833333333333333 a6 + 0.0833333333333333 a7 + 0.0833333333333333 a8 + 0.0833333333333333 a9 + 0.0833333333333333 a10 + 0.0833333333333333 a11 + 0.0833333333333333 a12 = 0
#>
#> Model 1: restricted model
#> Model 2: y ~ a
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 13 6085.81
#> 2 12 18.37 1 6067.44 3963.488 < 2.22e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> ------
#---BLUP Estimates (Table 9.1)
coef <- unlist(ranef(Exam9.1lmer))
BLUPa <- NULL
for( i in 1:length(coef)) {
BLUPa[i] <- (mean(DataSet9.1$y)+coef[i])
}
print(BLUPa)
#> [1] 16.11951 17.17318 12.60729 16.29513 14.58292 17.61221 19.54393 18.09514
#> [9] 17.30489 13.00241 13.26583 15.19755