library(readr)
library(car)
library(dplyr)
library(tidyr)
library(lme4)
library(epitools)
library(lattice)
library(blmeco)
library(ggplot2)
library(lmerTest)
library(performance)
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 (MY2) is a data subset with information about the association between lameness and test-day milk yield (kg/day) on 16 farms and 5301 cows with available data. Variables associated with milk yield were included in the analyses to verify whether lameness remained significant when these variables were taken into account.

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)

MY2 = read_csv("Cowlevel.csv")

names(MY2)

##  [1] "ID"                           "farm_ID"                     
##  [3] "Cow"                          "Breed"                       
##  [5] "Housing system"               "Group_cat"                   
##  [7] "Pen"                          "Scored cows per farm"        
##  [9] "MilkProduction"               "Lameness"                    
## [11] "Gait score"                   "BCS_dic"                     
## [13] "BCS"                          "Parity"                      
## [15] "Parity_"                      "Lactation"                   
## [17] "DIM"                          "DIM_Cat"                     
## [19] "Report_date"                  "Clinical lameness prevalence"
## [21] "Clin_lam"                     "Severe lameness prevalence"  
## [23] "Sev_lam"

Description of the variables:

“farm_ID”: identification for each visited farm

“Cow”: identification of each cow that was examined

“Breed”: breed of each cow that was examined

“Housing system”: housing type for each visited farm

“Group_cat”: identification of each group of cows that was clustered by pens. The pens were nested within each farm

“Pen”: categories for groups of cows that were housed within each pen

“MilkProduction”: test-day milk yield (kg/day) obtained from farmers and from a regional dairy herd improvement association database

“Gait score”: a five-point scale visual score to classify the cows as sound (score 1 and 2), clinically lame (score ≥ 3) or severely lame (scores 4 and 5)

“BCS”: body condition score for each cow. A five-point categorical scale with 0.25 increments

“BCS_dic”: a category for BCS (low or high)

“Lameness”: a category for gait score defined as follows: a) non-lame cows (gait score = 1 and 2) b) moderately lame cows (gait score = 3) c) severely lame cows (gait score 4 and 5)

“Scored cows per farm”: number of cows that were examined in each farm (gait score and body condition score (BCS))

“Lactation”: lactation number for each cow

“Parity”: a category for lactation number

“Parity_”: another category for lactation number (character)

“DIM”: days in milk for each cow

“DIM_cat”: a category for days in milk

“Report_date”: date of the test-day milk report

“Clinical lameness prevalence”: herd-level lameness prevalence (gait score equal or greater than 3)

“Severe lameness prevalence”: herd-level lameness prevalence (gait score equal or greater than 4)

“Sev_lam”: a colum indicating which cow was classified as severely lame (gait score 4 and 5 = 1) or non-severely lame (gait score 1,2 or 3 = 0)

“Clin_lam”: a colum indicating which cow was classified as clinically lame (gait score 3, 4 and 5 = 1) or non-lame (gait score 1 and 2 = 0)

Data description

Farms

Sixteen farms (“farm_ID”) and cows (n) with available data on cow-level variables considered in the analyses.

There are missing data for some of the scored cows (cows assessed for gait score and BCS)

count(MY2, farm_ID, `Scored cows per farm`)

farm_ID	Scored cows per farm	n
7	42	39
13	432	368
16	422	410
18	680	623
20	599	562
23	836	736
26	147	113
28	901	759
30	339	286
31	222	200
32	371	344
34	193	154
39	109	87
40	335	328
45	294	271
47	348	249

Test-day milk yield (kg/day)

Description

summary(MY2$MilkProduction)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3.50   28.10   35.60   35.71   43.20   84.20

Test-day milk yield by parity

MY2%>%
  group_by(Parity_) %>%
  summarise(`Mean test-day milk` = mean(MilkProduction),
            SD = sd(MilkProduction),
            Min = min(MilkProduction),
            Max = max(MilkProduction),
            Median = median(MilkProduction),
            Observations = n())

Parity_	Mean test-day milk	SD	Min	Max	Median	Observations
First	33.65360	8.539462	4.2	56.8	34.1	2084
Second	36.57227	11.187714	4.0	69.9	37.0	1482
Third or greater	37.25186	13.067101	3.5	84.2	37.3	1963

Test-day milk yield and lameness status

MY2%>%
  group_by(Lameness) %>%
  summarise(`Mean test-day milk` = mean(MilkProduction),
            SD = sd(MilkProduction),
            Min = min(MilkProduction),
            Max = max(MilkProduction),
            Median = median(MilkProduction),
            Observations = n())

Lameness	Mean test-day milk	SD	Min	Max	Median	Observations
Moderate	36.46029	12.06568	3.5	84.2	36.00	1103
Non-lame	35.88852	10.46191	4.2	76.2	35.95	3372
Severe	34.37173	12.12671	4.0	69.9	34.40	1054

Test-day milk yield and categories for days in milk

MY2%>%
  group_by(DIM_Cat) %>%
  summarise(`Mean test-day milk` = mean(MilkProduction),
            SD = sd(MilkProduction),
            Min = min(MilkProduction),
            Max = max(MilkProduction),
            Median = median(MilkProduction),
            Observations = n())

DIM_Cat	Mean test-day milk	SD	Min	Max	Median	Observations
(1,120]	41.38000	10.825074	8.6	84.2	41.3	1700
(120,175]	39.29850	9.348385	7.0	61.8	39.2	802
(175,230]	36.81962	9.065520	8.3	62.7	36.8	780
(230,280]	33.31505	8.725969	7.8	57.8	33.2	651
(280,335]	30.26259	8.830857	4.2	63.7	29.8	532
(335,600]	26.14007	8.842411	3.5	51.8	26.4	836
NA	31.73684	11.859811	5.6	65.6	32.0	228

Test-day milk yield and body condition score

MY2%>%
  group_by(BCS_dic) %>%
  summarise(`Mean test-day milk` = mean(MilkProduction),
            SD = sd(MilkProduction),
            Min = min(MilkProduction),
            Max = max(MilkProduction),
            Median = median(MilkProduction),
            Observations = n())

BCS_dic	Mean test-day milk	SD	Min	Max	Median	Observations
High BCS (3-4.5)	34.72883	11.09245	3.5	84.2	34.7	3059
Low BCS (2-2.75)	36.93283	11.09544	5.6	73.4	37.2	2470

Test-day milk yield by farm and group within farm

Milk_per_farm = MY2%>%
  group_by(farm_ID, Group_cat)%>%
  summarise(`Mean test-day milk yield` = mean(MilkProduction),
            SD = sd(MilkProduction),
            Min = min(MilkProduction),
            Max = max(MilkProduction),
            Median = median(MilkProduction),
            Cows = n())
Milk_per_farm

farm_ID	Group_cat	Mean test-day milk yield	SD	Min	Max	Median	Cows
7	B	32.25128	9.3145860	14.0	57.0	32.60	39
13	B	35.55699	6.4073904	18.8	47.6	36.60	93
13	C	38.44754	8.6105856	14.4	59.2	38.40	122
13	D	31.85902	9.9687073	4.6	49.0	31.20	61
13	E	24.45161	10.4305312	6.2	47.2	23.40	62
13	F	25.49412	10.7350169	4.2	39.8	30.20	17
13	G	25.80000	10.5337553	10.0	44.6	28.60	13
16	B	28.76250	6.9910300	20.6	41.6	28.45	16
16	C	29.70667	6.1181468	14.3	41.9	30.05	120
16	D	40.58873	6.1168981	26.6	57.1	40.00	71
16	E	33.25634	6.1286851	11.4	57.3	32.80	71
16	F	18.83214	6.5567820	3.5	33.7	19.55	112
16	G	18.58000	6.2674262	7.9	32.8	18.10	20
18	B	40.20556	7.7218408	23.9	53.9	40.55	36
18	C	40.58820	5.9440661	25.0	55.6	40.80	161
18	D	46.62767	8.4075758	10.7	64.4	47.80	159
18	E	34.11772	6.8219533	13.2	52.8	34.50	158
18	F	23.67538	7.7442323	12.4	55.7	22.60	65
18	G	32.74815	11.2152153	12.4	58.1	29.30	27
18	H	30.08824	6.3122780	16.8	46.7	29.40	17
20	B	40.19778	5.6485054	24.7	50.1	40.20	90
20	C	34.69167	7.8679648	14.3	58.3	35.55	84
20	D	46.05341	8.5508076	20.1	65.6	46.35	88
20	E	47.25000	8.6355322	24.0	67.1	46.85	90
20	F	42.12500	8.3057060	19.9	60.6	41.15	92
20	G	25.55750	6.8093713	12.9	59.3	25.40	80
20	H	40.33684	9.9195527	20.6	60.9	38.60	38
23	A	40.17778	9.4035602	27.3	58.2	35.70	9
23	B	38.88291	5.5663442	20.8	52.4	39.60	117
23	C	47.80676	8.3929411	29.3	66.3	47.60	74
23	D	50.11053	8.7801227	31.5	66.6	50.25	76
23	E	49.17965	8.8609916	28.0	76.2	48.60	113
23	F	35.12569	6.0153262	20.9	52.6	35.20	109
23	G	34.39224	7.2176613	10.7	64.0	35.15	116
23	H	45.79375	12.7630499	22.5	69.9	48.20	32
23	I	32.50833	9.3113282	18.0	46.9	31.40	12
23	J	33.39130	11.0046235	10.0	59.8	32.50	23
23	K	25.05091	6.1786332	10.7	38.9	24.90	55
26	B	33.07619	7.0947096	23.3	48.3	32.60	21
26	C	35.99722	7.8720658	19.2	53.0	37.20	36
26	D	23.24545	8.6358976	17.5	42.0	19.10	11
26	E	35.78800	9.1705289	19.2	53.6	36.00	25
26	F	29.28000	6.0477573	19.5	42.3	29.60	20
28	B	39.28852	7.9512192	13.6	56.4	40.00	209
28	C	43.70993	11.4145040	12.7	69.0	44.80	141
28	D	41.14671	10.1233971	10.4	61.1	42.05	152
28	E	32.13721	9.8496278	4.2	58.8	31.40	129
28	F	35.40952	15.7029781	7.0	84.2	32.35	84
28	G	46.69310	9.8001797	19.7	62.5	43.90	29
28	H	37.80667	10.2668305	25.2	61.1	34.90	15
30	B	29.34865	6.9754019	18.8	55.8	28.40	74
30	C	33.83699	7.6950472	18.8	52.0	33.00	73
30	D	19.76835	4.9831651	7.4	31.8	20.40	79
30	E	22.41351	7.8776111	8.2	40.2	22.40	37
30	F	14.28696	4.8926015	5.8	25.8	14.60	23
31	B	42.48000	7.2182485	24.4	58.8	42.20	55
31	C	23.10698	4.9706470	10.0	33.6	24.00	43
31	D	33.48315	7.5067732	8.6	52.2	33.60	89
31	E	31.20000	11.1424115	7.0	49.2	34.40	13
32	B	27.95000	0.7778175	27.4	28.5	27.95	2
32	C	38.21148	7.2513240	19.9	51.9	39.20	61
32	D	32.61071	6.9720658	18.7	48.1	32.95	56
32	E	48.32407	7.8129533	31.3	67.5	48.10	54
32	F	47.29848	8.3123355	22.0	61.0	48.60	66
32	G	30.07500	7.6505835	14.1	51.5	29.50	64
32	H	31.31071	9.0403761	10.9	45.5	31.40	28
32	I	37.03846	13.5444047	12.6	65.7	36.90	13
34	A	36.56250	6.9021606	27.2	45.5	35.85	8
34	B	25.93636	4.4638127	16.5	37.3	25.40	44
34	C	32.06444	5.3079338	22.4	50.1	31.10	45
34	D	31.75152	8.5817802	14.7	56.6	31.80	33
34	E	23.94118	7.1179403	14.9	42.6	23.50	17
34	F	28.20000	7.9688979	18.8	41.8	27.00	7
39	B	34.82143	6.5130000	15.4	50.4	34.60	56
39	C	31.68387	8.3702687	17.9	52.9	29.90	31
40	B	36.00299	6.1283240	15.9	52.2	36.60	67
40	C	46.37727	6.1570019	36.7	63.4	44.85	88
40	D	40.01791	11.5642333	9.6	61.4	40.70	67
40	E	28.90000	6.0289461	16.2	47.0	28.95	106
45	B	27.76226	5.4790239	12.4	37.2	27.80	53
45	C	33.62679	5.5199635	20.5	46.6	34.10	56
45	D	39.50759	7.4290969	18.7	60.3	39.70	79
45	E	24.88154	9.7281147	5.4	50.6	25.30	65
45	F	34.29231	8.7132525	20.3	50.9	35.70	13
45	G	38.84000	11.6528108	23.1	50.7	42.00	5
47	B	29.76585	6.4644559	11.0	47.2	30.50	82
47	C	29.56792	5.3372340	15.1	40.8	29.60	53
47	D	27.18077	11.3215574	12.9	63.7	23.80	52
47	E	33.24400	11.1649261	9.3	50.2	35.40	50
47	F	30.52500	6.3183895	18.0	41.0	31.50	12

Visualizing the data structure: cows housed in pens within farms

ggplot(Milk_per_farm, aes(x=Group_cat, y = `Mean test-day milk yield`, col= Cows))+
  geom_point()+
  facet_wrap( ~ farm_ID)+
  scale_colour_gradientn(colours = c("yellow", "orange", "red", "darkred","brown"))+
  labs(x = "Groups within farms", y = "Mean test-day milk yield (kg/d)")

Days in milk

The variable was previously modified: typing errors, abnormal data, or cows with days in milk equal to zero were recoded as “missing data” (n = 169)

summary(MY2$DIM)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     6.0    99.0   188.0   207.1   287.0   782.0     169

The variable was recoded again in order to rule out extreme values (59): a superior limit of 600 days was set. A category with reference to the lactation curve was defined

count(MY2, DIM_Cat)

DIM_Cat	n
(1,120]	1700
(120,175]	802
(175,230]	780
(230,280]	651
(280,335]	532
(335,600]	836
NA	228

Correlations between variables

MY2%>%
  select_if(is.numeric)%>%
  cor(method = "spearman", use = "complete.obs")

##                                       ID      farm_ID         Cow
## ID                            1.00000000 -0.395973598  0.99833232
## farm_ID                      -0.39597360  1.000000000 -0.39654545
## Cow                           0.99833232 -0.396545455  1.00000000
## Scored cows per farm          0.71258574 -0.426868882  0.71161345
## MilkProduction                0.23503908 -0.066905996  0.23633851
## Gait score                   -0.14516510  0.005461752 -0.14410029
## BCS                           0.07365849 -0.024794338  0.07222020
## Lactation                    -0.22608927  0.037227360 -0.22540424
## DIM                          -0.04107659 -0.008075889 -0.04177529
## Clinical lameness prevalence -0.07028257  0.147638527 -0.07016966
## Clin_lam                     -0.13128360 -0.003359945 -0.13026869
## Severe lameness prevalence   -0.24082300  0.247531355 -0.24039268
## Sev_lam                      -0.13639718  0.029254132 -0.13541876
##                              Scored cows per farm MilkProduction
## ID                                    0.712585740     0.23503908
## farm_ID                              -0.426868882    -0.06690600
## Cow                                   0.711613453     0.23633851
## Scored cows per farm                  1.000000000     0.27507198
## MilkProduction                        0.275071983     1.00000000
## Gait score                           -0.013344723    -0.03474112
## BCS                                   0.120441122    -0.12594529
## Lactation                             0.046055896     0.11736535
## DIM                                   0.064426083    -0.48160060
## Clinical lameness prevalence         -0.054214985    -0.02486358
## Clin_lam                             -0.001017819    -0.02026971
## Severe lameness prevalence           -0.207740460    -0.08280501
## Sev_lam                              -0.038847637    -0.05531124
##                                Gait score         BCS   Lactation
## ID                           -0.145165097  0.07365849 -0.22608927
## farm_ID                       0.005461752 -0.02479434  0.03722736
## Cow                          -0.144100293  0.07222020 -0.22540424
## Scored cows per farm         -0.013344723  0.12044112  0.04605590
## MilkProduction               -0.034741117 -0.12594529  0.11736535
## Gait score                    1.000000000 -0.20629079  0.47689512
## BCS                          -0.206290789  1.00000000 -0.05949765
## Lactation                     0.476895119 -0.05949765  1.00000000
## DIM                           0.063013828  0.22087987 -0.02567865
## Clinical lameness prevalence  0.124601976 -0.12535933  0.02758264
## Clin_lam                      0.968561772 -0.17871548  0.46301216
## Severe lameness prevalence    0.116017769 -0.07061631 -0.03809471
## Sev_lam                       0.782203241 -0.20285212  0.37557027
##                                       DIM Clinical lameness prevalence
## ID                           -0.041076585                  -0.07028257
## farm_ID                      -0.008075889                   0.14763853
## Cow                          -0.041775288                  -0.07016966
## Scored cows per farm          0.064426083                  -0.05421499
## MilkProduction               -0.481600597                  -0.02486358
## Gait score                    0.063013828                   0.12460198
## BCS                           0.220879872                  -0.12535933
## Lactation                    -0.025678648                   0.02758264
## DIM                           1.000000000                   0.01617390
## Clinical lameness prevalence  0.016173901                   1.00000000
## Clin_lam                      0.055889493                   0.11959973
## Severe lameness prevalence    0.012653271                   0.62899381
## Sev_lam                       0.064381615                   0.10245986
##                                  Clin_lam Severe lameness prevalence
## ID                           -0.131283595                -0.24082300
## farm_ID                      -0.003359945                 0.24753135
## Cow                          -0.130268690                -0.24039268
## Scored cows per farm         -0.001017819                -0.20774046
## MilkProduction               -0.020269710                -0.08280501
## Gait score                    0.968561772                 0.11601777
## BCS                          -0.178715484                -0.07061631
## Lactation                     0.463012156                -0.03809471
## DIM                           0.055889493                 0.01265327
## Clinical lameness prevalence  0.119599730                 0.62899381
## Clin_lam                      1.000000000                 0.09700784
## Severe lameness prevalence    0.097007840                 1.00000000
## Sev_lam                       0.609936534                 0.13245978
##                                  Sev_lam
## ID                           -0.13639718
## farm_ID                       0.02925413
## Cow                          -0.13541876
## Scored cows per farm         -0.03884764
## MilkProduction               -0.05531124
## Gait score                    0.78220324
## BCS                          -0.20285212
## Lactation                     0.37557027
## DIM                           0.06438161
## Clinical lameness prevalence  0.10245986
## Clin_lam                      0.60993653
## Severe lameness prevalence    0.13245978
## Sev_lam                       1.00000000

MY2 %>% 
  select_if(is.numeric)%>% 
  cor(method = "spearman", use = "complete.obs")%>% 
  corrplot(method = "circle", tl.cex = 0.8, tl.srt = 45, tl.col = "black", type = "upper", order = "hclust", diag = FALSE)

Building a mixed-effects linear model to investigate associations between lameness and test-day milk yield

Mixed-effects linear regressions were built to investigate the associations between test-day milk yield and lameness. Additional variables were considered given their influence on milk production. Farm and group within farm were used as random effects.

Univariable associations were tested previously.

The multivariable model was tested using all the available data (incuding complete observations for days in milk) and also biologically plausible interactions between variables were checked. The model including complete observations for days in milk was similar to the model we reported (fmm). We observed an interaction between days in milk and parity. The interaction did not affect the coefficients for the association between lameness and milk yield, thus we reported the model without interaction.

Influence measures were checked. Removing the points with higher influence (cooks’ distance, leverage) did not change the coefficients in the model.

Contrasts for lameness status:

MY2$Lameness <- as.factor(MY2$Lameness)
contrasts(MY2$Lameness) <- contr.Treatment(levels(MY2$Lameness), base = 2)
contrasts(MY2$Lameness)

##          [T.Moderate] [T.Severe]
## Moderate            1          0
## Non-lame            0          0
## Severe              0          1

First, excluding missing data

MY2 = subset(MY2, !is.na(DIM_Cat))

Final multivariable model

summary(fmm<-lmer(MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1|farm_ID/Group_cat), data = MY2))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1 |  
##     farm_ID/Group_cat)
##    Data: MY2
## 
## REML criterion at convergence: 36577.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9193 -0.5829  0.0275  0.6208  5.5802 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  Group_cat:farm_ID (Intercept) 24.68    4.968   
##  farm_ID           (Intercept) 14.45    3.802   
##  Residual                      55.04    7.419   
## Number of obs: 5301, groups:  Group_cat:farm_ID, 90; farm_ID, 16
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)               36.3021     1.1608   18.5990  31.272  < 2e-16
## Lameness[T.Moderate]      -0.3587     0.2856 5216.9551  -1.256  0.20912
## Lameness[T.Severe]        -1.3251     0.3231 5241.3409  -4.101 4.17e-05
## BCS_dicLow BCS (2-2.75)    0.4697     0.2276 5238.3477   2.063  0.03912
## Parity(1,2]                0.9569     0.3644 5091.9684   2.626  0.00868
## Parity(2,3]                2.6586     0.4110 5135.6108   6.468 1.08e-10
## Parity(3,11]               1.3359     0.4397 5153.1700   3.038  0.00239
## DIM_Cat(120,175]          -1.5061     0.3390 5258.6851  -4.443 9.07e-06
## DIM_Cat(175,230]          -2.9757     0.3483 5259.1964  -8.542  < 2e-16
## DIM_Cat(230,280]          -5.3194     0.3774 5266.1063 -14.094  < 2e-16
## DIM_Cat(280,335]          -7.2198     0.4073 5267.8019 -17.724  < 2e-16
## DIM_Cat(335,600]         -10.0601     0.3814 5278.9236 -26.377  < 2e-16
##                            
## (Intercept)             ***
## Lameness[T.Moderate]       
## Lameness[T.Severe]      ***
## BCS_dicLow BCS (2-2.75) *  
## Parity(1,2]             ** 
## Parity(2,3]             ***
## Parity(3,11]            ** 
## DIM_Cat(120,175]        ***
## DIM_Cat(175,230]        ***
## DIM_Cat(230,280]        ***
## DIM_Cat(280,335]        ***
## DIM_Cat(335,600]        ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) L[T.M] L[T.S] BCS_B( P(1,2] P(2,3] P(3,11 DIM_C(12
## Lmnss[T.Md] -0.022                                                   
## Lmnss[T.Sv]  0.003  0.374                                            
## BCS_LBCS(2- -0.096 -0.108 -0.217                                     
## Parity(1,2] -0.175 -0.087 -0.116  0.005                              
## Parity(2,3] -0.160 -0.159 -0.215  0.017  0.663                       
## Party(3,11] -0.153 -0.188 -0.348  0.017  0.636  0.649                
## DIM_C(120,1 -0.112 -0.037 -0.051  0.054  0.026  0.025  0.034         
## DIM_C(175,2 -0.119  0.007 -0.036  0.052  0.044  0.044  0.045  0.381  
## DIM_C(230,2 -0.120 -0.020 -0.054  0.072  0.059  0.063  0.065  0.354  
## DIM_C(280,3 -0.117 -0.017 -0.038  0.070  0.064  0.058  0.063  0.319  
## DIM_C(335,6 -0.140 -0.059 -0.116  0.123  0.101  0.102  0.108  0.349  
##             DIM_C(17 DIM_C(23 DIM_C(28
## Lmnss[T.Md]                           
## Lmnss[T.Sv]                           
## BCS_LBCS(2-                           
## Parity(1,2]                           
## Parity(2,3]                           
## Party(3,11]                           
## DIM_C(120,1                           
## DIM_C(175,2                           
## DIM_C(230,2  0.392                    
## DIM_C(280,3  0.359    0.375           
## DIM_C(335,6  0.394    0.421    0.421

confint(fmm, method = "Wald")

##                                2.5 %     97.5 %
## .sig01                            NA         NA
## .sig02                            NA         NA
## .sigma                            NA         NA
## (Intercept)              34.02693528 38.5773559
## Lameness[T.Moderate]     -0.91843162  0.2009893
## Lameness[T.Severe]       -1.95839319 -0.6918439
## BCS_dicLow BCS (2-2.75)   0.02356384  0.9159038
## Parity(1,2]               0.24258222  1.6711629
## Parity(2,3]               1.85303795  3.4641617
## Parity(3,11]              0.47412782  2.1976258
## DIM_Cat(120,175]         -2.17053346 -0.8416276
## DIM_Cat(175,230]         -3.65844458 -2.2929531
## DIM_Cat(230,280]         -6.05918384 -4.5796857
## DIM_Cat(280,335]         -8.01814105 -6.4213848
## DIM_Cat(335,600]        -10.80758923 -9.3125631

Intercept_by_group<-coef(fmm)$`Group_cat:farm_ID`

Intercept_by_farm<-coef(fmm)$`farm_ID`

Intra Class Correlation coefficient: pen within farm

24.68/(24.68 + 14.45 + 55.04)*100

## [1] 26.20792

Intra Class Correlation coefficient: farm

14.45/(24.68 + 14.45 + 55.04)*100

## [1] 15.34459

Checking collinearity

vif(fmm)

##              GVIF Df GVIF^(1/(2*Df))
## Lameness 1.248251  2        1.057001
## BCS_dic  1.068291  1        1.033581
## Parity   1.194675  3        1.030090
## DIM_Cat  1.039160  5        1.003849

Including residuals, influence and fitted values in the table:

MY2$Residuals <-resid(fmm)
MY2$Pearson_residuals<-resid(fmm, type = "pearson")
MY2$Fitted_values <-fitted(fmm)
MY2$Cook_dist <- cooks.distance(fmm)
MY2$Hat_values <-hatvalues(fmm)

Residuals: distribution

ggplot(MY2, aes(Residuals))+
  geom_histogram(aes(y =..density..), colour="black", fill="white")+
  geom_density(alpha = 0.3, fill="#FF4499")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Boxplot(resid(fmm))

##  [1] 3490 3292 2580 3932 1837 3380 4501 2313 5199 2490 4077  379 1981 3895
## [15] 1336 4023 3611 2367 5169 4329

There are some large residuals, however, given the high number of points these values may have little influence on the coefficients.

ggplot(MY2, aes(Fitted_values, Pearson_residuals))+
  geom_jitter(size = 0.7)+
  geom_smooth(method = "lm")+
  theme_bw()

scatter.smooth(resid(fmm, type = "pearson")~fitted(fmm), col = "grey", las = 1)

Random effect: farm

qqPlot(ranef(fmm)$farm_ID[,1])

## [1] 9 8

Random effect: group

qqPlot(ranef(fmm)$Group_cat[,1])

## [1] 70 36

library(sjPlot)

plot_model(fmm, type = c("re"))

## [[1]]

## 
## [[2]]

Influence

which.max(cooks.distance(fmm))

## 886 
## 886

sort(cooks.distance(fmm), decreasing = T)[1:10]

##        886       2367       2707        555        540        499 
## 0.08313184 0.08247000 0.08143741 0.05819039 0.05383365 0.05260745 
##        379       3570        609       4077 
## 0.05162768 0.04842686 0.04483035 0.04362315

which(MY2$Cook_dist >0.03559153)

##  [1]  379  499  540  555  609  886 1337 2367 2707 3570 4077

plot(cooks.distance(fmm), type= "h", main = "Cook's distance",
ylab="D", xlab = "Observation number", las=1 )

Leverage

sort(hatvalues(fmm), decreasing = TRUE)[1:11]

##      2639      2790        76       499       632       555       592 
## 0.2515117 0.2513768 0.1468069 0.1466749 0.1465257 0.1464997 0.1455317 
##      1722       858      1475      1324 
## 0.1154274 0.1151666 0.1146265 0.1139543

which(hatvalues(fmm) > 0.2)

## 2639 2790 
## 2639 2790

plot(hatvalues(fmm), type= "h", main = "Hat",
ylab="", xlab = "Observation number", las=1 )

influencePlot(fmm)

	StudRes	Hat	CookD
379	5.2924885	0.0216393	0.0516277
886	-3.6561694	0.0694446	0.0831318
2367	3.6799800	0.0681012	0.0824700
2639	-1.0211746	0.2515117	0.0292006
2790	-0.7765471	0.2513768	0.0168739
4077	5.6261834	0.0162685	0.0436231

summary(fmm_lev1<-lmer(MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1|farm_ID/Group_cat), data = subset(MY2 [-379,])))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1 |  
##     farm_ID/Group_cat)
##    Data: subset(MY2[-379, ])
## 
## REML criterion at convergence: 36543.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9357 -0.5853  0.0289  0.6217  5.5952 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  Group_cat:farm_ID (Intercept) 24.73    4.972   
##  farm_ID           (Intercept) 14.49    3.807   
##  Residual                      54.76    7.400   
## Number of obs: 5300, groups:  Group_cat:farm_ID, 90; farm_ID, 16
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)               36.3171     1.1620   18.5730  31.255  < 2e-16
## Lameness[T.Moderate]      -0.3578     0.2848 5215.8648  -1.256  0.20911
## Lameness[T.Severe]        -1.3528     0.3223 5240.1846  -4.197 2.75e-05
## BCS_dicLow BCS (2-2.75)    0.4552     0.2271 5237.1585   2.005  0.04503
## Parity(1,2]                0.9593     0.3635 5094.6656   2.639  0.00834
## Parity(2,3]                2.6670     0.4100 5137.6195   6.505 8.50e-11
## Parity(3,11]               1.3098     0.4386 5154.9707   2.986  0.00284
## DIM_Cat(120,175]          -1.5153     0.3381 5257.4574  -4.481 7.58e-06
## DIM_Cat(175,230]          -2.9835     0.3474 5257.9497  -8.587  < 2e-16
## DIM_Cat(230,280]          -5.3272     0.3765 5264.8552 -14.151  < 2e-16
## DIM_Cat(280,335]          -7.2950     0.4065 5266.4978 -17.944  < 2e-16
## DIM_Cat(335,600]         -10.0548     0.3804 5278.0693 -26.431  < 2e-16
##                            
## (Intercept)             ***
## Lameness[T.Moderate]       
## Lameness[T.Severe]      ***
## BCS_dicLow BCS (2-2.75) *  
## Parity(1,2]             ** 
## Parity(2,3]             ***
## Parity(3,11]            ** 
## DIM_Cat(120,175]        ***
## DIM_Cat(175,230]        ***
## DIM_Cat(230,280]        ***
## DIM_Cat(280,335]        ***
## DIM_Cat(335,600]        ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) L[T.M] L[T.S] BCS_B( P(1,2] P(2,3] P(3,11 DIM_C(12
## Lmnss[T.Md] -0.022                                                   
## Lmnss[T.Sv]  0.003  0.374                                            
## BCS_LBCS(2- -0.096 -0.108 -0.217                                     
## Parity(1,2] -0.174 -0.087 -0.116  0.005                              
## Parity(2,3] -0.160 -0.159 -0.215  0.017  0.663                       
## Party(3,11] -0.152 -0.188 -0.347  0.017  0.636  0.649                
## DIM_C(120,1 -0.111 -0.037 -0.051  0.054  0.026  0.025  0.034         
## DIM_C(175,2 -0.119  0.007 -0.036  0.052  0.044  0.044  0.045  0.381  
## DIM_C(230,2 -0.119 -0.020 -0.054  0.072  0.059  0.063  0.065  0.354  
## DIM_C(280,3 -0.117 -0.017 -0.037  0.070  0.064  0.058  0.063  0.319  
## DIM_C(335,6 -0.140 -0.059 -0.116  0.123  0.101  0.102  0.108  0.349  
##             DIM_C(17 DIM_C(23 DIM_C(28
## Lmnss[T.Md]                           
## Lmnss[T.Sv]                           
## BCS_LBCS(2-                           
## Parity(1,2]                           
## Parity(2,3]                           
## Party(3,11]                           
## DIM_C(120,1                           
## DIM_C(175,2                           
## DIM_C(230,2  0.392                    
## DIM_C(280,3  0.359    0.374           
## DIM_C(335,6  0.394    0.421    0.420

summary(fmm_lev2<-lmer(MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1|farm_ID/Group_cat), data = subset(MY2, Hat_values < 0.10)))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MilkProduction ~ Lameness + BCS_dic + Parity + DIM_Cat + (1 |  
##     farm_ID/Group_cat)
##    Data: subset(MY2, Hat_values < 0.1)
## 
## REML criterion at convergence: 36424.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9154 -0.5825  0.0312  0.6213  5.5807 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  Group_cat:farm_ID (Intercept) 24.81    4.981   
##  farm_ID           (Intercept) 15.29    3.911   
##  Residual                      55.06    7.420   
## Number of obs: 5279, groups:  Group_cat:farm_ID, 86; farm_ID, 16
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)               36.3081     1.1896   18.5039  30.522  < 2e-16
## Lameness[T.Moderate]      -0.3765     0.2858 5196.0878  -1.317  0.18788
## Lameness[T.Severe]        -1.3317     0.3237 5214.4688  -4.114 3.95e-05
## BCS_dicLow BCS (2-2.75)    0.4744     0.2281 5216.1055   2.080  0.03756
## Parity(1,2]                0.9176     0.3655 5077.2906   2.510  0.01210
## Parity(2,3]                2.6192     0.4121 5113.2587   6.355 2.26e-10
## Parity(3,11]               1.3177     0.4406 5130.4288   2.991  0.00279
## DIM_Cat(120,175]          -1.5167     0.3404 5237.1293  -4.455 8.55e-06
## DIM_Cat(175,230]          -2.9469     0.3492 5238.4712  -8.438  < 2e-16
## DIM_Cat(230,280]          -5.3003     0.3780 5245.0500 -14.022  < 2e-16
## DIM_Cat(280,335]          -7.1983     0.4081 5246.8226 -17.640  < 2e-16
## DIM_Cat(335,600]         -10.0433     0.3820 5255.8737 -26.290  < 2e-16
##                            
## (Intercept)             ***
## Lameness[T.Moderate]       
## Lameness[T.Severe]      ***
## BCS_dicLow BCS (2-2.75) *  
## Parity(1,2]             *  
## Parity(2,3]             ***
## Parity(3,11]            ** 
## DIM_Cat(120,175]        ***
## DIM_Cat(175,230]        ***
## DIM_Cat(230,280]        ***
## DIM_Cat(280,335]        ***
## DIM_Cat(335,600]        ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) L[T.M] L[T.S] BCS_B( P(1,2] P(2,3] P(3,11 DIM_C(12
## Lmnss[T.Md] -0.022                                                   
## Lmnss[T.Sv]  0.003  0.374                                            
## BCS_LBCS(2- -0.093 -0.108 -0.217                                     
## Parity(1,2] -0.170 -0.087 -0.115  0.005                              
## Parity(2,3] -0.156 -0.159 -0.213  0.018  0.663                       
## Party(3,11] -0.149 -0.188 -0.347  0.017  0.636  0.650                
## DIM_C(120,1 -0.109 -0.036 -0.050  0.054  0.028  0.026  0.035         
## DIM_C(175,2 -0.117  0.006 -0.037  0.052  0.045  0.044  0.045  0.381  
## DIM_C(230,2 -0.118 -0.019 -0.054  0.071  0.060  0.064  0.066  0.354  
## DIM_C(280,3 -0.116 -0.017 -0.038  0.070  0.064  0.060  0.064  0.319  
## DIM_C(335,6 -0.139 -0.059 -0.116  0.123  0.102  0.103  0.110  0.349  
##             DIM_C(17 DIM_C(23 DIM_C(28
## Lmnss[T.Md]                           
## Lmnss[T.Sv]                           
## BCS_LBCS(2-                           
## Parity(1,2]                           
## Parity(2,3]                           
## Party(3,11]                           
## DIM_C(120,1                           
## DIM_C(175,2                           
## DIM_C(230,2  0.393                    
## DIM_C(280,3  0.360    0.375           
## DIM_C(335,6  0.395    0.422    0.422

Residuals and explanatory variables considered in the model

ggplot(MY2, aes(Pearson_residuals , MilkProduction))+
  geom_point()+
  geom_smooth(method = "lm")

ggplot(MY2, aes(Lameness, Pearson_residuals))+
  geom_boxplot()

ggplot(MY2, aes(Lameness, Pearson_residuals))+
  geom_boxplot()

ggplot(MY2, aes(Parity, Pearson_residuals))+
  geom_boxplot()

Lameness prevalence in lactating cows housed in freestall and compost-bedded pack dairy farms in southern Brazil: associations between milk yield and lameness

José A. Bran, Joao H.C. Costa, Marina A.G. von Keyserlingk, Maria José Hötzel

Sep 2019