6.2 - Single Random Factor Model: Battery Life Example

Example
Example

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:

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)	.	.	.	.

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

^[1]	Link
↥	Has Tooltip/Popover
	Toggleable Visibility