Herd-level factors associated with lameness in lactating dairy cows housed in freestall and compost-bedded pack farms in southern Brazil, 2016
José A. Bran, Joao H. C. Costa, Marina A.G. von Keyserlingk, Maria José Hötzel
Sep 2019
library(sjstats)
library(readr)
library(car)
library(dplyr)
library(lme4)
library(epitools)
library(lattice)
library(blmeco)
library(ggplot2)
library(tidyr)
library(performance)
library(arm)
library(corrplot)Data
The data correspond to a cross-sectional study designed to investigate factors associated with lameness in dairy cows on intensive farms in southern Brazil.
Farms: 38 freestall and 12 compost-bedded pack visited once in 2016 (March to October). All lactating cows (n = 13,716) were examined and body condition score (BCS) and gait score were assessed.
The following dataframe (Herd) is a data subset with information about the associations between lameness and herd-level explanatory variables. The variables presented in the dataframe “Herd” were collected through inspection of facilities and using data from an interview with farmers on routine management practices applied on each farm.
Additional information is provided in the published paper (“Factors associated with lameness prevalence in lactating cows housed in freestall and compost-bedded pack dairy farms in southern Brazil” https://doi.org/10.1016/j.prevetmed.2019.104773)
Herd<- read_csv("Herd.csv")Variables included in the dataframe (a description is available below: final lines on the script)
names(Herd)## [1] "farm_ID" "Scored cows"
## [3] "LS_1" "LS_2"
## [5] "LS_3" "LS_4"
## [7] "LS_5" "Clinical lameness"
## [9] "Severe lameness" "Average milk yield per cow"
## [11] "Total cows" "Workers number"
## [13] "Farm type" "Bedding type"
## [15] "Dry period length" "Alley scrapping type"
## [17] "Farm type herd size" "Barn building year"
## [19] "Feeding type" "Heat control system type"
## [21] "Years of milk activity" "Feeding times per day"
## [23] "Pasture access to milking cows" "Pasture access to cows"
## [25] "Milkings per day" "Stall cleaning"
## [27] "Footbath" "Footbath solution"
## [29] "Preventive hoof trimming" "Hoof trimming freq"
## [31] "Track slope" "Track covering"
## [33] "Track pavement type" "Track moisture"
## [35] "Feed bunk flooring" "Feed bunk type"
## [37] "Holding area water trough" "Holding area flooring"
## [39] "Holding area roof cover" "Holding area Area"
## [41] "Do cows fall on your farm" "Cows falls category"
## [43] "Feed bunk flooring category" "Slipperines feed bunk"
## [45] "Herd size"Exploring the variable “farm type and herd size” used as random effect in the models
A category of farm type (i.e., compost barn and freestall farms) by herd size (quartiles) was built to be used as random effect
Lameness prevalence distribution by farm type and size
Herd%>%
group_by( `Farm type herd size`)%>%
summarise(Prevalence = (mean(`Clinical lameness`)),
Observations = n())| Farm type herd size | Prevalence | Observations |
|---|---|---|
| Compost_(140,204] | 29.57529 | 3 |
| Compost_(204,359] | 32.23141 | 1 |
| Compost_(40,140] | 33.16685 | 8 |
| Freestall_(140,204] | 47.34190 | 9 |
| Freestall_(204,359] | 44.41726 | 11 |
| Freestall_(359,901] | 45.52086 | 13 |
| Freestall_(40,140] | 32.82182 | 5 |
Ordering the variable “Herd size”
Herd$`Herd size` = factor(Herd$`Herd size`, levels = c("(40,140]", "(140,204]", "(204,359]", "(359,901]"))Plot: lameness prevalence distribution by farm type and herd size
ggplot(Herd, aes(x = `Herd size` , y = `Clinical lameness`, color = `Farm type`))+
geom_violin()+
geom_dotplot(binaxis = "y",
stackdir = "center",
dotsize = 0.6,
stackratio = 1.1,
binwidth = 1.3)+
facet_wrap(~`Farm type`)+
labs( x = "Lactating herd size (n)",
y ="Herd clinical lameness prevalence (%)")+
scale_color_discrete(guide = FALSE)
Within-herd lameness prevalence
summary(Herd$`Clinical lameness`)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.76 34.38 42.47 41.14 47.93 64.46sd(Herd$`Clinical lameness`)## [1] 11.31161Exploring correlations between variables
Herd%>%
select_if(is.numeric)%>%
cor(method = "spearman", use = "complete.obs")## farm_ID Scored cows LS_1
## farm_ID 1.00000000 0.37864684 0.139784587
## Scored cows 0.37864684 1.00000000 0.264044719
## LS_1 0.13978459 0.26404472 1.000000000
## LS_2 0.31821513 0.95662560 0.220578948
## LS_3 0.40571135 0.97103537 0.257887901
## LS_4 0.37794849 0.91710099 0.245446240
## LS_5 0.23238246 0.77186853 0.118494612
## Clinical lameness 0.36810682 0.34396287 0.009318972
## Severe lameness 0.29418150 0.29652346 -0.027956917
## Average milk yield per cow 0.08193936 0.40600802 0.196172817
## Total cows 0.39487611 0.99603756 0.276470353
## Workers number 0.34930642 0.89608899 0.252539186
## Barn building year -0.52003252 -0.47481216 -0.352230724
## Years of milk activity 0.23000058 0.43261378 0.236221714
## Feeding times per day -0.16037319 0.08849234 0.254575753
## Milkings per day 0.02344600 0.48670474 0.099923955
## Holding area Area 0.07984801 0.47182716 0.093197313
## LS_2 LS_3 LS_4 LS_5
## farm_ID 0.31821513 0.4057114 0.3779485 0.2323825
## Scored cows 0.95662560 0.9710354 0.9171010 0.7718685
## LS_1 0.22057895 0.2578879 0.2454462 0.1184946
## LS_2 1.00000000 0.8977304 0.8158323 0.6515493
## LS_3 0.89773036 1.0000000 0.9113814 0.7582206
## LS_4 0.81583234 0.9113814 1.0000000 0.8509594
## LS_5 0.65154935 0.7582206 0.8509594 1.0000000
## Clinical lameness 0.12735119 0.4446918 0.5985723 0.5978195
## Severe lameness 0.10504031 0.3224843 0.5962922 0.6297483
## Average milk yield per cow 0.47352793 0.3320280 0.3603894 0.2630949
## Total cows 0.95608274 0.9680493 0.9168295 0.7646217
## Workers number 0.90305132 0.8427469 0.8207427 0.6590288
## Barn building year -0.46208992 -0.4596341 -0.4383637 -0.3402406
## Years of milk activity 0.40091257 0.4152701 0.3333334 0.4053449
## Feeding times per day 0.05285204 0.1136353 0.1114343 0.1416704
## Milkings per day 0.52065288 0.4031701 0.4528869 0.4354203
## Holding area Area 0.44579262 0.4654414 0.3794451 0.3948644
## Clinical lameness Severe lameness
## farm_ID 0.368106817 0.29418150
## Scored cows 0.343962874 0.29652346
## LS_1 0.009318972 -0.02795692
## LS_2 0.127351194 0.10504031
## LS_3 0.444691780 0.32248433
## LS_4 0.598572254 0.59629219
## LS_5 0.597819541 0.62974829
## Clinical lameness 1.000000000 0.90914025
## Severe lameness 0.909140252 1.00000000
## Average milk yield per cow -0.073342798 -0.06061118
## Total cows 0.337937960 0.28729612
## Workers number 0.240352454 0.24274944
## Barn building year -0.201357376 -0.15685941
## Years of milk activity 0.151360479 0.10373108
## Feeding times per day 0.194832072 0.12158470
## Milkings per day 0.112927401 0.17035763
## Holding area Area 0.072357172 0.05259872
## Average milk yield per cow Total cows
## farm_ID 0.081939357 0.3948761
## Scored cows 0.406008019 0.9960376
## LS_1 0.196172817 0.2764704
## LS_2 0.473527928 0.9560827
## LS_3 0.332027998 0.9680493
## LS_4 0.360389447 0.9168295
## LS_5 0.263094918 0.7646217
## Clinical lameness -0.073342798 0.3379380
## Severe lameness -0.060611184 0.2872961
## Average milk yield per cow 1.000000000 0.4147680
## Total cows 0.414768042 1.0000000
## Workers number 0.476218156 0.8967700
## Barn building year -0.300229469 -0.4815388
## Years of milk activity -0.002123204 0.4378819
## Feeding times per day 0.256244085 0.1061671
## Milkings per day 0.708866728 0.4764468
## Holding area Area 0.064479445 0.4656389
## Workers number Barn building year
## farm_ID 0.3493064 -0.52003252
## Scored cows 0.8960890 -0.47481216
## LS_1 0.2525392 -0.35223072
## LS_2 0.9030513 -0.46208992
## LS_3 0.8427469 -0.45963410
## LS_4 0.8207427 -0.43836373
## LS_5 0.6590288 -0.34024059
## Clinical lameness 0.2403525 -0.20135738
## Severe lameness 0.2427494 -0.15685941
## Average milk yield per cow 0.4762182 -0.30022947
## Total cows 0.8967700 -0.48153878
## Workers number 1.0000000 -0.46684525
## Barn building year -0.4668453 1.00000000
## Years of milk activity 0.3960663 -0.31976479
## Feeding times per day 0.1601008 0.09549196
## Milkings per day 0.6162540 -0.27315987
## Holding area Area 0.3685123 0.06067903
## Years of milk activity Feeding times per day
## farm_ID 0.230000584 -0.16037319
## Scored cows 0.432613784 0.08849234
## LS_1 0.236221714 0.25457575
## LS_2 0.400912565 0.05285204
## LS_3 0.415270146 0.11363532
## LS_4 0.333333361 0.11143428
## LS_5 0.405344916 0.14167037
## Clinical lameness 0.151360479 0.19483207
## Severe lameness 0.103731078 0.12158470
## Average milk yield per cow -0.002123204 0.25624409
## Total cows 0.437881945 0.10616708
## Workers number 0.396066308 0.16010084
## Barn building year -0.319764790 0.09549196
## Years of milk activity 1.000000000 0.11260724
## Feeding times per day 0.112607239 1.00000000
## Milkings per day 0.167620518 0.43339225
## Holding area Area 0.274285130 0.15661978
## Milkings per day Holding area Area
## farm_ID 0.02344600 0.07984801
## Scored cows 0.48670474 0.47182716
## LS_1 0.09992395 0.09319731
## LS_2 0.52065288 0.44579262
## LS_3 0.40317012 0.46544142
## LS_4 0.45288688 0.37944514
## LS_5 0.43542032 0.39486439
## Clinical lameness 0.11292740 0.07235717
## Severe lameness 0.17035763 0.05259872
## Average milk yield per cow 0.70886673 0.06447945
## Total cows 0.47644684 0.46563891
## Workers number 0.61625400 0.36851233
## Barn building year -0.27315987 0.06067903
## Years of milk activity 0.16762052 0.27428513
## Feeding times per day 0.43339225 0.15661978
## Milkings per day 1.00000000 0.12026407
## Holding area Area 0.12026407 1.00000000Herd%>%
select_if(is.numeric)%>%
cor(method = "spearman", use = "complete.obs")%>%
corrplot()
Explanatory variables considered in the multivariable model
a) Stall bedding type on each farm
Distribution of lameness prevalence by bedding type:
Herd%>%
group_by(`Bedding type`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Bedding type | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| Compost | 32.19101 | 31.96936 | 6.812866 | 12 |
| Deep-bedding | 39.33045 | 42.81609 | 12.081737 | 11 |
| Mattress | 47.92168 | 47.56933 | 11.579573 | 16 |
| Other | 42.83323 | 42.12963 | 6.203443 | 11 |
Mixed-effects linear model
summary(fma<-lmer(`Clinical lameness` ~ `Bedding type` + (1|`Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula:
## `Clinical lameness` ~ `Bedding type` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 349.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.02910 -0.61269 -0.00059 0.48509 1.99377
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 22.51 4.745
## Residual 85.42 9.242
## Number of obs: 50, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 31.826 4.117 7.731
## `Bedding type`Deep-bedding 7.087 5.522 1.283
## `Bedding type`Mattress 15.949 5.308 3.005
## `Bedding type`Other 7.853 5.731 1.370
##
## Correlation of Fixed Effects:
## (Intr) `Bt`D- `Bty`M
## `Bddtyp`Dp- -0.745
## `Bddngtyp`M -0.776 0.773
## `Bddngtyp`O -0.718 0.719 0.729Anova(fma)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Bedding type |
12.13249 | 3 | 0.0069428 |
b) Feed bunk alley floor slipperiness score
Distribution of lameness prevalence by floor slipperiness score:
Herd%>%
group_by(`Slipperines feed bunk`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Slipperines feed bunk | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| Non-slippery | 33.63538 | 33.84057 | 12.963852 | 12 |
| Slightly slippery | 41.40932 | 42.30357 | 12.270597 | 20 |
| Very slippery | 47.90648 | 47.67442 | 8.196792 | 15 |
| NA | 35.47705 | 31.52174 | 6.295290 | 3 |
Distribution of lameness prevalence by farmers answer to the following question:
Do the cows usually fall at the barn?
No = cows never fall, Other = cows fall often, somethimes or rare times
Herd%>%
group_by(`Cows falls category`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Cows falls category | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| No | 34.42170 | 31.61453 | 10.06720 | 14 |
| Other | 43.74821 | 43.66131 | 10.63536 | 36 |
Mixed-effects linear model
summary(fmb<-lmer(`Clinical lameness` ~ `Slipperines feed bunk` + (1|`Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula:
## `Clinical lameness` ~ `Slipperines feed bunk` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 330.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78617 -0.62470 -0.05037 0.48615 2.21003
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 44.94 6.704
## Residual 72.88 8.537
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 32.054 3.588 8.934
## `Slipperines feed bunk`Slightly slippery 6.379 3.406 1.873
## `Slipperines feed bunk`Very slippery 14.141 3.367 4.200
##
## Correlation of Fixed Effects:
## (Intr) `Sfb`Ss
## `Slpfbnk`Ss -0.511
## `Slpfbnk`Vs -0.487 0.559Anova(fmb)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Slipperines feed bunk |
17.9674 | 2 | 0.0001254 |
c) Presence of a roof covering the holding area
Distribution of lameness prevalence filtered by levels of the variable:
Herd%>%
group_by(`Holding area roof cover`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Holding area roof cover | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| Covered | 39.85330 | 41.68981 | 15.163049 | 40 |
| Uncovered | 46.27075 | 45.50234 | 6.974189 | 10 |
Mixed-effects linear model
summary(fmc<-lmer(`Clinical lameness` ~ `Holding area roof cover` + (1|`Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula:
## `Clinical lameness` ~ `Holding area roof cover` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 363.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.84102 -0.71619 0.01489 0.56614 1.80405
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 52.69 7.259
## Residual 82.80 9.099
## Number of obs: 50, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 36.683 3.259 11.26
## `Holding area roof cover`Uncovered 10.198 3.642 2.80
##
## Correlation of Fixed Effects:
## (Intr)
## `Hldarcvr`U -0.239Anova(fmc)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Holding area roof cover |
7.84033 | 1 | 0.0051093 |
d) Dry period length criteria defined by farmers in each herd
Distribution of lameness prevalence by dry period length categories:
Herd%>%
group_by(`Dry period length`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Dry period length | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| Less than sixty days | 44.77446 | 45.50234 | 16.30728 | 18 |
| Sixty days | 39.09060 | 41.26200 | 14.35805 | 32 |
changing the reference category for this variable
Herd$`Dry period length`<- as.factor(Herd$`Dry period length`)
contrasts(Herd$`Dry period length`) <- contr.Treatment(levels(Herd$`Dry period length`), base = 2)
contrasts(Herd$`Dry period length`)## [T.Less than sixty days]
## Less than sixty days 1
## Sixty days 0Mixed-effects linear model
summary(fmd<-lmer(`Clinical lameness` ~ `Dry period length` + (1|`Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula:
## `Clinical lameness` ~ `Dry period length` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 364.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1787 -0.6007 0.1024 0.6224 1.9539
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 45.27 6.728
## Residual 85.52 9.248
## Number of obs: 50, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 36.591 3.132 11.683
## `Dry period length`[T.Less than sixty days] 7.613 2.926 2.602
##
## Correlation of Fixed Effects:
## (Intr)
## `Dpl`[T.tsd -0.293Anova(fmd)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Dry period length |
6.770323 | 1 | 0.0092686 |
e) Alley manure handling system
Distribution of lameness prevalence by categories of this variable:
Herd%>%
group_by(`Alley scrapping type`)%>%
summarise(`Mean prevalence` = mean(`Clinical lameness`),
`Median prevalence` = median(`Clinical lameness`),
Iqr = IQR(`Clinical lameness`),
Farms = n())| Alley scrapping type | Mean prevalence | Median prevalence | Iqr | Farms |
|---|---|---|---|---|
| Scrapper | 44.42257 | 46.29630 | 12.06680 | 17 |
| Tractor Bob Cat | 39.44411 | 40.90909 | 10.92948 | 33 |
Mixed-effects linear model
summary(fme<-lmer(`Clinical lameness` ~ `Alley scrapping type` + (1|`Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula:
## `Clinical lameness` ~ `Alley scrapping type` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 367.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.97773 -0.67592 -0.02163 0.61746 2.13113
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 41.29 6.426
## Residual 92.53 9.619
## Number of obs: 50, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 42.614 3.509 12.144
## `Alley scrapping type`Tractor Bob Cat -5.236 2.887 -1.813
##
## Correlation of Fixed Effects:
## (Intr)
## `Alstyp`TBC -0.554Anova(fme)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Alley scrapping type |
3.288649 | 1 | 0.0697604 |
Herd-level linear analyses: multivariable model
Mixed-effects linear regressions were fitted. Significant explanatory variables identified at the univariable models (P<0.2) were included in the multivariable model
Full model
summary(fm_a <- lmer(`Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` + `Holding area roof cover` + `Dry period length` +`Alley scrapping type` + (1| `Farm type herd size`), data = Herd))## Linear mixed model fit by REML ['lmerMod']
## Formula: `Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` +
## `Holding area roof cover` + `Dry period length` + `Alley scrapping type` +
## (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 283.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.99010 -0.59177 0.00729 0.62826 1.90547
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 40.39 6.356
## Residual 49.46 7.033
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 25.438 5.161 4.929
## `Bedding type`Deep-bedding 6.940 5.997 1.157
## `Bedding type`Mattress 14.503 5.908 2.455
## `Bedding type`Other 5.169 6.186 0.836
## `Slipperines feed bunk`Slightly slippery 5.879 2.968 1.981
## `Slipperines feed bunk`Very slippery 11.336 2.971 3.816
## `Holding area roof cover`Uncovered 4.765 3.389 1.406
## `Dry period length`[T.Less than sixty days] 4.909 2.476 1.983
## `Alley scrapping type`Tractor Bob Cat -1.896 2.477 -0.765
##
## Correlation of Fixed Effects:
## (Intr) `Bt`D- `Bty`M `Bty`O `Sfb`Ss `Sfb`Vs `Harc` `pltsd
## `Bddtyp`Dp- -0.609
## `Bddngtyp`M -0.700 0.863
## `Bddngtyp`O -0.620 0.843 0.844
## `Slpfbnk`Ss -0.265 -0.084 0.003 0.032
## `Slpfbnk`Vs -0.223 -0.017 -0.023 0.003 0.571
## `Hldarcvr`U 0.076 -0.079 -0.109 -0.169 -0.235 -0.320
## `Dpl`[T.tsd -0.241 0.019 0.055 0.023 0.014 -0.091 -0.133
## `Alstyp`TBC -0.367 -0.086 0.102 -0.056 0.000 -0.011 -0.022 0.281Anova(fm_a)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Bedding type |
12.4037229 | 3 | 0.0061207 |
Slipperines feed bunk |
14.6182602 | 2 | 0.0006694 |
Holding area roof cover |
1.9764095 | 1 | 0.1597692 |
Dry period length |
3.9314782 | 1 | 0.0473903 |
Alley scrapping type |
0.5856873 | 1 | 0.4440917 |
vif(fm_a)## GVIF Df GVIF^(1/(2*Df))
## `Bedding type` 1.319540 3 1.047298
## `Slipperines feed bunk` 1.236991 2 1.054609
## `Holding area roof cover` 1.202281 1 1.096486
## `Dry period length` 1.139710 1 1.067572
## `Alley scrapping type` 1.275237 1 1.129264We removed the variable Alley scrapping type
summary(fm_b<-update(fm_a, .~.- `Alley scrapping type`))## Linear mixed model fit by REML ['lmerMod']
## Formula: `Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` +
## `Holding area roof cover` + `Dry period length` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 287.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.14892 -0.69650 0.06796 0.56668 1.93930
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 41.08 6.410
## Residual 48.79 6.985
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 24.014 4.813 4.989
## `Bedding type`Deep-bedding 6.516 5.997 1.086
## `Bedding type`Mattress 14.926 5.901 2.529
## `Bedding type`Other 4.842 6.197 0.781
## `Slipperines feed bunk`Slightly slippery 5.862 2.951 1.986
## `Slipperines feed bunk`Very slippery 11.301 2.952 3.828
## `Holding area roof cover`Uncovered 4.750 3.371 1.409
## `Dry period length`[T.Less than sixty days] 5.443 2.361 2.305
##
## Correlation of Fixed Effects:
## (Intr) `Bt`D- `Bty`M `Bty`O `Sfb`Ss `Sfb`Vs `Harc`
## `Bddtyp`Dp- -0.693
## `Bddngtyp`M -0.717 0.882
## `Bddngtyp`O -0.691 0.845 0.859
## `Slpfbnk`Ss -0.282 -0.083 0.003 0.032
## `Slpfbnk`Vs -0.241 -0.018 -0.022 0.002 0.572
## `Hldarcvr`U 0.072 -0.081 -0.106 -0.170 -0.236 -0.321
## `Dpl`[T.tsd -0.152 0.044 0.027 0.040 0.016 -0.091 -0.133Anova(fm_b)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Bedding type |
16.134362 | 3 | 0.0010643 |
Slipperines feed bunk |
14.717090 | 2 | 0.0006371 |
Holding area roof cover |
1.986333 | 1 | 0.1587249 |
Dry period length |
5.313939 | 1 | 0.0211555 |
We excluded the variable: presence of a roof covering the holding area
Final model
summary(fm_c<-update(fm_b, .~.-`Holding area roof cover`))## Linear mixed model fit by REML ['lmerMod']
## Formula: `Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` +
## `Dry period length` + (1 | `Farm type herd size`)
## Data: Herd
##
## REML criterion at convergence: 293.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25175 -0.69175 0.09233 0.62219 1.87320
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 22.71 4.765
## Residual 52.48 7.244
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 23.078 4.163 5.544
## `Bedding type`Deep-bedding 7.729 5.030 1.537
## `Bedding type`Mattress 16.459 4.887 3.368
## `Bedding type`Other 7.354 5.171 1.422
## `Slipperines feed bunk`Slightly slippery 6.934 2.927 2.369
## `Slipperines feed bunk`Very slippery 12.556 2.889 4.346
## `Dry period length`[T.Less than sixty days] 5.732 2.396 2.393
##
## Correlation of Fixed Effects:
## (Intr) `Bt`D- `Bty`M `Bty`O `Sfb`Ss `Sfb`Vs
## `Bddtyp`Dp- -0.657
## `Bddngtyp`M -0.693 0.819
## `Bddngtyp`O -0.659 0.773 0.790
## `Slpfbnk`Ss -0.337 -0.118 -0.019 -0.002
## `Slpfbnk`Vs -0.286 -0.047 -0.063 -0.061 0.541
## `Dpl`[T.tsd -0.182 0.056 0.031 0.033 -0.032 -0.150Anova(fm_c)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
Bedding type |
17.658599 | 3 | 0.0005172 |
Slipperines feed bunk |
18.892468 | 2 | 0.0000790 |
Dry period length |
5.724872 | 1 | 0.0167263 |
confint(fm_c, method ="Wald")## 2.5 % 97.5 %
## .sig01 NA NA
## .sigma NA NA
## (Intercept) 14.918885 31.23700
## `Bedding type`Deep-bedding -2.129215 17.58627
## `Bedding type`Mattress 6.881337 26.03676
## `Bedding type`Other -2.781136 17.48815
## `Slipperines feed bunk`Slightly slippery 1.198007 12.66999
## `Slipperines feed bunk`Very slippery 6.894061 18.21771
## `Dry period length`[T.Less than sixty days] 1.036599 10.42724vif(fm_c)## GVIF Df GVIF^(1/(2*Df))
## `Bedding type` 1.077740 3 1.012556
## `Slipperines feed bunk` 1.102008 2 1.024581
## `Dry period length` 1.032383 1 1.016063plot(fm_c)
scatter.smooth(resid(fm_c, type = "pearson")~fitted(fm_c), col = "grey", las = 1)
r2(fm_c)## # R2 for mixed models
##
## Conditional R2: 0.637
## Marginal R2: 0.481icc(fm_c)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.302
## Conditional ICC: 0.157coef(fm_c)## $`Farm type herd size`
## (Intercept) `Bedding type`Deep-bedding
## Compost_(140,204] 23.07835 7.728526
## Compost_(204,359] 25.84215 7.728526
## Compost_(40,140] 20.31332 7.728526
## Freestall_(140,204] 26.25905 7.728526
## Freestall_(204,359] 23.06450 7.728526
## Freestall_(359,901] 26.02315 7.728526
## Freestall_(40,140] 16.96507 7.728526
## `Bedding type`Mattress `Bedding type`Other
## Compost_(140,204] 16.45905 7.353506
## Compost_(204,359] 16.45905 7.353506
## Compost_(40,140] 16.45905 7.353506
## Freestall_(140,204] 16.45905 7.353506
## Freestall_(204,359] 16.45905 7.353506
## Freestall_(359,901] 16.45905 7.353506
## Freestall_(40,140] 16.45905 7.353506
## `Slipperines feed bunk`Slightly slippery
## Compost_(140,204] 6.934
## Compost_(204,359] 6.934
## Compost_(40,140] 6.934
## Freestall_(140,204] 6.934
## Freestall_(204,359] 6.934
## Freestall_(359,901] 6.934
## Freestall_(40,140] 6.934
## `Slipperines feed bunk`Very slippery
## Compost_(140,204] 12.55588
## Compost_(204,359] 12.55588
## Compost_(40,140] 12.55588
## Freestall_(140,204] 12.55588
## Freestall_(204,359] 12.55588
## Freestall_(359,901] 12.55588
## Freestall_(40,140] 12.55588
## `Dry period length`[T.Less than sixty days]
## Compost_(140,204] 5.731921
## Compost_(204,359] 5.731921
## Compost_(40,140] 5.731921
## Freestall_(140,204] 5.731921
## Freestall_(204,359] 5.731921
## Freestall_(359,901] 5.731921
## Freestall_(40,140] 5.731921
##
## attr(,"class")
## [1] "coef.mer"qqPlot(ranef(fm_c)$`Farm type herd size`[,1])
## [1] 7 4ranef(fm_c)## $`Farm type herd size`
## (Intercept)
## Compost_(140,204] 0.0004099261
## Compost_(204,359] 2.7642117367
## Compost_(40,140] -2.7646216627
## Freestall_(140,204] 3.1811080891
## Freestall_(204,359] -0.0134435372
## Freestall_(359,901] 2.9452114335
## Freestall_(40,140] -6.1128759855
##
## with conditional variances for "Farm type herd size"qqmath(ranef(fm_c, condVar = T))## $`Farm type herd size`
dotplot(ranef(fm_c, condVar = TRUE))## $`Farm type herd size`
Influence and leverage
plot(cooks.distance(fm_c), type= "h", main = "Cook's distance",
ylab="D", xlab = "Observation number", las=1 )
Leverage
sort(hatvalues(fm_c), decreasing = TRUE)[1:11]## 6 3 12 8 42 22 39
## 0.4628679 0.3012295 0.2927468 0.2919880 0.2896449 0.2840125 0.2745458
## 36 46 33 21
## 0.2738912 0.2669690 0.2644980 0.2561675which(hatvalues(fm_c) > 0.2)## 1 3 6 8 10 12 13 19 21 22 25 26 27 33 36 37 39 41 42 43 45 46
## 1 3 5 7 8 10 11 17 19 20 22 23 24 30 33 34 36 38 39 40 42 43plot(hatvalues(fm_c), type= "h", main = "Hat",
ylab="", xlab = "Observation number", las=1 )
influencePlot(fm_c)
| StudRes | Hat | CookD | |
|---|---|---|---|
| 3 | -0.1360865 | 0.3012295 | 0.0011405 |
| 6 | 1.2033933 | 0.4628679 | 0.1782761 |
| 7 | -2.4975468 | 0.1871433 | 0.2051585 |
| 22 | 2.2137653 | 0.2840125 | 0.2777135 |
Likelihood ratio test
Excluding three farms with missing data on “Slipperiness…” variable
Herd1 = subset(Herd, !is.na(`Slipperines feed bunk`))Null model
summary(fmn<-lmer(`Clinical lameness` ~ 1 + (1|`Farm type herd size`), data = Herd1))## Linear mixed model fit by REML ['lmerMod']
## Formula: `Clinical lameness` ~ 1 + (1 | `Farm type herd size`)
## Data: Herd1
##
## REML criterion at convergence: 353.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15523 -0.55807 -0.04257 0.62041 1.87396
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 41.63 6.452
## Residual 101.37 10.068
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 39.115 3.001 13.04r2(fmn)## # R2 for mixed models
##
## Conditional R2: 0.291
## Marginal R2: 0.000summary(fmn2<- lmer(`Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` + `Dry period length` + (1 | `Farm type herd size`), data = Herd1))## Linear mixed model fit by REML ['lmerMod']
## Formula: `Clinical lameness` ~ `Bedding type` + `Slipperines feed bunk` +
## `Dry period length` + (1 | `Farm type herd size`)
## Data: Herd1
##
## REML criterion at convergence: 293.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.25175 -0.69175 0.09233 0.62219 1.87320
##
## Random effects:
## Groups Name Variance Std.Dev.
## Farm type herd size (Intercept) 22.71 4.765
## Residual 52.48 7.244
## Number of obs: 47, groups: Farm type herd size, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 23.078 4.163 5.544
## `Bedding type`Deep-bedding 7.729 5.030 1.537
## `Bedding type`Mattress 16.459 4.887 3.368
## `Bedding type`Other 7.354 5.171 1.422
## `Slipperines feed bunk`Slightly slippery 6.934 2.927 2.369
## `Slipperines feed bunk`Very slippery 12.556 2.889 4.346
## `Dry period length`[T.Less than sixty days] 5.732 2.396 2.393
##
## Correlation of Fixed Effects:
## (Intr) `Bt`D- `Bty`M `Bty`O `Sfb`Ss `Sfb`Vs
## `Bddtyp`Dp- -0.657
## `Bddngtyp`M -0.693 0.819
## `Bddngtyp`O -0.659 0.773 0.790
## `Slpfbnk`Ss -0.337 -0.118 -0.019 -0.002
## `Slpfbnk`Vs -0.286 -0.047 -0.063 -0.061 0.541
## `Dpl`[T.tsd -0.182 0.056 0.031 0.033 -0.032 -0.150likelihood ratio test
anova(fmn2, fmn)## refitting model(s) with ML (instead of REML)| Df | AIC | BIC | logLik | deviance | Chisq | Chi Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| fmn | 3 | 363.8399 | 369.3903 | -178.9200 | 357.8399 | NA | NA | NA |
| fmn2 | 9 | 338.3377 | 354.9891 | -160.1689 | 320.3377 | 37.50216 | 6 | 1.4e-06 |
Standard error for fixed and random effecs
library(arm)
se.coef(fmn2)## $fixef
## [1] 4.162860 5.029552 4.886679 5.170831 2.926581 2.888738 2.395617
##
## $`Farm type herd size`
## (Intercept)
## Compost_(140,204] 3.143378
## Compost_(204,359] 3.981024
## Compost_(40,140] 2.512837
## Freestall_(140,204] 2.153988
## Freestall_(204,359] 1.985590
## Freestall_(359,901] 1.851374
## Freestall_(40,140] 2.883622seran=se.ranef(fmn2)
seran## $`Farm type herd size`
## (Intercept)
## Compost_(140,204] 3.143378
## Compost_(204,359] 3.981024
## Compost_(40,140] 2.512837
## Freestall_(140,204] 2.153988
## Freestall_(204,359] 1.985590
## Freestall_(359,901] 1.851374
## Freestall_(40,140] 2.883622ggplot(Herd , aes(farm_ID, `Clinical lameness`))+
geom_col()+
facet_wrap(~`Farm type herd size`)
Confidence interval for random effect
coef(fmn2)## $`Farm type herd size`
## (Intercept) `Bedding type`Deep-bedding
## Compost_(140,204] 23.07835 7.728526
## Compost_(204,359] 25.84215 7.728526
## Compost_(40,140] 20.31332 7.728526
## Freestall_(140,204] 26.25905 7.728526
## Freestall_(204,359] 23.06450 7.728526
## Freestall_(359,901] 26.02315 7.728526
## Freestall_(40,140] 16.96507 7.728526
## `Bedding type`Mattress `Bedding type`Other
## Compost_(140,204] 16.45905 7.353506
## Compost_(204,359] 16.45905 7.353506
## Compost_(40,140] 16.45905 7.353506
## Freestall_(140,204] 16.45905 7.353506
## Freestall_(204,359] 16.45905 7.353506
## Freestall_(359,901] 16.45905 7.353506
## Freestall_(40,140] 16.45905 7.353506
## `Slipperines feed bunk`Slightly slippery
## Compost_(140,204] 6.934
## Compost_(204,359] 6.934
## Compost_(40,140] 6.934
## Freestall_(140,204] 6.934
## Freestall_(204,359] 6.934
## Freestall_(359,901] 6.934
## Freestall_(40,140] 6.934
## `Slipperines feed bunk`Very slippery
## Compost_(140,204] 12.55588
## Compost_(204,359] 12.55588
## Compost_(40,140] 12.55588
## Freestall_(140,204] 12.55588
## Freestall_(204,359] 12.55588
## Freestall_(359,901] 12.55588
## Freestall_(40,140] 12.55588
## `Dry period length`[T.Less than sixty days]
## Compost_(140,204] 5.731921
## Compost_(204,359] 5.731921
## Compost_(40,140] 5.731921
## Freestall_(140,204] 5.731921
## Freestall_(204,359] 5.731921
## Freestall_(359,901] 5.731921
## Freestall_(40,140] 5.731921
##
## attr(,"class")
## [1] "coef.mer"Variance ratio is around 50%
22.71/52.48## [1] 0.43273631/0.4327363## [1] 2.310876Description of the variables included in the dataframe:
- Scored cows = all inspected cows by farm (milking cows)
- LS_1, 2, 3, 4 and 5 = number of inspectec cows classified by each locomotion scoring category (1-2: non-lame cows; 3-4-5 = clinically lame; 4-5 severely lame)
- Clinical lameness = prevalence of clinical lameness (%)
- Severe lameness = prevalence of severe lameness (%)
- Average milk yield per cow = farmer report on average milk yield per cow (kg)
- Total cows = total number of cows within each herd
- Workers number = number of persons working on each farm
- Farm type = housing system (freestall or compost-bedded pack farms)
- Bedding type = type of bedding predominantly used on each farm
- Dry period length = category of dry period length criteria established by each farmer
- Alley scrapping type = manure handling system
- Farm type herd size = category including each farm housing type (compost and freestall farms) by lactating herd size (categorized by quartils)
- Barn building year = year in which the barn was built
- Feeding type = feeding practices adopted by each farm (TMR, top-dressing)
- Heat control system type = barn ventilation - farm heat abatement devices
- Yearsof milk activity = years that farmers have been engaged in dairy activity
- Feeding times per day = number of times that feed was delivered to the cows (daily)
- Pasture access to milking cows = farmers’ report on whether the cows are allowed access to pasture
- Pasture access to cows = category of cows generally allowed access to pasture
- Milkings per day = number of milkings per day
- Stall cleaning = number of times the barn is cleaned (daily)
- Footbath = farmers’ report on the regular use of footbatHerd
- Footbath solution = solution used in the footbath
- Preventive hoof trimming = farmers’ report on the regular application of foot trimming in cows
- Hoof trimming freq = farmers’ report on the frequency of foot trimming application in cows
- Track slope = slope of alley traffic (category; track = alley traffic to milking parlour)
- Track covering = cover (roof) of alley traffic (category)
- Track pavement type = floor type of alley traffic (category)
- Track moisture = alley traffic floor moisture
- Feed bunk flooring = floor type of feed bunk
- Feed bunk type = type of feed bunk
- Holding area water trough = presence of water trough at the holding area
- Holding area flooring = floor type of holding area
- Holding area roof cover = presence of a roof on holding area
- Holding area Area = area (squared meters) of holding area
- Do cows fall on your farm = farmers’ answer to question “Do the cows usually fall at the barn?”
- Cows falls category = farmers’ answer to question “Do the cows usually fall at the barn?” (categorized variable)
- Feed bunk flooring category = floor type of feed bunk (categorized)
- Slipperines feed bunk = feed bunk alley floor slipperiness score
- Herd size = lactating herd size. Categorized by quartils