REML and Mixed-Effects Model: Milk Yield Analysis

Introduction

This project aims to analyze the relationship between feed intake and milk yield in dairy cows using a mixed-effects model with REML (Restricted Maximum Likelihood) for random effects and BLUE (Best Linear Unbiased Estimators) for fixed effects. We use synthetic data to simulate cow milk production, feed intake, and cow-specific characteristics like breed, age, and lactation stage.

Key Objectives:

• Model milk yield based on feed intake using a linear mixed-effects model.
• Estimate variance due to individual cow differences using REML.
• Visualize relationships between variables and interpret model results.

Step 1: Load Libraries and Prepare Synthetic Data

# Load required libraries

Warning message:
package ‘RMySQL’ was built under R version 4.2.3 

library(lme4)  # For mixed-effects models and REML

Loading required package: Matrix

library(ggplot2)

Warning: package ‘ggplot2’ was built under R version 4.2.3

library(dplyr)

Warning: package ‘dplyr’ was built under R version 4.2.3
Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyr)

Attaching package: ‘tidyr’

The following objects are masked from ‘package:Matrix’:

    expand, pack, unpack

# 1.1 Load open-source or synthetic data
# Here we generate synthetic data as an example
set.seed(123)

# Generate synthetic data for 100 cows over 365 days
n_cows <- 100
n_days <- 365
cow_ids <- rep(1:n_cows, each = n_days)

# Generate synthetic cow factors (breed, age, lactation stage)
breed <- sample(c("Holstein", "Jersey", "Guernsey"), n_cows, replace = TRUE)
age <- sample(3:8, n_cows, replace = TRUE)
lactation_stage <- sample(1:3, n_cows, replace = TRUE)

# Generate daily feed intake (kg), milk yield (liters), and individual variability
feed_protein <- rnorm(n_cows * n_days, mean = 16, sd = 2)
feed_fat <- rnorm(n_cows * n_days, mean = 4, sd = 1)
feed_energy <- rnorm(n_cows * n_days, mean = 2800, sd = 500)  # kcal

# Random effect: individual cow variability
cow_random_effect <- rnorm(n_cows, mean = 0, sd = 2)

# Milk yield (liters) with random effect for cow
milk_yield <- 25 + 0.5 * feed_protein + 0.3 * feed_fat + cow_random_effect[rep(1:n_cows, each = n_days)] + 
              rnorm(n_cows * n_days, mean = 0, sd = 5)

# Combine into a data frame
synthetic_data <- data.frame(
  cow_id = cow_ids,
  breed = rep(breed, each = n_days),
  age = rep(age, each = n_days),
  lactation_stage = rep(lactation_stage, each = n_days),
  feed_protein = feed_protein,
  feed_fat = feed_fat,
  feed_energy = feed_energy,
  milk_yield = milk_yield
)

# View the first few rows of the data
head(synthetic_data)

NA

Step 2: Exploratory Data Analysis (EDA)

# 2.1 Summary statistics
summary(synthetic_data)

     cow_id          breed                age       lactation_stage
 Min.   :  1.00   Length:36500       Min.   :3.00   Min.   :1.00   
 1st Qu.: 25.75   Class :character   1st Qu.:4.00   1st Qu.:1.00   
 Median : 50.50   Mode  :character   Median :6.00   Median :1.50   
 Mean   : 50.50                      Mean   :5.77   Mean   :1.71   
 3rd Qu.: 75.25                      3rd Qu.:7.00   3rd Qu.:2.00   
 Max.   :100.00                      Max.   :8.00   Max.   :3.00   
  feed_protein      feed_fat       feed_energy       milk_yield   
 Min.   : 7.89   Min.   :-0.166   Min.   : 552.2   Min.   :13.08  
 1st Qu.:14.64   1st Qu.: 3.329   1st Qu.:2458.0   1st Qu.:30.38  
 Median :15.97   Median : 3.999   Median :2798.9   Median :34.08  
 Mean   :15.99   Mean   : 4.002   Mean   :2798.5   Mean   :34.06  
 3rd Qu.:17.33   3rd Qu.: 4.679   3rd Qu.:3135.7   3rd Qu.:37.75  
 Max.   :23.72   Max.   : 8.090   Max.   :5423.4   Max.   :55.99  

# 2.2 Visualize milk yield by breed
ggplot(synthetic_data, aes(x = breed, y = milk_yield)) +
  geom_boxplot() +
  labs(title = "Milk Yield by Breed", x = "Breed", y = "Milk Yield (liters)")

# 2.3 Visualize the relationship between feed protein and milk yield
ggplot(synthetic_data, aes(x = feed_protein, y = milk_yield)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", color = "blue") +
  labs(title = "Milk Yield vs Feed Protein Intake", x = "Feed Protein (kg)", y = "Milk Yield (liters)")

Step 3: Fit Mixed-Effects Model Using REML

We will use REML to fit a linear mixed-effects model that accounts for the random effect of individual cow variability and fixed effects like feed intake.

# 3.1 Fit a mixed-effects model using REML
# Random effect: cow-specific variation
# Fixed effects: feed_protein, feed_fat
model_reml <- lmer(milk_yield ~ feed_protein + feed_fat + (1 | cow_id), data = synthetic_data, REML = TRUE)

# 3.2 Model summary
summary(model_reml)

Linear mixed model fit by REML ['lmerMod']
Formula: milk_yield ~ feed_protein + feed_fat + (1 | cow_id)
   Data: synthetic_data

REML criterion at convergence: 221699.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9714 -0.6788  0.0002  0.6776  4.1166 

Random effects:
 Groups   Name        Variance Std.Dev.
 cow_id   (Intercept)  3.099   1.760   
 Residual             25.161   5.016   
Number of obs: 36500, groups:  cow_id, 100

Fixed effects:
             Estimate Std. Error t value
(Intercept)  24.85800    0.29480   84.32
feed_protein  0.49356    0.01317   37.47
feed_fat      0.32652    0.02632   12.40

Correlation of Fixed Effects:
            (Intr) fd_prt
feed_proten -0.713       
feed_fat    -0.354 -0.005

# 3.3 Extract variance components (random effect)
VarCorr(model_reml)

 Groups   Name        Std.Dev.
 cow_id   (Intercept) 1.7603  
 Residual             5.0161  

# 3.4 Visualize fitted vs actual milk yield
synthetic_data$fitted_values <- fitted(model_reml)

ggplot(synthetic_data, aes(x = milk_yield, y = fitted_values)) +
  geom_point(alpha = 0.5, color = "green") +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  labs(title = "Actual vs Fitted Milk Yield", x = "Actual Milk Yield (liters)", y = "Fitted Milk Yield (liters)")

Step 4: Calculate BLUE for Fixed Effects

We will calculate the Best Linear Unbiased Estimator (BLUE) for the fixed effects in the model, which will give us unbiased estimates of the effects of feed intake on milk yield.

# 4.1 Fixed effects are already estimated by the model
fixed_effects <- fixef(model_reml)
print(fixed_effects)

 (Intercept) feed_protein     feed_fat 
  24.8579981    0.4935608    0.3265210 

# 4.2 BLUE values (fixed effect estimates)
# The fixed effect estimates are the BLUE for this mixed-effects model:
# These values represent the unbiased estimates of feed_protein and feed_fat effects

Step 5: Model Diagnostics and Interpretation

# 5.1 Residual plot to check model assumptions
plot(residuals(model_reml), main = "Residuals of Mixed-Effects Model")

# 5.2 Interpret the results
# The model shows how much feed protein and feed fat influence milk yield on average,
# while accounting for the random variation between cows.

# 5.3 Calculate random effect variance
random_effect_variance <- as.data.frame(VarCorr(model_reml))$vcov[2]  # Cow-specific variance
print(paste("Random Effect Variance (Cow):", round(random_effect_variance, 2)))

[1] "Random Effect Variance (Cow): 25.16"