Sires should be modeled as a random effect because they are a random sample from a larger population. The individual differences of the sires are not fixed. Instead they are random because they depend on which sample is pulled from the larger population.
yij = mu (overall mean) + alphai (random effect) + eij (random error) In this model, sires are modeled as alphai (random effect) and yij is the weight gain of the sires.
library(lme4)
## Loading required package: Matrix
Weight <- c(1.46, 1.23, 1.12, 1.23, 1.02, 1.15, 1.17, 1.08, 1.20, 1.08, 1.01, 0.86, 0.98, 1.06, 1.15, 1.11, 0.83, 0.86, 0.95, 1.10, 1.07, 1.11, 0.89, 1.12)
Sired <- c("A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D")
Calves <- data.frame(Sired, Weight)
head(Calves)
## Sired Weight
## 1 A 1.46
## 2 A 1.23
## 3 A 1.12
## 4 A 1.23
## 5 A 1.02
## 6 A 1.15
modC <- lmer(Weight ~ 1 + (1|Sired), data = Calves)
summary(modC)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Weight ~ 1 + (1 | Sired)
## Data: Calves
##
## REML criterion at convergence: -23.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6639 -0.5601 0.1081 0.5645 2.3859
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sired (Intercept) 0.005078 0.07126
## Residual 0.015945 0.12627
## Number of obs: 24, groups: Sired, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.07667 0.04397 24.48
modD <- lm(Weight ~ 1, data = Calves)
summary(modD)
##
## Call:
## lm(formula = Weight ~ 1, data = Calves)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.24667 -0.07417 0.01333 0.07333 0.38333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.07667 0.02881 37.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1411 on 23 degrees of freedom
anova(modC, modD)
## refitting model(s) with ML (instead of REML)
## Data: Calves
## Models:
## modD: Weight ~ 1
## modC: Weight ~ 1 + (1 | Sired)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modD 2 -22.898 -20.542 13.449 -26.898
## modC 3 -22.095 -18.561 14.047 -28.095 1.1962 1 0.2741
The results indicate that there is no significant difference in the model with the random effect vs. the one without it. Therefore, it is not worth including the random effect for sire, because adding complexity to the model does not improve it.
modC <- lmer(Weight ~ 1 + (1|Sired), data = Calves)
summary(modC)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Weight ~ 1 + (1 | Sired)
## Data: Calves
##
## REML criterion at convergence: -23.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6639 -0.5601 0.1081 0.5645 2.3859
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sired (Intercept) 0.005078 0.07126
## Residual 0.015945 0.12627
## Number of obs: 24, groups: Sired, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.07667 0.04397 24.48
It is estimated that 0.005078 of the variance is explained by intercept differences in sire (random effects), and 0.015945 is residual variance (variance unaccounted for).
confint(modC)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.00000000 0.1795777
## .sigma 0.09534633 0.1783855
## (Intercept) 0.97994356 1.1733897
The CI for the standard deviation of the error is CI95 [0.095, 0.178] and for the standard deviation sire to sire is CI95 [0.000, 0.180] (CI’s rounded to 100th place).