Epidemiology of bovine digital dermatitis in Australian dairy cows

Acknowledgements

Authors and affiliations

Sam Rowe1 Bill Tranter,2,3 Ash Phipps,4 Andrew McPherson,1 Richard Laven,5 John House,1 Ruth Zadoks,1

1Sydney School of Veterinary Science, The University of Sydney, Camden, New South Wales 2570, Australia

2College of Public Health, Medical and Veterinary Sciences, James Cook University, Townsville, Queensland, Australia

3Tableland Veterinary Service, Malanda, Queensland, Australia

4Rochester Veterinary Practice, Rochester, Victoria, Australia

5Institute of Veterinary, Animal and Biomedical Sciences, Massey University, Palmerston North 4474, New Zealand

Thank you

Thank you to the farmers who volunteered to be involved in this study. Thank you to Aaron Yang (Nanjing Agricultural University) for sharing his expertise in Bayesian methods. Thanks to Luke Ingenhoff (University of Sydney), Peter Alexander (Bega and Cobargo Vet Hospitals) and Joey Rheinberger (Ironmines Veterinary Clinic) for helping with herd recruitment.

Import data and run packages

Show the code 
library(tidyverse)
library(janitor)
library(readr)
library(prevalence)
library(purrr)
library(ggpubr)
library(ggplot2)
library(readxl)
library(psych)
library(irr)
library(DescTools)
library(readxl)
library(epiR)
library(broom)
library(emmeans)

herd <- read_excel("bdd herd data entry.xlsx", 
    sheet = "entry wide", col_types = c("text", 
        "skip", "date", "date", "skip", "text", 
        "numeric", "text", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "text", "text", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "text", "text", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "text", "numeric", "numeric", "text", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))

data <- readRDS("~/Dropbox/R backup/2021/BDD2021/herds.rds")

data <- merge(data,herd,
              by="herd")

data$herd <- NULL

Training and evaluation of assessors

Sam Rowe (SR), Bill Tranter (WPT) and Ash Phipps (AP) visited farms to score cows feet. We used images from a study Vanhoudt et al (2019) for training and validation. 40 Images were used for training and 41 images were used for evaluating agreement between us. We were blinded to each others scores for images in the testing dataset. We compared our scores against each other and to the most common score for authors in the Vanhoudt et al (2019) study (referred to as ‘experts’ below). We used Cohen’s kappa and Fleiss’ kappa.

Vanhoudt, A., et al. “Interobserver agreement of digital dermatitis M-scores for photographs of the hind feet of standing dairy cattle.” Journal of dairy science 102.6 (2019): 5466-5474.

Show the code 
Mscore <- read_excel("M-score_IOA-raw_data.xlsx")
Mscore <- Mscore %>% select(!scorer)
MscoreT <- t(Mscore) %>% as.data.frame

rownam <- rownames(MscoreT)
MscoreT$image <- rownam
MscoreT$image <- parse_number(MscoreT$image)
rownames(MscoreT) <- NULL

colnames(MscoreT) <- c("S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","Mode","Image")

Mscore <- MscoreT %>% select(Image,Mode,S1:S11)
Mscore$Mode <- Mscore$Mode %>% as.numeric
MscoreTrain <- Mscore[1:40,]


#Import data with SR responses

test <- read_excel("Train.xlsx", col_types = c("numeric", 
                                                "numeric", "numeric", "text", "skip", 
                                                "text", "text", "skip", "skip"))

test <- test %>% filter(Image > 40)
test <- test %>% filter(Exclude <1)
test <- test[!is.na(test$SR),]

test$SR1 <- 0
test$SR1[test$SR!="0"] <-1

test$WPT1 <- 0
test$WPT1[test$WPT!="0"] <-1

test$AP1 <- 0
test$AP1[test$AP!="0"] <-1

test$Mode1 <- 0
test$Mode1[test$Mode!="0"] <-1

test$SR <- as.numeric(test$SR) %>% round(1)
test$WPT <- as.numeric(test$WPT) %>% round(1)
test$AP <- as.numeric(test$AP) %>% round(1)

test <- test[!is.na(test$SR),]

Show agreement dataset

Show the code 
table_test <- test %>% select(Mode,SR,WPT,AP,Mode1,SR1,WPT1,AP1)
colnames(table_test) <- c("Experts 6 point","SR 6 point","WPT 6 point","AP 6 point","Experts (0/1)","SR (0/1)","WPT (0/1)","WPT (0/1")
DT::datatable(table_test)

Agreement (Fleiss’ Kappa) between SR, WPT and AP when classifying an image as healthy (M0) or BDD (M1, M2, M3, M4, M4.1)

Show the code 
kappam.fleiss
function (ratings, exact = FALSE, detail = FALSE) 
{
    ratings <- as.matrix(na.omit(ratings))
    ns <- nrow(ratings)
    nr <- ncol(ratings)
    lev <- levels(as.factor(ratings))
    for (i in 1:ns) {
        frow <- factor(ratings[i, ], levels = lev)
        if (i == 1) 
            ttab <- as.numeric(table(frow))
        else ttab <- rbind(ttab, as.numeric(table(frow)))
    }
    ttab <- matrix(ttab, nrow = ns)
    agreeP <- sum((apply(ttab^2, 1, sum) - nr)/(nr * (nr - 1))/ns)
    if (!exact) {
        method <- "Fleiss' Kappa for m Raters"
        chanceP <- sum(apply(ttab, 2, sum)^2)/(ns * nr)^2
    }
    else {
        method <- "Fleiss' Kappa for m Raters (exact value)"
        for (i in 1:nr) {
            rcol <- factor(ratings[, i], levels = lev)
            if (i == 1) 
                rtab <- as.numeric(table(rcol))
            else rtab <- rbind(rtab, as.numeric(table(rcol)))
        }
        rtab <- rtab/ns
        chanceP <- sum(apply(ttab, 2, sum)^2)/(ns * nr)^2 - sum(apply(rtab, 
            2, var) * (nr - 1)/nr)/(nr - 1)
    }
    value <- (agreeP - chanceP)/(1 - chanceP)
    if (!exact) {
        pj <- apply(ttab, 2, sum)/(ns * nr)
        qj <- 1 - pj
        varkappa <- (2/(sum(pj * qj)^2 * (ns * nr * (nr - 1)))) * 
            (sum(pj * qj)^2 - sum(pj * qj * (qj - pj)))
        SEkappa <- sqrt(varkappa)
        u <- value/SEkappa
        p.value <- 2 * (1 - pnorm(abs(u)))
        if (detail) {
            pj <- apply(ttab, 2, sum)/(ns * nr)
            pjk <- (apply(ttab^2, 2, sum) - ns * nr * pj)/(ns * 
                nr * (nr - 1) * pj)
            kappaK <- (pjk - pj)/(1 - pj)
            varkappaK <- 2/(ns * nr * (nr - 1))
            SEkappaK <- sqrt(varkappaK)
            uK <- kappaK/SEkappaK
            p.valueK <- 2 * (1 - pnorm(abs(uK)))
            tableK <- as.table(round(cbind(kappaK, uK, p.valueK), 
                digits = 3))
            rownames(tableK) <- lev
            colnames(tableK) <- c("Kappa", "z", "p.value")
        }
    }
    if (!exact) {
        if (!detail) {
            rval <- list(method = method, subjects = ns, raters = nr, 
                irr.name = "Kappa", value = value)
        }
        else {
            rval <- list(method = method, subjects = ns, raters = nr, 
                irr.name = "Kappa", value = value, detail = tableK)
        }
        rval <- c(rval, stat.name = "z", statistic = u, p.value = p.value)
    }
    else {
        rval <- list(method = method, subjects = ns, raters = nr, 
            irr.name = "Kappa", value = value)
    }
    class(rval) <- "irrlist"
    return(rval)
}
<bytecode: 0x7f7ae4974938>
<environment: namespace:irr>
Show the code 
all <- kappam.fleiss(cbind(test$SR1,test$WPT1,test$AP1),detail=TRUE)

Agreement (Fleiss’ Kappa) between SR, WPT and AP when using a 6-point scoring system (M0, M1, M2, M3, M4, M4.1)

Show the code 
kappam.fleiss(cbind(test$SR,test$WPT,test$AP))
 Fleiss' Kappa for m Raters

 Subjects = 41 
   Raters = 3 
    Kappa = 0.499 

        z = 9.64 
  p-value = 0

Agreement between each scorer and experts from Vanhoudt et al (2019) when using a 6-point scoring system (M0, M1, M2, M3, M4, M4.1)

Sam Rowe

Show the code 
cohen.kappa(x=cbind(test$Mode,test$SR))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa  0.15     0.35  0.55
weighted kappa    0.63     0.63  0.63

 Number of subjects = 41

Bill Tranter

Show the code 
cohen.kappa(x=cbind(test$Mode,test$WPT))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa  0.28     0.47  0.65
weighted kappa    0.40     0.60  0.79

 Number of subjects = 41

Ash Phipps

Show the code 
cohen.kappa(x=cbind(test$Mode,test$AP))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa  0.41     0.60  0.78
weighted kappa    0.76     0.76  0.76

 Number of subjects = 41

Agreement between each scorer and experts from Vanhoudt et al (2019) when classifying an image as healthy (M0) or BDD (M1, M2, M3, M4, M4.1)

Sam Rowe

Show the code 
cohen.kappa(x=cbind(test$Mode1,test$SR1))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa  0.64     0.88     1
weighted kappa    0.64     0.88     1

 Number of subjects = 41

Bill Tranter

Show the code 
cohen.kappa(x=cbind(test$Mode1,test$WPT1))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa     1        1     1
weighted kappa       1        1     1

 Number of subjects = 41

Ash Phipps

Show the code 
cohen.kappa(x=cbind(test$Mode1,test$AP1))
Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)

Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
                 lower estimate upper
unweighted kappa     1        1     1
weighted kappa       1        1     1

 Number of subjects = 41

Analysis to inform sampling strategy

Visit analysis log (published on Rpubs) to see how the sampling strategy was determined.

Description of BDD prevalence

Show the code 
data <- data %>% mutate(
  region = case_when(
    region == "bega valley" ~ "NSW: Bega Valley",
    region == "central west" ~ "NSW: Central West",
    region == "FNQ" ~ "QLD: Atherton Tablelands",
    region == "north coast" ~ "NSW: North Coast",
    region == "south coast" ~ "NSW: South Coast",
    region == "sydney" ~ "NSW: Greater Western Sydney",
    region == "vic" ~ "Victoria"
      ),
  
  closed_herd_cow = case_when(
    biosecurity_most_recent_introduction_lactating_cows == 99 ~ 1,
    biosecurity_most_recent_introduction_lactating_cows != 99 ~ 0,
    ),
  
  closed_herd_bulls = case_when(
    biosecurity_most_recent_introduction_bulls == 99 ~ 1,
    biosecurity_most_recent_introduction_bulls != 99 ~ 0
  ),
  
  closed_herd = case_when(
    closed_herd_cow == 1 & closed_herd_bulls == 1 ~ 1,
    closed_herd_cow != 1 | closed_herd_bulls != 1 ~ 0
  ),
  
  who_does_most_trimming = case_when(
    hoof_care_who_trims_farmer > 49 ~ "farmer",
    hoof_care_who_trims_trimmer > 49 ~ "trimmer",
    hoof_care_who_trims_vet > 49 ~ "vet"
  ),
  
  farmer_trims = case_when(
    who_does_most_trimming == "farmer" ~ 1,
    who_does_most_trimming != "farmer" ~ 0,
  ),
  
  vet_trims = case_when(
    who_does_most_trimming == "vet" ~ 1,
    who_does_most_trimming != "vet" ~ 0,
  ),
  
  trimmer_trims = case_when(
    who_does_most_trimming == "trimmer" ~ 1,
    who_does_most_trimming != "trimmer" ~ 0,
  ),
  
  mgmt = case_when(
    housing_lactating_indoor == 1 ~ "barn",
    housing_lactating_feed_pad == 1 & housing_lactating_indoor == 0 ~ "feed pad",
    !is.na(housing_lactating_feed_pad) ~ "grazing"
  ),
  
  mgmt_barn = case_when(
    mgmt == "barn" ~ 1,
    mgmt != "barn" ~ 0
  ),
  
  mgmt_feed_pad = case_when(
    mgmt == "feed pad" ~ 1,
    mgmt != "feed pad" ~ 0
  ),
  
mgmt_grazing = case_when(
    mgmt == "grazing" ~ 1,
    mgmt != "grazubg" ~ 0
  )
)

data <- data %>% mutate(
  mgmt = relevel(factor(mgmt), ref = "grazing"),
  prev_cat = relevel(factor(prev_cat), ref ="(0,0.1]"),
  bdd_lame = case_when(
    seen_bdd_cause_lameness == "definitely yes" |
      seen_bdd_cause_lameness == "probably yes" ~ 1,
    !is.na(seen_bdd) ~ 0
    
  )
)

# data %>% tabyl(seen_bdd_cause_lameness,
#                bdd_lame)

# data %>% tabyl(who_does_most_trimming, farmer_trims)
# data %>% tabyl(who_does_most_trimming, vet_trims)
# data %>% tabyl(who_does_most_trimming, trimmer_trims)


# data %>% tabyl(biosecurity_most_recent_introduction_lactating_cows,closed_herd_cow)
# data %>% tabyl(closed_herd)
# data %>% tabyl(closed_herd_cow)
# data %>% tabyl(closed_herd_bulls)


#write.csv(data,'herds.csv')
print(summarytools::dfSummary(data,valid.col=FALSE, graph.magnif=0.8, style="grid"),method = "render")

Data Frame Summary

data

Dimensions: 71 x 71
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 region [character]
1. NSW: Bega Valley
2. NSW: Central West
3. NSW: Greater Western Sydn
4. NSW: North Coast
5. NSW: South Coast
6. QLD: Atherton Tablelands
7. Victoria
5(7.0%)
6(8.5%)
5(7.0%)
5(7.0%)
12(16.9%)
18(25.4%)
20(28.2%)
 0 (0.0%)
2 herd_size [numeric]
Mean (sd) : 440.6 (953.1)
min ≤ med ≤ max:
58 ≤ 243 ≤ 7992
IQR (CV) : 209.5 (2.2)
 63 distinct values 0 (0.0%)
3 cows_scored [integer]
Mean (sd) : 222.7 (156.5)
min ≤ med ≤ max:
57 ≤ 172 ≤ 893
IQR (CV) : 129.5 (0.7)
 66 distinct values 0 (0.0%)
4 prev_cow [numeric]
Mean (sd) : 0.1 (0.1)
min ≤ med ≤ max:
0 ≤ 0.1 ≤ 0.5
IQR (CV) : 0.1 (0.9)
 68 distinct values 0 (0.0%)
5 prev_left [numeric]
Mean (sd) : 0.1 (0.1)
min ≤ med ≤ max:
0 ≤ 0.1 ≤ 0.4
IQR (CV) : 0.1 (1)
 65 distinct values 0 (0.0%)
6 prev_right [numeric]
Mean (sd) : 0.1 (0.1)
min ≤ med ≤ max:
0 ≤ 0.1 ≤ 0.5
IQR (CV) : 0.1 (1)
 67 distinct values 0 (0.0%)
7 prev_m1 [numeric]
Mean (sd) : 0 (0)
min ≤ med ≤ max:
0 ≤ 0 ≤ 0
IQR (CV) : 0 (3.9)
 6 distinct values 0 (0.0%)
8 prev_m2 [numeric]
Mean (sd) : 0 (0)
min ≤ med ≤ max:
0 ≤ 0 ≤ 0
IQR (CV) : 0 (3.5)
 8 distinct values 0 (0.0%)
9 prev_m3 [numeric]
Mean (sd) : 0 (0)
min ≤ med ≤ max:
0 ≤ 0 ≤ 0
IQR (CV) : 0 (4)
 7 distinct values 0 (0.0%)
10 prev_m4 [numeric]
Mean (sd) : 0.1 (0.1)
min ≤ med ≤ max:
0 ≤ 0.1 ≤ 0.5
IQR (CV) : 0.1 (0.9)
 68 distinct values 0 (0.0%)
11 prev_m4.1 [numeric]
Mean (sd) : 0 (0)
min ≤ med ≤ max:
0 ≤ 0 ≤ 0.1
IQR (CV) : 0 (2.3)
 28 distinct values 0 (0.0%)
12 prev_cat [factor]
1. (0,0.1]
2. (-Inf,0]
3. (0.1,0.2]
4. (0.2,0.3]
5. (0.3,1]
35(49.3%)
3(4.2%)
17(23.9%)
9(12.7%)
7(9.9%)
 0 (0.0%)
13 prev_m4_cat [factor]
1. (-Inf,0.1]
2. (0.1,0.2]
3. (0.2,0.3]
4. (0.3,1]
39(54.9%)
17(23.9%)
8(11.3%)
7(9.9%)
 0 (0.0%)
14 m2_any [numeric]
Min : 0
Mean : 0.1
Max : 1
0:64(90.1%)
1:7(9.9%)
 0 (0.0%)
15 date_milking_visit [POSIXct, POSIXt]
min : 2021-06-07
med : 2022-04-11 12:00:00
max : 2022-07-19
range : 1y 1m 12d
 50 distinct values 1 (1.4%)
16 date_biopsy_1 [POSIXct, POSIXt]
1. 2021-06-07
2. 2021-06-17
3. 2021-07-14
4. 2021-07-30
5. 2021-08-13
6. 2021-11-01
7. 2022-02-15
8. 2022-02-23
9. 2022-04-29
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
1(11.1%)
 62 (87.3%)
17 breed [character]
1. ayrshire
2. brown swiss
3. hf
4. hf x j
5. illawarra
6. jersey
7. other
1(1.4%)
1(1.4%)
57(80.3%)
6(8.5%)
2(2.8%)
2(2.8%)
2(2.8%)
 0 (0.0%)
18 daily_milk [numeric]
Mean (sd) : 23.7 (5.3)
min ≤ med ≤ max:
15 ≤ 23 ≤ 42
IQR (CV) : 6 (0.2)
 22 distinct values 0 (0.0%)
19 calving_pattern [character]
1. seasonal
2. split
3. year round
4(5.6%)
15(21.1%)
52(73.2%)
 0 (0.0%)
20 lactating_diet_pasture [numeric]
Min : 0
Mean : 0.9
Max : 1
0:8(11.4%)
1:62(88.6%)
 1 (1.4%)
21 lactating_diet_grain [numeric]
Min : 0
Mean : 0.9
Max : 1
0:8(11.3%)
1:63(88.7%)
 0 (0.0%)
22 lactating_diet_hay [numeric]
Min : 0
Mean : 0.5
Max : 1
0:33(46.5%)
1:38(53.5%)
 0 (0.0%)
23 lactating_diet_silage [numeric]
Min : 0
Mean : 0.7
Max : 1
0:22(31.0%)
1:49(69.0%)
 0 (0.0%)
24 lactating_diet_mixed_ration [numeric]
Min : 0
Mean : 0.4
Max : 1
0:42(60.9%)
1:27(39.1%)
 2 (2.8%)
25 dry_cows_n [numeric]
Mean (sd) : 83.1 (116.3)
min ≤ med ≤ max:
1 ≤ 56.5 ≤ 900
IQR (CV) : 58 (1.4)
 41 distinct values 1 (1.4%)
26 youngstock_n [numeric]
Mean (sd) : 455.2 (1128)
min ≤ med ≤ max:
8 ≤ 200 ≤ 9000
IQR (CV) : 219 (2.5)
 45 distinct values 3 (4.2%)
27 seen_bdd [character]
1. 0
2. 1
38(53.5%)
33(46.5%)
 0 (0.0%)
28 seen_year [character]
1. not sure
2. 2016
3. 2017
4. 2015
5. 2018
6. 2019
7. 2020
8. 2021
9. 0
10. 2000
[ 3 others ]
7(21.9%)
5(15.6%)
5(15.6%)
2(6.2%)
2(6.2%)
2(6.2%)
2(6.2%)
2(6.2%)
1(3.1%)
1(3.1%)
3(9.4%)
 39 (54.9%)
29 seen_one_case [numeric]
Min : 0
Mean : 0
Max : 1
0:32(97.0%)
1:1(3.0%)
 38 (53.5%)
30 seen_odd_case [numeric]
Min : 0
Mean : 0.7
Max : 1
0:9(27.3%)
1:24(72.7%)
 38 (53.5%)
31 seen_increasing [numeric]
Min : 0
Mean : 0.2
Max : 1
0:25(75.8%)
1:8(24.2%)
 38 (53.5%)
32 seen_seasonal_any_time [numeric]
Min : 0
Mean : 0.4
Max : 1
0:20(62.5%)
1:12(37.5%)
 39 (54.9%)
33 seen_season_wet [numeric]
Min : 0
Mean : 0.6
Max : 1
0:12(37.5%)
1:20(62.5%)
 39 (54.9%)
34 seen_seasonal_season [character]
1. hot
2. wet
3. winter
1(5.0%)
14(70.0%)
5(25.0%)
 51 (71.8%)
35 seen_bdd_cause_lameness [character]
1. definitely yes
2. probably not
3. probably yes
4. unsure
19(59.4%)
4(12.5%)
8(25.0%)
1(3.1%)
 39 (54.9%)
36 hoof_care_preventative_trimming [numeric]
Min : 0
Mean : 0.2
Max : 1
0:57(80.3%)
1:14(19.7%)
 0 (0.0%)
37 hoof_care_therapeutic_trimming [numeric]
Min : 0
Mean : 1
Max : 1
0:1(1.4%)
1:70(98.6%)
 0 (0.0%)
38 hoof_care_who_trims_farmer [numeric]
Mean (sd) : 42 (44.7)
min ≤ med ≤ max:
0 ≤ 10 ≤ 100
IQR (CV) : 95 (1.1)
 12 distinct values 0 (0.0%)
39 hoof_care_who_trims_trimmer [numeric]
Mean (sd) : 8.9 (25)
min ≤ med ≤ max:
0 ≤ 0 ≤ 100
IQR (CV) : 0 (2.8)
0:60(85.7%)
10:1(1.4%)
20:1(1.4%)
30:1(1.4%)
50:1(1.4%)
70:1(1.4%)
80:2(2.9%)
90:1(1.4%)
95:1(1.4%)
100:1(1.4%)
 1 (1.4%)
40 hoof_care_who_trims_vet [numeric]
Mean (sd) : 49.9 (44.7)
min ≤ med ≤ max:
0 ≤ 30 ≤ 100
IQR (CV) : 95 (0.9)
0:14(20.0%)
5:8(11.4%)
10:6(8.6%)
20:4(5.7%)
25:1(1.4%)
30:3(4.3%)
50:2(2.9%)
80:3(4.3%)
90:4(5.7%)
100:25(35.7%)
 1 (1.4%)
41 hoof_care_outsider_vet_ever [numeric]
Min : 0
Mean : 0.9
Max : 1
0:8(11.3%)
1:63(88.7%)
 0 (0.0%)
42 hoof_care_outsider_trimmer_ever [numeric]
Min : 0
Mean : 0.3
Max : 1
0:47(66.2%)
1:24(33.8%)
 0 (0.0%)
43 hoof_care_foot_bath_intermittent [numeric]
Min : 0
Mean : 0.1
Max : 1
0:65(91.5%)
1:6(8.5%)
 0 (0.0%)
44 hoof_care_foot_bath_regular [numeric]
Min : 0
Mean : 0.1
Max : 1
0:63(88.7%)
1:8(11.3%)
 0 (0.0%)
45 hoof_care_foot_mats_intermittent [numeric]
Min : 0
Mean : 0.2
Max : 1
0:58(81.7%)
1:13(18.3%)
 0 (0.0%)
46 hoof_care_foot_mats_regular [numeric]
Min : 0
Mean : 0
Max : 1
0:68(95.8%)
1:3(4.2%)
 0 (0.0%)
47 hoof_care_foot_bath_times_per_week [numeric]
Mean (sd) : 4 (2.5)
min ≤ med ≤ max:
1 ≤ 3 ≤ 7
IQR (CV) : 5 (0.6)
1:2(14.3%)
2:4(28.6%)
3:2(14.3%)
5:1(7.1%)
7:5(35.7%)
 57 (80.3%)
48 hoof_care_foot_bath_chemical [character]
1. copper
2. formalin
3. zinc
5(35.7%)
7(50.0%)
2(14.3%)
 57 (80.3%)
49 hoof_care_foot_bath_concentration [numeric]
Mean (sd) : 0 (0)
min ≤ med ≤ max:
0 ≤ 0 ≤ 0.1
IQR (CV) : 0 (0.6)
 5 distinct values 64 (90.1%)
50 hoof_care_foot_mats_times_per_week [numeric]
Mean (sd) : 5.5 (2.7)
min ≤ med ≤ max:
1 ≤ 7 ≤ 8
IQR (CV) : 2.5 (0.5)
1:2(18.2%)
2:1(9.1%)
7:7(63.6%)
8:1(9.1%)
 60 (84.5%)
51 hoof_care_foot_mats_chemical [character]
1. copper
2. formalin
3. zinc
7(46.7%)
1(6.7%)
7(46.7%)
 56 (78.9%)
52 biosecurity_most_recent_introduction_lactating_cows [numeric]
Mean (sd) : 18.6 (33.1)
min ≤ med ≤ max:
1 ≤ 5 ≤ 99
IQR (CV) : 9 (1.8)
1:22(31.0%)
5:26(36.6%)
10:8(11.3%)
20:5(7.0%)
99:10(14.1%)
 0 (0.0%)
53 biosecurity_oldest_introduction_lactating_cows [numeric]
Mean (sd) : 22 (32)
min ≤ med ≤ max:
1 ≤ 10 ≤ 99
IQR (CV) : 15 (1.5)
1:10(14.1%)
5:21(29.6%)
10:14(19.7%)
20:16(22.5%)
99:10(14.1%)
 0 (0.0%)
54 biosecurity_most_recent_introduction_bulls [numeric]
Mean (sd) : 20.9 (37.4)
min ≤ med ≤ max:
1 ≤ 5 ≤ 99
IQR (CV) : 4 (1.8)
1:29(40.8%)
5:27(38.0%)
10:1(1.4%)
20:1(1.4%)
99:13(18.3%)
 0 (0.0%)
55 biosecurity_oldest_introduction_bulls [numeric]
Mean (sd) : 22.9 (36.5)
min ≤ med ≤ max:
0 ≤ 5 ≤ 99
IQR (CV) : 10 (1.6)
0:1(1.4%)
1:12(16.9%)
5:34(47.9%)
10:6(8.5%)
20:5(7.0%)
99:13(18.3%)
 0 (0.0%)
56 housing_lactating_pasture [numeric]
Min : 0
Mean : 0.9
Max : 1
0:9(12.7%)
1:62(87.3%)
 0 (0.0%)
57 housing_lactating_indoor [numeric]
Min : 0
Mean : 0.1
Max : 1
0:62(87.3%)
1:9(12.7%)
 0 (0.0%)
58 housing_lactating_feed_pad [numeric]
Min : 0
Mean : 0.7
Max : 1
0:19(26.8%)
1:52(73.2%)
 0 (0.0%)
59 housing_feed_pad_scrape_per_4_weeks [numeric]
Mean (sd) : 12.1 (19.2)
min ≤ med ≤ max:
0 ≤ 3 ≤ 84
IQR (CV) : 12.5 (1.6)
 13 distinct values 20 (28.2%)
60 closed_herd_cow [numeric]
Min : 0
Mean : 0.1
Max : 1
0:61(85.9%)
1:10(14.1%)
 0 (0.0%)
61 closed_herd_bulls [numeric]
Min : 0
Mean : 0.2
Max : 1
0:58(81.7%)
1:13(18.3%)
 0 (0.0%)
62 closed_herd [numeric]
Min : 0
Mean : 0
Max : 1
0:70(98.6%)
1:1(1.4%)
 0 (0.0%)
63 who_does_most_trimming [character]
1. farmer
2. trimmer
3. vet
32(45.1%)
7(9.9%)
32(45.1%)
 0 (0.0%)
64 farmer_trims [numeric]
Min : 0
Mean : 0.5
Max : 1
0:39(54.9%)
1:32(45.1%)
 0 (0.0%)
65 vet_trims [numeric]
Min : 0
Mean : 0.5
Max : 1
0:39(54.9%)
1:32(45.1%)
 0 (0.0%)
66 trimmer_trims [numeric]
Min : 0
Mean : 0.1
Max : 1
0:64(90.1%)
1:7(9.9%)
 0 (0.0%)
67 mgmt [factor]
1. grazing
2. barn
3. feed pad
19(26.8%)
9(12.7%)
43(60.6%)
 0 (0.0%)
68 mgmt_barn [numeric]
Min : 0
Mean : 0.1
Max : 1
0:62(87.3%)
1:9(12.7%)
 0 (0.0%)
69 mgmt_feed_pad [numeric]
Min : 0
Mean : 0.6
Max : 1
0:28(39.4%)
1:43(60.6%)
 0 (0.0%)
70 mgmt_grazing [numeric]
Min : 0
Mean : 0.3
Max : 1
0:52(73.2%)
1:19(26.8%)
 0 (0.0%)
71 bdd_lame [numeric]
Min : 0
Mean : 0.4
Max : 1
0:44(62.0%)
1:27(38.0%)
 0 (0.0%)

Generated by summarytools 1.0.1 (R version 4.2.1)
2023-10-29

Description of herds

Show the code 
summary <- data %>%
  mutate(across(c(region,breed,calving_pattern, hoof_care_preventative_trimming, hoof_care_foot_bath_regular, hoof_care_foot_bath_intermittent, hoof_care_foot_mats_regular, hoof_care_foot_mats_intermittent, farmer_trims, trimmer_trims, vet_trims, mgmt_grazing, mgmt_barn, mgmt_feed_pad, closed_herd,seen_bdd,seen_seasonal_season,seen_bdd_cause_lameness), as.factor)
  )

table1_output_1 <- table1::table1(~ region + herd_size + daily_milk + calving_pattern +
                                    mgmt_grazing + mgmt_barn + mgmt_feed_pad + closed_herd, data=summary) %>% as_tibble

table1_output_2 <- table1::table1(~ hoof_care_preventative_trimming+ hoof_care_foot_bath_regular+ hoof_care_foot_bath_intermittent+ hoof_care_foot_mats_regular+ hoof_care_foot_mats_intermittent+ farmer_trims+ trimmer_trims+ vet_trims+ seen_bdd+seen_seasonal_season+seen_bdd_cause_lameness , data=summary) %>% as_tibble


library(openxlsx)
wb <- createWorkbook()
# add worksheet
addWorksheet(wb, "Sheet1")
addWorksheet(wb, "Sheet2")
# write data to worksheet
writeData(wb, "Sheet1", table1_output_1)
writeData(wb, "Sheet2", table1_output_2)

# save workbook to file
#saveWorkbook(wb, "table1_output.xlsx", overwrite = TRUE)

Number of herds and cows examined

Show the code 
data$lesion <- (data$cows_scored * data$prev_cow) %>% round(0)

data %>% summarise(
    herds_visited = n(),
    cows_scored = sum(cows_scored),
    total_cows_in_herds_visited = sum(herd_size),
    cows_lesion = sum(lesion),
    apparent_prevalence_overall = cows_lesion / cows_scored

)
herds_visited cows_scored total_cows_in_herds_visited cows_lesion apparent_prevalence_overall
71 15813 31284 1817 0.1149055
Show the code 
data %>% group_by(region) %>% summarise(
    herds_visited = n(),
    cows_scored = sum(cows_scored),
    total_cows_in_herds_visited = sum(herd_size)
)
region herds_visited cows_scored total_cows_in_herds_visited
NSW: Bega Valley 5 987 2721
NSW: Central West 6 995 10091
NSW: Greater Western Sydney 5 971 2696
NSW: North Coast 5 582 1958
NSW: South Coast 12 2092 2630
QLD: Atherton Tablelands 18 3344 4236
Victoria 20 6842 6952

Apparent within-herd prevalence of lesions

All lesions (M1,M2,M3,M4,M4.1)

Show the code 
hist(data$prev_cow,breaks="FD")

Show the code 
data %>% tabyl(prev_cat)
prev_cat n percent
(0,0.1] 35 0.4929577
(-Inf,0] 3 0.0422535
(0.1,0.2] 17 0.2394366
(0.2,0.3] 9 0.1267606
(0.3,1] 7 0.0985915

Proportion of herds with at least 1 lesion observed

Show the code 
rbind(data %>% summarise(
  region = "All herds",
  herds_with_lesions_observed = sum(prev_cow != 0),
  herds_with_m2 = sum(prev_m2 != 0),
  herds_visited = n(),
  proportion_endemic_herds = herds_with_lesions_observed / herds_visited
),

data %>% group_by(region) %>% summarise(
  herds_with_lesions_observed = sum(prev_cow != 0),
  herds_with_m2 = sum(prev_m2 != 0),
  herds_visited = n(),
  proportion_endemic_herds = herds_with_lesions_observed / herds_visited
)
)
region herds_with_lesions_observed herds_with_m2 herds_visited proportion_endemic_herds
All herds 68 7 71 0.9577465
NSW: Bega Valley 5 1 5 1.0000000
NSW: Central West 6 3 6 1.0000000
NSW: Greater Western Sydney 5 1 5 1.0000000
NSW: North Coast 5 0 5 1.0000000
NSW: South Coast 12 0 12 1.0000000
QLD: Atherton Tablelands 17 1 18 0.9444444
Victoria 18 1 20 0.9000000
Show the code 
prev_all_scores <- rbind(data %>% summarise(
  region = "All herds",
  herd_prevalence = round(median(prev_cow*100),1),
  range = paste0(round(min(prev_cow*100),1)," to ",round(max(prev_cow*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence = round(median(prev_cow*100),1),
  range = paste0(round(min(prev_cow*100),1)," to ",round(max(prev_cow*100),1))
)
)
prev_all_scores
region herd_prevalence range
All herds 8.5 0 to 53
NSW: Bega Valley 8.7 2.3 to 27.6
NSW: Central West 25.1 10.3 to 53
NSW: Greater Western Sydney 33.0 11.5 to 35.6
NSW: North Coast 12.7 4.8 to 32.8
NSW: South Coast 14.6 2 to 35.7
QLD: Atherton Tablelands 3.5 0 to 17.4
Victoria 5.2 0 to 26.9

M1 lesions

Show the code 
prev_m1 <- rbind(
data %>% summarise(
  region = "All herds",
  herd_prevalence_m1 = round(median(prev_m1*100),1),
  range_m1 = paste0(round(min(prev_m1*100),1)," to ",round(max(prev_m1*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence_m1 = round(median(prev_m1*100),1),
  range_m1 = paste0(round(min(prev_m1*100),1)," to ",round(max(prev_m1*100),1))
)
)

prev_m1
region herd_prevalence_m1 range_m1
All herds 0 0 to 0.7
NSW: Bega Valley 0 0 to 0.4
NSW: Central West 0 0 to 0
NSW: Greater Western Sydney 0 0 to 0
NSW: North Coast 0 0 to 0
NSW: South Coast 0 0 to 0
QLD: Atherton Tablelands 0 0 to 0.7
Victoria 0 0 to 0.3

M2 lesions

Show the code 
prev_m2 <- rbind(
data %>% summarise(
  region = "All herds",
  herd_prevalence_m2 = round(median(prev_m2*100),1),
  range_m2 = paste0(round(min(prev_m2*100),1)," to ",round(max(prev_m2*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence_m2 = round(median(prev_m2*100),1),
  range_m2 = paste0(round(min(prev_m2*100),1)," to ",round(max(prev_m2*100),1))
)
)

prev_m2
region herd_prevalence_m2 range_m2
All herds 0.0 0 to 2.2
NSW: Bega Valley 0.0 0 to 1.2
NSW: Central West 0.3 0 to 2.2
NSW: Greater Western Sydney 0.0 0 to 0.2
NSW: North Coast 0.0 0 to 0
NSW: South Coast 0.0 0 to 0
QLD: Atherton Tablelands 0.0 0 to 1.3
Victoria 0.0 0 to 0.8

M2’s were very uncommon on endemic herds.

Show the code 
data <- data %>% mutate(
  m2_proportion_of_lesions = case_when(
    prev_cow == 0 ~ NA,
    prev_cow != 0 ~ prev_m2 / prev_cow
))

data$m2_proportion_of_lesions %>% summary
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
0.000000 0.000000 0.000000 0.005772 0.000000 0.120000        3

M3 lesions

Show the code 
prev_m3 <- rbind(
data %>% summarise(
  region = "All herds",
  herd_prevalence_m3 = round(median(prev_m3*100),1),
  range_m3 = paste0(round(min(prev_m3*100),1)," to ",round(max(prev_m3*100),1)),
),

data %>% group_by(region) %>% summarise(
  herd_prevalence_m3 = round(median(prev_m3*100),1),
  range_m3 = paste0(round(min(prev_m3*100),1)," to ",round(max(prev_m3*100),1)),
)
)

prev_m3
region herd_prevalence_m3 range_m3
All herds 0.0 0 to 1.6
NSW: Bega Valley 0.0 0 to 0.8
NSW: Central West 0.2 0 to 1.6
NSW: Greater Western Sydney 0.0 0 to 0
NSW: North Coast 0.0 0 to 0
NSW: South Coast 0.0 0 to 0
QLD: Atherton Tablelands 0.0 0 to 0
Victoria 0.0 0 to 0.4

M4 lesions

Show the code 
prev_m4 <- rbind(
data %>% summarise(
  region = "All herds",
  herd_prevalence_m4 = round(median(prev_m4*100),1),
  range_m4 = paste0(round(min(prev_m4*100),1)," to ",round(max(prev_m4*100),1))
),
  
data %>% group_by(region) %>% summarise(
  herd_prevalence_m4 = round(median(prev_m4*100),1),
  range_m4 = paste0(round(min(prev_m4*100),1)," to ",round(max(prev_m4*100),1))
))

prev_m4
region herd_prevalence_m4 range_m4
All herds 8.5 0 to 45.7
NSW: Bega Valley 8.5 2.3 to 26
NSW: Central West 22.5 10.3 to 45.7
NSW: Greater Western Sydney 32.1 10.8 to 35.6
NSW: North Coast 11.9 4.8 to 31
NSW: South Coast 14.6 2 to 35.7
QLD: Atherton Tablelands 3.5 0 to 16.8
Victoria 4.9 0 to 25.7

M4’s are the predominant lesion on endemic farms

Show the code 
data <- data %>% mutate(
  m4_proportion_of_lesions = case_when(
    prev_cow == 0 ~ NA,
    prev_cow != 0 ~ prev_m4 / prev_cow
))

data$m4_proportion_of_lesions %>% summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.7000  0.9514  1.0000  0.9627  1.0000  1.0000       3

M4.1 lesions

Show the code 
prev_m4.1  <- rbind(
data %>% summarise(
  region = "All herds",
  herd_prevalence_m4.1 = round(median(prev_m4.1*100),1),
  range_m4.1 = paste0(round(min(prev_m4.1*100),1)," to ",round(max(prev_m4.1*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence_m4.1 = round(median(prev_m4.1*100),1),
  range_m4.1 = paste0(round(min(prev_m4.1*100),1)," to ",round(max(prev_m4.1*100),1))
)
)

prev_m4.1
region herd_prevalence_m4.1 range_m4.1
All herds 0.0 0 to 8.6
NSW: Bega Valley 0.8 0 to 1.6
NSW: Central West 1.5 0 to 6
NSW: Greater Western Sydney 0.7 0 to 3.8
NSW: North Coast 0.0 0 to 8.6
NSW: South Coast 0.0 0 to 3.1
QLD: Atherton Tablelands 0.0 0 to 0.4
Victoria 0.0 0 to 1.2
Show the code 
table_prevalence <- cbind(
  prev_all_scores,
  prev_m1,
  prev_m2,
  prev_m3,
  prev_m4,
  prev_m4.1
)


write.csv(table_prevalence,'table_prev.csv')

Estimating true prevalence

Method 1: Probabilistic bias analysis (Monte Carlo simulation)

Assumptions based on:

Yang, D. A., & Laven, R. A. (2019). Detecting bovine digital dermatitis in the milking parlour: to wash or not to wash, a Bayesian superpopulation approach. The Veterinary Journal, 247, 38-43.

Se = 0.63 (95%CrI: 0.46–0.78)

Sp = 0.99 (95%CrI: 0.97-1.0)

True prevalence = (Apparent prevalence + Sp - 1)/(Se + Sp - 1)

Show the code 
#Initial attempt using fixed values. Will not use this. 
# calculate_true_prev <- function(){
# sens <- 0.63
# spec <- 0.99
# 
# true_prev <- function(prev){
#   (prev + spec - 1)/(sens + spec - 1)
# }
# 
# true_prev(data$prev_cow)
# }


#Using a distribution of values. 
set.seed(1)

calculate_true_prev <- function(){
sens <- rbeta(1, 17.996, 10.152, ncp = 0)
spec <- rbeta(1, 100, 2, ncp = 0)

true_prev <- function(prev){
  (prev + spec - 1)/(sens + spec - 1)
}

output <- true_prev(data$prev_cow)
ifelse(output < 0,0,output)
}

mc <- replicate(
  10000,
  calculate_true_prev(),
  simplify = T
) %>% t %>% as.data.frame()

mc <- mc %>% pivot_longer(
  cols = everything(),
  names_to = "herd",
  values_to = "true_prev"
)

mc$herd <- parse_number(mc$herd)

mc <- mc %>% group_by(herd) %>% summarise(
  median_mc = median(true_prev),
  mean_mc = mean(true_prev),
  q5_mc = quantile(true_prev,probs = 0.025),
  q95_mc = quantile(true_prev,probs = 0.975)
) %>% arrange(herd)

data$true_prev_mc_median <- mc$median_mc
data$true_prev_mc <- paste0(round(mc$median_mc*100,1),
                         " 95% CrI: ",
                         round(mc$q5_mc*100,1),
                         " to ",
                         round(mc$q95_mc*100,1))

Method 2: Bayesian estimation using priors

Settings:

  • 2 markov chains

  • Burn-in: 10,000

  • Update: 10,000

  • Priors - Yang, D. A., & Laven, R. A. (2019). Detecting bovine digital dermatitis in the milking parlour: to wash or not to wash, a Bayesian superpopulation approach. The Veterinary Journal, 247, 38-43.

    • Se = 0.63 (95%CrI: 0.46-0.78) = beta(17.996, 10.152)

    • Sp = 0.99 (95%CrI: 0.97-1.0) = dbeta(100, 2)

    • Prevalence = uniform(0,1)

      Show priors for Se and Sp

Show the code 
N = 10000
prior <- data.frame(Value = c(rbeta(N,17.996, 10.152),rbeta(N,100, 2)),
           Distribution = rep(c("Sensitivity","Specificity"),each=N))

ggplot(data = prior, aes(x = Value, group = Distribution)) + 
  geom_density(col = NA, alpha = 0.5,aes(fill = Distribution, color = Distribution)) + theme_classic()+ theme(axis.title.y=element_blank(),
                                            axis.text.y=element_blank(),
                                   legend.position = "bottom",
                                   text=element_text(family="Times")) + ggtitle("Priors") + xlim(0,1)

Show the code 
# tiff("plot1.tiff", units="in", width=4, height=4, res=200)
# ggplot(data = prior, aes(x = Value, group = Distribution)) +
#   geom_density(col = NA, alpha = 0.5,aes(fill = Distribution, color = Distribution)) +
#   theme_classic()+ theme(axis.title.y=element_blank(), axis.text.y=element_blank(),
#                                    legend.position =  c(0.2, 0.6),legend.title = element_blank(),
#                                    text=element_text(family="Times")) + ggtitle("") + xlim(0,1)
# dev.off()

Estimated true within-herd, cow-level BDD prevalence based on probabilistic bias analysis

Show the code 
data$positive_n <- as.integer(round(data$prev_cow * data$cows_scored,0))

# truePrev(x = as.integer(round(data$prev_cow * data$cows_scored,0)), n = data$cows_scored,
#          SE = ~dbeta(17.996, 10.152), SP = ~dbeta(100, 2))


check <- truePrev(x = 1, n = 1000,
         SE = ~dbeta(17.996, 10.152), SP = ~dbeta(100, 2))


true_prev_all <- function(positive,scored){
mcmc <- truePrev(x = positive, n = scored,
         SE = ~dbeta(17.996, 10.152), SP = ~dbeta(100, 2),
         nchains = 2, burnin = 10000, update = 10000)
chains <- c(mcmc@mcmc[["TP"]][[1]],mcmc@mcmc[["TP"]][[2]])
median_mcmc <- median(chains)
lower_mcmc <- quantile(chains,probs= 0.025)
upper_mcmc <- quantile(chains,probs= 0.975)
tibble(median_mcmc,lower_mcmc,upper_mcmc)
}

true_prev <- map2(.x=data$cows_scored,
                  .y=data$positive_n,
                  ~true_prev_all(
                    scored = .x,
                    positive = .y
                  ))



data <- cbind(data,do.call(rbind.data.frame, true_prev))

data$true_prev_mcmc <- paste0(round(data$median_mcmc*100,1),
                         " 95% CrI: ",
                         round(data$lower_mcmc*100,1),
                         " to ",
                         round(data$upper_mcmc*100,1))


data$apparent_prev <- data$prev_cow
data$true_prev_probabilistic_bias_analysis <- data$true_prev_mc_median
data$true_prev_bayesian <- data$median_mcmc

prev <- data %>% select(region,apparent_prev,true_prev_probabilistic_bias_analysis,true_prev_bayesian) %>% pivot_longer(
  cols = c(apparent_prev,
           true_prev_probabilistic_bias_analysis,
           true_prev_bayesian),
  names_to = "Measure of Prevalence",
  values_to = "Prevalence"
)


ggplot(data = prev, aes(x = Prevalence, group = `Measure of Prevalence`)) + geom_density(col = NA, alpha = 0.5,aes(fill = `Measure of Prevalence`, color = `Measure of Prevalence`)) + theme_classic()+ theme(axis.title.y=element_blank(),
                                            axis.text.y=element_blank(),
                                   legend.position = "bottom",
                                   text=element_text(family="Times")) + ggtitle("Distribution of herd prevalence of Bovine Digital Dermatitis for 71 Australian Dairy herds") + xlim(0,1)

Show the code 
prev1 <- prev %>% filter(`Measure of Prevalence` != "true_prev_probabilistic_bias_analysis")
prev1 <- prev1 %>% mutate(
  `Measure of Prevalence` = case_when(
    `Measure of Prevalence` == "apparent_prev" ~ "Apparent Prevalence",
    `Measure of Prevalence` == "true_prev_bayesian" ~ "Estimated True Prevalence (Bayesian)",
  )
)

tiff("plot2.tiff", units="in", width=4, height=4, res=200)
ggplot(data = prev1, aes(x = Prevalence, group = `Measure of Prevalence`)) + geom_density(col = NA, alpha = 0.5,aes(fill = `Measure of Prevalence`, color = `Measure of Prevalence`)) + theme_classic()+ theme(axis.title.y=element_blank(),
                                            axis.text.y=element_blank(),
                                   legend.position =  c(0.6, 0.6),legend.title = element_blank(),
                                   text=element_text(family="Times")) + ggtitle("") + xlim(0,1)
dev.off()
quartz_off_screen 
                2
Show the code 
rbind(data %>% summarise(
  region = "All herds",
  herd_prevalence = round(mean(true_prev_probabilistic_bias_analysis*100),1),
  range = paste0(round(min(true_prev_probabilistic_bias_analysis*100),1)," to ",round(max(true_prev_probabilistic_bias_analysis*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence = round(mean(true_prev_probabilistic_bias_analysis*100),1),
  range = paste0(round(min(true_prev_probabilistic_bias_analysis*100),1)," to ",round(max(true_prev_probabilistic_bias_analysis*100),1))
)
)
region herd_prevalence range
All herds 17.0 0 to 81.9
NSW: Bega Valley 16.8 1.1 to 41.1
NSW: Central West 41.2 13.6 to 81.9
NSW: Greater Western Sydney 37.5 15.4 to 53.9
NSW: North Coast 21.7 4.9 to 49.4
NSW: South Coast 21.3 0.5 to 54.2
QLD: Atherton Tablelands 6.9 0 to 24.7
Victoria 10.0 0 to 40

Estimated true within-herd, cow-level BDD prevalence based on Bayesian estimation

Show the code 
prev_bayes <- rbind(data %>% summarise(
  region = "All herds",
  herd_prevalence = round(mean(true_prev_bayesian*100),1),
  range = paste0(round(min(true_prev_bayesian*100),1)," to ",round(max(true_prev_bayesian*100),1))
),

data %>% group_by(region) %>% summarise(
  herd_prevalence = round(mean(true_prev_bayesian*100),1),
  range = paste0(round(min(true_prev_bayesian*100),1)," to ",round(max(true_prev_bayesian*100),1))
)
)

prev_bayes
region herd_prevalence range
All herds 18.1 0.4 to 80.8
NSW: Bega Valley 17.8 2.8 to 42.3
NSW: Central West 42.1 14.7 to 80.8
NSW: Greater Western Sydney 38.7 16.2 to 55.3
NSW: North Coast 23.0 5.4 to 51.1
NSW: South Coast 22.8 2.8 to 55.9
QLD: Atherton Tablelands 8.0 0.5 to 25.9
Victoria 10.8 0.4 to 41.1
Show the code 
write.csv(merge(prev_all_scores,prev_bayes,
      by="region"),"prev_table.csv")

Investigation of possible risk factors for BDD prevalence

Research questions:

  • Does preventative trimming increase the spread of BDD?

  • Do hoof trimmers increase the spread BDD?

  • Do barns or feed pads increase the spread of BDD?

  • Can farmers identify it in their herds?

  • Does BDD cause lameness?

Show the code 
save <- data
data$prev_cow <- data$true_prev_bayesian

Univariable assocations for prevalence odds ratio

Method 1 is to create a cow-level dataset and then use mixed effects logistic regression

Show the code 
odds <- data %>% mutate(
  farm_number = seq(1:71)
)

# 
# farm <- odds %>% filter(farm_number==1)
# pos_cows <- round(farm$true_prev_bayesian * farm$cows_scored,0)
# neg_cows <- farm$cows_scored - pos_cows
# 
# farm <- farm %>% slice(rep(1:n(), each = cows_scored))
# farm <- farm %>% mutate(
#   lesion_status = c(rep(x=1,times = pos_cows),rep(x=0,times=neg_cows))
# )


# Create an empty list to store dataframes
result_list <- list()

# Loop through farm_number from 1 to 71
for (i in 1:71) {
  farm <- odds %>% filter(farm_number == i)
  pos_cows <- round(farm$true_prev_bayesian * farm$cows_scored, 0)
  neg_cows <- farm$cows_scored - pos_cows
  farm <- farm %>% slice(rep(1:n(), each = farm$cows_scored))

  # pos_cows <- round(100*farm$true_prev_bayesian, 0)
  # neg_cows <- 100 - pos_cows
  # farm <- farm %>% slice(rep(1:n(), each = 100))

  farm <- farm %>% mutate(
    lesion_status = c(rep(x = 1, times = pos_cows), rep(x = 0, times = neg_cows))
  )
  
  # Append the dataframe to the list
  result_list[[i]] <- farm
}

# Combine all dataframes in the list into one large dataframe
odds_df <- bind_rows(result_list)

# odds_df %>% tabyl(lesion_status)

library(lme4)
library(broom.mixed)

# prev_odds ratio function

glmm <- glmer(lesion_status ~ trimmer_trims + (1 | farm_number),
              data = odds_df,family = binomial(link="logit"))

# glmm <- glmer(lesion_status ~ trimmer_trims + (1 | farm_number),
#               data = odds_df,family = binomial(link="log"))
# 
# glmm <- glmer(lesion_status ~ trimmer_trims + (1 | farm_number),
#               data = odds_df,family = binomial(link="identity"))
summary(glmm)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: lesion_status ~ trimmer_trims + (1 | farm_number)
   Data: odds_df

     AIC      BIC   logLik deviance df.resid 
 11641.2  11664.2  -5817.6  11635.2    15810 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9616 -0.4319 -0.2541 -0.1474 10.1491 

Random effects:
 Groups      Name        Variance Std.Dev.
 farm_number (Intercept) 1.658    1.288   
Number of obs: 15813, groups:  farm_number, 71

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -2.2120     0.1663 -13.302  < 2e-16 ***
trimmer_trims   1.7875     0.5199   3.438 0.000586 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
trimmr_trms -0.321
Show the code 
extract_por <- function(glmm,term_ref){
tidy(glmm,conf.int=TRUE,exponentiate=T,effects="fixed") %>% 
  slice(2:term_ref) %>%
  mutate(
    prev_odds_ratio = round(estimate,1),
    conf = paste0(round(conf.low,1)," to ",round(conf.high,1))
    ) %>% 
  select(term,prev_odds_ratio,conf)
}



# data$odds <- data$true_prev_bayesian / (1 - data$true_prev_bayesian)
# 
# lm1 <- glm(true_prev_bayesian ~ trimmer_trims,data=data,
#            family = gaussian(link="log"))
# 
# lm1 <- glm(odds ~ trimmer_trims,data=data,
#            family = gaussian(link="log"))
# 
# data %>% group_by(trimmer_trims) %>% summarise(
#   odds = mean(odds),
#   prev = mean(true_prev_bayesian)
# )
#
Show the code 
odds_df <- odds_df %>% mutate(
  foot_bath = case_when(
    hoof_care_foot_bath_intermittent == 1 ~ "intermittent",
    hoof_care_foot_bath_regular == 1 ~ "regular",
    T ~ "never"
  ) %>% as.factor %>% relevel("never"),
  foot_mat = case_when(
    hoof_care_foot_mats_intermittent == 1 ~ "intermittent",
    hoof_care_foot_mats_regular == 1 ~ "regular",
    T ~ "never"
  ) %>% as.factor %>% relevel("never"),
  
  outsider = case_when(
    hoof_care_outsider_trimmer_ever == 1 & hoof_care_outsider_vet_ever==1 ~ "Trimmer and Vet",
    hoof_care_outsider_trimmer_ever == 1 & hoof_care_outsider_vet_ever==0 ~ "Trimmer only",
    hoof_care_outsider_trimmer_ever == 0 & hoof_care_outsider_vet_ever==1 ~ "Vet only",
    hoof_care_outsider_trimmer_ever == 0 & hoof_care_outsider_vet_ever==1 ~ "None"
  ) %>% as.factor %>% relevel("Vet only"),
  
  awareness = case_when(
    seen_bdd_cause_lameness == "definitely yes" ~ "seen_lame",
    seen_bdd == 1 ~ "seen_no_lame",
    T ~ "not"
  ),
  
  open_herd_cow = case_when(
    closed_herd_cow == 1 ~ 0,
    closed_herd_cow == 0 ~ 1
  ),
  
  open_herd_bull = case_when(
    closed_herd_bulls == 1 ~ 0,
    closed_herd_bulls == 0 ~ 1
  )
)

data <- data %>% mutate(
  foot_bath = case_when(
    hoof_care_foot_bath_intermittent == 1 ~ "intermittent",
    hoof_care_foot_bath_regular == 1 ~ "regular",
    T ~ "never"
  ) %>% as.factor %>% relevel("never"),
  foot_mat = case_when(
    hoof_care_foot_mats_intermittent == 1 ~ "intermittent",
    hoof_care_foot_mats_regular == 1 ~ "regular",
    T ~ "never"
  ) %>% as.factor %>% relevel("never"),
  
  outsider = case_when(
    hoof_care_outsider_trimmer_ever == 1 & hoof_care_outsider_vet_ever==1 ~ "Trimmer and Vet",
    hoof_care_outsider_trimmer_ever == 1 & hoof_care_outsider_vet_ever==0 ~ "Trimmer only",
    hoof_care_outsider_trimmer_ever == 0 & hoof_care_outsider_vet_ever==1 ~ "Vet only",
    hoof_care_outsider_trimmer_ever == 0 & hoof_care_outsider_vet_ever==1 ~ "None"
  ) %>% as.factor %>% relevel("Vet only"),
  
  awareness = case_when(
    seen_bdd_cause_lameness == "definitely yes" ~ "seen_lame",
    seen_bdd == 1 ~ "seen_no_lame",
    T ~ "not"
  ),
  
  open_herd_cow = case_when(
    closed_herd_cow == 1 ~ 0,
    closed_herd_cow == 0 ~ 1
  ),
  
  open_herd_bull = case_when(
    closed_herd_bulls == 1 ~ 0,
    closed_herd_bulls == 0 ~ 1
  )
)

univar1 <- rbind(
  glmer(lesion_status ~ mgmt + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 3),
  glmer(lesion_status ~ hoof_care_preventative_trimming + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 2),
  glmer(lesion_status ~ foot_mat + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 3),
  glmer(lesion_status ~ foot_bath + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 3),
  glmer(lesion_status ~ outsider + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 4),
  glmer(lesion_status ~ who_does_most_trimming + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 3),
  glmer(lesion_status ~ awareness + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 3),
  glmer(lesion_status ~ open_herd_cow + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 2),
  glmer(lesion_status ~ open_herd_bull + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 2),
  #glmer(lesion_status ~ herd_size + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 2),
  glmer(lesion_status ~ daily_milk + (1 | farm_number),data = odds_df,family = binomial(link="logit")) %>% extract_por(term_ref = 2)
  
)
#Notes - it looks as though using the mixed model approach gives different effect estimates to the population averaged approach (glm), which gives very similar estimates of OR to what we get when collapsing each herd into a single row. It doesn't seem right to me to use glm, as the confidence intervals will be wrong, and I think that the effect estimates might be more appropriate as we want subject (herd) specific estimates.

univar1
term prev_odds_ratio conf
mgmtbarn 10.5 4 to 27.9
mgmtfeed pad 3.2 1.6 to 6.3
hoof_care_preventative_trimming 2.8 1.3 to 6.2
foot_matintermittent 1.7 0.7 to 3.9
foot_matregular 2.6 0.5 to 13.1
foot_bathintermittent 1.4 0.4 to 4.7
foot_bathregular 1.3 0.4 to 3.7
outsiderTrimmer and Vet 3.1 1.5 to 6.3
outsiderTrimmer only 6.3 2.5 to 15.8
who_does_most_trimmingtrimmer 5.2 1.8 to 15
who_does_most_trimmingvet 0.8 0.4 to 1.4
awarenessseen_lame 3.7 1.8 to 7.7
awarenessseen_no_lame 1.6 0.7 to 3.6
open_herd_cow 2.0 0.8 to 5.1
open_herd_bull 0.8 0.3 to 1.8
daily_milk 1.1 1.1 to 1.2
Show the code 
# odds_df %>% tabyl(trimmer_trims,
#                   lesion_status)


# glm1 <- glm(lesion_status ~ trimmer_trims,
#             data = odds_df,family=binomial)
# 
# 
# glm1 <- glmer(lesion_status ~ trimmer_trims + (1|farm_number),
#             data = odds_df,family=binomial)
# 
# library(geepack)
# gee <- geeglm(lesion_status ~ trimmer_trims, id = farm_number,
#               data= odds_df, family = binomial)
# 
# tidy(gee,exp=T)
# 
# tidy(glm1,exp=T)

univar1 %>% write.csv("univar prev odds ratio.csv")

#data %>% tabyl(hoof_care_foot_bath_regular)

# glmm <- glmer(lesion_status ~ who_does_most_trimming + (1 | farm_number),
#               data = odds_df,family = binomial(link="logit")) %>% extract_por
# 
# summary(glmm)
# data %>% tabyl(who_does_most_trimming)

Diagnostics from the MCMC

These were not uploaded to the analysis log as the file was too large to publish to RPUBS.

Show the code 
true_prev_all_vals <- function(positive,scored){
mcmc <- truePrev(x = positive, n = scored,
         SE = ~dbeta(17.996, 10.152), SP = ~dbeta(100, 2),
         nchains = 2, burnin = 10000, update = 10000)
# chains <- c(mcmc@mcmc[["TP"]][[1]],mcmc@mcmc[["TP"]][[2]])
# chains
plot(mcmc, ask=FALSE)
}

library(coda)

# true_prev_all_vals(positive=10,scored=1000)

true_prev_vals <- map2(.x=data$cows_scored,
                  .y=data$positive_n,
                  ~true_prev_all_vals(
                    scored = .x,
                    positive = .y
                  ))

plot(mcmc)