6.5.2: Using R

Last updated
Save as PDF

Page ID: 33669

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

R Fully Nested Random Effects Model

Load the data.
Obtain the ANOVA for the fully nested random effects.

1. Load the data by using the following commands:

setwd("~/path-to-folder/")
fullnest_data <- read.table("fullnest_data.txt",header=T)
attach(fullnest_data)

2. Obtain the ANOVA for the fully nested random effects by using the following commands:

library(lmerTest)
library(lme4)
random_fullnest<-lmer(Temp ~ (1 | Plant) + (1 | Plant:Operator) +
(1 | Plant:(Operator:Shift)) ,fullnest_data)
summary(random_fullnest)
Linear mixed model fit by REML. t-tests use Satterthwaites method ['lmerModLmerTest']
Formula: Temp ~ (1 | Plant) + (1 | Plant:Operator) + (1 | Plant:(Operator:Shift))
   Data: fullnest_data
REML criterion at convergence: 1097.2

#Scaled residuals:
#     Min       1Q   Median       3Q      Max
#-2.78620 -0.61163  0.00414  0.56721  1.99397

#Random effects:
# Groups                 Name        Variance Std.Dev.
# Plant:(Operator:Shift) (Intercept)  6.5237  2.5542  
# Plant:Operator         (Intercept)  0.8061  0.8979  
# Plant                  (Intercept)  4.2123  2.0524  
# Residual                           12.4063  3.5223  
# Number of obs: 192, groups:  Plant:(Operator:Shift), 64; Plant:Operator, 16; Plant, 4

#Fixed effects:
#            Estimate Std. Error      df t value Pr(>|t|)    
#(Intercept)  474.880      1.127   3.000   421.4 2.95e-08 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
confint(random_fullnest)

 #                2.5 %     97.5 %
#.sig01        1.7251242   3.487550
#.sig02        0.0000000   2.475048
#.sig03        0.1192372   4.695585
#.sigma        3.1311707   4.002066
#(Intercept) 472.4015615 477.358858

Note that the command lmer() gives the ANOVA table only for the fixed effects. Therefore, in this example, since there are no fixed effects, we won’t get the ANOVA table. In the "Random effects" section of the output, under the column variance, we get the estimates for \(\sigma_{\gamma}^{2}\), \(\sigma_{\beta}^{2}\), \(\sigma_{\alpha}^{2}\), and \(\sigma^{2}\) which are equal to 6.5237, 0.8061, 4.2123, and 12.4063 respectively. In the "Fixed effects" section under the column estimate, we get the estimate of \(\mu\) for the overall mean, which is equal to 474.880.

With the command confint() we will get confidence intervals for the standard deviations and the overall mean. If you take the square of the lower and upper bounds, you will get a confidence interval for the model variances.

Alternatively, we can use the command aov() which gives a partial ANOVA table.

random_fullnest1<-aov(Temp ~ Error(factor(Plant) + factor(Plant)/factor(Operator) + factor(Plant)/(factor(Operator)/factor(Shift))) ,fullnest_data)
summary(random_fullnest1)

#Error: factor(Plant)
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals  3  731.5   243.8               

#Error: factor(Plant):factor(Operator)
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 12  499.8   41.65               

#Error: factor(Plant):factor(Operator):factor(Shift)
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 48   1535   31.98               

#Error: Within
#            Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 128   1588   12.41               
detach(fullnest_data)

Note that both commands in R don’t give the \(F\)-values and the \(p\)-values for the tests. Therefore, these must be done manually.