qPCR Dataset Questions

Input Data

q1 <- read_csv("Q1_data_set_6-11-24.csv", 
    col_types = cols(treatment = col_factor(levels = c("3", 
        "4")), block = col_factor(levels = c("1", 
        "4", "6", "7", "8", "9", "10", "11", 
        "12")), round = col_factor(levels = c("1", 
        "2", "3")), inoc_binary = col_factor(levels = c("1", 
        "0"))))
q1$colony <- as.factor(q1$colony)

rounds <- read_csv("rounds.csv")

q3 <- read_csv("q3_non_inoc_bees.csv",
      col_types = cols(treatment = col_factor(levels = c("3", 
        "4")), block = col_factor(levels = c("1", 
        "4", "6", "7", "8", "9", "10", "11", 
        "12")), round = col_factor(levels = c("1", 
        "2", "3"))))

## Warning: The following named parsers don't match the column names: round

q3$colony <- as.factor(q3$colony)
q3$bee <- as.factor(q3$bee)

q3 <- merge(q3, rounds, by = "colony", all = FALSE)

q4 <- read_csv("q4_all_bees_per_bee.csv", 
    col_types = cols(treatment = col_factor(levels = c("3", 
        "4")), block = col_factor(levels = c("1", 
        "4", "6", "7", "8", "9", "10", "11", 
        "12")), inoc_binary = col_factor(levels = c("1", 
        "0")), round = col_factor(levels = c("1", 
        "2", "3"))))

q4$colony <- as.factor(q4$colony)

Input additional data for merge

all_bees <- read_csv("qpcr_inoc_bees.csv", col_types = cols(treatment = col_factor(levels = c("1","2", "3","4")),replicate = col_factor(levels = c("1","4", "6", "7", "8", "9", "10", "11", "12")), start = col_date(format = "%m/%d/%Y"), Innoculation_date = col_date(format = "%m/%d/%Y"), date = col_date(format = "%m/%d/%Y"), censor_status = col_factor(levels = c("1","2"))))

pollen <- read_csv("pollen.csv")
pollen$colony <- as.factor(pollen$colony)
pollen$treatment <- as.factor(pollen$treatment)
pollen$block <- as.factor(pollen$block)

avg.pol <- pollen %>%
  group_by(colony) %>%
  summarise(avg.pol = mean(whole_dif))


avg.pol <- as.data.frame(avg.pol)

all_bees <- merge(all_bees, avg.pol, by="colony", all = FALSE)
q1 <- merge(q1, avg.pol, by = "colony", all = FALSE)
q3 <- merge(q3, avg.pol, by = "colony", all = FALSE)
q4 <- merge(q4, avg.pol, by = "colony", all = FALSE)

How long since inoculation does it take for an infection to be detected in the inoculated bee?

days_detec <- glm.nb(day_first_detection ~ treatment + avg.pol + round, data = q1)

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

drop1(days_detec, test = "Chisq")

## Single term deletions
## 
## Model:
## day_first_detection ~ treatment + avg.pol + round
##           Df Deviance    AIC    LRT Pr(>Chi)
## <none>         2.6217 48.107                
## treatment  1   3.5071 46.992 0.8854   0.3467
## avg.pol    1   4.3335 47.819 1.7118   0.1908
## round      2   4.9317 46.417 2.3100   0.3151

ddet1 <- update(days_detec, .~. -round)

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

drop1(ddet1, test = "Chisq")

## Single term deletions
## 
## Model:
## day_first_detection ~ treatment + avg.pol
##           Df Deviance    AIC     LRT Pr(>Chi)
## <none>         4.9316 46.417                 
## treatment  1   5.8877 45.373 0.95607   0.3282
## avg.pol    1   5.4377 44.923 0.50612   0.4768

ddet2 <- update(ddet1, .~. -avg.pol)

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached

drop1(ddet2, test = "Chisq")

## Single term deletions
## 
## Model:
## day_first_detection ~ treatment
##           Df Deviance    AIC    LRT Pr(>Chi)
## <none>         5.4377 44.923                
## treatment  1   6.5334 44.019 1.0957   0.2952

Anova(ddet2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: day_first_detection
##           LR Chisq Df Pr(>Chisq)
## treatment   1.0957  1     0.2952

plot(q1$treatment, q1$day_first_detection)

days <- q1 %>%
  group_by(treatment) %>%
  summarise(m = mean(day_first_detection))

days

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          1.75
## 2 4          1.12

How long since inoculation does it take for an infection to be detected in non-inoculated bees?

days_detec.ni <- glmer.nb(days_first_det ~ treatment + avg.pol + (1|colony) + round, data = q3)

## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached

drop1(days_detec.ni, test = "Chisq")

## Single term deletions
## 
## Model:
## days_first_det ~ treatment + avg.pol + (1 | colony) + round
##           npar    AIC    LRT  Pr(Chi)   
## <none>         195.39                   
## treatment    1 194.25 0.8555 0.354993   
## avg.pol      1 194.09 0.6999 0.402828   
## round        1 201.21 7.8167 0.005177 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dd1 <- update(days_detec.ni, .~. -avg.pol)
drop1(dd1, test = "Chisq")

## Single term deletions
## 
## Model:
## days_first_det ~ treatment + (1 | colony) + round
##           npar    AIC    LRT  Pr(Chi)   
## <none>         194.09                   
## treatment    1 193.05 0.9535 0.328833   
## round        1 199.22 7.1251 0.007601 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q3$treatment, q3$days_first_det)

days.ni <- q3 %>%
  group_by(treatment) %>%
  summarise(m = mean(days_first_det))

days.ni

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          2.13
## 2 4          1.73

How long does it take for infection to reach ADL in inoculated bees?

days_detec_adl <- glm.nb(day_first_ADL ~ treatment + avg.pol + round, data = q1)
drop1(days_detec_adl, test = "Chisq")

## Single term deletions
## 
## Model:
## day_first_ADL ~ treatment + avg.pol + round
##           Df Deviance    AIC     LRT Pr(>Chi)
## <none>         11.664 65.652                 
## treatment  1   13.884 65.872 2.22045   0.1362
## avg.pol    1   11.765 63.754 0.10181   0.7497
## round      2   14.566 64.554 2.90186   0.2344

dadl <- update(days_detec_adl, .~. -avg.pol)
drop1(dadl, test = "Chisq")

## Single term deletions
## 
## Model:
## day_first_ADL ~ treatment + round
##           Df Deviance    AIC    LRT Pr(>Chi)
## <none>         11.775 63.754                
## treatment  1   13.996 63.975 2.2210   0.1361
## round      2   14.584 62.562 2.8084   0.2456

dadl2 <- update(dadl, .~. -round)
Anova(dadl2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: day_first_ADL
##           LR Chisq Df Pr(>Chisq)  
## treatment    2.993  1    0.08363 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q1$treatment, q1$day_first_ADL)

days.adl <- q1 %>%
  group_by(treatment) %>%
  summarise(m = mean(day_first_ADL))

days.adl

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3           3.5
## 2 4           2

How long does it take for infection to reach ADL in non-inoculated bees?

q3.adl <- q3[complete.cases(q3$days_first_adl), ]

days_detec_adl.ni <- glmer.nb(days_first_adl ~ treatment + avg.pol + (1|colony) + round, data = q3.adl)
drop1(days_detec_adl.ni, test = "Chisq")

## Single term deletions
## 
## Model:
## days_first_adl ~ treatment + avg.pol + (1 | colony) + round
##           npar    AIC     LRT   Pr(Chi)    
## <none>         221.34                      
## treatment    1 220.84  1.5064 0.2196816    
## avg.pol      1 224.33  4.9936 0.0254419 *  
## round        1 232.69 13.3492 0.0002585 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(days_detec_adl.ni)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: days_first_adl
##             Chisq Df Pr(>Chisq)    
## treatment  1.6428  1     0.1999    
## avg.pol    5.8933  1     0.0152 *  
## round     16.8066  1  4.139e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q3$treatment, q3$days_first_adl)

days.adl.ni <- q3.adl %>%
  group_by(treatment) %>%
  summarise(m = mean(days_first_adl))

days.adl.ni

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          4.91
## 2 4          6

How many days out of 14 did an inoc bee test above adl for infection?

total_adl <- glm.nb(num_days_adl ~ treatment + avg.pol + round, data = q1)
drop1(total_adl, test = "Chisq")

## Single term deletions
## 
## Model:
## num_days_adl ~ treatment + avg.pol + round
##           Df Deviance    AIC    LRT Pr(>Chi)
## <none>         17.667 90.571                
## treatment  1   17.981 88.885 0.3138   0.5753
## avg.pol    1   18.681 89.586 1.0145   0.3138
## round      2   21.194 90.099 3.5272   0.1714

tot_ad <- update(total_adl, .~. -avg.pol)
drop1(tot_ad, test = "Chisq")

## Single term deletions
## 
## Model:
## num_days_adl ~ treatment + round
##           Df Deviance    AIC    LRT Pr(>Chi)  
## <none>         17.505 89.549                  
## treatment  1   17.781 87.826 0.2763  0.59913  
## round      2   23.247 91.292 5.7430  0.05661 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

totad2 <- update(tot_ad, .~. -round)
Anova(totad2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: num_days_adl
##           LR Chisq Df Pr(>Chisq)
## treatment  0.24634  1     0.6197

plot(q1$treatment, q1$num_days_adl)

num_adl <- q1 %>%
  group_by(treatment) %>%
  summarise(m = mean(num_days_adl))
num_adl

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3           5.5
## 2 4           6.5

How many days out of 14 did a non-inoc bee test above adl for infection?

total_adl.ni <- glmer.nb(total_days_adl ~ treatment + avg.pol + (1|colony) + round, data = q3)
drop1(total_adl.ni, test = "Chisq")

## Single term deletions
## 
## Model:
## total_days_adl ~ treatment + avg.pol + (1 | colony) + round
##           npar    AIC    LRT Pr(Chi)  
## <none>         234.80                 
## treatment    1 233.68 0.8801 0.34817  
## avg.pol      1 235.43 2.6287 0.10495  
## round        1 238.09 5.2826 0.02154 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tot_ad.ni <- update(total_adl.ni, .~. -avg.pol)
drop1(tot_ad.ni, test = "Chisq")

## Single term deletions
## 
## Model:
## total_days_adl ~ treatment + (1 | colony) + round
##           npar    AIC    LRT  Pr(Chi)   
## <none>         235.43                   
## treatment    1 234.85 1.4154 0.234158   
## round        1 240.81 7.3743 0.006616 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(tot_ad.ni)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: total_days_adl
##            Chisq Df Pr(>Chisq)   
## treatment 1.5796  1    0.20882   
## round     8.8756  1    0.00289 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q3$treatment, q3$total_days_adl)

num_adl.ni <- q3 %>%
  group_by(treatment) %>%
  summarise(m = mean(total_days_adl))
num_adl.ni

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          2.8 
## 2 4          1.83

How many days out of 14 did an inoc bee test positive for infection?

total_det <- glm(num_days_det ~ treatment + avg.pol + round, data = q1)
drop1(total_det, test = "Chisq")

## Single term deletions
## 
## Model:
## num_days_det ~ treatment + avg.pol + round
##           Df Deviance    AIC scaled dev. Pr(>Chi)  
## <none>         86.433 84.395                       
## treatment  1   95.525 83.995      1.6003  0.20586  
## avg.pol    1  123.034 88.044      5.6495  0.01746 *
## round      2  142.224 88.363      7.9685  0.01861 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tot_det1 <- update(total_det, .~. -avg.pol)
drop1(tot_det1, test = "Chisq")

## Single term deletions
## 
## Model:
## num_days_det ~ treatment + round
##           Df Deviance    AIC scaled dev. Pr(>Chi)  
## <none>         123.03 88.044                       
## treatment  1   134.30 87.446      1.4018  0.23643  
## round      2   203.50 92.095      8.0512  0.01785 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(tot_det1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: num_days_det
##           LR Chisq Df Pr(>Chisq)  
## treatment   1.0988  1    0.29454  
## round       7.8481  2    0.01976 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q1$treatment, q1$num_days_det)

num_det <- q1 %>%
  group_by(treatment) %>%
  summarise(m = mean(num_days_det))
num_det

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          9.25
## 2 4         10.5

How many days out of 14 did a non-inoc bee test positive for infection?

total_det.ni <- glmer.nb(total_days_detected ~ treatment + avg.pol + (1|colony) + round, data = q3)

## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached

drop1(total_det.ni, test = "Chisq")

## Single term deletions
## 
## Model:
## total_days_detected ~ treatment + avg.pol + (1 | colony) + round
##           npar    AIC      LRT Pr(Chi)
## <none>         298.53                 
## treatment    1 296.69 0.154873  0.6939
## avg.pol      1 296.57 0.036981  0.8475
## round        1 296.78 0.245371  0.6204

totdetni <- update(total_det.ni, .~. -avg.pol)
drop1(totdetni, test = "Chisq")

## Single term deletions
## 
## Model:
## total_days_detected ~ treatment + (1 | colony) + round
##           npar    AIC     LRT Pr(Chi)
## <none>         296.57                
## treatment    1 294.71 0.13935  0.7089
## round        1 294.90 0.32909  0.5662

totdetni1 <- update(totdetni, .~. -round)
Anova(totdetni1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: total_days_detected
##            Chisq Df Pr(>Chisq)
## treatment 0.1693  1     0.6807

plot(q3$treatment, q3$total_days_detected)

num_detni <- q3 %>%
  group_by(treatment) %>%
  summarise(m = mean(total_days_detected))
num_detni

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          7.67
## 2 4          7.97

How many days for an inoculated bee to reach max spores?

days_max <- glm.nb(days_max_spores ~ treatment + avg.pol + round, data = q1)
drop1(days_max, test = "Chisq")

## Single term deletions
## 
## Model:
## days_max_spores ~ treatment + avg.pol + round
##           Df Deviance    AIC    LRT Pr(>Chi)
## <none>         17.045 95.278                
## treatment  1   17.046 93.279 0.0013   0.9710
## avg.pol    1   17.123 93.356 0.0781   0.7798
## round      2   20.548 94.781 3.5035   0.1735

dm1 <- update(days_max, .~. -avg.pol)
drop1(dm1, test = "Chisq")

## Single term deletions
## 
## Model:
## days_max_spores ~ treatment + round
##           Df Deviance    AIC    LRT Pr(>Chi)  
## <none>         16.970 93.355                  
## treatment  1   16.973 91.358 0.0030  0.95629  
## round      2   21.720 94.106 4.7506  0.09298 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dm2 <- update(dm1, .~. -round)
Anova(dm2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: days_max_spores
##           LR Chisq Df Pr(>Chisq)
## treatment 0.012681  1     0.9103

plot(q1$treatment, q1$days_max_spores)

num_det <- q1 %>%
  group_by(treatment) %>%
  summarise(m = mean(days_max_spores))
num_det

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3          6.62
## 2 4          6.88

How many days for a non-inoculated bee to reach max spores?

non.max <- glmer.nb(days_max_spores ~ treatment + round + avg.pol + (1|colony), data = q3)

## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached

drop1(non.max, test = "Chisq")

## Single term deletions
## 
## Model:
## days_max_spores ~ treatment + round + avg.pol + (1 | colony)
##           npar    AIC    LRT  Pr(Chi)   
## <none>         325.35                   
## treatment    1 323.36 0.0012 0.972704   
## round        1 332.56 9.2065 0.002412 **
## avg.pol      1 323.37 0.0127 0.910135   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nm1 <- update(non.max, .~. -avg.pol)
drop1(nm1, test = "Chisq")

## Single term deletions
## 
## Model:
## days_max_spores ~ treatment + round + (1 | colony)
##           npar    AIC    LRT  Pr(Chi)   
## <none>         323.37                   
## treatment    1 321.37 0.0024 0.960551   
## round        1 331.26 9.8911 0.001661 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(nm1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: days_max_spores
##             Chisq Df Pr(>Chisq)    
## treatment  0.0025  1  0.9604486    
## round     12.1176  1  0.0004995 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(q3$treatment, q3$days_max_spores)

num_max_ni <- q3 %>%
  group_by(treatment) %>%
  summarise(m = mean(days_max_spores))
num_max_ni

## # A tibble: 2 × 2
##   treatment     m
##   <fct>     <dbl>
## 1 3           7.5
## 2 4           7.1

What is the max spore count per treatment?

max.spore <- lmer(max_spores ~ treatment + round + avg.pol + (1|colony), data = q4)
Anova(max.spore)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: max_spores
##            Chisq Df Pr(>Chisq)
## treatment 0.4283  1     0.5128
## round     4.1195  2     0.1275
## avg.pol   0.4208  1     0.5165

plot(q4$treatment, q4$max_spores)

How many peaks of infection occur in IB vs. NIB per treatment?