Bumble bee disease dynamics

input data

custom_labels <- c("Control", "Fungicide",  "Fungicide + Crithidia", "Crithidia")

p <- read_csv("pollen.dates.csv", 
    col_types = cols(treatment = col_factor(levels = c("1", 
        "2", "3", "4")), pollen.start = col_date(format = "%m/%d/%Y"), 
        pollen.end = col_date(format = "%m/%d/%Y"), 
        colony.start = col_date(format = "%m/%d/%Y")))

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

p <- p %>%
  mutate(
    fungicide = factor(fungicide),
    crithidia = factor(crithidia),
    block = factor(block),
    qro = factor(qro),
    id = factor(id)
  )

shapiro.test(p$whole_dif)

## 
##  Shapiro-Wilk normality test
## 
## data:  p$whole_dif
## W = 0.76908, p-value < 2.2e-16

hist(p$whole_dif)

p$box <- bcPower(p$whole_dif, -3, gamma=1.1)
shapiro.test(p$box)

## 
##  Shapiro-Wilk normality test
## 
## data:  p$box
## W = 0.92686, p-value < 2.2e-16

hist(p$box)

p$log <- log(p$whole_dif)
shapiro.test(p$log)

## 
##  Shapiro-Wilk normality test
## 
## data:  p$log
## W = 0.9363, p-value < 2.2e-16

hist(p$log)

pol.mod1 <- lmer(box ~ crithidia + fungicide + pollen.time + workers_alive + block + (1|colony), data = p)
pol.mod.int <- lmer(box ~ crithidia*fungicide + pollen.time + workers_alive + block + (1|colony), data = p)
drop1(pol.mod.int, test = "Chisq")

## Single term deletions
## 
## Model:
## box ~ crithidia * fungicide + pollen.time + workers_alive + block + 
##     (1 | colony)
##                     npar     AIC     LRT   Pr(Chi)    
## <none>                   -3275.0                      
## pollen.time            1 -3005.9 271.188 < 2.2e-16 ***
## workers_alive          1 -3182.7  94.362 < 2.2e-16 ***
## block                  8 -3252.8  38.285 6.669e-06 ***
## crithidia:fungicide    1 -3277.0   0.042    0.8385    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

drop1(pol.mod1)

## Single term deletions
## 
## Model:
## box ~ crithidia + fungicide + pollen.time + workers_alive + block + 
##     (1 | colony)
##               npar     AIC
## <none>             -3277.0
## crithidia        1 -3274.7
## fungicide        1 -3277.0
## pollen.time      1 -3007.8
## workers_alive    1 -3184.7
## block            8 -3254.8

Anova(pol.mod1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##                  Chisq Df Pr(>Chisq)    
## crithidia       3.1483  1    0.07601 .  
## fungicide       1.4167  1    0.23395    
## pollen.time   322.3806  1  < 2.2e-16 ***
## workers_alive  96.6335  1  < 2.2e-16 ***
## block          47.2314  8  1.385e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(pol.mod1, pol.mod.int)

## Data: p
## Models:
## pol.mod1: box ~ crithidia + fungicide + pollen.time + workers_alive + block + (1 | colony)
## pol.mod.int: box ~ crithidia * fungicide + pollen.time + workers_alive + block + (1 | colony)
##             npar   AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## pol.mod1      15 -3277 -3208.0 1653.5    -3307                     
## pol.mod.int   16 -3275 -3201.5 1653.5    -3307 0.0415  1     0.8385

AIC(pol.mod1, pol.mod.int)

##             df       AIC
## pol.mod1    15 -3164.936
## pol.mod.int 16 -3156.316

drop1(pol.mod1, test = "Chisq")

## Single term deletions
## 
## Model:
## box ~ crithidia + fungicide + pollen.time + workers_alive + block + 
##     (1 | colony)
##               npar     AIC     LRT   Pr(Chi)    
## <none>             -3277.0                      
## crithidia        1 -3274.7   4.280   0.03856 *  
## fungicide        1 -3277.0   1.989   0.15847    
## pollen.time      1 -3007.8 271.190 < 2.2e-16 ***
## workers_alive    1 -3184.7  94.341 < 2.2e-16 ***
## block            8 -3254.8  38.246 6.781e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(pol.mod1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##                  Chisq Df Pr(>Chisq)    
## crithidia       3.1483  1    0.07601 .  
## fungicide       1.4167  1    0.23395    
## pollen.time   322.3806  1  < 2.2e-16 ***
## workers_alive  96.6335  1  < 2.2e-16 ***
## block          47.2314  8  1.385e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(pol.mod1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: box ~ crithidia + fungicide + pollen.time + workers_alive + block +  
##     (1 | colony)
##    Data: p
## 
## REML criterion at convergence: -3194.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.73866 -0.62142 -0.01707  0.61692  2.80644 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  colony   (Intercept) 0.0004090 0.02022 
##  Residual             0.0005754 0.02399 
## Number of obs: 733, groups:  colony, 35
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)    1.201e-01  1.447e-02   8.298
## crithidiaTRUE -1.266e-02  7.135e-03  -1.774
## fungicideTRUE -8.462e-03  7.110e-03  -1.190
## pollen.time    1.606e-03  8.944e-05  17.955
## workers_alive  1.473e-02  1.498e-03   9.830
## block4         5.313e-02  1.480e-02   3.590
## block6        -3.596e-02  1.475e-02  -2.437
## block7         1.160e-02  1.481e-02   0.783
## block8         1.948e-02  1.478e-02   1.318
## block9         7.993e-03  1.478e-02   0.541
## block10        3.679e-02  1.606e-02   2.291
## block11       -1.202e-02  1.480e-02  -0.813
## block12        8.685e-03  1.477e-02   0.588

qqnorm(resid(pol.mod1));qqline(resid(pol.mod1))

Anova(pol.mod1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##                  Chisq Df Pr(>Chisq)    
## crithidia       3.1483  1    0.07601 .  
## fungicide       1.4167  1    0.23395    
## pollen.time   322.3806  1  < 2.2e-16 ***
## workers_alive  96.6335  1  < 2.2e-16 ***
## block          47.2314  8  1.385e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

residuals <- resid(pol.mod1)
plot(residuals)

ps <- p[p$days >= 4 & p$days <=25, ]

range(ps$pollen.time)

## [1]  6 48

shapiro.test(ps$box)

## 
##  Shapiro-Wilk normality test
## 
## data:  ps$box
## W = 0.92339, p-value < 2.2e-16

hist(ps$box)

ps1 <- lmer(box ~ crithidia + fungicide + days + workers_alive + block +  (1|colony), data = ps)
Anova(ps1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##                  Chisq Df Pr(>Chisq)    
## crithidia       3.2185  1    0.07281 .  
## fungicide       1.2300  1    0.26741    
## days          229.4236  1  < 2.2e-16 ***
## workers_alive  73.0433  1  < 2.2e-16 ***
## block          43.3831  8  7.439e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pe <- emmeans(ps1, pairwise ~ crithidia, type = "response")
pe

## $emmeans
##  crithidia emmean      SE   df lower.CL upper.CL
##  FALSE      0.229 0.00526 24.2    0.219    0.240
##  TRUE       0.216 0.00546 24.0    0.205    0.227
## 
## Results are averaged over the levels of: fungicide, block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate     SE   df t.ratio p.value
##  FALSE - TRUE   0.0136 0.0076 24.3   1.794  0.0853
## 
## Results are averaged over the levels of: fungicide, block 
## Degrees-of-freedom method: kenward-roger

qqnorm(resid(ps1));qqline(resid(ps1))

residuals <- resid(ps1)
plot(residuals)

Anova(ps1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##                  Chisq Df Pr(>Chisq)    
## crithidia       3.2185  1    0.07281 .  
## fungicide       1.2300  1    0.26741    
## days          229.4236  1  < 2.2e-16 ***
## workers_alive  73.0433  1  < 2.2e-16 ***
## block          43.3831  8  7.439e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ps_sum <- ps %>%
  group_by(treatment) %>%
  summarise(mean = mean(whole_dif),
            sd = sd(whole_dif),
            n = length(whole_dif)) %>%
  mutate(se = sd/sqrt(n))
ps_sum

## # A tibble: 4 × 5
##   treatment  mean    sd     n     se
##   <fct>     <dbl> <dbl> <int>  <dbl>
## 1 1         0.514 0.341   167 0.0264
## 2 2         0.442 0.322   170 0.0247
## 3 3         0.337 0.257   177 0.0193
## 4 4         0.363 0.270   149 0.0221

range(p$days)

## [1]  2 26

pollen_subset <- p[p$pollen.time >= 2 & p$pollen.time <=42, ]
pol.mod1 <- lmer(box ~ block + crithidia + fungicide + days + (1|colony), data = pollen_subset)
drop1(pol.mod1, test = "Chisq")

## Single term deletions
## 
## Model:
## box ~ block + crithidia + fungicide + days + (1 | colony)
##           npar     AIC     LRT   Pr(Chi)    
## <none>         -3061.4                      
## block        8 -3043.1  34.324 3.548e-05 ***
## crithidia    1 -3056.9   6.502   0.01078 *  
## fungicide    1 -3060.8   2.599   0.10693    
## days         1 -2843.6 219.852 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pm1 <- update(pol.mod1, .~. -days)
drop1(pm1, test = "Chisq")

## Single term deletions
## 
## Model:
## box ~ block + crithidia + fungicide + (1 | colony)
##           npar     AIC    LRT   Pr(Chi)    
## <none>         -2843.6                     
## block        8 -2825.7 33.935 4.174e-05 ***
## crithidia    1 -2839.3  6.275   0.01225 *  
## fungicide    1 -2843.1  2.454   0.11721    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(pm1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##             Chisq Df Pr(>Chisq)    
## block     39.2939  8  4.335e-06 ***
## crithidia  4.7072  1    0.03004 *  
## fungicide  1.7459  1    0.18639    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

qqnorm(resid(pm1));qqline(resid(pm1))

Anova(pm1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: box
##             Chisq Df Pr(>Chisq)    
## block     39.2939  8  4.335e-06 ***
## crithidia  4.7072  1    0.03004 *  
## fungicide  1.7459  1    0.18639    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

residuals <- resid(pm1)
plot(residuals)

summary(pm1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: box ~ block + crithidia + fungicide + (1 | colony)
##    Data: pollen_subset
## 
## REML criterion at convergence: -2789.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.51660 -0.62354  0.01291  0.66121  2.27344 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  colony   (Intercept) 0.0005472 0.02339 
##  Residual             0.0008630 0.02938 
## Number of obs: 699, groups:  colony, 35
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    0.230527   0.013433  17.162
## block4         0.052682   0.017198   3.063
## block6        -0.040380   0.017174  -2.351
## block7        -0.001490   0.017191  -0.087
## block8         0.013249   0.017190   0.771
## block9         0.006216   0.017186   0.362
## block10        0.029247   0.018668   1.567
## block11       -0.025224   0.017179  -1.468
## block12        0.001247   0.017186   0.073
## crithidiaTRUE -0.017923   0.008261  -2.170
## fungicideTRUE -0.010916   0.008261  -1.321
## 
## Correlation of Fixed Effects:
##             (Intr) block4 block6 block7 block8 block9 blck10 blck11 blck12
## block4      -0.639                                                        
## block6      -0.640  0.500                                                 
## block7      -0.640  0.500  0.500                                          
## block8      -0.640  0.500  0.500  0.500                                   
## block9      -0.640  0.500  0.501  0.500  0.500                            
## block10     -0.589  0.460  0.461  0.460  0.460  0.460                     
## block11     -0.640  0.500  0.501  0.500  0.500  0.500  0.461              
## block12     -0.639  0.500  0.501  0.500  0.500  0.500  0.461  0.500       
## crithidTRUE -0.295 -0.001 -0.001  0.000 -0.001 -0.001  0.076 -0.001 -0.001
## fungicdTRUE -0.295  0.000  0.000  0.000  0.000  0.000 -0.077 -0.001 -0.001
##             crTRUE
## block4            
## block6            
## block7            
## block8            
## block9            
## block10           
## block11           
## block12           
## crithidTRUE       
## fungicdTRUE -0.039

range(pollen_subset$pollen.time)

## [1]  2 42

ggplot(data = p, aes(x = pollen.time, y = whole_dif, fill = treatment)) +
  geom_smooth() +
  labs(x = "Time (days)", y = "Mean Pollen Consumed (g)") +
  theme_cowplot()+
  scale_x_continuous(breaks = seq(min(p$pollen.time), max(ps$pollen.time), by = 1))

ggplot(data = pollen_subset, aes(x = pollen.time, y = whole_dif, fill = treatment)) +
  geom_smooth() +
  labs(x = "Time (days)", y = "Mean Pollen Consumed (g)") +
  theme_cowplot() +
  scale_x_continuous(breaks = seq(min(pollen_subset$pollen.time), max(pollen_subset$pollen.time), by = 1)) +
 scale_fill_manual(values = c("lightgreen", "lightblue", "sienna", "goldenrod"), 
                   labels = c("Control", "Fungicide", "Crithidia + Fungicide", "Crithidia"))

ggplot(data = ps_sum, aes(x = treatment, y = mean, fill = treatment)) +
  geom_col_pattern(
    aes(pattern = treatment),
    pattern_density = c(0, 0, 0.4, 0),  # Add density for the fourth column
    pattern_spacing = 0.03,
    position = position_dodge(0.9)
  ) +
  coord_cartesian(ylim = c(0, 0.65)) +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2, position = position_dodge(0.9)) +
  labs(x = "Treatment", y = "Average Pollen Consumed (g)") +
  annotate(
    geom = "text",
    x = 3,
    y = 0.6,
    label = "P = 0.05",
    size = 8
  ) + 
  theme_classic(base_size = 20) +
  scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
  scale_pattern_manual(values = c("none", "none", "stripe", "none")) +  # Add stripes to the fourth column
  scale_x_discrete(labels = custom_labels) +
  theme(legend.position = "none") +
  geom_segment(x = 1, xend = 2, y = 0.65, yend = 0.65, 
               lineend = "round", linejoin = "round") +
  geom_segment(x = 1, xend = 1, y = 0.64, yend = 0.66, 
               lineend = "round", linejoin = "round") +
  geom_segment(x = 2, xend = 2, y = 0.64, yend = 0.66, 
               lineend = "round", linejoin = "round") +
  geom_segment(x = 3, xend = 4, y = 0.47, yend = 0.47, 
               lineend = "round", linejoin = "round") +
  geom_segment(x = 3, xend = 3, y = 0.46, yend = 0.48, 
               lineend = "round", linejoin = "round") +
  geom_segment(x = 4, xend = 4, y = 0.46, yend = 0.48, 
               lineend = "round", linejoin = "round") +
  geom_text(x = 1.5, y = 0.65, label = "a", size = 6, vjust = -0.5) +
  geom_text(x = 3.5, y = 0.48, label = "b", size = 6, vjust = -0.5)