1. Why should we model sires as a random effect?

we treat sires as a random effect becasue we know that the possible mates for these bulls are part of a larger population, and were randomly selected as mates. This means if we were to perform the experiment again we would get different mates for these bulls.

  1. Write the factor effects model with a random effect for sire. (See slides for 11A) Hint: we are not interested in modeling trend over time, so there is only one effect in the model.

Yij= Mu + alphai +eij

Where Mu is the mean, alphai is the random effect and eij is the random error. In this model, the sires would be the alphai, or random effect.

  1. Fit a mixed model with a random intercept for sire using the lmer function in R.
library(Matrix)
library(lme4)

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

We should not include the random effect as there is no signifigant difference between the two models. Therefore, by adding sires we do not create a better model, nor do we better explain variance, instead we just complicate it.

  1. What are the variance component estimates for the model with a random effect for sire? (Model in part C)
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

We see that the variance component estimates for sired is 0.005078, whereas the residual variance is 0.015945

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

The 95 percent CI for the error standard deviation is [0.09534633, 0.1783855] and for sire to sire SD it is [0.00000000, 0.1795777]