Analysis of pig Growth records
Merk Project
Setup Code
Read file and check input
number of records, number of pigs, records per pig
# A tibble: 11,592 × 8
Pig Pen Trt Sex Parity age day weight
<dbl> <chr> <chr> <chr> <dbl> <chr> <dbl> <chr>
1 1 A02 F Male 2 0 0 1.821
2 1 A02 F Male 2 17 1 5.3659999999999997
3 1 A02 F Male 2 54 2 24.420999999999999
4 1 A02 F Male 2 110 3 85.474999999999994
5 1 A02 F Male 2 131 4 98.286000000000001
6 1 A02 F Male 2 151 5 110.149
7 1 A02 F Male 2 178 6 142.97800000000001
8 2 A02 F Male 2 0 0 1.2889999999999999
9 2 A02 F Male 2 20 1 7.1550000000000002
10 2 A02 F Male 2 57 2 26.123000000000001
# ℹ 11,582 more rows
[1] 1656
7
1656
# A tibble: 11,500 × 8
Pig Pen Trt Sex Parity age day weight
<dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 A02 F Male 2 0 0 1.82
2 1 A02 F Male 2 17 1 5.37
3 1 A02 F Male 2 54 2 24.4
4 1 A02 F Male 2 110 3 85.5
5 1 A02 F Male 2 131 4 98.3
6 1 A02 F Male 2 151 5 110.
7 1 A02 F Male 2 178 6 143.
8 2 A02 F Male 2 0 0 1.29
9 2 A02 F Male 2 20 1 7.16
10 2 A02 F Male 2 57 2 26.1
# ℹ 11,490 more rows
Over 11K growth records in over 1600 pigs
Summary of most important variables
Graphical summaries and tables to understand variable distributions
A B C D E F G
1653 1631 1653 1645 1637 1645 1636
Female Male
5729 5771
1 2 3 4 5
2864 3815 1063 2194 1564
Data are fairly “balanced” with respect to sex, number of records per pig, etc.
Weight data is clearly non-linear. A high (4th) order polinomial seems to capture the growth pattern.
To confirm this, we first look at the lowest order polynomial that models the mean trend.
model the mean of weights over time using polynomial regressions
Here we fit five models (m1…m5) corresponding to increasing polynomial degrees.
Models are compared using a likelihood ratio test.
Analysis of Variance Table
Model 1: weight ~ age
Model 2: weight ~ age + I(age^2)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11498 1786490
2 11497 1315940 1 470550 4111.1 < 0.00000000000000022 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table
Model 1: weight ~ age + I(age^2)
Model 2: weight ~ age + I(age^2) + I(age^3)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11497 1315940
2 11496 1239125 1 76814 712.65 < 0.00000000000000022 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table
Model 1: weight ~ age + I(age^2) + I(age^3)
Model 2: weight ~ age + I(age^2) + I(age^3) + I(age^4)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11496 1239125
2 11495 1238345 1 780.73 7.2472 0.007112 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table
Model 1: weight ~ age + I(age^2) + I(age^3) + I(age^4)
Model 2: weight ~ age + I(age^2) + I(age^3) + I(age^4) + I(age^5)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11495 1238345
2 11494 1238342 1 2.8583 0.0265 0.8706
Call:
lm(formula = weight ~ age + I(age^2) + I(age^3) + I(age^4), data = dts)
Residuals:
Min 1Q Median 3Q Max
-49.528 -4.011 -0.179 3.747 47.008
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.03688068790 0.24735261227 8.235 < 0.0000000000000002 ***
age -0.03871002658 0.02456354100 -1.576 0.11507
I(age^2) 0.00862460484 0.00061755249 13.966 < 0.0000000000000002 ***
I(age^3) -0.00003126580 0.00000525521 -5.949 0.00000000277 ***
I(age^4) 0.00000003825 0.00000001421 2.692 0.00711 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.38 on 11495 degrees of freedom
Multiple R-squared: 0.9539, Adjusted R-squared: 0.9539
F-statistic: 5.948e+04 on 4 and 11495 DF, p-value: < 0.00000000000000022
Call:
lm(formula = weight ~ Sex * poly(age, 4), data = dts)
Residuals:
Min 1Q Median 3Q Max
-47.547 -3.884 -0.219 3.709 49.483
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.7297 0.1348 435.607 < 0.0000000000000002 ***
SexMale 2.5002 0.1903 13.137 < 0.0000000000000002 ***
poly(age, 4)1 4868.7504 14.4329 337.336 < 0.0000000000000002 ***
poly(age, 4)2 620.8709 14.4498 42.967 < 0.0000000000000002 ***
poly(age, 4)3 -269.4672 14.4685 -18.624 < 0.0000000000000002 ***
poly(age, 4)4 37.1649 14.4045 2.580 0.00989 **
SexMale:poly(age, 4)1 279.3990 20.4093 13.690 < 0.0000000000000002 ***
SexMale:poly(age, 4)2 131.1193 20.4106 6.424 0.000000000138 ***
SexMale:poly(age, 4)3 -13.6857 20.4109 -0.671 0.50255
SexMale:poly(age, 4)4 -16.9606 20.4094 -0.831 0.40598
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.2 on 11490 degrees of freedom
Multiple R-squared: 0.9555, Adjusted R-squared: 0.9554
F-statistic: 2.739e+04 on 9 and 11490 DF, p-value: < 0.00000000000000022
Analysis of Variance Table
Model 1: weight ~ age + I(age^2) + I(age^3) + I(age^4)
Model 2: weight ~ Sex * age + Sex * I(age^2) + Sex * I(age^3) + Sex *
I(age^4)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 11495 1238345
2 11490 1196442 5 41902 80.482 < 0.00000000000000022 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
It looks like a 4th order polynomial models the mean trend appropriately. Then interaction with sex is significant.
Model variances and individual growth curves using random regression
Now that we have a good model for the mean trend, we can focus on random regressions of increasing order to model individual departures as well as covariance structures between repeated weights on the same animal.
Again, we use polynomials. But in this case, we use orthogonal polynomials, which are numerically more stable.
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Sex:t_1 + Sex:t_2 + Sex:t_2 + Sex:t_3 + Sex:t_4 + (1 |
Pen) + (t_1 - 1 | Pen) + (t_2 - 1 | Pen) + (1 | Pig) + (t_1 - 1 | Pig)
Data: dts
REML criterion at convergence: 77152.8
Scaled residuals:
Min 1Q Median 3Q Max
-3.7912 -0.5955 -0.0135 0.6733 3.5284
Random effects:
Groups Name Variance Std.Dev.
Pig t_1 17.307 4.160
Pig.1 (Intercept) 36.217 6.018
Pen t_2 7.104 2.665
Pen.1 t_1 18.197 4.266
Pen.2 (Intercept) 15.490 3.936
Residual 23.812 4.880
Number of obs: 11500, groups: Pig, 1656; Pen, 140
Fixed effects:
Estimate Std. Error t value
(Intercept) 60.27346 0.36694 164.262
SexFemale:t_1 46.06088 0.39394 116.922
SexMale:t_1 48.64063 0.39388 123.490
SexFemale:t_2 7.30294 0.23510 31.064
SexMale:t_2 8.51969 0.23504 36.247
SexFemale:t_3 -0.47055 0.06918 -6.801
SexMale:t_3 -0.62217 0.06869 -9.057
SexFemale:t_4 2.74739 0.07126 38.557
SexMale:t_4 2.62520 0.07143 36.753
Correlation of Fixed Effects:
(Intr) SxF:_1 SxM:_1 SxF:_2 SxM:_2 SxF:_3 SxM:_3 SxF:_4
SexFeml:t_1 0.000
SexMale:t_1 0.000 0.838
SexFeml:t_2 0.001 0.002 0.001
SexMale:t_2 0.001 0.001 0.002 0.922
SexFeml:t_3 0.004 0.006 0.005 0.024 0.017
SexMale:t_3 0.003 0.005 0.008 0.017 0.031 0.070
SexFeml:t_4 0.004 0.008 0.005 0.031 0.018 0.148 0.079
SexMale:t_4 0.004 0.006 0.008 0.019 0.033 0.079 0.154 0.091
integer(0)
We are fitting of two random regressions: one for the pen and one for the animal within the pen.
This means that the growth curve of the animal is broken down into three components:
the mean growth curve of the population: it’s a sex-specific 4th order polynomial. This predicts the mean of the whole… barn?
the pen deviation curve: it’s a random polinomial representing the pen-specific deviation from the mean curve.
the individual deviation: also a random polynomial, but in this case, it represents the deviation of a single animal.
How about predicting future weights?
let’s delete EVERY last observation, fit the same model, predict forward (forecast) weights and compare results
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Sex:t_1 + Sex:t_2 + Sex:t_2 + Sex:t_3 + Sex:t_4 + (1 |
Pen) + (t_1 - 1 | Pen) + (t_2 - 1 | Pen) + (1 | Pig) + (t_1 - 1 | Pig)
Data: dtmis
REML criterion at convergence: 65817.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.1157 -0.4585 0.0259 0.5322 3.0957
Random effects:
Groups Name Variance Std.Dev.
Pig t_1 14.46 3.803
Pig.1 (Intercept) 32.18 5.673
Pen t_2 10.01 3.164
Pen.1 t_1 22.35 4.728
Pen.2 (Intercept) 18.38 4.287
Residual 22.22 4.714
Number of obs: 9868, groups: Pig, 1656; Pen, 140
Fixed effects:
Estimate Std. Error t value
(Intercept) 57.2842 0.4009 142.879
SexFemale:t_1 42.5391 0.4557 93.341
SexMale:t_1 43.8927 0.4553 96.402
SexFemale:t_2 2.8545 0.3440 8.297
SexMale:t_2 2.4668 0.3433 7.186
SexFemale:t_3 -4.1758 0.1866 -22.383
SexMale:t_3 -5.6896 0.1841 -30.906
SexFemale:t_4 0.3756 0.1353 2.776
SexMale:t_4 -0.7071 0.1341 -5.273
Correlation of Fixed Effects:
(Intr) SxF:_1 SxM:_1 SxF:_2 SxM:_2 SxF:_3 SxM:_3 SxF:_4
SexFeml:t_1 0.061
SexMale:t_1 0.061 0.790
SexFeml:t_2 0.103 0.215 0.034
SexMale:t_2 0.102 0.034 0.213 0.663
SexFeml:t_3 0.159 0.333 0.053 0.565 0.088
SexMale:t_3 0.157 0.053 0.332 0.088 0.564 0.137
SexFeml:t_4 0.142 0.299 0.047 0.514 0.079 0.821 0.122
SexMale:t_4 0.140 0.047 0.297 0.079 0.511 0.122 0.818 0.108
integer(0)