- Why should we model sires as a random effect?
# Sires is a random effect because the sires come from a larger population, and we wouldn’t use the same sires again if we did the study again.
- Write the factor effects model with a random effect for sire.
# Standard factors effects model --> Yij = mew + alphai + epsilonij
# Only one effect --> Yij = mew + alphai (w/ alphai being the random effect for sire)
- Fit a mixed model with a random intercept for sire using the lmer function in R.
# install.packages("lme4")
library(lme4)
## Loading required package: Matrix
weight<-c(1.46, 1.17, 0.98, 0.95, 1.23, 1.08, 1.06, 1.10, 1.12, 1.20, 1.15, 1.07, 1.23, 1.08, 1.11, 1.11, 1.02, 1.01, 0.83, 0.89, 1.15, 0.86, 0.86, 1.12)
sired<-c("A","B","C","D")
calves<-data.frame(weight,sired)
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
- Fit a linear model with only an intercept using the lm function.
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
- In order to test the null hypothesis that there is no sire variability in the response, compare the models in part C and D with the anova function. Report the results. Is it worth including a random effect for sire?
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
# Based on the anova results' AICs, it would appear that there may be a slight worth to including a random effoect for sire. However, the chisq value is not significant, which would suggest there is not a significant differenc in considering a random effect for sire.
- What are the variance component estimates for the model with a random effect for sire? (Model in part C)
#Variance:
# Sired 0.005078
# Residual: 0.07126
- Find 95% confidence intervals for the error standard deviation and the sire to sire standard deviation using the confint function.
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