6.2 - Single Random Factor Model: Battery Life Example

various types of batteries

Consider a study of Battery Life, measured in hours, where 4 brands of batteries are evaluated using 4 replications in a completely randomized design (Battery Data):

Brand A Brand B Brand C Brand D
110 118 108 117
113 116 107 112
108 112 112 115
115 117 108 119

A reasonable question to ask in this study would be, should the brand of the battery be considered a fixed effect or a random effect?

If the researchers were interested in comparing the performance of the specific brands they chose for the study, then we have a fixed effect. However, if the researchers were actually interested in studying the overall variation in battery life, so that the results would be applicable to all brands of batteries, then they may have chosen (presumably with a random sampling process) a sample of 4 of the many brands available. In this latter case, the battery brand would add a dimension of variability to battery life and can be considered a random effect.

Now, let us use SAS proc mixed to compare the results of battery brand as a fixed vs. random effect:

  1. Fixed Effect model:
     
    Type 3 Analysis of Variance
    Source DF Sum of Squares Mean Square Expected Mean Square Error Term Error DF F Value Pr > F
    Brand 3 141.687500 47.229167 Var(Residual) + Q(Brand) MS(Residual) 12 6.21 0.0086
    Residual 12 91.250000 7.604167 Var(Residual) . . . .
  2. Random Effect model:
     
    Type 3 Analysis of Variance
    Source DF Sum of Squares Mean Square Expected Mean Square Error Term Error DF F Value Pr > F
    Brand 3 141.687500 47.229167 Var(Residual) + 4 Var(Brand) MS(Residual) 12 6.21 0.0086
    Residual 12 91.250000 7.604167 Var(Residual) . . . .
     
    Covariance Parameter Estimates
    Cov Parm Estimate
    Brand 9.9063
    Residual 7.6042

We can verify the estimated variance component for the random treatment effect as:

\(s_{\text{among trts}}^{2}=\dfrac{MS_{trt}-MS_{error}}{n}=\dfrac{47.229-7.604}{4}=9.9063\)

With this, we can calculate the ICC as

\(ICC= \dfrac{9.9063}{9.9063+7.604}=0.5657\)

The key points in comparing these two ANOVAs are 1) the scope of inference and 2) the hypothesis being tested. For a fixed effect, the scope of inference is restricted to only 4 brands chosen for comparison and the null hypothesis is a statement of equality of means. With brand as a random effect, the scope of inference is the larger population of battery brands and the null hypothesis is a statement that the variance due to battery brand is 0.

First, load and attach the battery data in R.

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

We then apply the fixed effect model and obtain the ANOVA table.

options(contrasts=c("contr.sum","contr.poly"))
lm_fix = lm(lifetime ~ trt)
aov3_fix = car::Anova(lm_fix, type=3) 
aov3_fix
Anova Table (Type III tests)

Response: lifetime
            Sum Sq Df   F value    Pr(>F)     
(Intercept) 204078  1 26837.663 < 2.2e-16 ***
trt            142  3     6.211  0.008632 ** 
Residuals       91 12                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Random and mixed effect models can be applied in R using the lmer function from the lme4 package. We load the package and fit the model with a random treatment effect.

library(lme4)
lm_rand = lmer(lifetime ~ (1|trt))

We can obtain the ANOVA table for the fixed effects with the following command. The resulting table is blank since there are no fixed effects in this model.

anova(lm_rand)
Type III Analysis of Variance Table with Satterthwaite's method
     Sum Sq Mean Sq NumDF DenDF F value Pr(>F)

R does not compute F-tests for random effects. However, you can test the random effect components by using a Likelihood Ratio Test (LRT) via the ranova function from the lmerTest package

ranova(lm_rand)
ANOVA-like table for random-effects: Single term deletions
Model:
lifetime ~ (1 | trt)
               npar  logLik    AIC    LRT Df Pr(>Chisq)  
                  3 -40.625 87.250                       
(1 | trt)         2 -43.241 90.482 5.2314  1    0.02218 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

To extract the variance estimates from the model, we can do the following.

vc = VarCorr(lm_rand)
print(vc,comp=c("Variance"))
Groups   Name        Variance
trt      (Intercept) 9.9063  
Residual             7.6042  

As an alternative to testing the random effects using the LRT above, we can also examine the CIs of the variance estimates using the extractor function confint. It automatically computes the standard deviations, so we square and rename the outputs before displaying them. The resulting 95% CIs do not include zero, which is further evidence that the random effects are significant.

ci_sig = confint(lm_rand, parm="theta_", oldNames = FALSE)
ci_sig2 = ci_sig^2
row.names(ci_sig2) = sapply(row.names(ci_sig2),function(x){gsub("sd_","var_",x)},USE.NAMES=FALSE)
row.names(ci_sig2)[length(row.names(ci_sig2))] = "sigma2"
ci_sig2
                        2.5 %   97.5 %
var_(Intercept)|trt 0.4265072 51.36464
sigma2              3.7525970 19.13200