Problem 1: Random Cows

library(tidyverse)
## ── Attaching packages ─────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   1.0.0     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
  1. Why should we model sires as a random effect?

Sired is modeled as a random effect because the four bulls (labeled A, B, C, D) were randomly selected from a population of bulls and if this experiment were to be repeated, the next random sample would use different levels, different bulls from the overall population.

calves <- read.csv("calves.csv", header = TRUE)
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
## 7      B   1.17
## 8      B   1.08
## 9      B   1.20
## 10     B   1.08
## 11     B   1.01
## 12     B   0.86
## 13     C   0.98
## 14     C   1.06
## 15     C   1.15
## 16     C   1.11
## 17     C   0.83
## 18     C   0.86
## 19     D   0.95
## 20     D   1.10
## 21     D   1.07
## 22     D   1.11
## 23     D   0.89
## 24     D   1.12
  1. Write the factor effects model with a random effect for sire.

Yij = mu + ai + Eij

mu - overall mean weight

ai - random effect for ith sire

Eij - random error

  1. Fit a mixed model with a random intercept for sire using the lmer function in R.
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

With a p-value of 0.2741 (much greater than a = 0.05) being insignificant, it doesn’t seem worth it to include a random effect for sire.

  1. What are the variance component estimates for the model with a random effect for sire? (Model in part C)

variance estimate for sired (random effect) = 0.005078

variance estimate for the residual = 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