The sires should be labeled as random effects because in this data they were mated with random female cows, therefore their effect on the offspring would be random due to how the data was collected.
Factor effects model: \(Y_{ij}\) = \(\mu + \alpha_i +\epsilon_{ij}\) Where, \(Y_{ij}\) = jth observation of average daily weight gain of calves for the ith sire \(\mu\) = overall mean average of daily weight gain of calves, \(\alpha_i\) = the fixed effect of the ith sire, \(\epsilon_{ij}\) = error of the data
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 0.8.3 ✔ 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 object is masked from 'package:tidyr':
##
## expand
sired <- c("A", "B", "C", "D")
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)
calves <- data.frame(sired, weight)
ggplot(calves, aes(sired, weight))+
geom_point()
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
According to the ANOVA test it is reasonable to include sired as a random affect because of the large p value, meaning that the null hypothesis is true.
Variance components for modC
var_random_effect <- as.numeric(VarCorr(modC))
var_random_effect
## [1] 0.00507769
var_residual <- attr(VarCorr(modC), "sc")^2
var_residual
## [1] 0.015945
confint(modC, level = 0.95)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.00000000 0.1795777
## .sigma 0.09534633 0.1783855
## (Intercept) 0.97994356 1.1733897