Exploring the phylogenetic diversity and antimicrobial activity of non-aureus staphylococci and mammaliicocci isolated from teat apices of organic dairy cows. Part 1 - Associations with in vitro antimicrobial activity
Felipe Pena Mosca
2025-08-06
1 Counts of NASM on teat apex
# NASM presence by IMI status
table1(~ as.factor(nas_present_sample) + as.factor(chromogenes_present_sample) + as.factor(cohnii_present_sample) + as.factor(devriesei_present_sample) + as.factor(equorum_present_sample) + as.factor(haemolyticus_present_sample) + as.factor(pseudoxylosus_present_sample) + as.factor(sciuri_present_sample) + as.factor(succinus_present_sample) + as.factor(xylosus_present_sample) + as.factor(xylosus_pseudoxylosus_present_sample) + as.factor(other_nasm_present_sample)| case_control,data=df_sau_per_cow_2)| Case (N=35) |
Control (N=36) |
Overall (N=71) |
|
|---|---|---|---|
| as.factor(nas_present_sample) | |||
| 0 | 13 (37.1%) | 12 (33.3%) | 25 (35.2%) |
| 1 | 22 (62.9%) | 24 (66.7%) | 46 (64.8%) |
| as.factor(chromogenes_present_sample) | |||
| 0 | 27 (77.1%) | 27 (75.0%) | 54 (76.1%) |
| 1 | 8 (22.9%) | 9 (25.0%) | 17 (23.9%) |
| as.factor(cohnii_present_sample) | |||
| 0 | 34 (97.1%) | 36 (100%) | 70 (98.6%) |
| 1 | 1 (2.9%) | 0 (0%) | 1 (1.4%) |
| as.factor(devriesei_present_sample) | |||
| 0 | 35 (100%) | 34 (94.4%) | 69 (97.2%) |
| 1 | 0 (0%) | 2 (5.6%) | 2 (2.8%) |
| as.factor(equorum_present_sample) | |||
| 0 | 35 (100%) | 35 (97.2%) | 70 (98.6%) |
| 1 | 0 (0%) | 1 (2.8%) | 1 (1.4%) |
| as.factor(haemolyticus_present_sample) | |||
| 0 | 27 (77.1%) | 27 (75.0%) | 54 (76.1%) |
| 1 | 8 (22.9%) | 9 (25.0%) | 17 (23.9%) |
| as.factor(pseudoxylosus_present_sample) | |||
| 0 | 33 (94.3%) | 34 (94.4%) | 67 (94.4%) |
| 1 | 2 (5.7%) | 2 (5.6%) | 4 (5.6%) |
| as.factor(sciuri_present_sample) | |||
| 0 | 33 (94.3%) | 33 (91.7%) | 66 (93.0%) |
| 1 | 2 (5.7%) | 3 (8.3%) | 5 (7.0%) |
| as.factor(succinus_present_sample) | |||
| 0 | 31 (88.6%) | 32 (88.9%) | 63 (88.7%) |
| 1 | 4 (11.4%) | 4 (11.1%) | 8 (11.3%) |
| as.factor(xylosus_present_sample) | |||
| 0 | 30 (85.7%) | 33 (91.7%) | 63 (88.7%) |
| 1 | 5 (14.3%) | 3 (8.3%) | 8 (11.3%) |
| as.factor(xylosus_pseudoxylosus_present_sample) | |||
| 0 | 29 (82.9%) | 32 (88.9%) | 61 (85.9%) |
| 1 | 6 (17.1%) | 4 (11.1%) | 10 (14.1%) |
| as.factor(other_nasm_present_sample) | |||
| 0 | 34 (97.1%) | 33 (91.7%) | 67 (94.4%) |
| 1 | 1 (2.9%) | 3 (8.3%) | 4 (5.6%) |
# NAS presence by farm
table1(~ as.factor(nas_present_sample) + as.factor(chromogenes_present_sample) + as.factor(cohnii_present_sample) + as.factor(devriesei_present_sample) + as.factor(equorum_present_sample) + as.factor(haemolyticus_present_sample) + as.factor(pseudoxylosus_present_sample) + as.factor(sciuri_present_sample) + as.factor(succinus_present_sample) + as.factor(xylosus_present_sample) + as.factor(xylosus_pseudoxylosus_present_sample) + as.factor(other_nasm_present_sample)| farm,data=df_sau_per_cow_2)| A (N=57) |
B (N=14) |
Overall (N=71) |
|
|---|---|---|---|
| as.factor(nas_present_sample) | |||
| 0 | 20 (35.1%) | 5 (35.7%) | 25 (35.2%) |
| 1 | 37 (64.9%) | 9 (64.3%) | 46 (64.8%) |
| as.factor(chromogenes_present_sample) | |||
| 0 | 41 (71.9%) | 13 (92.9%) | 54 (76.1%) |
| 1 | 16 (28.1%) | 1 (7.1%) | 17 (23.9%) |
| as.factor(cohnii_present_sample) | |||
| 0 | 57 (100%) | 13 (92.9%) | 70 (98.6%) |
| 1 | 0 (0%) | 1 (7.1%) | 1 (1.4%) |
| as.factor(devriesei_present_sample) | |||
| 0 | 55 (96.5%) | 14 (100%) | 69 (97.2%) |
| 1 | 2 (3.5%) | 0 (0%) | 2 (2.8%) |
| as.factor(equorum_present_sample) | |||
| 0 | 57 (100%) | 13 (92.9%) | 70 (98.6%) |
| 1 | 0 (0%) | 1 (7.1%) | 1 (1.4%) |
| as.factor(haemolyticus_present_sample) | |||
| 0 | 41 (71.9%) | 13 (92.9%) | 54 (76.1%) |
| 1 | 16 (28.1%) | 1 (7.1%) | 17 (23.9%) |
| as.factor(pseudoxylosus_present_sample) | |||
| 0 | 56 (98.2%) | 11 (78.6%) | 67 (94.4%) |
| 1 | 1 (1.8%) | 3 (21.4%) | 4 (5.6%) |
| as.factor(sciuri_present_sample) | |||
| 0 | 54 (94.7%) | 12 (85.7%) | 66 (93.0%) |
| 1 | 3 (5.3%) | 2 (14.3%) | 5 (7.0%) |
| as.factor(succinus_present_sample) | |||
| 0 | 49 (86.0%) | 14 (100%) | 63 (88.7%) |
| 1 | 8 (14.0%) | 0 (0%) | 8 (11.3%) |
| as.factor(xylosus_present_sample) | |||
| 0 | 56 (98.2%) | 7 (50.0%) | 63 (88.7%) |
| 1 | 1 (1.8%) | 7 (50.0%) | 8 (11.3%) |
| as.factor(xylosus_pseudoxylosus_present_sample) | |||
| 0 | 55 (96.5%) | 6 (42.9%) | 61 (85.9%) |
| 1 | 2 (3.5%) | 8 (57.1%) | 10 (14.1%) |
| as.factor(other_nasm_present_sample) | |||
| 0 | 55 (96.5%) | 12 (85.7%) | 67 (94.4%) |
| 1 | 2 (3.5%) | 2 (14.3%) | 4 (5.6%) |
# Counts in log scale in cases and controls:
table1(~ log_nas + log_chromogenes + log_cohnii + log_devriesei + log_equorum + log_haemolyticus + log_pseudoxylosus + log_sciuri + log_succinus + log_xylosus| case_control,data=df_sau_per_cow_2)| Case (N=35) |
Control (N=36) |
Overall (N=71) |
|
|---|---|---|---|
| log_nas | |||
| Mean (SD) | 0.944 (0.950) | 1.03 (0.871) | 0.987 (0.905) |
| Median [Min, Max] | 0.954 [0, 2.81] | 1.03 [0, 2.75] | 1.00 [0, 2.81] |
| log_chromogenes | |||
| Mean (SD) | 0.292 (0.567) | 0.350 (0.649) | 0.321 (0.606) |
| Median [Min, Max] | 0 [0, 2.00] | 0 [0, 1.90] | 0 [0, 2.00] |
| log_cohnii | |||
| Mean (SD) | 0.0458 (0.271) | 0 (0) | 0.0226 (0.190) |
| Median [Min, Max] | 0 [0, 1.60] | 0 [0, 0] | 0 [0, 1.60] |
| log_devriesei | |||
| Mean (SD) | 0 (0) | 0.0625 (0.266) | 0.0317 (0.191) |
| Median [Min, Max] | 0 [0, 0] | 0 [0, 1.32] | 0 [0, 1.32] |
| log_equorum | |||
| Mean (SD) | 0 (0) | 0.0661 (0.397) | 0.0335 (0.282) |
| Median [Min, Max] | 0 [0, 0] | 0 [0, 2.38] | 0 [0, 2.38] |
| log_haemolyticus | |||
| Mean (SD) | 0.303 (0.726) | 0.296 (0.609) | 0.300 (0.665) |
| Median [Min, Max] | 0 [0, 2.37] | 0 [0, 2.35] | 0 [0, 2.37] |
| log_pseudoxylosus | |||
| Mean (SD) | 0.104 (0.436) | 0.0737 (0.318) | 0.0886 (0.379) |
| Median [Min, Max] | 0 [0, 2.16] | 0 [0, 1.65] | 0 [0, 2.16] |
| log_sciuri | |||
| Mean (SD) | 0.141 (0.582) | 0.126 (0.438) | 0.133 (0.511) |
| Median [Min, Max] | 0 [0, 2.66] | 0 [0, 1.93] | 0 [0, 2.66] |
| log_succinus | |||
| Mean (SD) | 0.132 (0.450) | 0.169 (0.525) | 0.151 (0.486) |
| Median [Min, Max] | 0 [0, 1.85] | 0 [0, 2.16] | 0 [0, 2.16] |
| log_xylosus | |||
| Mean (SD) | 0.200 (0.540) | 0.157 (0.553) | 0.179 (0.543) |
| Median [Min, Max] | 0 [0, 2.20] | 0 [0, 2.51] | 0 [0, 2.51] |
# --- NAS ----
# Presence
model <- glm(nas_present_sample~case+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 2.0101629 0.9851076 4.333702
## case 0.8462728 0.3161614 2.249185
## farmB 0.9743576 0.2936964 3.536861# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.665 0.0894 Inf 0.475 0.813
## 1 0.627 0.0927 Inf 0.436 0.785
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_nas~case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_nas ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22176 -0.89566 0.01813 0.78528 1.67736
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9819 0.1610 6.098 5.72e-08 ***
## case -0.0862 0.2166 -0.398 0.692
## farmB 0.2399 0.2722 0.881 0.381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9124 on 68 degrees of freedom
## Multiple R-squared: 0.0135, Adjusted R-squared: -0.01552
## F-statistic: 0.4653 on 2 and 68 DF, p-value: 0.63Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.98186685 0.6605812 1.3031525
## case -0.08620335 -0.5183922 0.3459855
## farmB 0.23988875 -0.3031832 0.7829607# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 1.10 0.173 68 0.756 1.45
## 1 1.02 0.174 68 0.667 1.36
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95# --- Staph chromogenes ----
# Presence
model <- glm(chromogenes_present_sample~ case+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
# OR and 95% CI
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.4127831 0.17856064 0.8778753
## case 0.8906506 0.28708487 2.7259299
## farmB 0.1971774 0.01042697 1.1197874# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.155 0.0790 Inf 0.0532 0.374
## 1 0.140 0.0743 Inf 0.0466 0.353
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_chromogenes~ case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_chromogenes ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.40895 -0.40895 -0.35306 -0.04782 1.64694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40895 0.10612 3.854 0.000261 ***
## case -0.05589 0.14276 -0.391 0.696652
## farmB -0.30524 0.17938 -1.702 0.093394 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6014 on 68 degrees of freedom
## Multiple R-squared: 0.04304, Adjusted R-squared: 0.01489
## F-statistic: 1.529 on 2 and 68 DF, p-value: 0.2241Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.40894909 0.1971829 0.62071529
## case -0.05588883 -0.3407536 0.22897599
## farmB -0.30524054 -0.6631908 0.05270968# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.256 0.114 68 0.0284 0.484
## 1 0.200 0.115 68 -0.0291 0.430
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95# --- haemolyticus ----
# Presence
model <- glm(haemolyticus_present_sample~ case+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.4127831 0.17856064 0.8778753
## case 0.8906506 0.28708487 2.7259299
## farmB 0.1971774 0.01042697 1.1197874# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.155 0.0790 Inf 0.0532 0.374
## 1 0.140 0.0743 Inf 0.0466 0.353
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_haemolyticus~ case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_haemolyticus ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36521 -0.36521 -0.35688 -0.04576 2.00586
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.356885 0.116878 3.053 0.00323 **
## case 0.008324 0.157223 0.053 0.95793
## farmB -0.311121 0.197560 -1.575 0.11994
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6623 on 68 degrees of freedom
## Multiple R-squared: 0.03521, Adjusted R-squared: 0.006836
## F-statistic: 1.241 on 2 and 68 DF, p-value: 0.2956Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.356884718 0.1236576 0.59011179
## case 0.008324475 -0.3054092 0.32205814
## farmB -0.311120527 -0.7053462 0.08310518# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.201 0.126 68 -0.0497 0.452
## 1 0.210 0.127 68 -0.0431 0.462
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95# --- sciuri ----
# Presence
model <- glm(sciuri_present_sample~ case+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.06700151 0.01270191 0.2193419
## case 0.65605763 0.08101237 4.2777039
## farmB 3.02285622 0.36704516 20.3698904# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.104 0.0583 Inf 0.0332 0.284
## 1 0.071 0.0489 Inf 0.0175 0.247
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_sciuri~ case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_sciuri ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31153 -0.09819 -0.09819 -0.08439 2.35123
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08439 0.09014 0.936 0.352
## case 0.01380 0.12126 0.114 0.910
## farmB 0.21334 0.15237 1.400 0.166
##
## Residual standard error: 0.5108 on 68 degrees of freedom
## Multiple R-squared: 0.02823, Adjusted R-squared: -0.0003477
## F-statistic: 0.9878 on 2 and 68 DF, p-value: 0.3777Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.08439152 -0.09548714 0.2642702
## case 0.01379549 -0.22817467 0.2557656
## farmB 0.21334057 -0.09070987 0.5173910# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.191 0.0970 68 -0.00257 0.385
## 1 0.205 0.0977 68 0.00991 0.400
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95# --- succinus ----
# Presence
model <- glm(succinus_present_sample~ case,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.125000 0.03722053 0.314947
## case 1.032258 0.22598700 4.715691# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.111 0.0524 Inf 0.0423 0.261
## 1 0.114 0.0538 Inf 0.0436 0.268
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_succinus~ case,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_succinus ~ case, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1694 -0.1694 -0.1316 -0.1316 1.9920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.16936 0.08155 2.077 0.0416 *
## case -0.03775 0.11615 -0.325 0.7462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4893 on 69 degrees of freedom
## Multiple R-squared: 0.001528, Adjusted R-squared: -0.01294
## F-statistic: 0.1056 on 1 and 69 DF, p-value: 0.7462Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.16936464 0.006670748 0.3320585
## case -0.03774607 -0.269467609 0.1939755# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.169 0.0816 69 0.00667 0.332
## 1 0.132 0.0827 69 -0.03338 0.297
##
## Confidence level used: 0.95# --- xylosus_pseudoxylosus ----
# Presence
model <- glm(xylosus_pseudoxylosus_present_sample~ case,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.125000 0.03722053 0.314947
## case 1.655172 0.42986656 7.028589# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.111 0.0524 Inf 0.0423 0.261
## 1 0.171 0.0637 Inf 0.0791 0.333
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_xylosus_pseudoxylosus~ case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_xylosus_pseudoxylosus ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97956 -0.08558 -0.01601 -0.01601 2.07579
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.01601 0.08799 0.182 0.856
## case 0.06957 0.11836 0.588 0.559
## farmB 0.89398 0.14873 6.011 8.14e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4986 on 68 degrees of freedom
## Multiple R-squared: 0.3495, Adjusted R-squared: 0.3303
## F-statistic: 18.26 on 2 and 68 DF, p-value: 4.483e-07Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.01600678 -0.1595714 0.1915850
## case 0.06956969 -0.1666155 0.3057549
## farmB 0.89397977 0.5971985 1.1907611# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.463 0.0947 68 0.274 0.652
## 1 0.533 0.0954 68 0.342 0.723
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95# --- other_nasm ----
# Presence
model <- glm(other_nasm_present_sample~ case+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.05566145 0.008197893 0.1999215
## case 0.30733264 0.014433101 2.6340789
## farmB 4.77952519 0.520106214 44.5598020# Adj risk
emmeans(model,~case,type="response")## case prob SE df asymp.LCL asymp.UCL
## 0 0.1085 0.0601 Inf 0.03473 0.292
## 1 0.0361 0.0356 Inf 0.00499 0.218
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Counts
model <- lm(log_other_nasm~ case+farm,data=df_sau_per_cow_2)
summary(model)##
## Call:
## lm(formula = log_other_nasm ~ case + farm, data = df_sau_per_cow_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30972 -0.08496 -0.08496 -0.03442 2.07050
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08496 0.07047 1.206 0.2321
## case -0.05054 0.09479 -0.533 0.5957
## farmB 0.22476 0.11911 1.887 0.0634 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3993 on 68 degrees of freedom
## Multiple R-squared: 0.05333, Adjusted R-squared: 0.02549
## F-statistic: 1.915 on 2 and 68 DF, p-value: 0.1551Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 0.08496000 -0.05565743 0.2255774
## case -0.05053538 -0.23969190 0.1386211
## farmB 0.22475563 -0.01293121 0.4624425# Adjusted means
emmeans(model,~case)## case emmean SE df lower.CL upper.CL
## 0 0.197 0.0759 68 0.0460 0.349
## 1 0.147 0.0764 68 -0.0056 0.299
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.952 Isolates identified on teat apex
# ---- Exploratory analysis ----
# Identified isolates from cases and controls
table1(~ as.factor(primary_result) | farm:case_control,data=ln_4)A |
B |
Overall |
||||
|---|---|---|---|---|---|---|
| Case (N=103) |
Control (N=101) |
Case (N=38) |
Control (N=31) |
Case (N=141) |
Control (N=132) |
|
| as.factor(primary_result) | ||||||
| Aerococcus sp. | 0 (0%) | 0 (0%) | 2 (5.3%) | 0 (0%) | 2 (1.4%) | 0 (0%) |
| Aerococcus viridans | 3 (2.9%) | 0 (0%) | 1 (2.6%) | 1 (3.2%) | 4 (2.8%) | 1 (0.8%) |
| Bacillus sp. | 61 (59.2%) | 65 (64.4%) | 26 (68.4%) | 18 (58.1%) | 87 (61.7%) | 83 (62.9%) |
| Gram Positive Rod | 7 (6.8%) | 6 (5.9%) | 0 (0%) | 0 (0%) | 7 (5.0%) | 6 (4.5%) |
| Staphylococcus aureus | 2 (1.9%) | 0 (0%) | 0 (0%) | 0 (0%) | 2 (1.4%) | 0 (0%) |
| Staphylococcus chromogenes | 8 (7.8%) | 11 (10.9%) | 0 (0%) | 1 (3.2%) | 8 (5.7%) | 12 (9.1%) |
| Staphylococcus haemolyticus | 7 (6.8%) | 7 (6.9%) | 0 (0%) | 1 (3.2%) | 7 (5.0%) | 8 (6.1%) |
| Staphylococcus sp. | 12 (11.7%) | 12 (11.9%) | 6 (15.8%) | 4 (12.9%) | 18 (12.8%) | 16 (12.1%) |
| Staphylococcus xylosus/saprophyticus | 3 (2.9%) | 0 (0%) | 3 (7.9%) | 6 (19.4%) | 6 (4.3%) | 6 (4.5%) |
table1(~ as.factor(wgs_species) | farm:case_control,data=ln_4)A |
B |
Overall |
||||
|---|---|---|---|---|---|---|
| Case (N=103) |
Control (N=101) |
Case (N=38) |
Control (N=31) |
Case (N=141) |
Control (N=132) |
|
| as.factor(wgs_species) | ||||||
| Staphylococcus aureus | 2 (1.9%) | 0 (0%) | 0 (0%) | 0 (0%) | 2 (1.4%) | 0 (0%) |
| Staphylococcus chromogenes | 10 (9.7%) | 11 (10.9%) | 0 (0%) | 1 (3.2%) | 10 (7.1%) | 12 (9.1%) |
| Staphylococcus cohnii | 0 (0%) | 0 (0%) | 1 (2.6%) | 0 (0%) | 1 (0.7%) | 0 (0%) |
| Staphylococcus devriesei | 0 (0%) | 2 (2.0%) | 0 (0%) | 0 (0%) | 0 (0%) | 2 (1.5%) |
| Staphylococcus equorum | 0 (0%) | 0 (0%) | 0 (0%) | 1 (3.2%) | 0 (0%) | 1 (0.8%) |
| Staphylococcus haemolyticus | 9 (8.7%) | 9 (8.9%) | 0 (0%) | 1 (3.2%) | 9 (6.4%) | 10 (7.6%) |
| Staphylococcus pseudoxylosus | 2 (1.9%) | 0 (0%) | 1 (2.6%) | 2 (6.5%) | 3 (2.1%) | 2 (1.5%) |
| Staphylococcus sciuri | 2 (1.9%) | 2 (2.0%) | 1 (2.6%) | 1 (3.2%) | 3 (2.1%) | 3 (2.3%) |
| Staphylococcus succinus | 6 (5.8%) | 6 (5.9%) | 0 (0%) | 0 (0%) | 6 (4.3%) | 6 (4.5%) |
| Staphylococcus xylosus | 1 (1.0%) | 0 (0%) | 6 (15.8%) | 6 (19.4%) | 7 (5.0%) | 6 (4.5%) |
| Missing | 71 (68.9%) | 71 (70.3%) | 29 (76.3%) | 19 (61.3%) | 100 (70.9%) | 90 (68.2%) |
table1(~ as.factor(primary_result) | case_control,data=ln_4)| Case (N=141) |
Control (N=132) |
Overall (N=273) |
|
|---|---|---|---|
| as.factor(primary_result) | |||
| Aerococcus sp. | 2 (1.4%) | 0 (0%) | 2 (0.7%) |
| Aerococcus viridans | 4 (2.8%) | 1 (0.8%) | 5 (1.8%) |
| Bacillus sp. | 87 (61.7%) | 83 (62.9%) | 170 (62.3%) |
| Gram Positive Rod | 7 (5.0%) | 6 (4.5%) | 13 (4.8%) |
| Staphylococcus aureus | 2 (1.4%) | 0 (0%) | 2 (0.7%) |
| Staphylococcus chromogenes | 8 (5.7%) | 12 (9.1%) | 20 (7.3%) |
| Staphylococcus haemolyticus | 7 (5.0%) | 8 (6.1%) | 15 (5.5%) |
| Staphylococcus sp. | 18 (12.8%) | 16 (12.1%) | 34 (12.5%) |
| Staphylococcus xylosus/saprophyticus | 6 (4.3%) | 6 (4.5%) | 12 (4.4%) |
table1(~ as.factor(wgs_species) | case_control,data=ln_4)| Case (N=141) |
Control (N=132) |
Overall (N=273) |
|
|---|---|---|---|
| as.factor(wgs_species) | |||
| Staphylococcus aureus | 2 (1.4%) | 0 (0%) | 2 (0.7%) |
| Staphylococcus chromogenes | 10 (7.1%) | 12 (9.1%) | 22 (8.1%) |
| Staphylococcus cohnii | 1 (0.7%) | 0 (0%) | 1 (0.4%) |
| Staphylococcus devriesei | 0 (0%) | 2 (1.5%) | 2 (0.7%) |
| Staphylococcus equorum | 0 (0%) | 1 (0.8%) | 1 (0.4%) |
| Staphylococcus haemolyticus | 9 (6.4%) | 10 (7.6%) | 19 (7.0%) |
| Staphylococcus pseudoxylosus | 3 (2.1%) | 2 (1.5%) | 5 (1.8%) |
| Staphylococcus sciuri | 3 (2.1%) | 3 (2.3%) | 6 (2.2%) |
| Staphylococcus succinus | 6 (4.3%) | 6 (4.5%) | 12 (4.4%) |
| Staphylococcus xylosus | 7 (5.0%) | 6 (4.5%) | 13 (4.8%) |
| Missing | 100 (70.9%) | 90 (68.2%) | 190 (69.6%) |
# Isolates with antimicrobial activity
table1(~ as.factor(inhibition_sau) | case_control,data=df_sau)| Control (N=40) |
Case (N=37) |
Overall (N=77) |
|
|---|---|---|---|
| as.factor(inhibition_sau) | |||
| 0 | 3 (7.5%) | 2 (5.4%) | 5 (6.5%) |
| 1 | 37 (92.5%) | 35 (94.6%) | 72 (93.5%) |
table1(~ as.factor(inhibition_sau) | wgs_species_2,data=df_sau)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| as.factor(inhibition_sau) | |||||||
| 0 | 0 (0%) | 4 (21.1%) | 0 (0%) | 0 (0%) | 0 (0%) | 1 (25.0%) | 5 (6.5%) |
| 1 | 19 (100%) | 15 (78.9%) | 6 (100%) | 11 (100%) | 18 (100%) | 3 (75.0%) | 72 (93.5%) |
table1(~ as.factor(inhibition_sub) | case_control,data=df_sub)| Control (N=40) |
Case (N=37) |
Overall (N=77) |
|
|---|---|---|---|
| as.factor(inhibition_sub) | |||
| 0 | 8 (20.0%) | 7 (18.9%) | 15 (19.5%) |
| 1 | 32 (80.0%) | 30 (81.1%) | 62 (80.5%) |
table1(~ as.factor(inhibition_sub) | wgs_species_2,data=df_sub)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| as.factor(inhibition_sub) | |||||||
| 0 | 0 (0%) | 13 (68.4%) | 0 (0%) | 0 (0%) | 0 (0%) | 2 (50.0%) | 15 (19.5%) |
| 1 | 19 (100%) | 6 (31.6%) | 6 (100%) | 11 (100%) | 18 (100%) | 2 (50.0%) | 62 (80.5%) |
# Isolates with antimicrobial activity
table1(~ as.factor(high_ama) | case_control,data=df_sau)| Control (N=40) |
Case (N=37) |
Overall (N=77) |
|
|---|---|---|---|
| as.factor(high_ama) | |||
| 0 | 32 (80.0%) | 35 (94.6%) | 67 (87.0%) |
| 1 | 8 (20.0%) | 2 (5.4%) | 10 (13.0%) |
table1(~ as.factor(high_ama) | wgs_species_2,data=df_sau)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| as.factor(high_ama) | |||||||
| 0 | 17 (89.5%) | 18 (94.7%) | 5 (83.3%) | 6 (54.5%) | 17 (94.4%) | 4 (100%) | 67 (87.0%) |
| 1 | 2 (10.5%) | 1 (5.3%) | 1 (16.7%) | 5 (45.5%) | 1 (5.6%) | 0 (0%) | 10 (13.0%) |
table1(~ as.factor(high_ama) | case_control,data=df_sub)| Control (N=40) |
Case (N=37) |
Overall (N=77) |
|
|---|---|---|---|
| as.factor(high_ama) | |||
| 0 | 36 (90.0%) | 31 (83.8%) | 67 (87.0%) |
| 1 | 4 (10.0%) | 6 (16.2%) | 10 (13.0%) |
table1(~ as.factor(high_ama) | wgs_species_2,data=df_sub)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| as.factor(high_ama) | |||||||
| 0 | 17 (89.5%) | 18 (94.7%) | 6 (100%) | 6 (54.5%) | 16 (88.9%) | 4 (100%) | 67 (87.0%) |
| 1 | 2 (10.5%) | 1 (5.3%) | 0 (0%) | 5 (45.5%) | 2 (11.1%) | 0 (0%) | 10 (13.0%) |
# MIC across different species
table1(~ log_nas_media | wgs_species_2,data=df_sau)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| log_nas_media | |||||||
| Mean (SD) | 3.32 (0.497) | 3.95 (0.603) | 3.32 (0.560) | 2.92 (0.690) | 3.05 (0.416) | 4.39 (0.220) | 3.41 (0.678) |
| Median [Min, Max] | 3.25 [2.59, 4.51] | 4.09 [2.16, 4.66] | 3.49 [2.60, 4.03] | 3.03 [1.92, 4.11] | 2.91 [2.61, 4.09] | 4.44 [4.08, 4.59] | 3.28 [1.92, 4.66] |
table1(~ log_nas_media | wgs_species_2,data=df_sub)| Staphylococcus chromogenes (N=19) |
Staphylococcus haemolyticus (N=19) |
Staphylococcus sciuri (N=6) |
Staphylococcus succinus (N=11) |
Staphylococcus xylosus/pseudoxylosus (N=18) |
Other (N=4) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|
| log_nas_media | |||||||
| Mean (SD) | 3.06 (0.604) | 4.06 (0.720) | 3.49 (0.480) | 2.92 (0.930) | 3.05 (0.637) | 3.89 (0.714) | 3.36 (0.812) |
| Median [Min, Max] | 3.00 [2.51, 5.18] | 4.04 [2.16, 5.14] | 3.58 [2.60, 4.03] | 2.58 [1.65, 4.14] | 2.97 [1.71, 4.09] | 3.94 [3.08, 4.59] | 3.24 [1.65, 5.18] |
# MIC across different SAU and SUB categories
table1(~ log_nas_media | inhibition_category,data=df_sau)| Top 10 (N=10) |
11-20 (N=10) |
21-30 (N=10) |
31-40 (N=10) |
41-50 (N=10) |
51-60 (N=10) |
61-70 (N=10) |
>70 (N=2) |
No inhibition (N=5) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|---|---|---|
| log_nas_media | ||||||||||
| Mean (SD) | 2.41 (0.263) | 2.76 (0.0489) | 3.02 (0.0786) | 3.25 (0.104) | 3.59 (0.0517) | 4.00 (0.121) | 4.34 (0.166) | 4.63 (0.0547) | 3.98 (0.453) | 3.41 (0.678) |
| Median [Min, Max] | 2.55 [1.92, 2.66] | 2.75 [2.70, 2.85] | 3.02 [2.90, 3.14] | 3.20 [3.16, 3.45] | 3.59 [3.51, 3.68] | 4.04 [3.75, 4.11] | 4.34 [4.14, 4.56] | 4.63 [4.59, 4.66] | 4.04 [3.24, 4.38] | 3.28 [1.92, 4.66] |
table1(~ log_nas_media | inhibition_category,data=df_sub)| Top 10 (N=10) |
11-20 (N=10) |
21-30 (N=10) |
31-40 (N=9) |
41-50 (N=9) |
51-60 (N=8) |
61-70 (N=5) |
>70 (N=1) |
No inhibition (N=15) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|---|---|---|
| log_nas_media | ||||||||||
| Mean (SD) | 2.16 (0.324) | 2.63 (0.0575) | 2.91 (0.132) | 3.17 (0.0494) | 3.49 (0.119) | 3.88 (0.121) | 4.13 (0.0851) | 5.18 (NA) | 4.33 (0.519) | 3.36 (0.812) |
| Median [Min, Max] | 2.17 [1.65, 2.56] | 2.61 [2.56, 2.74] | 2.91 [2.74, 3.08] | 3.17 [3.09, 3.25] | 3.50 [3.28, 3.68] | 3.89 [3.70, 4.03] | 4.11 [4.04, 4.27] | 5.18 [5.18, 5.18] | 4.38 [3.24, 5.14] | 3.24 [1.65, 5.18] |
# About ATCC
table1(~nas_count_hundred + log(nas_count_hundred,10)|isolate,data=sau_inhibition_atcc)| atcc_sau (N=1) |
atcc_sub (N=1) |
Overall (N=2) |
|
|---|---|---|---|
| nas_count_hundred | |||
| Mean (SD) | 2640 (NA) | 2640 (NA) | 2640 (0) |
| Median [Min, Max] | 2640 [2640, 2640] | 2640 [2640, 2640] | 2640 [2640, 2640] |
| log(nas_count_hundred, 10) | |||
| Mean (SD) | 3.42 (NA) | 3.42 (NA) | 3.42 (0) |
| Median [Min, Max] | 3.42 [3.42, 3.42] | 3.42 [3.42, 3.42] | 3.42 [3.42, 3.42] |
table(df_sau$log_nas_media<3.42)##
## FALSE TRUE
## 37 40prop.table(table(df_sau$log_nas_media<3.42))##
## FALSE TRUE
## 0.4805195 0.5194805table(df_sau$log_nas_media<3.42,df_sau$case_control)##
## Control Case
## FALSE 19 18
## TRUE 21 19table(df_sau$log_nas_media<3.42,df_sau$primary_result)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## FALSE 6 12
## TRUE 11 3
##
## Staphylococcus sp. Staphylococcus xylosus/saprophyticus
## FALSE 17 2
## TRUE 16 10table1(~nas_count_hundred + log(nas_count_hundred,10)|isolate,data=sub_inhibition_atcc)| atcc_sau (N=1) |
atcc_sub (N=1) |
Overall (N=2) |
|
|---|---|---|---|
| nas_count_hundred | |||
| Mean (SD) | 817 (NA) | 817 (NA) | 817 (0) |
| Median [Min, Max] | 817 [817, 817] | 817 [817, 817] | 817 [817, 817] |
| log(nas_count_hundred, 10) | |||
| Mean (SD) | 2.91 (NA) | 2.91 (NA) | 2.91 (0) |
| Median [Min, Max] | 2.91 [2.91, 2.91] | 2.91 [2.91, 2.91] | 2.91 [2.91, 2.91] |
table(df_sub$log_nas_media<2.91)##
## FALSE TRUE
## 52 25prop.table(table(df_sub$log_nas_media<2.91))##
## FALSE TRUE
## 0.6753247 0.3246753table(df_sub$log_nas_media<2.91,df_sub$case_control)##
## Control Case
## FALSE 26 26
## TRUE 14 11table(df_sub$log_nas_media<2.91,df_sub$primary_result)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## FALSE 9 14
## TRUE 8 1
##
## Staphylococcus sp. Staphylococcus xylosus/saprophyticus
## FALSE 23 6
## TRUE 10 63 In vitro inhibitory activity
3.1 In vitro inhibitory activity and presence of SAU-SSLO-IMI
# ---- Analysis of relationship between presence of SAU-SSLO-IMI and MIC ----
# MIC for sau not different between cases and controls
model <- lmer(log_nas_media~as.factor(case_control)+farm+(1|cow_id),data=df_sau_inhibition)
summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_nas_media ~ as.factor(case_control) + farm + (1 | cow_id)
## Data: df_sau_inhibition
##
## REML criterion at convergence: 148.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2277 -0.6211 -0.1400 0.8946 2.1059
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 0.0114 0.1068
## Residual 0.4281 0.6543
## Number of obs: 72, groups: cow_id, 43
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.4157 0.1221 36.4085 27.971 <2e-16 ***
## as.factor(case_control)Case 0.1244 0.1590 27.2563 0.782 0.4408
## farmB -0.3554 0.1782 19.5567 -1.994 0.0603 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(_)C
## as.fctr(_)C -0.647
## farmB -0.429 0.047Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 3.4156991 3.1763525 3.655045624
## as.factor(case_control)Case 0.1243566 -0.1871830 0.435896209
## farmB -0.3554036 -0.7047117 -0.006095431emmeans(model,pairwise~case_control)## $emmeans
## case_control emmean SE df lower.CL upper.CL
## Control 3.24 0.118 28.5 3.00 3.48
## Case 3.36 0.130 18.2 3.09 3.64
##
## Results are averaged over the levels of: farm
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control - Case -0.124 0.165 27.3 -0.753 0.4579
##
## Results are averaged over the levels of: farm
## Degrees-of-freedom method: kenward-rogerperformance::icc(model,by_group=TRUE)## # ICC by Group
##
## Group | ICC
## --------------
## cow_id | 0.026# Survival analysis
df_sau$survival <- Surv(event = df_sau$inhibition_sau, time = df_sau$log_nas_media, type="right")
kmsurv_sau_1 <- survfit(survival ~ case_control , data = df_sau)
survplots <- list()
km_1 <- ggsurvplot(kmsurv_sau_1, data = df_sau,
title = "Staphylococcus aureus",
conf.int = FALSE,
risk.table = FALSE,
surv.plot.height = 5,
legend.labs = c("No", "Yes"),
ggtheme = theme_classic(),
legend.title = "IMI",
ylim = c(0, 1),
xlim = c(2,5),
xlab = "NASM (log10 CFU/mL)",
ylab = expression(italic("Staphylococcus aureus") ~ " growth (%)"),
censor = FALSE,
palette = c("#009999", "#990000"),
surv.scale = "percent") # Converts proportions to percentages
km_1_plot <- km_1$plot + theme(plot.title = element_text(face = "italic"))
model <- coxph(survival~ case_control+farm+cluster(cow_id), data=df_sau)
summary(model)## Call:
## coxph(formula = survival ~ case_control + farm, data = df_sau,
## cluster = cow_id)
##
## n= 77, number of events= 72
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## case_controlCase -0.1072 0.8984 0.2436 0.2440 -0.439 0.661
## farmB 0.4310 1.5388 0.2689 0.3178 1.356 0.175
##
## exp(coef) exp(-coef) lower .95 upper .95
## case_controlCase 0.8984 1.1131 0.5569 1.449
## farmB 1.5388 0.6499 0.8254 2.869
##
## Concordance= 0.583 (se = 0.034 )
## Likelihood ratio test= 2.77 on 2 df, p=0.2
## Wald test = 2.72 on 2 df, p=0.3
## Score (logrank) test = 3 on 2 df, p=0.2, Robust = 2.97 p=0.2
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).cox.zph(model)## chisq df p
## case_control 0.0647 1 0.80
## farm 2.5777 1 0.11
## GLOBAL 2.6531 2 0.27df_sub$survival <- Surv(event = df_sub$inhibition_sub, time = df_sub$log_nas_media, type="right")
kmsurv_sub_1 <- survfit(survival ~ case_control , data = df_sub)
km_2 <- ggsurvplot(kmsurv_sub_1, data = df_sub,
title = "Streptococcus uberis",
conf.int = FALSE,
risk.table = FALSE,
surv.plot.height = 5,
legend.labs = c("No", "Yes"),
ggtheme = theme_classic(),
legend.title = "IMI",
ylim = c(0, 1),
xlim = c(0,5),
xlab = "NASM (log10 CFU/mL)",
ylab = expression(italic("Streptococcus uberis") ~ " growth (%)"),
censor = FALSE,
palette = c("#009999", "#990000"),
surv.scale = "percent") # Converts proportions to percentages
km_2_plot <- km_2$plot + theme(plot.title = element_text(face = "italic"))
# Cox regression model
model <- coxph(survival~ case_control+farm+cluster(cow_id), data=df_sub)
# Summary model
summary(model)## Call:
## coxph(formula = survival ~ case_control + farm, data = df_sub,
## cluster = cow_id)
##
## n= 77, number of events= 62
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## case_controlCase 0.1570 1.1700 0.2607 0.2658 0.591 0.555
## farmB 0.3582 1.4307 0.2820 0.2961 1.210 0.226
##
## exp(coef) exp(-coef) lower .95 upper .95
## case_controlCase 1.170 0.8547 0.6949 1.970
## farmB 1.431 0.6990 0.8008 2.556
##
## Concordance= 0.528 (se = 0.043 )
## Likelihood ratio test= 1.7 on 2 df, p=0.4
## Wald test = 1.65 on 2 df, p=0.4
## Score (logrank) test = 1.78 on 2 df, p=0.4, Robust = 1.34 p=0.5
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).# Proportional hazards assumption
cox.zph(model)## chisq df p
## case_control 0.237 1 0.63
## farm 0.111 1 0.74
## GLOBAL 0.400 2 0.82# ---- MIC for sub not different between cases and controls ----
model <- lmer(log_nas_media~as.factor(case_control)+farm+(1|cow_id),data=df_sub_inhibition)## boundary (singular) fit: see help('isSingular')summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_nas_media ~ as.factor(case_control) + farm + (1 | cow_id)
## Data: df_sub_inhibition
##
## REML criterion at convergence: 134.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18824 -0.72571 -0.06739 0.70970 2.88008
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 0.0000 0.0000
## Residual 0.4847 0.6962
## Number of obs: 62, groups: cow_id, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.17237 0.13973 59.00000 22.704 <2e-16 ***
## as.factor(case_control)Case -0.05991 0.17755 59.00000 -0.337 0.737
## farmB -0.05825 0.19246 59.00000 -0.303 0.763
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(_)C
## as.fctr(_)C -0.650
## farmB -0.473 0.084
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Confint(model)## Estimate 2.5 % 97.5 %
## (Intercept) 3.17237250 2.8985074 3.4462376
## as.factor(case_control)Case -0.05990544 -0.4079000 0.2880891
## farmB -0.05824505 -0.4354671 0.3189770emmeans(model,~case_control)## case_control emmean SE df lower.CL upper.CL
## Control 3.14 0.129 25.2 2.88 3.41
## Case 3.08 0.147 12.9 2.77 3.40
##
## Results are averaged over the levels of: farm
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95performance::icc(model,by_group=TRUE)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## [1] NAkm_plots_1 <- ggarrange(km_1_plot, km_2_plot, ncol = 2, nrow = 1, common.legend = TRUE,align = "h",labels = c("A", "B"))
km_plots_1
ggsave(plot = last_plot(),"./figures/km_1_annotated.png",width = 20, height = 10, units = "cm") 3.2 In vitro inhibitory activity across different NASM species
# ---- MIC for sau was different between different species ----
model <- lmer(log_nas_media~wgs_species_2+(1|cow_id),data=df_sau_inhibition)## boundary (singular) fit: see help('isSingular')summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_nas_media ~ wgs_species_2 + (1 | cow_id)
## Data: df_sau_inhibition
##
## REML criterion at convergence: 121.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3079 -0.5624 0.1201 0.4621 2.1815
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 0.0000 0.0000
## Residual 0.2995 0.5473
## Number of obs: 72, groups: cow_id, 43
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 3.320351 0.125551
## wgs_species_2Staphylococcus haemolyticus 0.652940 0.189022
## wgs_species_2Staphylococcus sciuri 0.002349 0.256279
## wgs_species_2Staphylococcus succinus -0.397583 0.207340
## wgs_species_2Staphylococcus xylosus/pseudoxylosus -0.266120 0.180005
## wgs_species_2Other 1.068612 0.339993
## df t value Pr(>|t|)
## (Intercept) 66.000000 26.446 < 2e-16
## wgs_species_2Staphylococcus haemolyticus 66.000000 3.454 0.000969
## wgs_species_2Staphylococcus sciuri 66.000000 0.009 0.992715
## wgs_species_2Staphylococcus succinus 66.000000 -1.918 0.059499
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 66.000000 -1.478 0.144056
## wgs_species_2Other 66.000000 3.143 0.002505
##
## (Intercept) ***
## wgs_species_2Staphylococcus haemolyticus ***
## wgs_species_2Staphylococcus sciuri
## wgs_species_2Staphylococcus succinus .
## wgs_species_2Staphylococcus xylosus/pseudoxylosus
## wgs_species_2Other **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) w__2Sh wgs_spcs_2Stphylcccsscr
## wgs_spc_2Sh -0.664
## wgs_spcs_2Stphylcccsscr -0.490 0.325
## wgs_spcs_2Stphylcccsscc -0.606 0.402 0.297
## wgs_sp_2Sx/ -0.697 0.463 0.342
## wgs_spcs_2O -0.369 0.245 0.181
## wgs_spcs_2Stphylcccsscc w__2Sx
## wgs_spc_2Sh
## wgs_spcs_2Stphylcccsscr
## wgs_spcs_2Stphylcccsscc
## wgs_sp_2Sx/ 0.422
## wgs_spcs_2O 0.224 0.258
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Anova(model,type="III")## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log_nas_media
## Chisq Df Pr(>Chisq)
## (Intercept) 699.406 1 < 2.2e-16 ***
## wgs_species_2 42.166 5 5.453e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model,pairwise~wgs_species_2)## $emmeans
## wgs_species_2 emmean SE df lower.CL upper.CL
## Staphylococcus chromogenes 3.32 0.128 50.2 3.06 3.58
## Staphylococcus haemolyticus 3.97 0.143 64.0 3.69 4.26
## Staphylococcus sciuri 3.32 0.239 58.0 2.84 3.80
## Staphylococcus succinus 2.92 0.174 41.9 2.57 3.27
## Staphylococcus xylosus/pseudoxylosus 3.05 0.137 23.3 2.77 3.34
## Other 4.39 0.322 65.9 3.75 5.03
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate
## Staphylococcus chromogenes - Staphylococcus haemolyticus -0.65294
## Staphylococcus chromogenes - Staphylococcus sciuri -0.00235
## Staphylococcus chromogenes - Staphylococcus succinus 0.39758
## Staphylococcus chromogenes - (Staphylococcus xylosus/pseudoxylosus) 0.26612
## Staphylococcus chromogenes - Other -1.06861
## Staphylococcus haemolyticus - Staphylococcus sciuri 0.65059
## Staphylococcus haemolyticus - Staphylococcus succinus 1.05052
## Staphylococcus haemolyticus - (Staphylococcus xylosus/pseudoxylosus) 0.91906
## Staphylococcus haemolyticus - Other -0.41567
## Staphylococcus sciuri - Staphylococcus succinus 0.39993
## Staphylococcus sciuri - (Staphylococcus xylosus/pseudoxylosus) 0.26847
## Staphylococcus sciuri - Other -1.06626
## Staphylococcus succinus - (Staphylococcus xylosus/pseudoxylosus) -0.13146
## Staphylococcus succinus - Other -1.46619
## (Staphylococcus xylosus/pseudoxylosus) - Other -1.33473
## SE df t.ratio p.value
## 0.192 60.9 -3.407 0.0142
## 0.272 59.2 -0.009 1.0000
## 0.217 53.2 1.835 0.4527
## 0.187 34.1 1.420 0.7147
## 0.346 65.1 -3.087 0.0337
## 0.279 60.0 2.333 0.1972
## 0.225 52.3 4.664 0.0003
## 0.199 48.9 4.625 0.0004
## 0.352 65.9 -1.180 0.8446
## 0.285 60.6 1.402 0.7258
## 0.273 65.0 0.983 0.9216
## 0.402 64.5 -2.655 0.0988
## 0.224 33.5 -0.586 0.9913
## 0.366 63.5 -4.008 0.0022
## 0.351 64.8 -3.804 0.0041
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesperformance::icc(model,by_group=TRUE)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## [1] NAkmsurv_sau_2 <- survfit(survival ~ wgs_species_2 , data = df_sau)
km_3 <- ggsurvplot(kmsurv_sau_2, data = df_sub, title = "Staphylococcus aureus", conf.int = F,risk.table = F, surv.plot.height = 5, legend.labs = c("SCH","SHAEM","MSC","SSUC","SXYL","Other"), ggtheme = theme_classic(),legend.title="NASM taxonomy",ylim=c(0,1),xlim=c(2,5),xlab="NASM (log10 CFU/mL)",ylab=expression(italic("Staphylococcus aureus") ~ " growth (%)"),censor=F,palette = "Set1",surv.scale="percent")
km_3_plot <- km_3$plot + theme(plot.title = element_text(face = "italic"))
# Cox regression model
table(df_sau$inhibition_sau,df_sau$wgs_species_2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 0 4
## 1 19 15
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 0 0
## 1 6 11
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 0 1
## 1 18 3prop.table(table(df_sau$inhibition_sau,df_sau$wgs_species_2),2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 0.0000000 0.2105263
## 1 1.0000000 0.7894737
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 0.0000000 0.0000000
## 1 1.0000000 1.0000000
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 0.0000000 0.2500000
## 1 1.0000000 0.7500000model <- coxph(survival~ wgs_species_2+cluster(cow_id), data=df_sau)
# Summary model
summary(model)## Call:
## coxph(formula = survival ~ wgs_species_2, data = df_sau, cluster = cow_id)
##
## n= 77, number of events= 72
##
## coef exp(coef) se(coef)
## wgs_species_2Staphylococcus haemolyticus -1.3364 0.2628 0.3717
## wgs_species_2Staphylococcus sciuri 0.2582 1.2946 0.4789
## wgs_species_2Staphylococcus succinus 0.5869 1.7985 0.3887
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.6757 1.9655 0.3444
## wgs_species_2Other -1.9025 0.1492 0.6512
## robust se z Pr(>|z|)
## wgs_species_2Staphylococcus haemolyticus 0.4156 -3.216 0.001301 **
## wgs_species_2Staphylococcus sciuri 0.4480 0.576 0.564405
## wgs_species_2Staphylococcus succinus 0.4421 1.328 0.184334
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.3555 1.901 0.057353 .
## wgs_species_2Other 0.5480 -3.472 0.000518 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef)
## wgs_species_2Staphylococcus haemolyticus 0.2628 3.8055
## wgs_species_2Staphylococcus sciuri 1.2946 0.7725
## wgs_species_2Staphylococcus succinus 1.7985 0.5560
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 1.9655 0.5088
## wgs_species_2Other 0.1492 6.7028
## lower .95 upper .95
## wgs_species_2Staphylococcus haemolyticus 0.11637 0.5934
## wgs_species_2Staphylococcus sciuri 0.53804 3.1148
## wgs_species_2Staphylococcus succinus 0.75607 4.2781
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.97912 3.9456
## wgs_species_2Other 0.05096 0.4368
##
## Concordance= 0.713 (se = 0.033 )
## Likelihood ratio test= 41.29 on 5 df, p=8e-08
## Wald test = 52.27 on 5 df, p=5e-10
## Score (logrank) test = 40.6 on 5 df, p=1e-07, Robust = 27.19 p=5e-05
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).# Type III p-value
Anova(model,type="III")## Warning in Anova.coxph(model, type = "III"): LR tests unavailable with robust variances
## Wald tests substituted
## Warning in Anova.coxph(model, type = "III"): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III tests)
##
## Response: survival
## Df Chisq Pr(>Chisq)
## wgs_species_2 5 52.271 4.748e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Pairwise comparisons
confint(emmeans(model,pairwise~wgs_species_2,type="response"))$contrasts## contrast ratio
## Staphylococcus chromogenes / Staphylococcus haemolyticus 3.805
## Staphylococcus chromogenes / Staphylococcus sciuri 0.772
## Staphylococcus chromogenes / Staphylococcus succinus 0.556
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus) 0.509
## Staphylococcus chromogenes / Other 6.703
## Staphylococcus haemolyticus / Staphylococcus sciuri 0.203
## Staphylococcus haemolyticus / Staphylococcus succinus 0.146
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus) 0.134
## Staphylococcus haemolyticus / Other 1.761
## Staphylococcus sciuri / Staphylococcus succinus 0.720
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus) 0.659
## Staphylococcus sciuri / Other 8.677
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus) 0.915
## Staphylococcus succinus / Other 12.055
## (Staphylococcus xylosus/pseudoxylosus) / Other 13.174
## SE df asymp.LCL asymp.UCL
## 1.5815 Inf 1.1643 12.438
## 0.3460 Inf 0.2155 2.769
## 0.2458 Inf 0.1577 1.960
## 0.1809 Inf 0.1847 1.401
## 3.6734 Inf 1.4060 31.953
## 0.0741 Inf 0.0717 0.575
## 0.0544 Inf 0.0506 0.422
## 0.0487 Inf 0.0473 0.378
## 0.8069 Inf 0.4773 6.499
## 0.2327 Inf 0.2865 1.809
## 0.2269 Inf 0.2468 1.758
## 4.5771 Inf 1.9300 39.012
## 0.3528 Inf 0.3049 2.746
## 6.3862 Inf 2.6639 54.550
## 6.2485 Inf 3.4098 50.900
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 6 estimates
## Intervals are back-transformed from the log scaleemmeans(model,pairwise~wgs_species_2,type="response")$contrasts## contrast ratio
## Staphylococcus chromogenes / Staphylococcus haemolyticus 3.805
## Staphylococcus chromogenes / Staphylococcus sciuri 0.772
## Staphylococcus chromogenes / Staphylococcus succinus 0.556
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus) 0.509
## Staphylococcus chromogenes / Other 6.703
## Staphylococcus haemolyticus / Staphylococcus sciuri 0.203
## Staphylococcus haemolyticus / Staphylococcus succinus 0.146
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus) 0.134
## Staphylococcus haemolyticus / Other 1.761
## Staphylococcus sciuri / Staphylococcus succinus 0.720
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus) 0.659
## Staphylococcus sciuri / Other 8.677
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus) 0.915
## Staphylococcus succinus / Other 12.055
## (Staphylococcus xylosus/pseudoxylosus) / Other 13.174
## SE df null z.ratio p.value
## 1.5815 Inf 1 3.216 0.0164
## 0.3460 Inf 1 -0.576 0.9926
## 0.2458 Inf 1 -1.328 0.7698
## 0.1809 Inf 1 -1.901 0.4017
## 3.6734 Inf 1 3.472 0.0069
## 0.0741 Inf 1 -4.366 0.0002
## 0.0544 Inf 1 -5.169 <.0001
## 0.0487 Inf 1 -5.520 <.0001
## 0.8069 Inf 1 1.236 0.8194
## 0.2327 Inf 1 -1.017 0.9125
## 0.2269 Inf 1 -1.212 0.8310
## 4.5771 Inf 1 4.096 0.0006
## 0.3528 Inf 1 -0.230 0.9999
## 6.3862 Inf 1 4.699 <.0001
## 6.2485 Inf 1 5.436 <.0001
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log scaleconfint(emmeans(model,revpairwise~wgs_species_2,type="response"))$contrasts## contrast ratio
## Staphylococcus haemolyticus / Staphylococcus chromogenes 0.2628
## Staphylococcus sciuri / Staphylococcus chromogenes 1.2946
## Staphylococcus sciuri / Staphylococcus haemolyticus 4.9264
## Staphylococcus succinus / Staphylococcus chromogenes 1.7985
## Staphylococcus succinus / Staphylococcus haemolyticus 6.8441
## Staphylococcus succinus / Staphylococcus sciuri 1.3893
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus chromogenes 1.9655
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus haemolyticus 7.4797
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus sciuri 1.5183
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus succinus 1.0929
## Other / Staphylococcus chromogenes 0.1492
## Other / Staphylococcus haemolyticus 0.5677
## Other / Staphylococcus sciuri 0.1152
## Other / Staphylococcus succinus 0.0830
## Other / (Staphylococcus xylosus/pseudoxylosus) 0.0759
## SE df asymp.LCL asymp.UCL
## 0.1092 Inf 0.0804 0.859
## 0.5799 Inf 0.3612 4.640
## 1.7994 Inf 1.7398 13.950
## 0.7952 Inf 0.5102 6.340
## 2.5468 Inf 2.3702 19.763
## 0.4492 Inf 0.5529 3.491
## 0.6988 Inf 0.7136 5.414
## 2.7267 Inf 2.6468 21.137
## 0.5229 Inf 0.5690 4.052
## 0.4214 Inf 0.3642 3.279
## 0.0818 Inf 0.0313 0.711
## 0.2601 Inf 0.1539 2.095
## 0.0608 Inf 0.0256 0.518
## 0.0439 Inf 0.0183 0.375
## 0.0360 Inf 0.0196 0.293
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 6 estimates
## Intervals are back-transformed from the log scaleemmeans(model,revpairwise~wgs_species_2,type="response")$contrasts## contrast ratio
## Staphylococcus haemolyticus / Staphylococcus chromogenes 0.2628
## Staphylococcus sciuri / Staphylococcus chromogenes 1.2946
## Staphylococcus sciuri / Staphylococcus haemolyticus 4.9264
## Staphylococcus succinus / Staphylococcus chromogenes 1.7985
## Staphylococcus succinus / Staphylococcus haemolyticus 6.8441
## Staphylococcus succinus / Staphylococcus sciuri 1.3893
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus chromogenes 1.9655
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus haemolyticus 7.4797
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus sciuri 1.5183
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus succinus 1.0929
## Other / Staphylococcus chromogenes 0.1492
## Other / Staphylococcus haemolyticus 0.5677
## Other / Staphylococcus sciuri 0.1152
## Other / Staphylococcus succinus 0.0830
## Other / (Staphylococcus xylosus/pseudoxylosus) 0.0759
## SE df null z.ratio p.value
## 0.1092 Inf 1 -3.216 0.0164
## 0.5799 Inf 1 0.576 0.9926
## 1.7994 Inf 1 4.366 0.0002
## 0.7952 Inf 1 1.328 0.7698
## 2.5468 Inf 1 5.169 <.0001
## 0.4492 Inf 1 1.017 0.9125
## 0.6988 Inf 1 1.901 0.4017
## 2.7267 Inf 1 5.520 <.0001
## 0.5229 Inf 1 1.212 0.8310
## 0.4214 Inf 1 0.230 0.9999
## 0.0818 Inf 1 -3.472 0.0069
## 0.2601 Inf 1 -1.236 0.8194
## 0.0608 Inf 1 -4.096 0.0006
## 0.0439 Inf 1 -4.699 <.0001
## 0.0360 Inf 1 -5.436 <.0001
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log scale# Proportional hazards assumption
cox.zph(model)## chisq df p
## wgs_species_2 4.6 5 0.47
## GLOBAL 4.6 5 0.47# MIC for sub was different between different species
model <- lmer(log_nas_media~wgs_species_2 + (1|cow_id),data=df_sub_inhibition)## boundary (singular) fit: see help('isSingular')summary(model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_nas_media ~ wgs_species_2 + (1 | cow_id)
## Data: df_sub_inhibition
##
## REML criterion at convergence: 129.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96300 -0.59520 -0.05971 0.45991 3.09240
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 0.0000 0.0000
## Residual 0.4701 0.6857
## Number of obs: 62, groups: cow_id, 36
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 3.057193 0.157301
## wgs_species_2Staphylococcus haemolyticus 0.451753 0.321089
## wgs_species_2Staphylococcus sciuri 0.432173 0.321089
## wgs_species_2Staphylococcus succinus -0.134425 0.259774
## wgs_species_2Staphylococcus xylosus/pseudoxylosus -0.002962 0.225525
## wgs_species_2Other 0.232879 0.509713
## df t value Pr(>|t|)
## (Intercept) 56.000000 19.435 <2e-16
## wgs_species_2Staphylococcus haemolyticus 56.000000 1.407 0.165
## wgs_species_2Staphylococcus sciuri 56.000000 1.346 0.184
## wgs_species_2Staphylococcus succinus 56.000000 -0.517 0.607
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 56.000000 -0.013 0.990
## wgs_species_2Other 56.000000 0.457 0.650
##
## (Intercept) ***
## wgs_species_2Staphylococcus haemolyticus
## wgs_species_2Staphylococcus sciuri
## wgs_species_2Staphylococcus succinus
## wgs_species_2Staphylococcus xylosus/pseudoxylosus
## wgs_species_2Other
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) w__2Sh wgs_spcs_2Stphylcccsscr
## wgs_spc_2Sh -0.490
## wgs_spcs_2Stphylcccsscr -0.490 0.240
## wgs_spcs_2Stphylcccsscc -0.606 0.297 0.297
## wgs_sp_2Sx/ -0.697 0.342 0.342
## wgs_spcs_2O -0.309 0.151 0.151
## wgs_spcs_2Stphylcccsscc w__2Sx
## wgs_spc_2Sh
## wgs_spcs_2Stphylcccsscr
## wgs_spcs_2Stphylcccsscc
## wgs_sp_2Sx/ 0.422
## wgs_spcs_2O 0.187 0.215
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Anova(model,type="III")## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log_nas_media
## Chisq Df Pr(>Chisq)
## (Intercept) 377.7319 1 <2e-16 ***
## wgs_species_2 5.0261 5 0.4127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model,pairwise~wgs_species_2)## $emmeans
## wgs_species_2 emmean SE df lower.CL upper.CL
## Staphylococcus chromogenes 3.06 0.160 43.2 2.73 3.38
## Staphylococcus haemolyticus 3.51 0.283 55.9 2.94 4.08
## Staphylococcus sciuri 3.49 0.301 49.5 2.88 4.09
## Staphylococcus succinus 2.92 0.219 36.4 2.48 3.37
## Staphylococcus xylosus/pseudoxylosus 3.05 0.172 20.7 2.70 3.41
## Other 3.29 0.491 55.9 2.31 4.27
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate
## Staphylococcus chromogenes - Staphylococcus haemolyticus -0.45175
## Staphylococcus chromogenes - Staphylococcus sciuri -0.43217
## Staphylococcus chromogenes - Staphylococcus succinus 0.13442
## Staphylococcus chromogenes - (Staphylococcus xylosus/pseudoxylosus) 0.00296
## Staphylococcus chromogenes - Other -0.23288
## Staphylococcus haemolyticus - Staphylococcus sciuri 0.01958
## Staphylococcus haemolyticus - Staphylococcus succinus 0.58618
## Staphylococcus haemolyticus - (Staphylococcus xylosus/pseudoxylosus) 0.45472
## Staphylococcus haemolyticus - Other 0.21887
## Staphylococcus sciuri - Staphylococcus succinus 0.56660
## Staphylococcus sciuri - (Staphylococcus xylosus/pseudoxylosus) 0.43514
## Staphylococcus sciuri - Other 0.19929
## Staphylococcus succinus - (Staphylococcus xylosus/pseudoxylosus) -0.13146
## Staphylococcus succinus - Other -0.36730
## (Staphylococcus xylosus/pseudoxylosus) - Other -0.23584
## SE df t.ratio p.value
## 0.325 54.5 -1.389 0.7334
## 0.342 50.5 -1.263 0.8034
## 0.272 45.6 0.494 0.9962
## 0.235 29.9 0.013 1.0000
## 0.517 55.5 -0.451 0.9975
## 0.415 53.7 0.047 1.0000
## 0.358 51.6 1.637 0.5788
## 0.334 53.7 1.363 0.7484
## 0.567 54.1 0.386 0.9988
## 0.358 51.9 1.581 0.6142
## 0.343 55.1 1.267 0.8015
## 0.576 55.1 0.346 0.9993
## 0.282 29.4 -0.466 0.9970
## 0.538 54.7 -0.683 0.9831
## 0.521 54.3 -0.453 0.9975
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 6 estimatesperformance::icc(model,by_group=TRUE)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## [1] NAkmsurv_sub_2 <- survfit(survival ~ wgs_species_2 , data = df_sub)
km_4 <- ggsurvplot(kmsurv_sub_2, data = df_sub, title = "Streptococcus uberis", conf.int = F,risk.table = F, surv.plot.height = 5, legend.labs = c("SCH","SHAEM","MSC","SSUC","SXYL","Other"), ggtheme = theme_classic(),legend.title="NASM taxonomy",ylim=c(0,1),xlim=c(2,5),xlab="NASM (log10 CFU/mL)",ylab=expression(italic("Streptococcus uberis") ~ " growth (%)"),censor=F,palette = "Set1",surv.scale="percent")
km_4_plot <- km_4$plot + theme(plot.title = element_text(face = "italic"))
# Cox regression model
table(df_sub$inhibition_sub,df_sub$wgs_species_2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 0 13
## 1 19 6
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 0 0
## 1 6 11
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 0 2
## 1 18 2prop.table(table(df_sub$inhibition_sub,df_sub$wgs_species_2),2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 0.0000000 0.6842105
## 1 1.0000000 0.3157895
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 0.0000000 0.0000000
## 1 1.0000000 1.0000000
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 0.0000000 0.5000000
## 1 1.0000000 0.5000000model <- coxph(survival~ wgs_species_2+cluster(cow_id), data=df_sub)
# Summary model
summary(model)## Call:
## coxph(formula = survival ~ wgs_species_2, data = df_sub, cluster = cow_id)
##
## n= 77, number of events= 62
##
## coef exp(coef) se(coef)
## wgs_species_2Staphylococcus haemolyticus -2.26405 0.10393 0.48838
## wgs_species_2Staphylococcus sciuri -0.51368 0.59829 0.48000
## wgs_species_2Staphylococcus succinus -0.21751 0.80452 0.39682
## wgs_species_2Staphylococcus xylosus/pseudoxylosus -0.06414 0.93788 0.33960
## wgs_species_2Other -1.70442 0.18188 0.75629
## robust se z Pr(>|z|)
## wgs_species_2Staphylococcus haemolyticus 0.54963 -4.119 3.8e-05 ***
## wgs_species_2Staphylococcus sciuri 0.43882 -1.171 0.2418
## wgs_species_2Staphylococcus succinus 0.58032 -0.375 0.7078
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.41398 -0.155 0.8769
## wgs_species_2Other 0.88350 -1.929 0.0537 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef)
## wgs_species_2Staphylococcus haemolyticus 0.1039 9.622
## wgs_species_2Staphylococcus sciuri 0.5983 1.671
## wgs_species_2Staphylococcus succinus 0.8045 1.243
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.9379 1.066
## wgs_species_2Other 0.1819 5.498
## lower .95 upper .95
## wgs_species_2Staphylococcus haemolyticus 0.03539 0.3052
## wgs_species_2Staphylococcus sciuri 0.25316 1.4140
## wgs_species_2Staphylococcus succinus 0.25797 2.5090
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.41664 2.1112
## wgs_species_2Other 0.03219 1.0276
##
## Concordance= 0.683 (se = 0.041 )
## Likelihood ratio test= 39.48 on 5 df, p=2e-07
## Wald test = 27.76 on 5 df, p=4e-05
## Score (logrank) test = 35.34 on 5 df, p=1e-06, Robust = 26.85 p=6e-05
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).# Type III p-value
Anova(model,type="III")## Warning in Anova.coxph(model, type = "III"): LR tests unavailable with robust variances
## Wald tests substituted
## Warning in Anova.coxph(model, type = "III"): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III tests)
##
## Response: survival
## Df Chisq Pr(>Chisq)
## wgs_species_2 5 27.758 4.059e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Pairwise comparisons
confint(emmeans(model,pairwise~wgs_species_2,type="response"))$contrasts## contrast ratio
## Staphylococcus chromogenes / Staphylococcus haemolyticus 9.622
## Staphylococcus chromogenes / Staphylococcus sciuri 1.671
## Staphylococcus chromogenes / Staphylococcus succinus 1.243
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus) 1.066
## Staphylococcus chromogenes / Other 5.498
## Staphylococcus haemolyticus / Staphylococcus sciuri 0.174
## Staphylococcus haemolyticus / Staphylococcus succinus 0.129
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus) 0.111
## Staphylococcus haemolyticus / Other 0.571
## Staphylococcus sciuri / Staphylococcus succinus 0.744
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus) 0.638
## Staphylococcus sciuri / Other 3.290
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus) 0.858
## Staphylococcus succinus / Other 4.423
## (Staphylococcus xylosus/pseudoxylosus) / Other 5.157
## SE df asymp.LCL asymp.UCL
## 5.2885 Inf 2.0093 46.078
## 0.7335 Inf 0.4786 5.837
## 0.7213 Inf 0.2378 6.496
## 0.4414 Inf 0.3277 3.469
## 4.8577 Inf 0.4434 68.180
## 0.0688 Inf 0.0562 0.537
## 0.0654 Inf 0.0305 0.547
## 0.0558 Inf 0.0264 0.466
## 0.4729 Inf 0.0540 6.041
## 0.3000 Inf 0.2356 2.348
## 0.2437 Inf 0.2147 1.895
## 2.7198 Inf 0.3118 34.705
## 0.4354 Inf 0.2020 3.644
## 3.8095 Inf 0.3801 51.477
## 3.8947 Inf 0.5993 44.373
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 6 estimates
## Intervals are back-transformed from the log scaleemmeans(model,pairwise~wgs_species_2,type="response")$contrasts## contrast ratio
## Staphylococcus chromogenes / Staphylococcus haemolyticus 9.622
## Staphylococcus chromogenes / Staphylococcus sciuri 1.671
## Staphylococcus chromogenes / Staphylococcus succinus 1.243
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus) 1.066
## Staphylococcus chromogenes / Other 5.498
## Staphylococcus haemolyticus / Staphylococcus sciuri 0.174
## Staphylococcus haemolyticus / Staphylococcus succinus 0.129
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus) 0.111
## Staphylococcus haemolyticus / Other 0.571
## Staphylococcus sciuri / Staphylococcus succinus 0.744
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus) 0.638
## Staphylococcus sciuri / Other 3.290
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus) 0.858
## Staphylococcus succinus / Other 4.423
## (Staphylococcus xylosus/pseudoxylosus) / Other 5.157
## SE df null z.ratio p.value
## 5.2885 Inf 1 4.119 0.0005
## 0.7335 Inf 1 1.171 0.8510
## 0.7213 Inf 1 0.375 0.9990
## 0.4414 Inf 1 0.155 1.0000
## 4.8577 Inf 1 1.929 0.3841
## 0.0688 Inf 1 -4.418 0.0001
## 0.0654 Inf 1 -4.043 0.0008
## 0.0558 Inf 1 -4.365 0.0002
## 0.4729 Inf 1 -0.676 0.9846
## 0.3000 Inf 1 -0.734 0.9777
## 0.2437 Inf 1 -1.177 0.8482
## 2.7198 Inf 1 1.440 0.7023
## 0.4354 Inf 1 -0.302 0.9997
## 3.8095 Inf 1 1.727 0.5141
## 3.8947 Inf 1 2.172 0.2510
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log scaleconfint(emmeans(model,revpairwise~wgs_species_2,type="response"))$contrasts## contrast ratio
## Staphylococcus haemolyticus / Staphylococcus chromogenes 0.104
## Staphylococcus sciuri / Staphylococcus chromogenes 0.598
## Staphylococcus sciuri / Staphylococcus haemolyticus 5.757
## Staphylococcus succinus / Staphylococcus chromogenes 0.805
## Staphylococcus succinus / Staphylococcus haemolyticus 7.741
## Staphylococcus succinus / Staphylococcus sciuri 1.345
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus chromogenes 0.938
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus haemolyticus 9.024
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus sciuri 1.568
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus succinus 1.166
## Other / Staphylococcus chromogenes 0.182
## Other / Staphylococcus haemolyticus 1.750
## Other / Staphylococcus sciuri 0.304
## Other / Staphylococcus succinus 0.226
## Other / (Staphylococcus xylosus/pseudoxylosus) 0.194
## SE df asymp.LCL asymp.UCL
## 0.0571 Inf 0.0217 0.498
## 0.2625 Inf 0.1713 2.089
## 2.2807 Inf 1.8615 17.803
## 0.4669 Inf 0.1539 4.205
## 3.9188 Inf 1.8292 32.759
## 0.5424 Inf 0.4260 4.245
## 0.3883 Inf 0.2883 3.051
## 4.5477 Inf 2.1464 37.940
## 0.5989 Inf 0.5277 4.657
## 0.5917 Inf 0.2745 4.952
## 0.1607 Inf 0.0147 2.255
## 1.4482 Inf 0.1655 18.502
## 0.2513 Inf 0.0288 3.207
## 0.1947 Inf 0.0194 2.631
## 0.1465 Inf 0.0225 1.669
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 6 estimates
## Intervals are back-transformed from the log scaleemmeans(model,revpairwise~wgs_species_2,type="response")$contrasts## contrast ratio
## Staphylococcus haemolyticus / Staphylococcus chromogenes 0.104
## Staphylococcus sciuri / Staphylococcus chromogenes 0.598
## Staphylococcus sciuri / Staphylococcus haemolyticus 5.757
## Staphylococcus succinus / Staphylococcus chromogenes 0.805
## Staphylococcus succinus / Staphylococcus haemolyticus 7.741
## Staphylococcus succinus / Staphylococcus sciuri 1.345
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus chromogenes 0.938
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus haemolyticus 9.024
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus sciuri 1.568
## (Staphylococcus xylosus/pseudoxylosus) / Staphylococcus succinus 1.166
## Other / Staphylococcus chromogenes 0.182
## Other / Staphylococcus haemolyticus 1.750
## Other / Staphylococcus sciuri 0.304
## Other / Staphylococcus succinus 0.226
## Other / (Staphylococcus xylosus/pseudoxylosus) 0.194
## SE df null z.ratio p.value
## 0.0571 Inf 1 -4.119 0.0005
## 0.2625 Inf 1 -1.171 0.8510
## 2.2807 Inf 1 4.418 0.0001
## 0.4669 Inf 1 -0.375 0.9990
## 3.9188 Inf 1 4.043 0.0008
## 0.5424 Inf 1 0.734 0.9777
## 0.3883 Inf 1 -0.155 1.0000
## 4.5477 Inf 1 4.365 0.0002
## 0.5989 Inf 1 1.177 0.8482
## 0.5917 Inf 1 0.302 0.9997
## 0.1607 Inf 1 -1.929 0.3841
## 1.4482 Inf 1 0.676 0.9846
## 0.2513 Inf 1 -1.440 0.7023
## 0.1947 Inf 1 -1.727 0.5141
## 0.1465 Inf 1 -2.172 0.2510
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log scale# Proportional hazards assumption
cox.zph(model)## chisq df p
## wgs_species_2 7.88 5 0.16
## GLOBAL 7.88 5 0.16# High ama different species
table(df_sau$high_ama,df_sau$wgs_species_2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 17 18
## 1 2 1
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 5 6
## 1 1 5
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 17 4
## 1 1 0prop.table(table(df_sau$high_ama,df_sau$wgs_species_2),2)##
## Staphylococcus chromogenes Staphylococcus haemolyticus
## 0 0.89473684 0.94736842
## 1 0.10526316 0.05263158
##
## Staphylococcus sciuri Staphylococcus succinus
## 0 0.83333333 0.54545455
## 1 0.16666667 0.45454545
##
## Staphylococcus xylosus/pseudoxylosus Other
## 0 0.94444444 1.00000000
## 1 0.05555556 0.00000000model <- glmer(high_ama~ wgs_species_2 + (1|cow_id),family=binomial(link="logit"), data=df_sau)## boundary (singular) fit: see help('isSingular')summary(model)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: high_ama ~ wgs_species_2 + (1 | cow_id)
## Data: df_sau
##
## AIC BIC logLik deviance df.resid
## 62.9 79.3 -24.5 48.9 70
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9129 -0.3430 -0.2425 -0.2357 4.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 5.263e-15 7.255e-08
## Number of obs: 77, groups: cow_id, 44
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.140e+00 7.475e-01 -2.863
## wgs_species_2Staphylococcus haemolyticus -7.503e-01 1.271e+00 -0.591
## wgs_species_2Staphylococcus sciuri 5.306e-01 1.326e+00 0.400
## wgs_species_2Staphylococcus succinus 1.958e+00 9.620e-01 2.035
## wgs_species_2Staphylococcus xylosus/pseudoxylosus -6.931e-01 1.272e+00 -0.545
## wgs_species_2Other -3.728e+01 3.355e+07 0.000
## Pr(>|z|)
## (Intercept) 0.0042 **
## wgs_species_2Staphylococcus haemolyticus 0.5548
## wgs_species_2Staphylococcus sciuri 0.6891
## wgs_species_2Staphylococcus succinus 0.0418 *
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.5858
## wgs_species_2Other 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) w__2Sh wgs_spcs_2Stphylcccsscr
## wgs_spc_2Sh -0.588
## wgs_spcs_2Stphylcccsscr -0.564 0.332
## wgs_spcs_2Stphylcccsscc -0.777 0.457 0.438
## wgs_sp_2Sx/ -0.588 0.346 0.331
## wgs_spcs_2O 0.000 0.000 0.000
## wgs_spcs_2Stphylcccsscc w__2Sx
## wgs_spc_2Sh
## wgs_spcs_2Stphylcccsscr
## wgs_spcs_2Stphylcccsscc
## wgs_sp_2Sx/ 0.457
## wgs_spcs_2O 0.000 0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')exp(Confint(model))## Warning in vcov.merMod(object, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Estimate 2.5 %
## (Intercept) 1.176471e-01 0.02718126
## wgs_species_2Staphylococcus haemolyticus 4.722222e-01 0.03914037
## wgs_species_2Staphylococcus sciuri 1.700000e+00 0.12635137
## wgs_species_2Staphylococcus succinus 7.083334e+00 1.07488349
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 5.000001e-01 0.04133848
## wgs_species_2Other 6.471762e-17 0.00000000
## 97.5 %
## (Intercept) 0.5092048
## wgs_species_2Staphylococcus haemolyticus 5.6972851
## wgs_species_2Staphylococcus sciuri 22.8727249
## wgs_species_2Staphylococcus succinus 46.6781941
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 6.0476361
## wgs_species_2Other InfAnova(model,type="III")## Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: high_ama
## Chisq Df Pr(>Chisq)
## (Intercept) 8.1956 1 0.004199 **
## wgs_species_2 8.9819 5 0.109790
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model,pairwise~wgs_species_2,type="response")## Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## $emmeans
## wgs_species_2 prob SE df asymp.LCL asymp.UCL
## Staphylococcus chromogenes 0.1053 7.04e-02 Inf 0.02646 0.337
## Staphylococcus haemolyticus 0.0526 5.12e-02 Inf 0.00736 0.294
## Staphylococcus sciuri 0.1667 1.52e-01 Inf 0.02283 0.631
## Staphylococcus succinus 0.4545 1.50e-01 Inf 0.20276 0.732
## Staphylococcus xylosus/pseudoxylosus 0.0556 5.40e-02 Inf 0.00777 0.307
## Other 0.0000 1.00e-08 Inf 0.00000 1.000
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast
## Staphylococcus chromogenes / Staphylococcus haemolyticus
## Staphylococcus chromogenes / Staphylococcus sciuri
## Staphylococcus chromogenes / Staphylococcus succinus
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus chromogenes / Other
## Staphylococcus haemolyticus / Staphylococcus sciuri
## Staphylococcus haemolyticus / Staphylococcus succinus
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus haemolyticus / Other
## Staphylococcus sciuri / Staphylococcus succinus
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus sciuri / Other
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus succinus / Other
## (Staphylococcus xylosus/pseudoxylosus) / Other
## odds.ratio SE df null z.ratio p.value
## 2.00e+00 3.00e+00 Inf 1 0.591 0.9917
## 1.00e+00 1.00e+00 Inf 1 -0.400 0.9987
## 0.00e+00 0.00e+00 Inf 1 -2.035 0.3223
## 2.00e+00 3.00e+00 Inf 1 0.545 0.9943
## 1.55e+16 5.18e+23 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 -0.853 0.9574
## 0.00e+00 0.00e+00 Inf 1 -2.271 0.2061
## 1.00e+00 1.00e+00 Inf 1 -0.039 1.0000
## 7.30e+15 2.45e+23 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 -1.140 0.8646
## 3.00e+00 5.00e+00 Inf 1 0.814 0.9650
## 2.63e+16 8.81e+23 Inf 1 0.000 1.0000
## 1.40e+01 1.70e+01 Inf 1 2.220 0.2283
## 1.09e+17 3.67e+24 Inf 1 0.000 1.0000
## 7.73e+15 2.59e+23 Inf 1 0.000 1.0000
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log odds ratio scaleperformance::icc(model,by_group=TRUE)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## [1] NA# High ama different species
table(df_sub$high_ama,df_sub$wgs_species)##
## Staphylococcus chromogenes Staphylococcus cohnii Staphylococcus devriesei
## 0 17 1 2
## 1 1 0 0
##
## Staphylococcus equorum Staphylococcus haemolyticus
## 0 1 18
## 1 0 1
##
## Staphylococcus pseudoxylosus Staphylococcus sciuri Staphylococcus succinus
## 0 5 6 5
## 1 0 0 5
##
## Staphylococcus xylosus Unconclusive for assigning CFUs
## 0 11 1
## 1 2 1prop.table(table(df_sub$high_ama,df_sub$wgs_species),2)##
## Staphylococcus chromogenes Staphylococcus cohnii Staphylococcus devriesei
## 0 0.94444444 1.00000000 1.00000000
## 1 0.05555556 0.00000000 0.00000000
##
## Staphylococcus equorum Staphylococcus haemolyticus
## 0 1.00000000 0.94736842
## 1 0.00000000 0.05263158
##
## Staphylococcus pseudoxylosus Staphylococcus sciuri Staphylococcus succinus
## 0 1.00000000 1.00000000 0.50000000
## 1 0.00000000 0.00000000 0.50000000
##
## Staphylococcus xylosus Unconclusive for assigning CFUs
## 0 0.84615385 0.50000000
## 1 0.15384615 0.50000000model <- glmer(high_ama~ wgs_species_2 + (1|cow_id),family=binomial(link="logit"), data=df_sub)## boundary (singular) fit: see help('isSingular')summary(model)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: high_ama ~ wgs_species_2 + (1 | cow_id)
## Data: df_sub
##
## AIC BIC logLik deviance df.resid
## 62.3 78.7 -24.2 48.3 70
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9129 -0.3536 -0.3430 0.0000 4.2426
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow_id (Intercept) 0 0
## Number of obs: 77, groups: cow_id, 44
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.140e+00 7.475e-01 -2.863
## wgs_species_2Staphylococcus haemolyticus -7.503e-01 1.271e+00 -0.591
## wgs_species_2Staphylococcus sciuri -4.767e+01 2.740e+07 0.000
## wgs_species_2Staphylococcus succinus 1.958e+00 9.620e-01 2.035
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 6.062e-02 1.059e+00 0.057
## wgs_species_2Other -5.276e+01 3.355e+07 0.000
## Pr(>|z|)
## (Intercept) 0.0042 **
## wgs_species_2Staphylococcus haemolyticus 0.5548
## wgs_species_2Staphylococcus sciuri 1.0000
## wgs_species_2Staphylococcus succinus 0.0418 *
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 0.9543
## wgs_species_2Other 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) w__2Sh wgs_spcs_2Stphylcccsscr
## wgs_spc_2Sh -0.588
## wgs_spcs_2Stphylcccsscr 0.000 0.000
## wgs_spcs_2Stphylcccsscc -0.777 0.457 0.000
## wgs_sp_2Sx/ -0.706 0.415 0.000
## wgs_spcs_2O 0.000 0.000 0.000
## wgs_spcs_2Stphylcccsscc w__2Sx
## wgs_spc_2Sh
## wgs_spcs_2Stphylcccsscr
## wgs_spcs_2Stphylcccsscc
## wgs_sp_2Sx/ 0.549
## wgs_spcs_2O 0.000 0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')exp(Confint(model))## Warning in vcov.merMod(object, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Estimate 2.5 %
## (Intercept) 1.176471e-01 0.02718126
## wgs_species_2Staphylococcus haemolyticus 4.722222e-01 0.03914037
## wgs_species_2Staphylococcus sciuri 1.974037e-21 0.00000000
## wgs_species_2Staphylococcus succinus 7.083333e+00 1.07488334
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 1.062500e+00 0.13334287
## wgs_species_2Other 1.222760e-23 0.00000000
## 97.5 %
## (Intercept) 0.5092049
## wgs_species_2Staphylococcus haemolyticus 5.6972846
## wgs_species_2Staphylococcus sciuri Inf
## wgs_species_2Staphylococcus succinus 46.6781855
## wgs_species_2Staphylococcus xylosus/pseudoxylosus 8.4661913
## wgs_species_2Other InfAnova(model,type="III")## Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
## arithmetic operators in their names;
## the printed representation of the hypothesis will be omitted## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: high_ama
## Chisq Df Pr(>Chisq)
## (Intercept) 8.1956 1 0.004199 **
## wgs_species_2 7.8953 5 0.162101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(model,pairwise~wgs_species_2,type="response")## Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## $emmeans
## wgs_species_2 prob SE df asymp.LCL asymp.UCL
## Staphylococcus chromogenes 0.1053 7.04e-02 Inf 0.02646 0.337
## Staphylococcus haemolyticus 0.0526 5.12e-02 Inf 0.00736 0.294
## Staphylococcus sciuri 0.0000 1.00e-08 Inf 0.00000 1.000
## Staphylococcus succinus 0.4545 1.50e-01 Inf 0.20276 0.732
## Staphylococcus xylosus/pseudoxylosus 0.1111 7.41e-02 Inf 0.02794 0.352
## Other 0.0000 1.00e-08 Inf 0.00000 1.000
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast
## Staphylococcus chromogenes / Staphylococcus haemolyticus
## Staphylococcus chromogenes / Staphylococcus sciuri
## Staphylococcus chromogenes / Staphylococcus succinus
## Staphylococcus chromogenes / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus chromogenes / Other
## Staphylococcus haemolyticus / Staphylococcus sciuri
## Staphylococcus haemolyticus / Staphylococcus succinus
## Staphylococcus haemolyticus / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus haemolyticus / Other
## Staphylococcus sciuri / Staphylococcus succinus
## Staphylococcus sciuri / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus sciuri / Other
## Staphylococcus succinus / (Staphylococcus xylosus/pseudoxylosus)
## Staphylococcus succinus / Other
## (Staphylococcus xylosus/pseudoxylosus) / Other
## odds.ratio SE df null z.ratio p.value
## 2.00e+00 3.00e+00 Inf 1 0.591 0.9917
## 5.07e+20 1.39e+28 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 -2.035 0.3223
## 1.00e+00 1.00e+00 Inf 1 -0.057 1.0000
## 8.18e+22 2.74e+30 Inf 1 0.000 1.0000
## 2.39e+20 6.55e+27 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 -2.271 0.2061
## 0.00e+00 1.00e+00 Inf 1 -0.638 0.9882
## 3.86e+22 1.30e+30 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 0.000 1.0000
## 0.00e+00 0.00e+00 Inf 1 0.000 1.0000
## 1.61e+02 6.99e+09 Inf 1 0.000 1.0000
## 7.00e+00 6.00e+00 Inf 1 1.968 0.3608
## 5.79e+23 1.94e+31 Inf 1 0.000 1.0000
## 8.69e+22 2.92e+30 Inf 1 0.000 1.0000
##
## P value adjustment: tukey method for comparing a family of 6 estimates
## Tests are performed on the log odds ratio scaleperformance::icc(model,by_group = TRUE)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## [1] NAkm_plots_2 <- ggarrange(km_3_plot, km_4_plot, ncol = 2, nrow = 1, common.legend = TRUE,align = "h",labels = c("C", "D"))
km_plots_2
ggsave(plot = last_plot(),"./figures/km_2_annotated.png",width = 20, height = 10, units = "cm")
km_plots_all <- ggarrange(km_plots_1, km_plots_2,
ncol = 1, nrow = 2, align = "v")
km_plots_all
ggsave(plot = last_plot(),"./figures/figure_2.png",width = 20, height = 20, units = "cm",dpi=600)
ggsave(plot = last_plot(),"./figures/figure_2.pdf",width = 20, height = 20, units = "cm",dpi=600)
ggsave(plot = last_plot(),"./figures/figure_2.jpeg",width = 20, height = 20, units = "cm",dpi=600) 4 High in vitro inhibitory activity
4.1 High in vitro inhibitory activity and presence of SAU-SSLO-IMI
# --- Presence of high in-vitro antimicrobial activity and SAU-SSLO-IMI ----
# If we sort isolates by antimicrobial activity against staph aureus, the top 10 isolates are mostly found in controls
table(df_sau$case_control ,df_sau$inhibition_category)##
## Top 10 11-20 21-30 31-40 41-50 51-60 61-70 >70 No inhibition
## Control 8 5 5 3 4 7 5 0 3
## Case 2 5 5 7 6 3 5 2 2table1(~nas_count_hundred+log_nas_media|inhibition_category,data=df_sau)| Top 10 (N=10) |
11-20 (N=10) |
21-30 (N=10) |
31-40 (N=10) |
41-50 (N=10) |
51-60 (N=10) |
61-70 (N=10) |
>70 (N=2) |
No inhibition (N=5) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|---|---|---|
| nas_count_hundred | ||||||||||
| Mean (SD) | 297 (142) | 578 (66.1) | 1070 (193) | 1810 (485) | 3920 (466) | 10300 (2460) | 23200 (8680) | 42400 (5320) | 13100 (9090) | 7290 (10300) |
| Median [Min, Max] | 355 [83.2, 455] | 564 [499, 703] | 1040 [802, 1400] | 1600 [1450, 2810] | 3940 [3230, 4750] | 10900 [5640, 12800] | 22500 [13700, 36600] | 42400 [38600, 46100] | 11100 [1750, 23900] | 1920 [83.2, 46100] |
| log_nas_media | ||||||||||
| Mean (SD) | 2.41 (0.263) | 2.76 (0.0489) | 3.02 (0.0786) | 3.25 (0.104) | 3.59 (0.0517) | 4.00 (0.121) | 4.34 (0.166) | 4.63 (0.0547) | 3.98 (0.453) | 3.41 (0.678) |
| Median [Min, Max] | 2.55 [1.92, 2.66] | 2.75 [2.70, 2.85] | 3.02 [2.90, 3.14] | 3.20 [3.16, 3.45] | 3.59 [3.51, 3.68] | 4.04 [3.75, 4.11] | 4.34 [4.14, 4.56] | 4.63 [4.59, 4.66] | 4.04 [3.24, 4.38] | 3.28 [1.92, 4.66] |
table1(~nas_count_hundred+log_nas_media|inhibition_order,data=subset(df_sau,df_sau$inhibition_category=="Top 10"))## Warning in table1.formula(~nas_count_hundred + log_nas_media |
## inhibition_order, : Terms to the right of '|' in formula 'x' define table
## columns and are expected to be factors with meaningful labels.| 1 (N=1) |
2 (N=1) |
3 (N=1) |
4 (N=1) |
5 (N=1) |
6 (N=1) |
7 (N=1) |
8 (N=1) |
9 (N=1) |
10 (N=1) |
Overall (N=10) |
|
|---|---|---|---|---|---|---|---|---|---|---|---|
| nas_count_hundred | |||||||||||
| Mean (SD) | 83.2 (NA) | 140 (NA) | 146 (NA) | 191 (NA) | 325 (NA) | 386 (NA) | 396 (NA) | 406 (NA) | 446 (NA) | 455 (NA) | 297 (142) |
| Median [Min, Max] | 83.2 [83.2, 83.2] | 140 [140, 140] | 146 [146, 146] | 191 [191, 191] | 325 [325, 325] | 386 [386, 386] | 396 [396, 396] | 406 [406, 406] | 446 [446, 446] | 455 [455, 455] | 355 [83.2, 455] |
| log_nas_media | |||||||||||
| Mean (SD) | 1.92 (NA) | 2.14 (NA) | 2.16 (NA) | 2.28 (NA) | 2.51 (NA) | 2.59 (NA) | 2.60 (NA) | 2.61 (NA) | 2.65 (NA) | 2.66 (NA) | 2.41 (0.263) |
| Median [Min, Max] | 1.92 [1.92, 1.92] | 2.14 [2.14, 2.14] | 2.16 [2.16, 2.16] | 2.28 [2.28, 2.28] | 2.51 [2.51, 2.51] | 2.59 [2.59, 2.59] | 2.60 [2.60, 2.60] | 2.61 [2.61, 2.61] | 2.65 [2.65, 2.65] | 2.66 [2.66, 2.66] | 2.55 [1.92, 2.66] |
data_plot <- as.data.frame(table(df_sau$imi ,df_sau$inhibition_category))
ranking_1 <- data_plot %>%
ggplot(aes(x=Var1, y =Freq,fill=Var2)) +
geom_bar(position="fill", stat="identity") +
theme_classic() +
labs(x="IMI",fill="Isolates categorization",title="Staphylococcus aureus") +
scale_x_discrete(limits = rev) +
scale_fill_viridis(discrete=TRUE) + coord_flip() +
theme(
plot.title = element_text(face = "italic"),
axis.title.x = element_blank() # Removes the y-axis label
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1))
# If we sort isolates by antimicrobial activity against strep uberis, there is no clear reltionship
table(df_sub$case_control ,df_sub$inhibition_category)##
## Top 10 11-20 21-30 31-40 41-50 51-60 61-70 >70 No inhibition
## Control 4 7 5 3 4 6 2 1 8
## Case 6 3 5 6 5 2 3 0 7table1(~nas_count_hundred+log_nas_media|inhibition_category,data=df_sub)| Top 10 (N=10) |
11-20 (N=10) |
21-30 (N=10) |
31-40 (N=9) |
41-50 (N=9) |
51-60 (N=8) |
61-70 (N=5) |
>70 (N=1) |
No inhibition (N=15) |
Overall (N=77) |
|
|---|---|---|---|---|---|---|---|---|---|---|
| nas_count_hundred | ||||||||||
| Mean (SD) | 182 (118) | 427 (60.1) | 838 (250) | 1490 (169) | 3210 (883) | 7900 (2150) | 13700 (2860) | 150000 (NA) | 38000 (41100) | 11800 (27800) |
| Median [Min, Max] | 147 [44.6, 364] | 406 [366, 554] | 822 [554, 1210] | 1480 [1240, 1760] | 3150 [1920, 4750] | 7720 [4990, 10600] | 12800 [11000, 18400] | 150000 [150000, 150000] | 23900 [1750, 137000] | 1750 [44.6, 150000] |
| log_nas_media | ||||||||||
| Mean (SD) | 2.16 (0.324) | 2.63 (0.0575) | 2.91 (0.132) | 3.17 (0.0494) | 3.49 (0.119) | 3.88 (0.121) | 4.13 (0.0851) | 5.18 (NA) | 4.33 (0.519) | 3.36 (0.812) |
| Median [Min, Max] | 2.17 [1.65, 2.56] | 2.61 [2.56, 2.74] | 2.91 [2.74, 3.08] | 3.17 [3.09, 3.25] | 3.50 [3.28, 3.68] | 3.89 [3.70, 4.03] | 4.11 [4.04, 4.27] | 5.18 [5.18, 5.18] | 4.38 [3.24, 5.14] | 3.24 [1.65, 5.18] |
table1(~nas_count_hundred+log_nas_media|inhibition_order,data=subset(df_sub,df_sub$inhibition_category=="Top 10"))## Warning in table1.formula(~nas_count_hundred + log_nas_media |
## inhibition_order, : Terms to the right of '|' in formula 'x' define table
## columns and are expected to be factors with meaningful labels.| 1 (N=1) |
2 (N=1) |
3 (N=1) |
4 (N=1) |
5 (N=1) |
6 (N=1) |
7 (N=1) |
8 (N=1) |
9 (N=1) |
10 (N=1) |
Overall (N=10) |
|
|---|---|---|---|---|---|---|---|---|---|---|---|
| nas_count_hundred | |||||||||||
| Mean (SD) | 44.6 (NA) | 51.5 (NA) | 83.2 (NA) | 140 (NA) | 146 (NA) | 149 (NA) | 191 (NA) | 325 (NA) | 327 (NA) | 364 (NA) | 182 (118) |
| Median [Min, Max] | 44.6 [44.6, 44.6] | 51.5 [51.5, 51.5] | 83.2 [83.2, 83.2] | 140 [140, 140] | 146 [146, 146] | 149 [149, 149] | 191 [191, 191] | 325 [325, 325] | 327 [327, 327] | 364 [364, 364] | 147 [44.6, 364] |
| log_nas_media | |||||||||||
| Mean (SD) | 1.65 (NA) | 1.71 (NA) | 1.92 (NA) | 2.14 (NA) | 2.16 (NA) | 2.17 (NA) | 2.28 (NA) | 2.51 (NA) | 2.51 (NA) | 2.56 (NA) | 2.16 (0.324) |
| Median [Min, Max] | 1.65 [1.65, 1.65] | 1.71 [1.71, 1.71] | 1.92 [1.92, 1.92] | 2.14 [2.14, 2.14] | 2.16 [2.16, 2.16] | 2.17 [2.17, 2.17] | 2.28 [2.28, 2.28] | 2.51 [2.51, 2.51] | 2.51 [2.51, 2.51] | 2.56 [2.56, 2.56] | 2.17 [1.65, 2.56] |
data_plot <- as.data.frame(table(df_sub$imi ,df_sub$inhibition_category))
ranking_2 <- data_plot %>%
ggplot(aes(x=Var1, y =Freq,fill=Var2)) +
geom_bar(position="fill", stat="identity") +
theme_classic() +
labs(x="IMI",fill="Isolates categorization",title="Streptococcus uberis") +
scale_x_discrete(limits = rev) +
scale_fill_viridis(discrete=TRUE) + coord_flip() +
theme(
plot.title = element_text(face = "italic"),
axis.title.x = element_blank() # Removes the y-axis label
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1))
# Run logisic regression model to compare the likelihood of exposure between cases and controls
table(df_sau_per_cow_2$imi,df_sau_per_cow_2$high_ama_cow)##
## 0 1
## No 29 7
## Yes 33 2model <- glm(high_ama_cow~imi+farm,family=binomial(link="logit"),data=df_sau_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.2320254 0.08545841 0.5375776
## imiYes 0.2505766 0.03539033 1.1347412
## farmB 1.2140638 0.16208626 6.1506383# Probability of being a case in cows with isolates with high antimicrobial activity (1) and those without high antimicrobial activity (0)
emmeans(model,~imi,type="response")## imi prob SE df asymp.LCL asymp.UCL
## No 0.2036 0.0798 Inf 0.0888 0.401
## Yes 0.0602 0.0434 Inf 0.0141 0.224
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Type III p-values
Anova(model,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: high_ama_cow
## LR Chisq Df Pr(>Chisq)
## imi 3.1932 1 0.07394 .
## farm 0.0469 1 0.82852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Run logistic regression model to compare the likelihood of exposure between cases and controls
table(df_sub_per_cow_2$imi,df_sub_per_cow_2$high_ama_cow)##
## 0 1
## No 32 4
## Yes 30 5model <- glm(high_ama_cow~imi+farm,family=binomial(link="logit"),data=df_sub_per_cow_2)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.1206605 0.0337768 0.3228133
## imiYes 1.3322226 0.3230293 5.8307269
## farmB 1.1879755 0.1629468 5.7186982# Probability of being a case in cows with isolates with high antimicrobial activity (1) and those without high antimicrobial activity (0)
emmeans(model,~imi,type="response")## imi prob SE df asymp.LCL asymp.UCL
## No 0.116 0.0600 Inf 0.0401 0.293
## Yes 0.149 0.0686 Inf 0.0573 0.336
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Type III p-values
Anova(model,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: high_ama_cow
## LR Chisq Df Pr(>Chisq)
## imi 0.160750 1 0.6885
## farm 0.038744 1 0.8440# ---- MIC below referent bacteria -----
table(df_sau_per_cow_2$high_ama_atcc_cow,df_sau_per_cow_2$case)##
## 0 1
## 0 24 28
## 1 12 7# Run logisic regression model to compare the likelihood of exposure between cases and controls
model <- glm(high_ama_atcc~imi+farm,family=binomial(link="logit"),data=df_sau)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.4050444 0.1821173 0.8401318
## imiYes 0.6631045 0.2373941 1.8081527
## farmB 4.8071177 1.6790957 14.6377125# Probability of being a case in cows with isolates with high antimicrobial activity (1) and those without high antimicrobial activity (0)
emmeans(model,~imi,type="response")## imi prob SE df asymp.LCL asymp.UCL
## No 0.470 0.0893 Inf 0.306 0.642
## Yes 0.371 0.0918 Inf 0.214 0.560
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Type III p-values
Anova(model,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: high_ama_atcc
## LR Chisq Df Pr(>Chisq)
## imi 0.6443 1 0.422170
## farm 8.6267 1 0.003313 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Run logistic regression model to compare the likelihood of exposure between cases and controls
model <- glm(high_ama_atcc~imi+farm,family=binomial(link="logit"),data=df_sub)
exp(Confint(model))## Estimate 2.5 % 97.5 %
## (Intercept) 0.2700152 0.1088310 0.5968939
## imiYes 0.8143358 0.2599106 2.4715003
## farmB 1.2604100 0.3515613 4.0784468# Probability of being a case in cows with isolates with high antimicrobial activity (1) and those without high antimicrobial activity (0)
emmeans(model,~imi,type="response")## imi prob SE df asymp.LCL asymp.UCL
## No 0.233 0.0703 Inf 0.1229 0.396
## Yes 0.198 0.0704 Inf 0.0938 0.370
##
## Results are averaged over the levels of: farm
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale# Type III p-values
Anova(model,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: high_ama_atcc
## LR Chisq Df Pr(>Chisq)
## imi 0.13181 1 0.7166
## farm 0.13947 1 0.70884.2 High in vitro inhibitory activity across different NASM species
# --- Staph aureus
table(df_sau$wgs_species_3 ,df_sau$inhibition_category)##
## Top 10 11-20 21-30 31-40 41-50 51-60 61-70 >70 No inhibition
## SCH 2 2 2 6 4 1 2 0 0
## SHAEM 1 0 1 0 1 4 7 1 4
## MSC 1 1 0 1 2 1 0 0 0
## SSUC 5 0 3 0 2 1 0 0 0
## SXYL 1 7 4 3 1 2 0 0 0
## Other 0 0 0 0 0 1 1 1 1table(df_sau$wgs_species_3 ,df_sau$high_ama)##
## 0 1
## SCH 17 2
## SHAEM 18 1
## MSC 5 1
## SSUC 6 5
## SXYL 17 1
## Other 4 0fisher.test(table(df_sau$wgs_species_3 ,df_sau$high_ama))##
## Fisher's Exact Test for Count Data
##
## data: table(df_sau$wgs_species_3, df_sau$high_ama)
## p-value = 0.04789
## alternative hypothesis: two.sideddata_plot <- as.data.frame(table(df_sau$wgs_species_3 ,df_sau$inhibition_category))
ranking_3 <- data_plot %>%
ggplot(aes(x=Var1, y =Freq,fill=Var2)) +
geom_bar(position="fill", stat="identity") +
theme_classic() +
labs(x="NASM taxonomy",fill="Isolates categorization",title="Staphylococcus aureus") +
scale_x_discrete(limits = rev) +
scale_fill_viridis(discrete=TRUE) + coord_flip() +
theme(
plot.title = element_text(face = "italic"),
axis.title.x = element_blank() # Removes the y-axis label
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1))
# Strep uberis
table(df_sub$wgs_species_3 ,df_sub$inhibition_category)##
## Top 10 11-20 21-30 31-40 41-50 51-60 61-70 >70 No inhibition
## SCH 2 5 3 5 3 0 0 1 0
## SHAEM 1 0 0 1 1 1 2 0 13
## MSC 0 1 0 0 3 2 0 0 0
## SSUC 5 1 0 1 1 1 2 0 0
## SXYL 2 3 6 2 0 4 1 0 0
## Other 0 0 1 0 1 0 0 0 2table(df_sub$wgs_species_3 ,df_sub$high_ama)##
## 0 1
## SCH 17 2
## SHAEM 18 1
## MSC 6 0
## SSUC 6 5
## SXYL 16 2
## Other 4 0fisher.test(table(df_sub$wgs_species_3 ,df_sub$high_ama))##
## Fisher's Exact Test for Count Data
##
## data: table(df_sub$wgs_species_3, df_sub$high_ama)
## p-value = 0.0691
## alternative hypothesis: two.sideddata_plot <- as.data.frame(table(df_sub$wgs_species_3 ,df_sub$inhibition_category))
ranking_4 <- data_plot %>%
ggplot(aes(x=Var1, y =Freq,fill=Var2)) +
geom_bar(position="fill", stat="identity") +
theme_classic() +
labs(x="NASM taxonomy",fill="Isolates categorization",title="Streptococcus uberis") +
scale_x_discrete(limits = rev) +
scale_fill_viridis(discrete=TRUE) + coord_flip() +
theme(
plot.title = element_text(face = "italic"),
axis.title.x = element_blank() # Removes the y-axis label
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1))
ranking_figure <- ggarrange(ranking_1 , ranking_2,ranking_3, ranking_4, ncol = 2, nrow = 2, common.legend = TRUE,align = "h",labels = c("A", "B","C","D"))
#plot_1 <- annotate_figure(ranking_figure,
# bottom = text_grob("%", color = "black",
# hjust = 0.75, x = 0.5))
#plot_1
ranking_figure
ggsave(plot = last_plot(),"./figures/figure_1.png",width = 20, height = 20, units = "cm",dpi=600)
ggsave(plot = last_plot(),"./figures/figure_1.pdf",width = 20, height = 20, units = "cm",dpi=600)
ggsave(plot = last_plot(),"./figures/figure_1.jpeg",width = 20, height = 20, units = "cm",dpi=600)