brood <- read_csv("brood.csv")
brood$colony <- as.factor(brood$colony)
brood$treatment <- as.factor(brood$treatment)
brood$block <- as.factor(brood$replicate)
pollen <- read_csv("pollen.csv")
pollen$colony <- as.factor(pollen$colony)
pollen$treatment <- as.factor(pollen$treatment)
pollen$block <- as.factor(pollen$block)
workers <- read_csv("workers.csv")
workers$colony <- as.factor(workers$colony)
workers$treatment <- as.factor(workers$treatment)
workers$block <- as.factor(workers$block)
workers$qro <- as.factor(workers$qro)
workers$inoculate <- as.logical(workers$inoculate)
workers_dry <- read_csv("workers_dry.csv")
workers_dry$colony <- as.factor(workers_dry$colony)
workers_dry$treatment <- as.factor(workers_dry$treatment)
workers_dry$block <- as.factor(workers_dry$block)
workers_dry$qro <- as.factor(workers_dry$qro)
workers_dry$inoculate <- as.logical(workers_dry$inoculate)
workers$bee_id <- as.factor(workers$sig_id)
duration <- read_csv("duration.csv")
duration$treatment <- as.factor(duration$treatment)
duration$block <- as.factor(duration$block)
duration$colony <- as.factor(duration$colony)
duration$qro <- as.factor(duration$qro)
drones <- read_csv("drones.csv")
drones$treatment <- as.factor(drones$treatment)
drones$block <- as.factor(drones$block)
drones$colony <- as.factor(drones$colony)
drones$id <- as.factor(drones$id)
drones$abdomen_post_ethyl <- as.numeric(drones$abdomen_post_ethyl)
drones_rf <- read_csv("drones_rf.csv")
drones_rf$treatment <- as.factor(drones_rf$treatment)
drones_rf$block <- as.factor(drones_rf$block)
drones_rf$colony <- as.factor(drones_rf$colony)
drones_rf$id <- as.factor(drones_rf$id)
drones_rf$abdomen_post_ethyl <- as.numeric(drones_rf$abdomen_post_ethyl)
qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)
qro$fungicide <- as.logical(qro$fungicide)
qro$crithidia <- as.logical(qro$crithidia)
qro_simple <- qro[c('colony', 'qro', 'fungicide', 'crithidia')]
qro_pol <- qro[c('colony', 'qro')]
pollen <- merge(pollen, qro_pol, by = "colony", all = FALSE)
brood1 <- merge(brood, duration, by = "colony", all = FALSE)
custom_labels <- c("Control", "Fungicide", "Fungicide + Crithidia", "Crithidia")
avg.pol <- pollen %>%
group_by(colony) %>%
summarise(avg.pol = mean(whole_dif))
duration <- merge(duration, avg.pol, by = "colony", all = FALSE)
new_dataframe <- duration[c('colony', 'days_active')]
brood <- merge(qro_simple, brood, by = "colony", all = FALSE)
workers <- merge(new_dataframe, workers, by = "colony", all = FALSE)
qpcr <- read_csv("qPCR results final.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"),
spores = col_number(), round = col_factor(levels = c("1",
"2", "3")), adl = col_logical(),
detected = col_logical()))
qpcr$colony <- as.factor(qpcr$colony)
qpcr$bee_id <- as.factor(qpcr$bee_id)
qpcr$plate <- as.factor(qpcr$plate)
all_bees <- read_csv("all_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"))))
## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)
all_bees$colony <- as.factor(all_bees$colony)
all_bees$bee_id <- as.factor(all_bees$bee_id)
workers_for_qpcr_merge <- read_csv("workers_for qpcr merge.csv",
col_types = cols(fungicide = col_logical(),
crithidia = col_logical(), inoculate_round = col_factor(levels = c("1",
"2", "3")), inoculate = col_logical(),
premature_death = col_logical(),
`end date` = col_date(format = "%m/%d/%Y")))
qpcr <- merge(workers_for_qpcr_merge, qpcr, by = "bee_id", all = FALSE)
all_bees <- merge(workers_for_qpcr_merge, all_bees, by = "bee_id", all = FALSE)
qpcr$inoculate <- as.logical(qpcr$inoculate)
qpcr$bee_id <- as.factor(qpcr$bee_id)
qpcr$fungicide <- as.logical(qpcr$fungicide)
qpcr$crithidia <- as.logical(qpcr$crithidia)
qpcr$qro <- as.factor(qpcr$qro)
qpcr$dry <- as.double(qpcr$dry)
## Warning: NAs introduced by coercion
# brood cells
brood.col <- lm(brood_cells~ treatment.x + block.x + workers_alive.x + qro + days_active + avg_pollen, data = brood1)
Anova(brood.col)
## Anova Table (Type II tests)
##
## Response: brood_cells
## Sum Sq Df F value Pr(>F)
## treatment.x 244.75 3 2.3913 0.09739 .
## block.x 29.69 4 0.2175 0.92569
## workers_alive.x 202.62 1 5.9392 0.02379 *
## qro 0
## days_active 0.53 1 0.0154 0.90241
## avg_pollen 1696.95 1 49.7397 5.839e-07 ***
## Residuals 716.45 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(brood.col, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ treatment.x + block.x + workers_alive.x + qro +
## days_active + avg_pollen
## Df Sum of Sq RSS AIC Pr(>Chi)
## <none> 716.45 137.67
## treatment.x 3 244.75 961.20 142.25 0.01423 *
## block.x 4 29.69 746.14 131.13 0.83343
## workers_alive.x 1 202.62 919.07 144.63 0.00275 **
## qro 0 0.00 716.45 137.67
## days_active 1 0.53 716.98 135.69 0.87094
## avg_pollen 1 1696.95 2413.40 179.39 3.786e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b1 <- update(brood.col, .~. -qro)
vif(b1)
## GVIF Df GVIF^(1/(2*Df))
## treatment.x 1.486219 3 1.068269
## block.x 8.358495 8 1.141913
## workers_alive.x 3.822801 1 1.955198
## days_active 3.191570 1 1.786497
## avg_pollen 5.733278 1 2.394426
b2 <- update(b1, .~. -block.x)
anova(brood.col, b1, b2)
## Analysis of Variance Table
##
## Model 1: brood_cells ~ treatment.x + block.x + workers_alive.x + qro +
## days_active + avg_pollen
## Model 2: brood_cells ~ treatment.x + block.x + workers_alive.x + days_active +
## avg_pollen
## Model 3: brood_cells ~ treatment.x + workers_alive.x + days_active + avg_pollen
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 716.45
## 2 21 716.45 0 0.0
## 3 29 2481.73 -8 -1765.3 6.4678 0.0002799 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(b1, b2)
## df AIC
## b1 16 241.8320
## b2 8 270.5585
vif(b2)
## GVIF Df GVIF^(1/(2*Df))
## treatment.x 1.238327 3 1.036269
## workers_alive.x 2.424156 1 1.556970
## days_active 1.734039 1 1.316829
## avg_pollen 2.116179 1 1.454709
drop1(b2, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ treatment.x + workers_alive.x + days_active + avg_pollen
## Df Sum of Sq RSS AIC Pr(>Chi)
## <none> 2481.7 166.40
## treatment.x 3 338.4 2820.2 165.00 0.20336
## workers_alive.x 1 263.1 2744.9 168.02 0.05682 .
## days_active 1 346.5 2828.3 169.10 0.03007 *
## avg_pollen 1 7460.0 9941.7 214.35 1.569e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(b1, b2)
## Analysis of Variance Table
##
## Model 1: brood_cells ~ treatment.x + block.x + workers_alive.x + days_active +
## avg_pollen
## Model 2: brood_cells ~ treatment.x + workers_alive.x + days_active + avg_pollen
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 716.45
## 2 29 2481.73 -8 -1765.3 6.4678 0.0002799 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(b1, brood.col)
## Analysis of Variance Table
##
## Model 1: brood_cells ~ treatment.x + block.x + workers_alive.x + days_active +
## avg_pollen
## Model 2: brood_cells ~ treatment.x + block.x + workers_alive.x + qro +
## days_active + avg_pollen
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 716.45
## 2 21 716.45 0 0
drop1(b1, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ treatment.x + block.x + workers_alive.x + days_active +
## avg_pollen
## Df Sum of Sq RSS AIC Pr(>Chi)
## <none> 716.45 137.67
## treatment.x 3 244.75 961.20 142.25 0.01423 *
## block.x 8 1765.28 2481.73 166.40 4.146e-07 ***
## workers_alive.x 1 202.62 919.07 144.63 0.00275 **
## days_active 1 0.53 716.98 135.69 0.87094
## avg_pollen 1 1696.95 2413.40 179.39 3.786e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b2 <- update(b1, .~. -days_active)
vif(b2)
## GVIF Df GVIF^(1/(2*Df))
## treatment.x 1.362270 3 1.052876
## block.x 4.541324 8 1.099193
## workers_alive.x 3.361139 1 1.833341
## avg_pollen 5.615940 1 2.369797
AIC(b1, b2)
## df AIC
## b1 16 241.8320
## b2 15 239.8584
anova(b1, b2)
## Analysis of Variance Table
##
## Model 1: brood_cells ~ treatment.x + block.x + workers_alive.x + days_active +
## avg_pollen
## Model 2: brood_cells ~ treatment.x + block.x + workers_alive.x + avg_pollen
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 716.45
## 2 22 716.98 -1 -0.52545 0.0154 0.9024
drop1(b2, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ treatment.x + block.x + workers_alive.x + avg_pollen
## Df Sum of Sq RSS AIC Pr(>Chi)
## <none> 716.98 135.69
## treatment.x 3 246.29 963.26 140.32 0.013902 *
## block.x 8 2111.28 2828.26 169.10 5.316e-08 ***
## workers_alive.x 1 222.37 939.35 143.42 0.001818 **
## avg_pollen 1 1723.70 2440.67 177.79 3.121e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#don't include qro
pol_consum_sum <- pollen %>%
group_by(colony) %>%
summarise(mean.pollen = mean(whole_dif))
pol_consum_sum <- as.data.frame(pol_consum_sum)
workers <- merge(workers, pol_consum_sum, by = "colony", all = FALSE)
#brood <- merge(brood, pol_consum_sum, by = "colony", all = FALSE)
duration <- merge(duration, pol_consum_sum, by = "colony", all = FALSE)
#drones <- na.omit(drones)
#brood <- na.omit(brood)
pollen$days <- pollen$`pollen ball id`
spores_sum.3.4 <- qpcr %>%
group_by(colony) %>%
summarise(mean.spores = mean(spores))
#write.csv(spores_sum, "C:/Users/runni/The Ohio State University/Sivakoff Lab - Runnion and Sivakoff - Runnion and Sivakoff/Synergism Experiment/Data analysis/Files for analysis/spores_sum.csv", row.names = FALSE)
spores_sum <- read_csv("spores_sum.csv")
brood <- merge(brood, spores_sum, by = "colony", all = FALSE)
duration <- merge(duration, spores_sum, by = "colony", all = FALSE)
pollen <- merge(pollen, spores_sum, by = "colony", all = FALSE)
pollen$fungicide <- as.logical(pollen$fungicide)
pollen$crithidia <- as.logical(pollen$crithidia)
pollen$id <- as.factor(pollen$`pollen ball id`)
spores_sum_workers.34 <- qpcr %>%
group_by(bee_id) %>%
summarise(mean.spores = mean(spores))
spores_sum_workers <- qpcr %>%
group_by(bee_id) %>%
summarise(mean.spores = mean(spores))
spores_sum_workers <- as.data.frame(spores_sum_workers)
#write.csv(spores_sum_workers, "C:/Users/runni/The Ohio State University/Sivakoff Lab - Runnion and Sivakoff - Runnion and Sivakoff/Synergism Experiment/Data analysis/Files for analysis/spores_sum_workers.csv", row.names = FALSE)
spores_sum_workers <- read_csv("spores_sum_workers.csv")
workers.34 <- merge(workers, spores_sum_workers.34, by = "bee_id", all = FALSE)
workers <- merge(workers, spores_sum_workers, by = "bee_id", all = FALSE)
shapiro.test(pollen$whole_dif)
##
## Shapiro-Wilk normality test
##
## data: pollen$whole_dif
## W = 0.76746, p-value < 2.2e-16
hist(pollen$whole_dif)
range(pollen$whole_dif)
## [1] 0.03316 1.39545
pollen$box <- bcPower(pollen$whole_dif, -5, gamma=1)
shapiro.test(pollen$box)
##
## Shapiro-Wilk normality test
##
## data: pollen$box
## W = 0.95847, p-value = 1.414e-13
hist(pollen$box)
pollen$log <- log(pollen$whole_dif)
shapiro.test(pollen$log)
##
## Shapiro-Wilk normality test
##
## data: pollen$log
## W = 0.93447, p-value < 2.2e-16
hist(pollen$log)
pollen$square <- pollen$whole_dif^2
shapiro.test(pollen$square)
##
## Shapiro-Wilk normality test
##
## data: pollen$square
## W = 0.60832, p-value < 2.2e-16
hist(pollen$square)
pollen$root <- sqrt(pollen$whole_dif)
shapiro.test(pollen$root)
##
## Shapiro-Wilk normality test
##
## data: pollen$root
## W = 0.85865, p-value < 2.2e-16
hist(pollen$root)
descdist(pollen$whole_dif, discrete = FALSE)
## summary statistics
## ------
## min: 0.03316 max: 1.39545
## median: 0.273935
## mean: 0.3984176
## estimated sd: 0.2985876
## estimated skewness: 1.595252
## estimated kurtosis: 4.524735
ggplot(pollen, aes(x = log, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.1, col = I("black")) +
scale_fill_viridis_d() + # Use viridis_d() for the color-blind friendly palette
ggtitle("Pollen Consumption(g)") +
labs(y = "Count", x = "Pollen (g)")
pol.mod <- lmer(box ~ fungicide*crithidia + id + block + days + workers_alive + (1|colony), data = pollen)
drop1(pol.mod, test = "Chisq")
## Single term deletions
##
## Model:
## box ~ fungicide * crithidia + id + block + days + workers_alive +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4009.9
## id 23 -3835.7 220.218 < 2.2e-16 ***
## block 8 -3989.3 36.664 1.327e-05 ***
## days 0 -4009.9 0.000
## workers_alive 1 -3953.7 58.246 2.314e-14 ***
## fungicide:crithidia 1 -4011.9 0.001 0.9811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pm1 <- update(pol.mod, .~. -days)
drop1(pm1, test = "Chisq")
## Single term deletions
##
## Model:
## box ~ fungicide + crithidia + id + block + workers_alive + (1 |
## colony) + fungicide:crithidia
## npar AIC LRT Pr(Chi)
## <none> -4009.9
## id 24 -3585.6 472.34 < 2.2e-16 ***
## block 8 -3989.3 36.66 1.327e-05 ***
## workers_alive 1 -3953.7 58.25 2.314e-14 ***
## fungicide:crithidia 1 -4011.9 0.00 0.9811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pol.mod1 <- lmer(box ~ block + crithidia + fungicide + workers_alive + mean.spores + days + (1|colony), data = pollen)
drop1(pol.mod1, test = "Chisq")
## Single term deletions
##
## Model:
## box ~ block + crithidia + fungicide + workers_alive + mean.spores +
## days + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3835.7
## block 8 -3813.7 37.976 7.605e-06 ***
## crithidia 1 -3835.2 2.553 0.1101
## fungicide 1 -3835.6 2.136 0.1439
## workers_alive 1 -3758.5 79.161 < 2.2e-16 ***
## mean.spores 1 -3837.7 0.031 0.8612
## days 1 -3585.6 252.147 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pm1 <- update(pol.mod1, .~. -mean.spores)
drop1(pm1, test = "Chisq")
## Single term deletions
##
## Model:
## box ~ block + crithidia + fungicide + workers_alive + days +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3837.7
## block 8 -3815.5 38.217 6.863e-06 ***
## crithidia 1 -3835.5 4.218 0.0400 *
## fungicide 1 -3837.6 2.116 0.1457
## workers_alive 1 -3760.5 79.130 < 2.2e-16 ***
## days 1 -3587.5 252.143 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pm2 <- update(pm1, .~. -fungicide)
drop1(pm2, test = "Chisq")
## Single term deletions
##
## Model:
## box ~ block + crithidia + workers_alive + days + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3837.6
## block 8 -3817.0 36.594 1.366e-05 ***
## crithidia 1 -3835.4 4.195 0.04054 *
## workers_alive 1 -3760.3 79.289 < 2.2e-16 ***
## days 1 -3587.5 252.050 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(pm2));qqline(resid(pm2))
Anova(pm2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: box
## Chisq Df Pr(>Chisq)
## block 45.7960 8 2.599e-07 ***
## crithidia 3.2042 1 0.07345 .
## workers_alive 80.4369 1 < 2.2e-16 ***
## days 294.2944 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
residuals <- resid(pm2)
plot(residuals)
wilcox.test(residuals ~ pollen$fungicide)
##
## Wilcoxon rank sum test with continuity correction
##
## data: residuals by pollen$fungicide
## W = 68072, p-value = 0.8699
## alternative hypothesis: true location shift is not equal to 0
pollen %>%
wilcox.test(whole_dif ~ crithidia, data = .)
##
## Wilcoxon rank sum test with continuity correction
##
## data: whole_dif by crithidia
## W = 86152, p-value = 1.585e-10
## alternative hypothesis: true location shift is not equal to 0
pollen %>%
ggplot(aes(x = factor(crithidia),
y = whole_dif)) +
geom_boxplot(aes(fill = factor(crithidia))) +
geom_jitter(alpha = 0.4) + # add data points
theme(legend.position = "none")
#this model says: average pollen consumed ~ yes/no Fung + yes/no Crit. + workers surviving when colony was frozen + time (id is pollen ball id, meaning it is the pollen ball number) + block + random effect of colony)
pe <- emmeans(pol.mod1, pairwise ~ crithidia, type = "response")
pe
## $emmeans
## crithidia emmean SE df lower.CL upper.CL
## FALSE 0.152 0.00431 23.1 0.143 0.160
## TRUE 0.142 0.00447 22.9 0.133 0.151
##
## Results are averaged over the levels of: block, fungicide
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## FALSE - TRUE 0.00948 0.00717 23 1.322 0.1991
##
## Results are averaged over the levels of: block, fungicide
## Degrees-of-freedom method: kenward-roger
kruskal.test(whole_dif ~ crithidia, data = pollen)
##
## Kruskal-Wallis rank sum test
##
## data: whole_dif by crithidia
## Kruskal-Wallis chi-squared = 40.923, df = 1, p-value = 1.583e-10
pairwise.wilcox.test(pollen$whole_dif, pollen$crithidia,
p.adjust.method = "BH")
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: pollen$whole_dif and pollen$crithidia
##
## FALSE
## TRUE 1.6e-10
##
## P value adjustment method: BH
ggplot(data = pollen, aes(x = days, y = whole_dif, fill = treatment)) +
geom_smooth() +
labs(x = "Time", y = "Mean Pollen Consumed (g)")
pollen_sum <- pollen %>%
group_by(treatment) %>%
summarise(mean = mean(whole_dif),
sd = sd(whole_dif),
n = length(whole_dif)) %>%
mutate(se = sd/sqrt(n))
pollen_box_sum <- pollen %>%
group_by(treatment) %>%
summarise(mean = mean(box),
sd = sd(box),
n = length(box)) %>%
mutate(se = sd/sqrt(n))
pollen_sum
## # A tibble: 4 × 5
## treatment mean sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.493 0.333 184 0.0245
## 2 2 0.426 0.317 188 0.0231
## 3 3 0.327 0.250 195 0.0179
## 4 4 0.348 0.258 169 0.0199
pollen_box_sum
## # A tibble: 4 × 5
## treatment mean sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.155 0.0313 184 0.00231
## 2 2 0.148 0.0314 188 0.00229
## 3 3 0.136 0.0285 195 0.00204
## 4 4 0.139 0.0295 169 0.00227
pollen_sum$plot <- pollen_sum$mean + pollen_sum$se
plot(pollen$id, pollen$whole_dif)
ggplot(data = pollen_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.55)) +
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.55,
label = "P = 0.04",
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.54, yend = 0.54,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 0.54, yend = 0.53,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 0.54, yend = 0.53,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 4, y = 0.42, yend = 0.42,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 0.42, yend = 0.41,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 0.42, yend = 0.41,
lineend = "round", linejoin = "round") +
geom_text(x = 1.5, y = 0.55, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 3.5, y = 0.43, label = "b", size = 6, vjust = -0.5)
ggplot(data = pollen_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.55)) +
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.55,
label = "P = 0.07",
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")
duration$fungicide <- as.logical(duration$fungicide)
duration$crithidia <- as.logical(duration$crithidia)
cbw1 <- glm(cbind(workers_alive, workers_dead) ~ fungicide*crithidia + mean.pollen + block + days_active + mean.spores, data = duration, family = binomial("logit"))
Anova(cbw1)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.7229 1 0.39521
## crithidia 0.3903 1 0.53213
## mean.pollen 18.6266 1 1.59e-05 ***
## block 16.9088 8 0.03107 *
## days_active 0.9681 1 0.32516
## mean.spores 0.5491 1 0.45867
## fungicide:crithidia 1.2877 1 0.25647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(cbw1, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(workers_alive, workers_dead) ~ fungicide * crithidia +
## mean.pollen + block + days_active + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 23.501 91.157
## mean.pollen 1 42.128 107.783 18.6266 1.59e-05 ***
## block 8 40.410 92.065 16.9088 0.03107 *
## days_active 1 24.469 90.125 0.9681 0.32516
## mean.spores 1 24.050 89.706 0.5491 0.45867
## fungicide:crithidia 1 24.789 90.444 1.2877 0.25647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cbw1 <- glm(cbind(workers_alive, workers_dead) ~ fungicide + crithidia + mean.pollen + block + days_active + mean.spores, data = duration, family = binomial("logit"))
Anova(cbw1)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.7229 1 0.39521
## crithidia 0.3903 1 0.53213
## mean.pollen 19.0696 1 1.26e-05 ***
## block 17.5505 8 0.02486 *
## days_active 0.6280 1 0.42809
## mean.spores 0.6865 1 0.40735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(cbw1, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(workers_alive, workers_dead) ~ fungicide + crithidia +
## mean.pollen + block + days_active + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 24.789 90.444
## fungicide 1 25.512 89.167 0.7229 0.39521
## crithidia 1 25.179 88.835 0.3903 0.53213
## mean.pollen 1 43.858 107.514 19.0696 1.26e-05 ***
## block 8 42.339 91.995 17.5505 0.02486 *
## days_active 1 25.417 89.072 0.6280 0.42809
## mean.spores 1 25.475 89.131 0.6865 0.40735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cbw2 <- update(cbw1, .~. -days_active)
drop1(cbw2, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(workers_alive, workers_dead) ~ fungicide + crithidia +
## mean.pollen + block + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 25.417 89.072
## fungicide 1 26.543 88.198 1.1260 0.28863
## crithidia 1 25.966 87.622 0.5492 0.45863
## mean.pollen 1 52.028 113.684 26.6116 2.487e-07 ***
## block 8 44.341 91.997 18.9246 0.01527 *
## mean.spores 1 26.298 87.953 0.8810 0.34793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cbw3 <- update(cbw2, .~. -mean.spores)
Anova(cbw3)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.9845 1 0.321099
## crithidia 3.8355 1 0.050179 .
## mean.pollen 25.9146 1 3.569e-07 ***
## block 20.3944 8 0.008943 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(cbw3));qqline(resid(cbw3))
Anova(cbw3)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.9845 1 0.321099
## crithidia 3.8355 1 0.050179 .
## mean.pollen 25.9146 1 3.569e-07 ***
## block 20.3944 8 0.008943 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(cbw3)
##
## Call:
## glm(formula = cbind(workers_alive, workers_dead) ~ fungicide +
## crithidia + mean.pollen + block, family = binomial("logit"),
## data = duration)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.90189 -0.55432 0.09474 0.65984 1.53335
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0986 1.0863 -1.932 0.0534 .
## fungicideTRUE 0.4434 0.4506 0.984 0.3251
## crithidiaTRUE -0.8905 0.4551 -1.957 0.0504 .
## mean.pollen 10.3678 2.5161 4.121 3.78e-05 ***
## block4 14.8725 3589.0517 0.004 0.9967
## block6 0.3725 0.7909 0.471 0.6377
## block7 -1.1756 0.7908 -1.487 0.1371
## block8 0.6140 0.8993 0.683 0.4948
## block9 0.9543 0.9529 1.002 0.3166
## block10 -2.4598 1.0361 -2.374 0.0176 *
## block11 -0.1592 0.7737 -0.206 0.8369
## block12 -1.9318 0.9305 -2.076 0.0379 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 103.807 on 34 degrees of freedom
## Residual deviance: 26.298 on 23 degrees of freedom
## AIC: 87.953
##
## Number of Fisher Scoring iterations: 18
worker_sum <-duration %>%
group_by(treatment) %>%
summarise(m = mean(workers_alive),
sd = sd(workers_alive))
worker_sum
## # A tibble: 4 × 3
## treatment m sd
## <fct> <dbl> <dbl>
## 1 1 4.33 0.866
## 2 2 3.78 1.20
## 3 3 2.78 1.86
## 4 4 2.75 2.05
worker_sum <-duration %>%
group_by(treatment) %>%
summarise(m = mean(workers_alive),
sd = sd(workers_alive),
l = length(workers_alive)) %>%
mutate(se = sd/sqrt(l))
worker_sum
## # A tibble: 4 × 5
## treatment m sd l se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 4.33 0.866 9 0.289
## 2 2 3.78 1.20 9 0.401
## 3 3 2.78 1.86 9 0.619
## 4 4 2.75 2.05 8 0.726
workers$prob <- workers$days_alive / workers$days_active
worker_prob_sum <-workers %>%
group_by(treatment) %>%
summarise(m = mean(prob),
sd = sd(prob),
l = length(prob)) %>%
mutate(se = sd/sqrt(l))
worker_prob_sum$plot <- worker_prob_sum$m + worker_prob_sum$se
worker_prob_sum$treatment <- as.factor(worker_prob_sum$treatment)
ggplot(data = worker_prob_sum, aes(x = treatment, y = m, fill = treatment)) +
geom_col_pattern(
aes(pattern = treatment),
pattern_density = c(0, 0, 0, 0.4), # Add density for the fourth column
pattern_spacing = 0.03,
position = position_dodge(0.9)) +
coord_cartesian(ylim = c(0.5, 1.05)) +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Probability") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 3.5,
y = 1.05,
label = "P = 0.05",
size = 7
) + # Add stripes to the fourth column
scale_fill_manual(values = c("lightgreen", "lightblue", "grey", "lightblue")) +
scale_pattern_manual(values = c("none", "none", "none", "stripe")) + # Add stripes to the fourth column
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none") +
geom_segment(x = 1, xend = 2, y = 1, yend = 1,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 0.98, yend = 1.02,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 0.98, yend = 1.02,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 4, y = 0.9, yend = 0.9,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 0.88, yend = 0.92,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 0.88, yend = 0.92,
lineend = "round", linejoin = "round") +
geom_text(x = 1.5, y = 1.01, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 3.5, y = 0.93, label = "b", size = 6, vjust = -0.5)
library(survival)
library(coxme)
## Warning: package 'coxme' was built under R version 4.2.3
library(survminer)
## Warning: package 'survminer' was built under R version 4.2.3
workers$censor_status <- ifelse(workers$premature_death == 0, 1, 2)
qpcr$censor_status <- ifelse(qpcr$premature_death == 0, 1, 2)
workers$fungicide <- as.logical(workers$fungicide)
workers$crithidia <- as.logical(workers$crithidia)
all_bees$bee_id <-as.factor(all_bees$bee_id)
workers$inoculate_round <- as.factor(workers$inoculate_round)
res.cox <- coxph(Surv(days_alive, censor_status) ~ crithidia + fungicide + block + qro + inoculate_round + avg_pollen, data = workers)
## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 3 ; coefficient may be infinite.
summary(res.cox)
## Call:
## coxph(formula = Surv(days_alive, censor_status) ~ crithidia +
## fungicide + block + qro + inoculate_round + avg_pollen, data = workers)
##
## n= 179, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## crithidiaTRUE 3.432e-01 1.410e+00 3.420e-01 1.004 0.315610
## fungicideTRUE -2.933e-01 7.458e-01 3.393e-01 -0.864 0.387346
## block4 -1.406e+01 7.850e-07 3.319e+03 -0.004 0.996621
## block6 -3.711e-01 6.899e-01 6.349e-01 -0.585 0.558854
## block7 3.893e-01 1.476e+00 9.518e-01 0.409 0.682562
## block8 -1.463e-01 8.639e-01 9.745e-01 -0.150 0.880672
## block9 -2.575e-01 7.730e-01 8.761e-01 -0.294 0.768859
## block10 2.322e+00 1.019e+01 8.600e-01 2.700 0.006935 **
## block11 9.888e-02 1.104e+00 8.432e-01 0.117 0.906653
## block12 5.961e-01 1.815e+00 6.372e-01 0.936 0.349515
## qro3 NA NA 0.000e+00 NA NA
## qro4 NA NA 0.000e+00 NA NA
## qro5 NA NA 0.000e+00 NA NA
## qro6 NA NA 0.000e+00 NA NA
## inoculate_round2 3.234e-01 1.382e+00 6.743e-01 0.480 0.631486
## inoculate_round3 1.026e+00 2.790e+00 7.284e-01 1.409 0.158982
## avg_pollen -7.009e+00 9.040e-04 1.981e+00 -3.537 0.000404 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## crithidiaTRUE 1.410e+00 7.095e-01 7.210e-01 2.75554
## fungicideTRUE 7.458e-01 1.341e+00 3.836e-01 1.45015
## block4 7.850e-07 1.274e+06 0.000e+00 Inf
## block6 6.899e-01 1.449e+00 1.988e-01 2.39476
## block7 1.476e+00 6.776e-01 2.285e-01 9.53293
## block8 8.639e-01 1.158e+00 1.279e-01 5.83398
## block9 7.730e-01 1.294e+00 1.388e-01 4.30459
## block10 1.019e+01 9.809e-02 1.890e+00 54.99924
## block11 1.104e+00 9.059e-01 2.114e-01 5.76375
## block12 1.815e+00 5.509e-01 5.206e-01 6.32845
## qro3 NA NA NA NA
## qro4 NA NA NA NA
## qro5 NA NA NA NA
## qro6 NA NA NA NA
## inoculate_round2 1.382e+00 7.237e-01 3.686e-01 5.18078
## inoculate_round3 2.790e+00 3.584e-01 6.692e-01 11.63115
## avg_pollen 9.040e-04 1.106e+03 1.861e-05 0.04392
##
## Concordance= 0.775 (se = 0.031 )
## Likelihood ratio test= 55.48 on 13 df, p=3e-07
## Wald test = 36.31 on 13 df, p=5e-04
## Score (logrank) test = 53.26 on 13 df, p=8e-07
fit <- survfit(res.cox, data = workers)
ggsurvplot(fit, data = workers, color = "#2E9FDF", ggtheme = theme_minimal())
## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
## instead of color = '#2E9FDF'
library(survminer)
all_bees$censor_status <- as.double(all_bees$censor_status)
res.cox <- coxme(Surv(days_since_innoculation, censor_status) ~ treatment + (1|bee_id), data = all_bees)
res.cox
## Cox mixed-effects model fit by maximum likelihood
## Data: all_bees
## events, n = 96, 2450
## Iterations= 50 310
## NULL Integrated Fitted
## Log-likelihood -626.0995 -491.9001 -410.6103
##
## Chisq df p AIC BIC
## Integrated loglik 268.40 4.00 0 260.40 250.14
## Penalized loglik 430.98 34.49 0 361.99 273.54
##
## Model: Surv(days_since_innoculation, censor_status) ~ treatment + (1 | bee_id)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## treatment2 1.155925 3.176960 0.5278473 2.19 0.0290
## treatment3 1.568494 4.799414 0.5122320 3.06 0.0022
## treatment4 1.422254 4.146455 0.5168415 2.75 0.0059
##
## Random effects
## Group Variable Std Dev Variance
## bee_id Intercept 0.8062065 0.6499689
Anova(res.cox)
## Analysis of Deviance Table (Type II tests)
##
## Response: Surv(days_since_innoculation, censor_status)
## Df Chisq Pr(>Chisq)
## treatment 3 9.9859 0.01869 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emm.cox <- emmeans(res.cox, pairwise ~ treatment, type = "response")
pairs(emm.cox)
## contrast ratio SE df null z.ratio p.value
## treatment1 / treatment2 0.315 0.166 Inf 1 -2.190 0.1260
## treatment1 / treatment3 0.208 0.107 Inf 1 -3.062 0.0118
## treatment1 / treatment4 0.241 0.125 Inf 1 -2.752 0.0302
## treatment2 / treatment3 0.662 0.251 Inf 1 -1.089 0.6963
## treatment2 / treatment4 0.766 0.294 Inf 1 -0.693 0.8998
## treatment3 / treatment4 1.157 0.418 Inf 1 0.404 0.9776
##
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
emmdf <- as.data.frame(emm.cox$contrasts)
emmdf
## contrast ratio SE df null z.ratio p.value
## treatment1 / treatment2 0.3147663 0.1661485 Inf 1 -2.190 0.1260
## treatment1 / treatment3 0.2083588 0.1067280 Inf 1 -3.062 0.0118
## treatment1 / treatment4 0.2411699 0.1246466 Inf 1 -2.752 0.0302
## treatment2 / treatment3 0.6619475 0.2507820 Inf 1 -1.089 0.6963
## treatment2 / treatment4 0.7661871 0.2944372 Inf 1 -0.693 0.8998
## treatment3 / treatment4 1.1574741 0.4184718 Inf 1 0.404 0.9776
##
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
emmdf <- setDT(emmdf)
workcld <- cld(object = emm.cox,
adjust = "Tukey",
alpha = 0.05,
Letters = letters)
workcld
## treatment response SE df asymp.LCL asymp.UCL .group
## 1 0.359 0.126 Inf 0.150 0.861 a
## 2 1.142 0.291 Inf 0.605 2.155 ab
## 4 1.490 0.367 Inf 0.807 2.753 b
## 3 1.725 0.419 Inf 0.942 3.161 b
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
require("survival")
ggsurvplot(survfit(Surv(days_since_innoculation, censor_status) ~ treatment, data = all_bees),
legend.title = "",
censor.shape = 124,
censor.size = 2.5)
ggsurvplot(
survfit(Surv(days_since_innoculation, censor_status) ~ treatment, data = all_bees),
legend.title = "",
censor.shape = 124,
censor.size = 2.5,
ylim = c(0.7, 1),
palette = c("green", "lightblue", "darkblue", "orange")
)
workers$fungicide <- as.logical(workers$fungicide)
workers$crithidia <- as.logical(workers$crithidia)
workers <- na.omit(workers)
dayswrk <- glmer.nb(days_alive ~ fungicide*crithidia + avg_pollen + inoculate + block + (1|colony), data = workers)
drop1(dayswrk, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide * crithidia + avg_pollen + inoculate +
## block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1314.5
## avg_pollen 1 1313.5 1.0424 0.30727
## inoculate 1 1312.9 0.3888 0.53291
## block 8 1314.1 15.5934 0.04858 *
## fungicide:crithidia 1 1315.1 2.6141 0.10592
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dayswrk <- glmer(days_alive ~ fungicide + crithidia + block + inoculate + avg_pollen + mean.spores + (1|colony), data = workers, family = "poisson")
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
summary(dayswrk)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: days_alive ~ fungicide + crithidia + block + inoculate + avg_pollen +
## mean.spores + (1 | colony)
## Data: workers
##
## AIC BIC logLik deviance df.resid
## 1346.4 1393.8 -658.2 1316.4 159
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8313 -0.0601 0.2587 0.5267 2.5613
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.004174 0.06461
## Number of obs: 174, groups: colony, 36
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.755458 0.076371 49.174 < 2e-16 ***
## fungicideTRUE 0.016932 0.033397 0.507 0.61217
## crithidiaTRUE 0.008048 0.041210 0.195 0.84516
## block4 -0.179672 0.089948 -1.998 0.04577 *
## block6 -0.018198 0.070036 -0.260 0.79499
## block7 -0.228252 0.071763 -3.181 0.00147 **
## block8 -0.165446 0.071174 -2.325 0.02010 *
## block9 -0.114500 0.067803 -1.689 0.09127 .
## block10 -0.223710 0.075007 -2.983 0.00286 **
## block11 -0.147517 0.069393 -2.126 0.03352 *
## block12 -0.109205 0.068282 -1.599 0.10975
## inoculateTRUE -0.022153 0.031750 -0.698 0.48535
## avg_pollen 0.143154 0.131200 1.091 0.27522
## mean.spores -0.001938 0.001824 -1.063 0.28780
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
dayswrk <- glmer.nb(days_alive ~ fungicide + crithidia + block + inoculate + avg_pollen + mean.spores + (1|colony), data = workers)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
drop1(dayswrk, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + inoculate + avg_pollen +
## mean.spores + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1315.8
## fungicide 1 1314.0 0.1770 0.67396
## crithidia 1 1314.0 0.0979 0.75442
## block 8 1313.7 13.7958 0.08725 .
## inoculate 1 1314.0 0.1692 0.68084
## avg_pollen 1 1314.8 0.9752 0.32338
## mean.spores 1 1315.1 1.2108 0.27118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dayswrk1 <- update(dayswrk, .~. -inoculate)
drop1(dayswrk1, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + avg_pollen + mean.spores +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1314.0
## fungicide 1 1312.2 0.1686 0.68133
## crithidia 1 1312.2 0.1254 0.72324
## block 8 1311.8 13.7660 0.08807 .
## avg_pollen 1 1313.0 1.0229 0.31184
## mean.spores 1 1313.4 1.3890 0.23857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dayswrk2 <- update(dayswrk1, .~. -avg_pollen)
drop1(dayswrk2, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + mean.spores + (1 |
## colony)
## npar AIC LRT Pr(Chi)
## <none> 1313.0
## fungicide 1 1311.1 0.0277 0.8677
## crithidia 1 1311.0 0.0030 0.9561
## block 8 1309.9 12.8615 0.1167
## mean.spores 1 1312.5 1.4633 0.2264
dayswrk3 <- update(dayswrk2, .~. -mean.spores)
drop1(dayswrk3, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1312.5
## fungicide 1 1310.5 0.0102 0.91945
## crithidia 1 1311.2 0.7047 0.40121
## block 8 1310.5 13.9981 0.08182 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dayswrk4 <- update(dayswrk3, .~. -fungicide)
drop1(dayswrk4, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ crithidia + block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1310.5
## crithidia 1 1309.2 0.7041 0.40142
## block 8 1308.5 13.9994 0.08178 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dayswrk1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: days_alive
## Chisq Df Pr(>Chisq)
## fungicide 0.1686 1 0.68133
## crithidia 0.1254 1 0.72325
## block 14.7567 8 0.06405 .
## avg_pollen 1.0227 1 0.31189
## mean.spores 1.3887 1 0.23863
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
worker_days_sum <-workers %>%
group_by(treatment) %>%
summarise(m = mean(days_alive),
sd = sd(days_alive),
l = length(days_alive)) %>%
mutate(se = sd/sqrt(l))
worker_days_sum$plot <- worker_days_sum$m + worker_days_sum$se
worker_days_sum$treatment <- as.factor(worker_days_sum$treatment)
ggplot(data = worker_days_sum, aes(x = treatment, y = m, fill = treatment, pattern = treatment)) +
geom_col_pattern(
aes(pattern = treatment),
pattern_density = c(0, 0, 0, 0.4), # Add density for the fourth column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
coord_cartesian(ylim = c(20, 45)) +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Days") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 45,
label = "P > 0.05",
size = 7
) + scale_fill_manual(values = c("lightgreen", "lightblue", "grey", "lightblue")) +
scale_pattern_manual(values = c("none", "none", "none", "stripe")) + # Add stripes to the fourth column
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none")
durmod <- glm.nb(days_active ~ fungicide + crithidia + mean.pollen + days_first_ov, data = duration)
## 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(durmod, test = "Chisq")
## Single term deletions
##
## Model:
## days_active ~ fungicide + crithidia + mean.pollen + days_first_ov
## Df Deviance AIC LRT Pr(>Chi)
## <none> 11.049 212.82
## fungicide 1 11.064 210.84 0.01539 0.9013
## crithidia 1 11.504 211.28 0.45466 0.5001
## mean.pollen 1 12.211 211.99 1.16261 0.2809
## days_first_ov 1 11.997 211.77 0.94845 0.3301
durmod <- glm.nb(days_active ~ fungicide + crithidia + mean.pollen + mean.spores, data = duration)
## 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(durmod, test = "Chisq")
## Single term deletions
##
## Model:
## days_active ~ fungicide + crithidia + mean.pollen + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 12.057 219.58
## fungicide 1 12.057 217.58 0.0000 0.99839
## crithidia 1 12.224 217.75 0.1667 0.68304
## mean.pollen 1 16.448 221.97 4.3914 0.03612 *
## mean.spores 1 12.071 217.60 0.0141 0.90547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dm1 <- update(durmod, .~. -mean.spores)
## 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(dm1 , test = "Chisq")
## Single term deletions
##
## Model:
## days_active ~ fungicide + crithidia + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 12.071 217.60
## fungicide 1 12.071 215.60 0.0000 0.99599
## crithidia 1 12.436 215.96 0.3644 0.54606
## mean.pollen 1 16.625 220.15 4.5536 0.03285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dm2 <- update(dm1, .~. -fungicide)
## 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(dm2, test = "Chisq")
## Single term deletions
##
## Model:
## days_active ~ crithidia + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 12.071 215.60
## crithidia 1 12.436 213.96 0.3644 0.54607
## mean.pollen 1 16.694 218.22 4.6232 0.03154 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(durmod)
##
## Call:
## glm.nb(formula = days_active ~ fungicide + crithidia + mean.pollen +
## mean.spores, data = duration, init.theta = 2516957.549, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2521 -0.2703 0.1053 0.2854 1.1096
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.8972984 0.0744619 52.339 <2e-16 ***
## fungicideTRUE 0.0001023 0.0508093 0.002 0.9984
## crithidiaTRUE 0.0269431 0.0659250 0.409 0.6828
## mean.pollen -0.2584826 0.1241103 -2.083 0.0373 *
## mean.spores 0.0003879 0.0032641 0.119 0.9054
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2516958) family taken to be 1)
##
## Null deviance: 18.284 on 34 degrees of freedom
## Residual deviance: 12.057 on 30 degrees of freedom
## AIC: 221.58
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2516958
## Std. Err.: 67035528
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -209.584
Anova(durmod)
## Analysis of Deviance Table (Type II tests)
##
## Response: days_active
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0000 1 0.99839
## crithidia 0.1667 1 0.68304
## mean.pollen 4.3914 1 0.03612 *
## mean.spores 0.0141 1 0.90547
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(duration$treatment, duration$days_active)
hist(workers_dry$dry)
shapiro.test(workers_dry$dry)
##
## Shapiro-Wilk normality test
##
## data: workers_dry$dry
## W = 0.96362, p-value = 0.0001443
workers_dry$logdry <- log(workers_dry$dry)
shapiro.test(workers_dry$logdry)
##
## Shapiro-Wilk normality test
##
## data: workers_dry$logdry
## W = 0.98991, p-value = 0.2444
hist(workers_dry$logdry)
wrkdry <- lmer(logdry ~ fungicide*crithidia + avg_pollen + inoculate +block + (1|colony), data = workers_dry)
drop1(wrkdry, test = "Chisq")
## Single term deletions
##
## Model:
## logdry ~ fungicide * crithidia + avg_pollen + inoculate + block +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 66.754
## avg_pollen 1 84.670 19.916 8.092e-06 ***
## inoculate 1 64.992 0.238 0.6257
## block 8 86.760 36.006 1.752e-05 ***
## fungicide:crithidia 1 65.151 0.397 0.5289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wrkdry <- lmer(logdry ~ fungicide + crithidia + avg_pollen + inoculate +block + (1|colony), data = workers_dry)
drop1(wrkdry, test = "Chisq")
## Single term deletions
##
## Model:
## logdry ~ fungicide + crithidia + avg_pollen + inoculate + block +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 65.151
## fungicide 1 65.866 2.715 0.099389 .
## crithidia 1 70.823 7.672 0.005608 **
## avg_pollen 1 82.788 19.638 9.360e-06 ***
## inoculate 1 63.400 0.250 0.617267
## block 8 84.831 35.680 2.009e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wrkdry1 <- update(wrkdry, .~. -inoculate)
drop1(wrkdry1, test = "Chisq")
## Single term deletions
##
## Model:
## logdry ~ fungicide + crithidia + avg_pollen + block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 63.400
## fungicide 1 64.165 2.765 0.096363 .
## crithidia 1 69.016 7.616 0.005786 **
## avg_pollen 1 81.172 19.772 8.725e-06 ***
## block 8 83.155 35.754 1.947e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wm2 <- update(wrkdry1, .~. -fungicide)
drop1(wm2, test = "Chisq")
## Single term deletions
##
## Model:
## logdry ~ crithidia + avg_pollen + block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 64.165
## crithidia 1 68.571 6.406 0.01137 *
## avg_pollen 1 84.273 22.108 2.578e-06 ***
## block 8 82.775 34.610 3.149e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrkdry, wm2, test = "Chisq")
## Data: workers_dry
## Models:
## wm2: logdry ~ crithidia + avg_pollen + block + (1 | colony)
## wrkdry: logdry ~ fungicide + crithidia + avg_pollen + inoculate + block + (1 | colony)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## wm2 13 64.165 105.45 -19.082 38.165
## wrkdry 15 65.151 112.79 -17.575 35.151 3.0145 2 0.2215
summary(wrkdry1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logdry ~ fungicide + crithidia + avg_pollen + block + (1 | colony)
## Data: workers_dry
##
## REML criterion at convergence: 77.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.02299 -0.56652 -0.04535 0.61291 2.16906
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.008674 0.09313
## Residual 0.070615 0.26574
## Number of obs: 177, groups: colony, 36
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.133210 0.120772 -25.943
## fungicide -0.072102 0.052330 -1.378
## crithidia -0.130746 0.055229 -2.367
## avg_pollen 0.849884 0.203664 4.173
## block4 -0.293828 0.140549 -2.091
## block6 0.122021 0.112098 1.089
## block7 -0.126237 0.107821 -1.171
## block8 0.004312 0.110376 0.039
## block9 -0.418020 0.107027 -3.906
## block10 -0.329383 0.117314 -2.808
## block11 -0.297229 0.109508 -2.714
## block12 -0.212682 0.107617 -1.976
##
## Correlation of Fixed Effects:
## (Intr) fungcd crithd avg_pl block4 block6 block7 block8 block9
## fungicide -0.386
## crithidia -0.481 0.105
## avg_pollen -0.722 0.253 0.400
## block4 0.134 -0.165 -0.260 -0.650
## block6 -0.621 0.063 0.118 0.278 0.181
## block7 -0.399 -0.020 -0.027 -0.049 0.408 0.458
## block8 -0.244 -0.064 -0.101 -0.254 0.532 0.390 0.491
## block9 -0.390 -0.018 -0.028 -0.070 0.424 0.456 0.497 0.500
## block10 -0.103 -0.105 -0.166 -0.414 0.615 0.318 0.471 0.545 0.483
## block11 -0.565 0.054 0.080 0.181 0.252 0.515 0.474 0.425 0.474
## block12 -0.348 -0.032 -0.050 -0.126 0.459 0.438 0.497 0.512 0.504
## blck10 blck11
## fungicide
## crithidia
## avg_pollen
## block4
## block6
## block7
## block8
## block9
## block10
## block11 0.368
## block12 0.504 0.461
Anova(wrkdry1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logdry
## Chisq Df Pr(>Chisq)
## fungicide 1.8984 1 0.16825
## crithidia 5.6043 1 0.01792 *
## avg_pollen 17.4136 1 3.007e-05 ***
## block 40.8700 8 2.205e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(wrkdry1));qqline(resid(wrkdry1))
pe <- emmeans(wrkdry1, pairwise ~ crithidia, type = "response")
pe
## $emmeans
## crithidia emmean SE df lower.CL upper.CL
## 0 -2.98 0.0372 23.8 -3.06 -2.91
## 1 -3.11 0.0377 24.1 -3.19 -3.04
##
## 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
## crithidia0 - crithidia1 0.131 0.0552 23.9 2.367 0.0264
##
## Results are averaged over the levels of: fungicide, block
## Degrees-of-freedom method: kenward-roger
pairs(pe)
## contrast estimate SE df t.ratio p.value
## crithidia0 - crithidia1 0.131 0.0552 23.9 2.367 0.0264
##
## Results are averaged over the levels of: fungicide, block
## Degrees-of-freedom method: kenward-roger
wrkdrysum <- workers_dry %>%
group_by(treatment) %>%
summarise(m = mean(dry),
sd = sd(dry),
n = length(dry)) %>%
mutate(se = sd/sqrt(n))
wrkdrysum
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.0586 0.0186 44 0.00280
## 2 2 0.0528 0.0169 45 0.00252
## 3 3 0.0416 0.0132 44 0.00199
## 4 4 0.0487 0.0182 44 0.00275
wtuk.means <- emmeans(object = wrkdry1,
specs = "crithidia",
adjust = "Tukey",
type = "response")
wtuk.means
## crithidia emmean SE df lower.CL upper.CL
## 0 -2.98 0.0372 23.8 -3.07 -2.89
## 1 -3.11 0.0377 24.1 -3.20 -3.02
##
## Results are averaged over the levels of: fungicide, block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
w.cld.model <- cld(object = wtuk.means,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
w.cld.model
## crithidia emmean SE df lower.CL upper.CL .group
## 1 -3.11 0.0377 24.1 -3.20 -3.02 a
## 0 -2.98 0.0372 23.8 -3.07 -2.89 b
##
## Results are averaged over the levels of: fungicide, block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
ggplot(data = wrkdrysum, aes(x = treatment, y = m, fill = treatment)) +
geom_col(col = "black") +
coord_cartesian(ylim = c(0, 0.08)) +
scale_fill_viridis_d() +
labs(x = "Treatment", y = "Worker Dry Weight (g)") +
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)
) +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
theme_classic(base_size = 20) +
theme(legend.position = "none") +
annotate(
geom = "text",
x = 4,
y = 0.079,
label = "P = 0.02",
size = 7
) +
theme(legend.position = "none") +
scale_x_discrete(labels = custom_labels) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_manual(values = c("none", "none", "stripe", "none")) +
geom_segment(x = 1, xend = 2, y = 0.07, yend = 0.07,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 0.07, yend = 0.069,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 0.07, yend = 0.069,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 4, y = 0.06, yend = 0.06,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 0.06, yend = 0.059,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 0.06, yend = 0.059,
lineend = "round", linejoin = "round") +
geom_text(x = 1.5, y = 0.071, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 3.5, y = 0.061, label = "b", size = 6, vjust = -0.5)
ggplot(data = wrkdrysum, aes(x = treatment, y = m, fill = treatment, pattern = treatment)) +
geom_col_pattern(
aes(pattern = treatment),
pattern_density = c(0, 0, 0, 0.4), # Add density for the fourth column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
coord_cartesian(ylim = c(0, 0.07)) +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Worker Dry Weight (g)") +
annotate(
geom = "text",
x = c(1, 2, 3, 4),
y = wrkdrysum$m + wrkdrysum$se + 0.01, # Adjust the y-position as needed
label = c("a", "a", "b", "ab"),
size = 8
) +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 4,
y = 0.07,
label = "P = 0.05",
size = 7
) +
scale_fill_manual(values = c("lightgreen", "lightblue", "grey", "lightblue")) +
scale_pattern_manual(values = c("none", "none", "none", "stripe")) + # Add stripes to the fourth column
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none")
duration$fungicide <- as.logical(duration$fungicide)
duration$crithidia <- as.logical(duration$crithidia)
ov <- glm.nb(days_first_ov ~ fungicide*crithidia + avg.pol + workers_alive + block, data = duration)
## 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(ov, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide * crithidia + avg.pol + workers_alive +
## block
## Df Deviance AIC LRT Pr(>Chi)
## <none> 16.973 189.49
## avg.pol 1 19.910 190.43 2.9371 0.08657 .
## workers_alive 1 17.675 188.19 0.7019 0.40215
## block 8 26.199 182.71 9.2259 0.32360
## fungicide:crithidia 1 17.548 188.06 0.5753 0.44816
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ov.pois <- glm(days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive + block, data = duration, family = "poisson")
summary(ov.pois)
##
## Call:
## glm(formula = days_first_ov ~ fungicide + crithidia + avg.pol +
## workers_alive + block, family = "poisson", data = duration)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.67973 -0.48815 -0.09477 0.52635 1.30610
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.84526 0.28315 10.049 <2e-16 ***
## fungicideTRUE -0.04245 0.10844 -0.391 0.6955
## crithidiaTRUE -0.13132 0.11964 -1.098 0.2723
## avg.pol -0.96857 0.60640 -1.597 0.1102
## workers_alive -0.05751 0.05788 -0.994 0.3204
## block4 0.48837 0.33516 1.457 0.1451
## block6 0.48388 0.21264 2.276 0.0229 *
## block7 0.22638 0.23090 0.980 0.3269
## block8 0.26358 0.27192 0.969 0.3324
## block9 0.19297 0.23071 0.836 0.4029
## block10 0.37820 0.28240 1.339 0.1805
## block11 0.21588 0.22163 0.974 0.3300
## block12 0.13299 0.24786 0.537 0.5916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 48.911 on 33 degrees of freedom
## Residual deviance: 17.549 on 21 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 188.06
##
## Number of Fisher Scoring iterations: 4
qqnorm(resid(ov.pois));qqline(resid(ov.pois))
ov <- glm.nb(days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive + block, data = duration)
## 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
anova(ov.pois, ov, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive +
## block
## Model 2: days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive +
## block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 21 17.549
## 2 21 17.548 0 0.00054669
qqnorm(resid(ov));qqline(resid(ov))
AIC(ov.pois, ov)
## df AIC
## ov.pois 13 188.0626
## ov 14 190.0630
drop1(ov.pois, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive +
## block
## Df Deviance AIC LRT Pr(>Chi)
## <none> 17.549 188.06
## fungicide 1 17.702 186.22 0.1531 0.6956
## crithidia 1 18.763 187.28 1.2140 0.2705
## avg.pol 1 20.184 188.70 2.6358 0.1045
## workers_alive 1 18.536 187.05 0.9873 0.3204
## block 8 26.710 181.22 9.1619 0.3288
ov1 <- update(ov.pois, .~. -block)
drop1(ov1, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 26.711 181.22
## fungicide 1 26.836 179.35 0.1252 0.72345
## crithidia 1 28.079 180.59 1.3682 0.24212
## avg.pol 1 32.343 184.86 5.6325 0.01763 *
## workers_alive 1 28.872 181.39 2.1614 0.14151
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ov2 <- update(ov1, .~. -workers_alive)
drop1(ov2, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide + crithidia + avg.pol
## Df Deviance AIC LRT Pr(>Chi)
## <none> 28.872 181.39
## fungicide 1 29.007 179.52 0.1348 0.7135
## crithidia 1 29.367 179.88 0.4948 0.4818
## avg.pol 1 48.597 199.11 19.7250 8.942e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(ov2)
## Analysis of Deviance Table (Type II tests)
##
## Response: days_first_ov
## LR Chisq Df Pr(>Chisq)
## fungicide 0.1348 1 0.7135
## crithidia 0.4948 1 0.4818
## avg.pol 19.7250 1 8.942e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ov2)
##
## Call:
## glm(formula = days_first_ov ~ fungicide + crithidia + avg.pol,
## family = "poisson", data = duration)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.71430 -0.66809 -0.00478 0.51505 1.84060
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.95558 0.14685 20.126 < 2e-16 ***
## fungicideTRUE -0.03751 0.10214 -0.367 0.713
## crithidiaTRUE -0.07313 0.10398 -0.703 0.482
## avg.pol -1.12781 0.26254 -4.296 1.74e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 48.911 on 33 degrees of freedom
## Residual deviance: 28.872 on 30 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 181.39
##
## Number of Fisher Scoring iterations: 4
plot(duration$treatment, duration$days_first_ov)
duration$fungicide <- as.factor(duration$fungicide)
duration$crithidia <- as.factor(duration$crithidia)
dm1 <- glm.nb(days_active ~ fungicide + crithidia + avg.pol, data = duration)
## 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(dm1, test = "Chisq")
## Single term deletions
##
## Model:
## days_active ~ fungicide + crithidia + avg.pol
## Df Deviance AIC LRT Pr(>Chi)
## <none> 12.071 217.60
## fungicide 1 12.071 215.60 0.0000 0.99599
## crithidia 1 12.436 215.96 0.3644 0.54606
## avg.pol 1 16.625 220.15 4.5536 0.03285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dm1)
##
## Call:
## glm.nb(formula = days_active ~ fungicide + crithidia + avg.pol,
## data = duration, init.theta = 2513383.999, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2556 -0.2870 0.1041 0.2764 1.1105
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.898150 0.074097 52.608 <2e-16 ***
## fungicideTRUE 0.000255 0.050792 0.005 0.9960
## crithidiaTRUE 0.031677 0.052474 0.604 0.5461
## avg.pol -0.260497 0.122924 -2.119 0.0341 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2513384) family taken to be 1)
##
## Null deviance: 18.284 on 34 degrees of freedom
## Residual deviance: 12.071 on 31 degrees of freedom
## AIC: 219.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2513384
## Std. Err.: 66918604
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -209.598
Anova(dm1)
## Analysis of Deviance Table (Type II tests)
##
## Response: days_active
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0000 1 0.99599
## crithidia 0.3644 1 0.54606
## avg.pol 4.5536 1 0.03285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
brood$fungicide <- as.factor(brood$fungicide)
brood$crithidia <- as.factor(brood$crithidia)
plot(brood$treatment, brood$brood_cells)
brood.mod <- glm.nb(brood_cells ~ fungicide*crithidia + block + workers_alive + duration + avg_pollen, data = brood)
## Warning in glm.nb(brood_cells ~ fungicide * crithidia + block + workers_alive +
## : alternation limit reached
Anova(brood.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 8.519 1 0.003515 **
## crithidia 0.407 1 0.523593
## block 79.894 8 5.134e-14 ***
## workers_alive 33.756 1 6.247e-09 ***
## duration 0.255 1 0.613372
## avg_pollen 17.581 1 2.754e-05 ***
## fungicide:crithidia 0.339 1 0.560394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(brood.mod, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide * crithidia + block + workers_alive +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 47.956 242.33
## block 8 127.850 306.22 79.894 5.134e-14 ***
## workers_alive 1 81.712 274.08 33.756 6.247e-09 ***
## duration 1 48.211 240.58 0.255 0.6134
## avg_pollen 1 65.537 257.91 17.581 2.754e-05 ***
## fungicide:crithidia 1 48.295 240.66 0.339 0.5604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
brood.mod <- glm(brood_cells ~ fungicide + crithidia + block + workers_alive + duration + avg_pollen + mean.spores, data = brood, family = "poisson")
Anova(brood.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 16.953 1 3.832e-05 ***
## crithidia 0.007 1 0.9316
## block 137.115 8 < 2.2e-16 ***
## workers_alive 47.044 1 6.940e-12 ***
## duration 0.052 1 0.8190
## avg_pollen 36.564 1 1.477e-09 ***
## mean.spores 0.513 1 0.4738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(brood.mod)
##
## Call:
## glm(formula = brood_cells ~ fungicide + crithidia + block + workers_alive +
## duration + avg_pollen + mean.spores, family = "poisson",
## data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2302 -1.4477 -0.1618 0.8937 2.5638
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2150353 0.8827013 1.376 0.169
## fungicideTRUE 0.3099288 0.0760404 4.076 4.58e-05 ***
## crithidiaTRUE 0.0092355 0.1075903 0.086 0.932
## block4 -0.7829479 0.1775619 -4.409 1.04e-05 ***
## block6 -1.5347587 0.3045493 -5.039 4.67e-07 ***
## block7 -0.7636176 0.1932362 -3.952 7.76e-05 ***
## block8 -1.1949215 0.1729063 -6.911 4.82e-12 ***
## block9 -0.8424385 0.1526241 -5.520 3.40e-08 ***
## block10 0.0008464 0.1650699 0.005 0.996
## block11 -1.3006112 0.2350764 -5.533 3.15e-08 ***
## block12 -0.3592129 0.1740248 -2.064 0.039 *
## workers_alive 0.3778677 0.0561759 6.727 1.74e-11 ***
## duration -0.0039645 0.0173020 -0.229 0.819
## avg_pollen 2.2288060 0.3730142 5.975 2.30e-09 ***
## mean.spores 0.0060080 0.0083686 0.718 0.473
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 802.894 on 34 degrees of freedom
## Residual deviance: 74.665 on 20 degrees of freedom
## AIC: 248.75
##
## Number of Fisher Scoring iterations: 5
brood.mod <- glm.nb(brood_cells ~ fungicide + crithidia + block + workers_alive + duration + avg_pollen, data = brood)
## Warning in glm.nb(brood_cells ~ fungicide + crithidia + block + workers_alive +
## : alternation limit reached
Anova(brood.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 8.647 1 0.003276 **
## crithidia 0.410 1 0.521935
## block 80.619 8 3.670e-14 ***
## workers_alive 34.070 1 5.316e-09 ***
## duration 0.150 1 0.698608
## avg_pollen 18.242 1 1.946e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(brood.mod)
##
## Call:
## glm.nb(formula = brood_cells ~ fungicide + crithidia + block +
## workers_alive + duration + avg_pollen, data = brood, init.theta = 27.48429804,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.92142 -0.95339 -0.06455 0.52130 2.04677
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.169617 1.175321 0.995 0.319665
## fungicideTRUE 0.335872 0.113681 2.955 0.003132 **
## crithidiaTRUE 0.076905 0.119019 0.646 0.518177
## block4 -0.935350 0.270173 -3.462 0.000536 ***
## block6 -1.475461 0.348371 -4.235 2.28e-05 ***
## block7 -0.791329 0.246527 -3.210 0.001328 **
## block8 -1.312512 0.242672 -5.409 6.35e-08 ***
## block9 -0.920666 0.214360 -4.295 1.75e-05 ***
## block10 0.028614 0.242710 0.118 0.906151
## block11 -1.394267 0.291455 -4.784 1.72e-06 ***
## block12 -0.379938 0.249180 -1.525 0.127321
## workers_alive 0.437542 0.077088 5.676 1.38e-08 ***
## duration -0.008819 0.022744 -0.388 0.698211
## avg_pollen 2.334388 0.552210 4.227 2.36e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(27.4843) family taken to be 1)
##
## Null deviance: 473.432 on 34 degrees of freedom
## Residual deviance: 48.711 on 21 degrees of freedom
## AIC: 242.66
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 27.5
## Std. Err.: 15.6
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -212.662
drop1(brood.mod, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + crithidia + block + workers_alive +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 48.711 240.66
## fungicide 1 57.358 247.31 8.647 0.003276 **
## crithidia 1 49.121 239.07 0.410 0.521935
## block 8 129.330 305.28 80.619 3.670e-14 ***
## workers_alive 1 82.781 272.73 34.070 5.316e-09 ***
## duration 1 48.861 238.81 0.150 0.698608
## avg_pollen 1 66.953 256.90 18.242 1.946e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bm1 <- update(brood.mod, .~. -duration)
## Warning in glm.nb(formula = brood_cells ~ fungicide + crithidia + block + :
## alternation limit reached
drop1(bm1, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + crithidia + block + workers_alive +
## avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.321 238.81
## fungicide 1 58.004 245.49 8.683 0.003212 **
## crithidia 1 49.654 237.14 0.334 0.563447
## block 8 138.670 312.16 89.349 6.303e-16 ***
## workers_alive 1 89.466 276.95 40.145 2.358e-10 ***
## avg_pollen 1 68.352 255.84 19.032 1.286e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(bm1)
##
## Call:
## glm.nb(formula = brood_cells ~ fungicide + crithidia + block +
## workers_alive + avg_pollen, data = brood, init.theta = 28.46029319,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9602 -0.9832 -0.1007 0.5530 1.9290
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.733699 0.307740 2.384 0.017119 *
## fungicideTRUE 0.332652 0.112427 2.959 0.003088 **
## crithidiaTRUE 0.067304 0.115745 0.581 0.560913
## block4 -0.904334 0.260709 -3.469 0.000523 ***
## block6 -1.503876 0.338113 -4.448 8.67e-06 ***
## block7 -0.751666 0.226423 -3.320 0.000901 ***
## block8 -1.287538 0.236208 -5.451 5.01e-08 ***
## block9 -0.897831 0.207652 -4.324 1.53e-05 ***
## block10 -0.001076 0.222935 -0.005 0.996150
## block11 -1.389201 0.289248 -4.803 1.56e-06 ***
## block12 -0.360657 0.242696 -1.486 0.137266
## workers_alive 0.445713 0.072512 6.147 7.91e-10 ***
## avg_pollen 2.353886 0.544220 4.325 1.52e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(28.4603) family taken to be 1)
##
## Null deviance: 479.581 on 34 degrees of freedom
## Residual deviance: 49.321 on 22 degrees of freedom
## AIC: 240.81
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 28.5
## Std. Err.: 16.4
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -212.81
bm2 <- update(bm1, .~. -mean.spores)
drop1(bm2, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + crithidia + block + workers_alive +
## avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.321 238.81
## fungicide 1 58.004 245.49 8.683 0.003212 **
## crithidia 1 49.655 237.14 0.334 0.563445
## block 8 138.671 312.16 89.350 6.301e-16 ***
## workers_alive 1 89.466 276.95 40.145 2.358e-10 ***
## avg_pollen 1 68.353 255.84 19.032 1.286e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bm3 <- update(bm2, .~. -crithidia)
## Warning in glm.nb(formula = brood_cells ~ fungicide + block + workers_alive + :
## alternation limit reached
drop1(bm3, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + block + workers_alive + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.219 237.14
## fungicide 1 57.672 243.59 8.453 0.003645 **
## block 8 137.453 309.38 88.234 1.061e-15 ***
## workers_alive 1 88.797 274.72 39.577 3.153e-10 ***
## avg_pollen 1 67.703 253.62 18.484 1.713e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(bm2, bm3)
## df AIC
## bm2 14 240.8095
## bm3 13 239.1408
Anova(bm2)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 8.683 1 0.003212 **
## crithidia 0.334 1 0.563445
## block 89.350 8 6.301e-16 ***
## workers_alive 40.145 1 2.358e-10 ***
## avg_pollen 19.032 1 1.286e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(bm3)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 8.453 1 0.003645 **
## block 88.234 8 1.061e-15 ***
## workers_alive 39.577 1 3.153e-10 ***
## avg_pollen 18.484 1 1.713e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(bm2)
##
## Call:
## glm.nb(formula = brood_cells ~ fungicide + crithidia + block +
## workers_alive + avg_pollen, data = brood, init.theta = 28.46069118,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9602 -0.9832 -0.1007 0.5530 1.9290
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.733701 0.307740 2.384 0.017118 *
## fungicideTRUE 0.332653 0.112426 2.959 0.003088 **
## crithidiaTRUE 0.067304 0.115745 0.581 0.560912
## block4 -0.904332 0.260709 -3.469 0.000523 ***
## block6 -1.503876 0.338113 -4.448 8.67e-06 ***
## block7 -0.751665 0.226423 -3.320 0.000901 ***
## block8 -1.287537 0.236207 -5.451 5.01e-08 ***
## block9 -0.897830 0.207652 -4.324 1.53e-05 ***
## block10 -0.001076 0.222935 -0.005 0.996148
## block11 -1.389199 0.289247 -4.803 1.56e-06 ***
## block12 -0.360657 0.242695 -1.486 0.137265
## workers_alive 0.445712 0.072512 6.147 7.91e-10 ***
## avg_pollen 2.353887 0.544218 4.325 1.52e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(28.4607) family taken to be 1)
##
## Null deviance: 479.584 on 34 degrees of freedom
## Residual deviance: 49.321 on 22 degrees of freedom
## AIC: 240.81
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 28.5
## Std. Err.: 16.4
##
## 2 x log-likelihood: -212.81
qqnorm(resid(bm3));qqline(resid(bm3))
broodem <- emmeans(bm2, pairwise ~ fungicide*crithidia, type = "response")
broodem
## $emmeans
## fungicide crithidia response SE df asymp.LCL asymp.UCL
## FALSE FALSE 11.4 1.28 Inf 9.15 14.2
## TRUE FALSE 15.9 1.55 Inf 13.15 19.2
## FALSE TRUE 12.2 1.41 Inf 9.73 15.3
## TRUE TRUE 17.0 1.77 Inf 13.87 20.9
##
## Results are averaged over the levels of: block
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## FALSE FALSE / TRUE FALSE 0.717 0.0806 Inf 1 -2.959 0.0163
## FALSE FALSE / FALSE TRUE 0.935 0.1082 Inf 1 -0.581 0.9377
## FALSE FALSE / TRUE TRUE 0.670 0.1096 Inf 1 -2.446 0.0687
## TRUE FALSE / FALSE TRUE 1.304 0.2075 Inf 1 1.667 0.3412
## TRUE FALSE / TRUE TRUE 0.935 0.1082 Inf 1 -0.581 0.9377
## FALSE TRUE / TRUE TRUE 0.717 0.0806 Inf 1 -2.959 0.0163
##
## Results are averaged over the levels of: block
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
broodem.df <- as.data.frame(broodem$emmeans)
broodem.df
## fungicide crithidia response SE df asymp.LCL asymp.UCL
## FALSE FALSE 11.40582 1.283451 Inf 9.148383 14.22030
## TRUE FALSE 15.90728 1.547015 Inf 13.146640 19.24762
## FALSE TRUE 12.19990 1.409618 Inf 9.727599 15.30055
## TRUE TRUE 17.01475 1.772253 Inf 13.872815 20.86828
##
## Results are averaged over the levels of: block
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
broodem.df$treatment <- c(1, 2, 4, 3)
broodem.df$treatment <-as.factor(broodem.df$treatment)
broodcld <- cld(object = broodem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
broodcld
## fungicide crithidia response SE df asymp.LCL asymp.UCL .group
## FALSE FALSE 11.4 1.28 Inf 8.62 15.1 ab
## FALSE TRUE 12.2 1.41 Inf 9.15 16.3 a c
## TRUE FALSE 15.9 1.55 Inf 12.49 20.3 cd
## TRUE TRUE 17.0 1.77 Inf 13.13 22.1 b d
##
## Results are averaged over the levels of: block
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
brood_sum <- brood %>%
group_by(fungicide) %>%
summarise(mb = mean(brood_cells),
nb = length(brood_cells),
sdb = sd(brood_cells)) %>%
mutate(seb = (sdb/sqrt(nb)))
brood_sum
## # A tibble: 2 × 5
## fungicide mb nb sdb seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 FALSE 25.4 17 23.6 5.72
## 2 TRUE 25.7 18 22.8 5.37
ggplot(data = broodem.df, aes(x = treatment, y = response, fill = treatment)) +
geom_col_pattern(
aes(pattern_density = treatment),
pattern = "stripe", # Set a common pattern type, but differentiate density
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, 25)) +
geom_errorbar(aes(ymin = response - SE, ymax = response + SE), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Brood Cells") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 20,
label = "P < 0.001",
size = 7
) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_density_manual(values = c("none", "none", "stripe", "none")) +
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none") +
geom_segment(x = 2, xend = 3, y = 21, yend = 21,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 20.5, yend = 21.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 20.5, yend = 21.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 4, y = 23, yend = 23,
lineend = "round", linejoin = "round") +
geom_segment(x =1, xend = 1, y = 22.5, yend = 23.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 22.5, yend = 23.5,
lineend = "round", linejoin = "round") +
geom_text(x = 2.5, y = 21.5, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 2.5, y = 23.5, label = "b", size = 6, vjust = -0.5)
ggplot(data = broodem.df, aes(x = treatment, y = response, fill = treatment)) +
geom_col_pattern(
aes(pattern_density = treatment),
pattern = "stripe", # Set a common pattern type, but differentiate density
pattern_density = c(0, 0, 0, 0.4), # Add density for the third column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
coord_cartesian(ylim = c(0, 25)) +
geom_errorbar(aes(ymin = response - SE, ymax = response + SE), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Brood Cells") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 20,
label = "P < 0.001",
size = 7
) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_density_manual(values = c(0, 0, 0, 0.4)) +
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none") +
geom_segment(x = 2, xend = 3, y = 21, yend = 21,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 20.5, yend = 21.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 20.5, yend = 21.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 4, y = 23, yend = 23,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 22.5, yend = 23.5,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 22.5, yend = 23.5,
lineend = "round", linejoin = "round") +
geom_text(x = 2.5, y = 21.5, label = "b", size = 6, vjust = -0.5) +
geom_text(x = 2.5, y = 23.5, label = "a", size = 6, vjust = -0.5)
#live Pupae
plot(brood$treatment, brood$live_pupae)
livepup.mod.int <- glm(live_pupae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
drop1(livepup.mod.int, test = "Chisq")
## Single term deletions
##
## Model:
## live_pupae ~ fungicide * crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 34.640 149.83
## workers_alive 1 45.775 158.96 11.1345 0.0008474 ***
## block 8 45.225 144.41 10.5851 0.2263347
## duration 1 38.178 151.36 3.5381 0.0599739 .
## avg_pollen 1 58.122 171.31 23.4822 1.261e-06 ***
## fungicide:crithidia 1 34.855 148.04 0.2143 0.6434311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
livepup.mod <- glm.nb(live_pupae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen, data = brood)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(live_pupae ~ fungicide * crithidia + workers_alive + block +
## : alternation limit reached
Anova(livepup.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.3066 1 0.5797439
## crithidia 2.1103 1 0.1463146
## workers_alive 11.1319 1 0.0008486 ***
## block 10.5838 8 0.2264115
## duration 3.5370 1 0.0600122 .
## avg_pollen 23.4723 1 1.267e-06 ***
## fungicide:crithidia 0.2135 1 0.6440189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(livepup.mod, test = "Chisq")
## Single term deletions
##
## Model:
## live_pupae ~ fungicide * crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 34.633 149.83
## workers_alive 1 45.765 158.96 11.1319 0.0008486 ***
## block 8 45.217 144.41 10.5838 0.2264115
## duration 1 38.170 151.37 3.5370 0.0600122 .
## avg_pollen 1 58.105 171.30 23.4723 1.267e-06 ***
## fungicide:crithidia 1 34.847 148.04 0.2135 0.6440189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
livepup.mod <- glm(live_pupae ~ fungicide + crithidia ++ workers_alive + mean.spores + block + duration + avg_pollen, data = brood, family = "poisson")
livepup.mod.nb <- glm.nb(live_pupae ~ fungicide + crithidia + mean.spores + workers_alive + block + duration + avg_pollen, data = brood)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(live_pupae ~ fungicide + crithidia + mean.spores +
## workers_alive + : alternation limit reached
summary(livepup.mod)
##
## Call:
## glm(formula = live_pupae ~ fungicide + crithidia + +workers_alive +
## mean.spores + block + duration + avg_pollen, family = "poisson",
## data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7907 -0.9821 -0.2128 0.4539 2.3285
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.91690 2.85862 -2.070 0.03847 *
## fungicideTRUE -0.09809 0.17411 -0.563 0.57317
## crithidiaTRUE -0.50947 0.25780 -1.976 0.04813 *
## workers_alive 0.51960 0.15991 3.249 0.00116 **
## mean.spores 0.02791 0.02024 1.379 0.16803
## block4 -0.56476 0.43995 -1.284 0.19925
## block6 -2.09899 1.05468 -1.990 0.04657 *
## block7 -0.45504 0.46780 -0.973 0.33069
## block8 -0.58431 0.41759 -1.399 0.16174
## block9 -0.27289 0.39305 -0.694 0.48751
## block10 -0.64187 0.46118 -1.392 0.16399
## block11 -0.73389 0.52754 -1.391 0.16418
## block12 -0.46293 0.49041 -0.944 0.34519
## duration 0.08893 0.05411 1.643 0.10029
## avg_pollen 4.01816 0.87986 4.567 4.95e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 252.740 on 34 degrees of freedom
## Residual deviance: 32.961 on 20 degrees of freedom
## AIC: 148.15
##
## Number of Fisher Scoring iterations: 6
qqnorm(resid(livepup.mod));qqline(resid(livepup.mod))
plot(livepup.mod)
AIC(livepup.mod, livepup.mod.nb)
## df AIC
## livepup.mod 15 148.1473
## livepup.mod.nb 16 150.1516
anova(livepup.mod, livepup.mod.nb, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: live_pupae ~ fungicide + crithidia + +workers_alive + mean.spores +
## block + duration + avg_pollen
## Model 2: live_pupae ~ fungicide + crithidia + mean.spores + workers_alive +
## block + duration + avg_pollen
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 20 32.961
## 2 20 32.954 0 0.007183
drop1(livepup.mod, test = "Chisq")
## Single term deletions
##
## Model:
## live_pupae ~ fungicide + crithidia + +workers_alive + mean.spores +
## block + duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 32.961 148.15
## fungicide 1 33.278 146.47 0.3172 0.5732658
## crithidia 1 36.964 150.15 4.0031 0.0454165 *
## workers_alive 1 44.568 157.75 11.6069 0.0006571 ***
## mean.spores 1 34.855 148.04 1.8933 0.1688305
## block 8 41.853 141.04 8.8920 0.3514918
## duration 1 35.811 149.00 2.8496 0.0913963 .
## avg_pollen 1 56.150 169.34 23.1885 1.469e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lp1 <- update(livepup.mod, .~. -block)
drop1(lp1, test = "Chisq")
## Single term deletions
##
## Model:
## live_pupae ~ fungicide + crithidia + workers_alive + mean.spores +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 41.853 141.04
## fungicide 1 42.264 139.45 0.411 0.52167
## crithidia 1 46.072 143.26 4.219 0.03997 *
## workers_alive 1 64.445 161.63 22.592 2.003e-06 ***
## mean.spores 1 45.812 143.00 3.959 0.04661 *
## duration 1 47.936 145.12 6.083 0.01365 *
## avg_pollen 1 102.098 199.28 60.244 8.378e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lp1)
##
## Call:
## glm(formula = live_pupae ~ fungicide + crithidia + workers_alive +
## mean.spores + duration + avg_pollen, family = "poisson",
## data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9620 -1.0652 -0.3355 0.4490 2.6177
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.37928 1.54690 -3.477 0.000506 ***
## fungicideTRUE -0.09883 0.15449 -0.640 0.522374
## crithidiaTRUE -0.44103 0.22170 -1.989 0.046663 *
## workers_alive 0.52980 0.12121 4.371 1.24e-05 ***
## mean.spores 0.02838 0.01383 2.052 0.040164 *
## duration 0.06793 0.02712 2.505 0.012252 *
## avg_pollen 3.65015 0.50856 7.177 7.10e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 252.740 on 34 degrees of freedom
## Residual deviance: 41.853 on 28 degrees of freedom
## AIC: 141.04
##
## Number of Fisher Scoring iterations: 6
Anova(lp1)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.411 1 0.52167
## crithidia 4.219 1 0.03997 *
## workers_alive 22.592 1 2.003e-06 ***
## mean.spores 3.959 1 0.04661 *
## duration 6.083 1 0.01365 *
## avg_pollen 60.244 1 8.378e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(lp1));qqline(resid(lp1))
anova(livepup.mod, lp1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: live_pupae ~ fungicide + crithidia + +workers_alive + mean.spores +
## block + duration + avg_pollen
## Model 2: live_pupae ~ fungicide + crithidia + workers_alive + mean.spores +
## duration + avg_pollen
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 20 32.961
## 2 28 41.853 -8 -8.892 0.3515
AIC(livepup.mod, lp1)
## df AIC
## livepup.mod 15 148.1473
## lp1 7 141.0393
be <- emmeans(lp1, "crithidia")
pairs(be)
## contrast estimate SE df z.ratio p.value
## FALSE - TRUE 0.441 0.222 Inf 1.989 0.0467
##
## Results are averaged over the levels of: fungicide
## Results are given on the log (not the response) scale.
broodem <- emmeans(lp1, pairwise ~ crithidia, type = "response")
broodem
## $emmeans
## crithidia rate SE df asymp.LCL asymp.UCL
## FALSE 3.06 0.505 Inf 2.21 4.23
## TRUE 1.97 0.378 Inf 1.35 2.87
##
## Results are averaged over the levels of: fungicide
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## FALSE / TRUE 1.55 0.345 Inf 1 1.989 0.0467
##
## Results are averaged over the levels of: fungicide
## Tests are performed on the log scale
broodem <- emmeans(lp1, pairwise ~ crithidia*fungicide, type = "response")
broodcld <- cld(object = broodem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
broodcld
## crithidia fungicide rate SE df asymp.LCL asymp.UCL .group
## TRUE TRUE 1.87 0.371 Inf 1.14 3.07 a
## TRUE FALSE 2.07 0.445 Inf 1.21 3.54 a
## FALSE TRUE 2.91 0.546 Inf 1.83 4.64 a
## FALSE FALSE 3.21 0.568 Inf 2.07 4.99 a
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
livepup_sum <- brood %>%
group_by(treatment) %>%
summarise(mb = mean(live_pupae),
nb = length(live_pupae),
sdb = sd(live_pupae)) %>%
mutate(seb = (sdb/sqrt(nb)))
livepup_sum
## # A tibble: 4 × 5
## treatment mb nb sdb seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 8.22 9 7.58 2.53
## 2 2 5.89 9 7.54 2.51
## 3 3 2.89 9 3.30 1.10
## 4 4 3.88 8 5.51 1.95
ggplot(data = livepup_sum, aes(x = treatment, y = mb, fill = treatment, pattern = 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, 13)) +
geom_errorbar(aes(ymin = mb - seb, ymax = mb + seb), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Live Pupae") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 4,
y = 12.5,
label = "P = 0.04",
size = 7
) +
scale_x_discrete(labels = custom_labels) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_manual(values = c("none", "none", "stripe", "none")) +
geom_segment(x = 1, xend = 2, y = 12, yend = 12,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 11.8, yend = 12.2,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 11.8, yend = 12.2,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 4, y = 7, yend = 7,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 6.8, yend = 7.2,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 6.8, yend = 7.2,
lineend = "round", linejoin = "round") +
geom_text(x = 1.5, y = 12, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 3.5, y = 7, label = "b", size = 6, vjust = -0.5) +
theme(legend.position = "none")
plot(brood$treatment, brood$live_larvae)
livelar.mod <- glm(live_larvae ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen + mean.spores, data = brood, family = "poisson")
summary(livelar.mod) #overdisp
##
## Call:
## glm(formula = live_larvae ~ fungicide + crithidia + workers_alive +
## block + duration + avg_pollen + mean.spores, family = "poisson",
## data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5438 -1.6464 -0.7831 1.0794 2.8360
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.730899 1.184374 1.461 0.143893
## fungicideTRUE 0.205617 0.100781 2.040 0.041327 *
## crithidiaTRUE 0.104241 0.145114 0.718 0.472547
## workers_alive 0.354239 0.074593 4.749 2.04e-06 ***
## block4 -0.323313 0.239673 -1.349 0.177345
## block6 -1.857634 0.608450 -3.053 0.002265 **
## block7 -0.592166 0.277349 -2.135 0.032753 *
## block8 -0.920914 0.244738 -3.763 0.000168 ***
## block9 -0.656017 0.227363 -2.885 0.003910 **
## block10 0.701209 0.238182 2.944 0.003240 **
## block11 -1.317346 0.386150 -3.411 0.000646 ***
## block12 0.329240 0.237823 1.384 0.166239
## duration -0.034907 0.023584 -1.480 0.138844
## avg_pollen 2.215138 0.486168 4.556 5.21e-06 ***
## mean.spores 0.005285 0.011887 0.445 0.656629
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 637.297 on 34 degrees of freedom
## Residual deviance: 91.016 on 20 degrees of freedom
## AIC: 232.97
##
## Number of Fisher Scoring iterations: 5
livelar.mod.int <- glm(live_larvae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen + mean.spores, data = brood, family = "poisson")
summary(livelar.mod.int) #overdisp
##
## Call:
## glm(formula = live_larvae ~ fungicide * crithidia + workers_alive +
## block + duration + avg_pollen + mean.spores, family = "poisson",
## data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7192 -1.6768 -0.3928 1.0847 2.6161
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.978218 1.181438 1.674 0.094049 .
## fungicideTRUE 0.108249 0.125521 0.862 0.388466
## crithidiaTRUE -0.040472 0.182830 -0.221 0.824811
## workers_alive 0.365091 0.075762 4.819 1.44e-06 ***
## block4 -0.334658 0.241567 -1.385 0.165941
## block6 -1.890510 0.608355 -3.108 0.001886 **
## block7 -0.563300 0.273068 -2.063 0.039126 *
## block8 -0.934116 0.243585 -3.835 0.000126 ***
## block9 -0.722523 0.234400 -3.082 0.002053 **
## block10 0.669531 0.235042 2.849 0.004392 **
## block11 -1.402393 0.391280 -3.584 0.000338 ***
## block12 0.300783 0.240316 1.252 0.210710
## duration -0.038075 0.023320 -1.633 0.102534
## avg_pollen 2.068682 0.498238 4.152 3.30e-05 ***
## mean.spores 0.004148 0.011662 0.356 0.722103
## fungicideTRUE:crithidiaTRUE 0.266514 0.204306 1.304 0.192069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 637.297 on 34 degrees of freedom
## Residual deviance: 89.303 on 19 degrees of freedom
## AIC: 233.25
##
## Number of Fisher Scoring iterations: 5
livelar.mod.nb.int <- glm.nb(live_larvae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen + mean.spores, data = brood)
## Warning in glm.nb(live_larvae ~ fungicide * crithidia + workers_alive + :
## alternation limit reached
summary(livelar.mod.nb.int)
##
## Call:
## glm.nb(formula = live_larvae ~ fungicide * crithidia + workers_alive +
## block + duration + avg_pollen + mean.spores, data = brood,
## init.theta = 7.784770105, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9108 -1.1908 -0.3695 0.5482 2.0425
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.244101 1.969836 0.632 0.527664
## fungicideTRUE 0.121084 0.239238 0.506 0.612769
## crithidiaTRUE -0.051580 0.342593 -0.151 0.880323
## workers_alive 0.440942 0.128031 3.444 0.000573 ***
## block4 -0.670663 0.465138 -1.442 0.149342
## block6 -1.766032 0.697272 -2.533 0.011316 *
## block7 -0.626983 0.451785 -1.388 0.165201
## block8 -1.226428 0.415884 -2.949 0.003188 **
## block9 -0.757775 0.376687 -2.012 0.044253 *
## block10 0.641467 0.418712 1.532 0.125522
## block11 -1.514036 0.529228 -2.861 0.004225 **
## block12 0.386750 0.410957 0.941 0.346656
## duration -0.033864 0.038102 -0.889 0.374124
## avg_pollen 2.731714 0.925783 2.951 0.003170 **
## mean.spores -0.002224 0.019461 -0.114 0.909003
## fungicideTRUE:crithidiaTRUE 0.336126 0.382670 0.878 0.379744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7.7848) family taken to be 1)
##
## Null deviance: 276.05 on 34 degrees of freedom
## Residual deviance: 48.72 on 19 degrees of freedom
## AIC: 223.73
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7.78
## Std. Err.: 4.31
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -189.734
drop1(livelar.mod.nb.int, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide * crithidia + workers_alive + block +
## duration + avg_pollen + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 48.720 221.73
## workers_alive 1 60.463 231.48 11.743 0.0006107 ***
## block 8 96.392 253.41 47.672 1.141e-07 ***
## duration 1 49.414 220.43 0.694 0.4048798
## avg_pollen 1 57.019 228.03 8.299 0.0039673 **
## mean.spores 1 48.732 219.75 0.012 0.9122416
## fungicide:crithidia 1 49.475 220.49 0.755 0.3848937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
livelar.mod.nb <- glm.nb(live_larvae ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen + mean.spores, data = brood) #start with this one
## Warning in glm.nb(live_larvae ~ fungicide + crithidia + workers_alive + :
## alternation limit reached
summary(livelar.mod.nb)
##
## Call:
## glm.nb(formula = live_larvae ~ fungicide + crithidia + workers_alive +
## block + duration + avg_pollen + mean.spores, data = brood,
## init.theta = 7.564716463, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7385 -1.1436 -0.5959 0.5551 2.1529
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.879108 1.956613 0.449 0.653215
## fungicideTRUE 0.244326 0.193382 1.263 0.206433
## crithidiaTRUE 0.100646 0.287746 0.350 0.726509
## workers_alive 0.434394 0.128058 3.392 0.000693 ***
## block4 -0.624050 0.466291 -1.338 0.180790
## block6 -1.782600 0.702001 -2.539 0.011107 *
## block7 -0.624096 0.458462 -1.361 0.173424
## block8 -1.181552 0.415829 -2.841 0.004491 **
## block9 -0.716770 0.372659 -1.923 0.054431 .
## block10 0.676809 0.423762 1.597 0.110234
## block11 -1.439430 0.521982 -2.758 0.005822 **
## block12 0.401871 0.411766 0.976 0.329080
## duration -0.027438 0.038099 -0.720 0.471427
## avg_pollen 2.783866 0.931899 2.987 0.002815 **
## mean.spores -0.000826 0.019723 -0.042 0.966593
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7.5647) family taken to be 1)
##
## Null deviance: 272.22 on 34 degrees of freedom
## Residual deviance: 49.01 on 20 degrees of freedom
## AIC: 222.49
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7.56
## Std. Err.: 4.20
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -190.487
drop1(livelar.mod.nb, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen + mean.spores
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.010 220.49
## fungicide 1 50.541 220.02 1.531 0.2159583
## crithidia 1 49.123 218.60 0.114 0.7361879
## workers_alive 1 60.326 229.80 11.317 0.0007682 ***
## block 8 95.413 250.89 46.404 1.992e-07 ***
## duration 1 49.476 218.95 0.466 0.4947161
## avg_pollen 1 57.498 226.98 8.489 0.0035738 **
## mean.spores 1 49.011 218.49 0.002 0.9675214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ll1 <- update(livelar.mod.nb, .~. -mean.spores)
drop1(ll1, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.070 218.49
## fungicide 1 50.610 218.03 1.540 0.2146048
## crithidia 1 49.270 216.69 0.200 0.6547016
## workers_alive 1 60.417 227.84 11.348 0.0007554 ***
## block 8 95.665 249.09 46.596 1.831e-07 ***
## duration 1 49.567 216.99 0.497 0.4806883
## avg_pollen 1 57.629 225.05 8.559 0.0034376 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ll1 <-update(ll1, .~. -duration)
drop1(ll1, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide + crithidia + workers_alive + block +
## avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 48.646 216.98
## fungicide 1 50.098 216.43 1.452 0.2282051
## crithidia 1 48.755 215.09 0.109 0.7415677
## workers_alive 1 63.119 229.45 14.473 0.0001422 ***
## block 8 95.177 247.51 46.532 1.883e-07 ***
## avg_pollen 1 57.480 223.81 8.834 0.0029565 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ll1)
##
## Call:
## glm.nb(formula = live_larvae ~ fungicide + crithidia + workers_alive +
## block + avg_pollen, data = brood, init.theta = 7.179765666,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8114 -1.0573 -0.4253 0.5787 2.1132
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.51697 0.53007 -0.975 0.329425
## fungicideTRUE 0.24082 0.19563 1.231 0.218323
## crithidiaTRUE 0.06821 0.20411 0.334 0.738250
## workers_alive 0.46773 0.12413 3.768 0.000164 ***
## block4 -0.54258 0.45839 -1.184 0.236545
## block6 -1.87282 0.69447 -2.697 0.007002 **
## block7 -0.52011 0.39823 -1.306 0.191529
## block8 -1.12332 0.41192 -2.727 0.006390 **
## block9 -0.63579 0.36575 -1.738 0.082154 .
## block10 0.59992 0.39818 1.507 0.131896
## block11 -1.42840 0.52356 -2.728 0.006367 **
## block12 0.45392 0.41203 1.102 0.270609
## avg_pollen 2.86018 0.93920 3.045 0.002324 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7.1798) family taken to be 1)
##
## Null deviance: 265.290 on 34 degrees of freedom
## Residual deviance: 48.646 on 22 degrees of freedom
## AIC: 218.98
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7.18
## Std. Err.: 3.91
##
## 2 x log-likelihood: -190.978
Anova(ll1)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 1.452 1 0.2282051
## crithidia 0.109 1 0.7415677
## workers_alive 14.473 1 0.0001422 ***
## block 46.532 8 1.883e-07 ***
## avg_pollen 8.834 1 0.0029565 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broodem <- emmeans(ll1, pairwise ~ crithidia, type = "response")
broodcld <- cld(object = broodem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
broodcld
## crithidia response SE df asymp.LCL asymp.UCL .group
## FALSE 6.13 0.973 Inf 4.30 8.74 a
## TRUE 6.56 1.092 Inf 4.52 9.52 a
##
## Results are averaged over the levels of: fungicide, block
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
## Intervals are back-transformed from the log scale
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
livelarv_sum <- brood %>%
group_by(treatment) %>%
summarise(mb = mean(live_larvae),
nb = length(live_larvae),
sdb = sd(live_larvae)) %>%
mutate(seb = (sdb/sqrt(nb)))
livelarv_sum
## # A tibble: 4 × 5
## treatment mb nb sdb seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 20.7 9 19.0 6.34
## 2 2 15.1 9 11.0 3.68
## 3 3 12.8 9 19.5 6.52
## 4 4 10.2 8 13.3 4.71
livelarv_sum$plot <- livelarv_sum$mb + livelarv_sum$seb
plot(brood$treatment, brood$live_larvae)
ggplot(data = livelarv_sum, aes(x = treatment, y = mb, fill = treatment)) +
geom_col_pattern(
aes(pattern = treatment),
pattern_density = c(0, 0, 0, 0.4), # Add density for the fourth column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
coord_cartesian(ylim = c(0, 35)) +
geom_errorbar(aes(ymin = mb - seb, ymax = mb + seb), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Live Pupae") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 34,
label = "P > 0.5",
size = 7
) +
scale_fill_manual(values = c("lightgreen", "lightblue", "grey", "lightblue")) +
scale_pattern_manual(values = c("none", "none", "none", "stripe")) + # Add stripes to the fourth column
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none")
dl.mod.pois <- glm(dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive + block + duration, data = brood, family = "poisson")
summary(dl.mod.pois) #overdisp
##
## Call:
## glm(formula = dead_larvae ~ fungicide + crithidia + avg_pollen +
## workers_alive + block + duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4732 -1.1508 -0.5425 0.3474 3.0203
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -13.899981 5.671894 -2.451 0.0143 *
## fungicideTRUE 0.193663 0.338075 0.573 0.5668
## crithidiaTRUE -0.339239 0.352774 -0.962 0.3362
## avg_pollen 4.550107 1.937007 2.349 0.0188 *
## workers_alive 0.311674 0.231644 1.345 0.1785
## block4 0.007266 0.850489 0.009 0.9932
## block6 -0.372933 0.783530 -0.476 0.6341
## block7 1.056425 0.716866 1.474 0.1406
## block8 -0.639657 0.867599 -0.737 0.4610
## block9 0.917755 0.749365 1.225 0.2207
## block10 -1.948176 0.923300 -2.110 0.0349 *
## block11 -1.139936 1.114845 -1.023 0.3065
## block12 -0.029739 0.806340 -0.037 0.9706
## duration 0.244663 0.105595 2.317 0.0205 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.190 on 34 degrees of freedom
## Residual deviance: 64.033 on 21 degrees of freedom
## AIC: 135.52
##
## Number of Fisher Scoring iterations: 6
dl.mod <- glm.nb(dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive + block + duration, data = brood)
drop1(dl.mod, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive +
## block + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 31.432 124.16
## fungicide 1 31.928 122.66 0.4962 0.48119
## crithidia 1 31.690 122.42 0.2579 0.61158
## avg_pollen 1 32.789 123.52 1.3567 0.24411
## workers_alive 1 32.658 123.39 1.2258 0.26823
## block 8 37.722 114.45 6.2900 0.61478
## duration 1 35.607 126.34 4.1749 0.04103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dl1 <- update(dl.mod, .~. -block)
drop1(dl1, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive +
## duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 32.107 113.94
## fungicide 1 32.672 112.50 0.56499 0.4523
## crithidia 1 32.121 111.95 0.01367 0.9069
## avg_pollen 1 34.166 113.99 2.05858 0.1514
## workers_alive 1 33.126 112.95 1.01838 0.3129
## duration 1 32.706 112.53 0.59828 0.4392
dl2 <- update(dl1, .~. -duration)
drop1(dl2, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 32.279 112.53
## fungicide 1 32.885 111.14 0.60569 0.4364
## crithidia 1 32.342 110.59 0.06231 0.8029
## avg_pollen 1 33.817 112.07 1.53714 0.2150
## workers_alive 1 32.857 111.11 0.57746 0.4473
dl3 <- update(dl2, .~. -workers_alive)
drop1(dl3, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 32.135 111.10
## fungicide 1 32.811 109.77 0.6759 0.4110
## crithidia 1 32.138 109.10 0.0028 0.9578
## avg_pollen 1 37.104 114.07 4.9691 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(brood$treatment, brood$dead_larvae)
Anova(dl3)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.6759 1 0.4110
## crithidia 0.0028 1 0.9578
## avg_pollen 4.9691 1 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dl3)
##
## Call:
## glm.nb(formula = dead_larvae ~ fungicide + crithidia + avg_pollen,
## data = brood, init.theta = 0.6377843997, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5607 -1.0379 -0.8436 0.2955 1.4610
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.27071 0.81370 -1.562 0.1184
## fungicideTRUE 0.47337 0.54520 0.868 0.3853
## crithidiaTRUE -0.02963 0.56091 -0.053 0.9579
## avg_pollen 2.74777 1.23052 2.233 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.6378) family taken to be 1)
##
## Null deviance: 37.455 on 34 degrees of freedom
## Residual deviance: 32.135 on 31 degrees of freedom
## AIC: 113.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.638
## Std. Err.: 0.305
##
## 2 x log-likelihood: -103.098
dp.mod.pois <- glm(dead_pupae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood, family = "poisson")
summary(dp.mod.pois)
##
## Call:
## glm(formula = dead_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive + duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.90219 -0.37359 -0.15693 -0.06104 1.56338
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.1769 13.6891 -1.255 0.210
## fungicideTRUE 0.5435 1.2911 0.421 0.674
## crithidiaTRUE 0.1942 1.4213 0.137 0.891
## avg_pollen 12.7579 7.7976 1.636 0.102
## workers_alive -1.1850 0.8174 -1.450 0.147
## duration 0.2606 0.2421 1.076 0.282
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17.3524 on 34 degrees of freedom
## Residual deviance: 9.3054 on 29 degrees of freedom
## AIC: 29.305
##
## Number of Fisher Scoring iterations: 7
dp.mod <- glm.nb(dead_pupae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood)
## 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(dp.mod, test = "Chisq")
## Single term deletions
##
## Model:
## dead_pupae ~ fungicide + crithidia + avg_pollen + workers_alive +
## duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 9.3051 29.306
## fungicide 1 9.4896 27.490 0.1845 0.66750
## crithidia 1 9.3237 27.324 0.0186 0.89145
## avg_pollen 1 15.6112 33.612 6.3061 0.01203 *
## workers_alive 1 12.2715 30.272 2.9664 0.08501 .
## duration 1 10.7087 28.709 1.4036 0.23612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dp1 <- update(dp.mod, .~. -duration)
## 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(dp1, test = "Chisq")
## Single term deletions
##
## Model:
## dead_pupae ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 10.709 28.709
## fungicide 1 11.759 27.759 1.0501 0.30549
## crithidia 1 10.781 26.781 0.0721 0.78830
## avg_pollen 1 15.824 31.825 5.1153 0.02372 *
## workers_alive 1 15.425 31.426 4.7165 0.02987 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dp1)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 1.0501 1 0.30549
## crithidia 0.0721 1 0.78830
## avg_pollen 5.1153 1 0.02372 *
## workers_alive 4.7165 1 0.02987 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dp1)
##
## Call:
## glm.nb(formula = dead_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive, data = brood, init.theta = 6567.642691, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0028 -0.3231 -0.1841 -0.1075 1.5299
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1910 1.9295 -1.654 0.0982 .
## fungicideTRUE 1.1184 1.1735 0.953 0.3405
## crithidiaTRUE -0.3406 1.2861 -0.265 0.7911
## avg_pollen 9.2015 5.2366 1.757 0.0789 .
## workers_alive -1.2795 0.7263 -1.762 0.0781 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(6567.643) family taken to be 1)
##
## Null deviance: 17.352 on 34 degrees of freedom
## Residual deviance: 10.709 on 30 degrees of freedom
## AIC: 30.709
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 6568
## Std. Err.: 266783
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -18.709
plot(brood$treatment, brood$dead_pupae)
tl.mod.pois <- glm(total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood, family = "poisson")
summary(tl.mod.pois)
##
## Call:
## glm(formula = total_larvae ~ fungicide + crithidia + avg_pollen +
## workers_alive + duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7993 -2.4311 -0.6997 1.0207 5.9745
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.26957 0.85010 -2.670 0.007591 **
## fungicideTRUE 0.01569 0.08606 0.182 0.855306
## crithidiaTRUE 0.07336 0.08899 0.824 0.409733
## avg_pollen 3.56189 0.28189 12.636 < 2e-16 ***
## workers_alive 0.22926 0.05740 3.994 6.49e-05 ***
## duration 0.05151 0.01550 3.323 0.000891 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 654.02 on 34 degrees of freedom
## Residual deviance: 201.89 on 29 degrees of freedom
## AIC: 332.85
##
## Number of Fisher Scoring iterations: 6
tl.mod <- glm.nb(total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood)
drop1(tl.mod, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive +
## duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.271 238.46
## fungicide 1 47.022 240.21 3.7505 0.05279 .
## crithidia 1 43.384 236.57 0.1131 0.73667
## avg_pollen 1 65.585 258.77 22.3136 2.316e-06 ***
## workers_alive 1 46.312 239.50 3.0415 0.08116 .
## duration 1 43.604 236.79 0.3330 0.56391
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tl.mod1 <- update(tl.mod, .~. -duration)
drop1(tl.mod1, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.853 236.79
## fungicide 1 47.028 238.96 4.1756 0.04101 *
## crithidia 1 43.005 234.94 0.1516 0.69700
## avg_pollen 1 64.511 256.44 21.6581 3.258e-06 ***
## workers_alive 1 45.525 237.46 2.6721 0.10212
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tl1 <- update(tl.mod1, .~. -workers_alive)
drop1(tl1, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.714 237.36
## fungicide 1 45.662 238.31 2.948 0.08596 .
## crithidia 1 42.767 235.41 0.053 0.81726
## avg_pollen 1 89.002 281.65 46.288 1.021e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tl2 <- update(tl1, .~. -crithidia)
drop1(tl2, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.713 235.41
## fungicide 1 45.663 236.36 2.951 0.08584 .
## avg_pollen 1 90.718 281.42 48.005 4.25e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(tl.mod, tl.mod1, tl1)
## df AIC
## tl.mod 7 240.4600
## tl.mod1 6 238.7858
## tl1 5 239.3597
anova(tl.mod1, tl1, test = "Chisq")
## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_larvae
## Model theta Resid. df
## 1 fungicide + crithidia + avg_pollen 1.387082 31
## 2 fungicide + crithidia + avg_pollen + workers_alive 1.543006 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -229.3597
## 2 -226.7858 1 vs 2 1 2.573877 0.1086412
Anova(tl1)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 2.948 1 0.08596 .
## crithidia 0.053 1 0.81726
## avg_pollen 46.288 1 1.021e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(tl1)
##
## Call:
## glm.nb(formula = total_larvae ~ fungicide + crithidia + avg_pollen,
## data = brood, init.theta = 1.387082107, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8569 -1.3615 -0.4205 0.4500 1.4367
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.44219 0.48239 -0.917 0.3593
## fungicideTRUE 0.57342 0.32289 1.776 0.0758 .
## crithidiaTRUE -0.07433 0.33183 -0.224 0.8228
## avg_pollen 5.82756 0.74711 7.800 6.19e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.3871) family taken to be 1)
##
## Null deviance: 90.992 on 34 degrees of freedom
## Residual deviance: 42.714 on 31 degrees of freedom
## AIC: 239.36
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.387
## Std. Err.: 0.463
##
## 2 x log-likelihood: -229.360
qqnorm(resid(tl.mod1));qqline(resid(tl.mod1))
plot(brood$treatment, brood$total_larvae)
tp.mod.pois <- glm(total_pupae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood, family = "poisson")
summary(tp.mod.pois)
##
## Call:
## glm(formula = total_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive + duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5037 -1.1070 -0.2122 0.3784 2.8710
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.87063 1.49720 -3.253 0.001141 **
## fungicideTRUE -0.13470 0.15032 -0.896 0.370213
## crithidiaTRUE -0.15252 0.16055 -0.950 0.342115
## avg_pollen 3.64523 0.49406 7.378 1.61e-13 ***
## workers_alive 0.43424 0.11376 3.817 0.000135 ***
## duration 0.06726 0.02666 2.523 0.011652 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 251.751 on 34 degrees of freedom
## Residual deviance: 49.396 on 29 degrees of freedom
## AIC: 147.77
##
## Number of Fisher Scoring iterations: 6
tp.mod <- glm.nb(total_pupae ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(total_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive + : alternation limit reached
drop1(tp.mod, test = "Chisq")
## Single term deletions
##
## Model:
## total_pupae ~ fungicide + crithidia + avg_pollen + workers_alive +
## duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.372 147.77
## fungicide 1 50.175 146.58 0.804 0.37002
## crithidia 1 50.283 146.69 0.912 0.33970
## avg_pollen 1 112.309 208.71 62.938 2.133e-15 ***
## workers_alive 1 65.909 162.31 16.537 4.771e-05 ***
## duration 1 55.584 151.99 6.213 0.01268 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(tp.mod));qqline(resid(tp.mod))
Anova(tp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.804 1 0.37002
## crithidia 0.912 1 0.33970
## avg_pollen 62.938 1 2.133e-15 ***
## workers_alive 16.537 1 4.771e-05 ***
## duration 6.213 1 0.01268 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(tp.mod)
##
## Call:
## glm.nb(formula = total_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive + duration, data = brood, init.theta = 7247.587331,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5030 -1.1068 -0.2124 0.3779 2.8698
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.87146 1.49850 -3.251 0.001150 **
## fungicideTRUE -0.13451 0.15044 -0.894 0.371256
## crithidiaTRUE -0.15218 0.16065 -0.947 0.343514
## avg_pollen 3.64561 0.49433 7.375 1.64e-13 ***
## workers_alive 0.43428 0.11381 3.816 0.000136 ***
## duration 0.06726 0.02669 2.520 0.011719 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7247.588) family taken to be 1)
##
## Null deviance: 251.560 on 34 degrees of freedom
## Residual deviance: 49.372 on 29 degrees of freedom
## AIC: 149.77
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7248
## Std. Err.: 115777
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -135.774
plot(brood$treatment, brood$total_pupae)
egg.mod.pois <- glm(eggs ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood, family = "poisson")
summary(egg.mod.pois)
##
## Call:
## glm(formula = eggs ~ fungicide + crithidia + avg_pollen + workers_alive +
## duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.3075 -3.2797 -0.5805 1.6062 9.8920
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44892 0.67512 0.665 0.5061
## fungicideTRUE 0.16831 0.07692 2.188 0.0287 *
## crithidiaTRUE 0.00982 0.07987 0.123 0.9022
## avg_pollen 1.62997 0.22926 7.110 1.16e-12 ***
## workers_alive 0.25396 0.04563 5.565 2.62e-08 ***
## duration 0.01561 0.01254 1.245 0.2132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 697.97 on 34 degrees of freedom
## Residual deviance: 440.07 on 29 degrees of freedom
## AIC: 583.88
##
## Number of Fisher Scoring iterations: 5
egg.mod <- glm.nb(eggs ~ fungicide + crithidia + avg_pollen + workers_alive + duration, data = brood)
drop1(egg.mod, test = "Chisq")
## Single term deletions
##
## Model:
## eggs ~ fungicide + crithidia + avg_pollen + workers_alive + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.558 273.31
## fungicide 1 44.110 272.86 1.5522 0.21281
## crithidia 1 42.648 271.40 0.0903 0.76378
## avg_pollen 1 44.550 273.30 1.9923 0.15810
## workers_alive 1 46.650 275.40 4.0923 0.04308 *
## duration 1 42.558 271.31 0.0001 0.99182
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em1 <- update(egg.mod, .~. -duration)
drop1(em1, test = "Chisq")
## Single term deletions
##
## Model:
## eggs ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.558 271.31
## fungicide 1 44.115 270.87 1.5568 0.21213
## crithidia 1 42.651 269.40 0.0925 0.76096
## avg_pollen 1 44.570 271.32 2.0121 0.15605
## workers_alive 1 47.181 273.93 4.6225 0.03155 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em2 <- update(em1, .~. -avg_pollen)
drop1(em2, test = "Chisq")
## Single term deletions
##
## Model:
## eggs ~ fungicide + crithidia + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.085 271.24
## fungicide 1 42.965 270.12 0.8804 0.3481013
## crithidia 1 42.397 269.55 0.3119 0.5765148
## workers_alive 1 57.092 284.25 15.0076 0.0001071 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(em2)
## Analysis of Deviance Table (Type II tests)
##
## Response: eggs
## LR Chisq Df Pr(>Chisq)
## fungicide 0.8804 1 0.3481013
## crithidia 0.3119 1 0.5765148
## workers_alive 15.0076 1 0.0001071 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(em2)
##
## Call:
## glm.nb(formula = eggs ~ fungicide + crithidia + workers_alive,
## data = brood, init.theta = 0.9090852056, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5070 -1.0515 -0.1384 0.3260 1.6643
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7034 0.6194 1.136 0.256
## fungicideTRUE 0.3526 0.3725 0.947 0.344
## crithidiaTRUE -0.2105 0.3952 -0.533 0.594
## workers_alive 0.5683 0.1283 4.431 9.4e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9091) family taken to be 1)
##
## Null deviance: 58.057 on 34 degrees of freedom
## Residual deviance: 42.085 on 31 degrees of freedom
## AIC: 273.24
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.909
## Std. Err.: 0.252
##
## 2 x log-likelihood: -263.241
plot(brood$treatment, brood$eggs)
hp.mod <- glm.nb(honey_pots ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
hp.mod.pois <- glm.nb(honey_pots ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
summary(hp.mod.pois)
##
## Call:
## glm.nb(formula = honey_pots ~ fungicide + crithidia + avg_pollen +
## workers_alive, data = brood, init.theta = 19.66709073, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9681 -0.8476 -0.2812 0.2471 2.1264
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4311 0.4339 -0.994 0.32044
## fungicideTRUE 0.4625 0.2213 2.090 0.03666 *
## crithidiaTRUE -0.2824 0.2340 -1.207 0.22752
## avg_pollen 1.6542 0.6204 2.667 0.00766 **
## workers_alive 0.1530 0.1137 1.346 0.17839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(19.6671) family taken to be 1)
##
## Null deviance: 66.900 on 34 degrees of freedom
## Residual deviance: 34.352 on 30 degrees of freedom
## AIC: 138.04
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 19.7
## Std. Err.: 31.1
##
## 2 x log-likelihood: -126.035
anova(hp.mod, hp.mod.pois, test = "Chisq")
## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pots
## Model theta Resid. df
## 1 fungicide + crithidia + avg_pollen + workers_alive 19.66709 30
## 2 fungicide + crithidia + avg_pollen + workers_alive 19.66709 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -126.0351
## 2 -126.0351 1 vs 2 0 0 1
AIC(hp.mod, hp.mod.pois)
## df AIC
## hp.mod 6 138.0351
## hp.mod.pois 6 138.0351
drop1(hp.mod.pois, test = "Chisq")
## Single term deletions
##
## Model:
## honey_pots ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 34.352 136.03
## fungicide 1 38.757 138.44 4.4048 0.035839 *
## crithidia 1 35.828 135.51 1.4753 0.224517
## avg_pollen 1 41.423 141.11 7.0710 0.007834 **
## workers_alive 1 36.234 135.92 1.8821 0.170097
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hp1 <- update(hp.mod.pois, .~. -workers_alive)
drop1(hp1, test = "Chisq")
## Single term deletions
##
## Model:
## honey_pots ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 35.690 135.91
## fungicide 1 39.550 137.77 3.8596 0.04946 *
## crithidia 1 37.849 136.07 2.1588 0.14176
## avg_pollen 1 56.441 154.66 20.7507 5.231e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(hp.mod, hp1)
## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pots
## Model theta Resid. df
## 1 fungicide + crithidia + avg_pollen 17.43233 31
## 2 fungicide + crithidia + avg_pollen + workers_alive 19.66709 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -127.9116
## 2 -126.0351 1 vs 2 1 1.8765 0.1707324
Anova(hp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pots
## LR Chisq Df Pr(>Chisq)
## fungicide 4.4048 1 0.035839 *
## crithidia 1.4753 1 0.224517
## avg_pollen 7.0710 1 0.007834 **
## workers_alive 1.8821 1 0.170097
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(hp.mod, hp1)
## df AIC
## hp.mod 6 138.0351
## hp1 5 137.9116
Anova(hp1)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pots
## LR Chisq Df Pr(>Chisq)
## fungicide 3.8596 1 0.04946 *
## crithidia 2.1588 1 0.14176
## avg_pollen 20.7507 1 5.231e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(hp1)
##
## Call:
## glm.nb(formula = honey_pots ~ fungicide + crithidia + avg_pollen,
## data = brood, init.theta = 17.43232521, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1038 -0.7175 -0.3289 0.4007 2.2807
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07772 0.33372 -0.233 0.8158
## fungicideTRUE 0.43470 0.22245 1.954 0.0507 .
## crithidiaTRUE -0.34171 0.23437 -1.458 0.1448
## avg_pollen 2.22358 0.47589 4.672 2.98e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(17.4323) family taken to be 1)
##
## Null deviance: 65.926 on 34 degrees of freedom
## Residual deviance: 35.690 on 31 degrees of freedom
## AIC: 137.91
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 17.4
## Std. Err.: 25.6
##
## 2 x log-likelihood: -127.912
hpem_contrast <- emmeans(hp1, pairwise ~ fungicide, type = "response")
hpem_contrast
## $emmeans
## fungicide response SE df asymp.LCL asymp.UCL
## FALSE 1.96 0.353 Inf 1.37 2.79
## TRUE 3.02 0.453 Inf 2.25 4.05
##
## Results are averaged over the levels of: crithidia
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## FALSE / TRUE 0.647 0.144 Inf 1 -1.954 0.0507
##
## Results are averaged over the levels of: crithidia
## Tests are performed on the log scale
hpem <- emmeans(hp1, pairwise ~ fungicide*crithidia, type = "response")
hpem
## $emmeans
## fungicide crithidia response SE df asymp.LCL asymp.UCL
## FALSE FALSE 2.32 0.467 Inf 1.56 3.44
## TRUE FALSE 3.58 0.640 Inf 2.53 5.09
## FALSE TRUE 1.65 0.376 Inf 1.05 2.58
## TRUE TRUE 2.55 0.513 Inf 1.72 3.78
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## FALSE FALSE / TRUE FALSE 0.647 0.144 Inf 1 -1.954 0.2057
## FALSE FALSE / FALSE TRUE 1.407 0.330 Inf 1 1.458 0.4631
## FALSE FALSE / TRUE TRUE 0.911 0.290 Inf 1 -0.292 0.9914
## TRUE FALSE / FALSE TRUE 2.174 0.712 Inf 1 2.371 0.0827
## TRUE FALSE / TRUE TRUE 1.407 0.330 Inf 1 1.458 0.4631
## FALSE TRUE / TRUE TRUE 0.647 0.144 Inf 1 -1.954 0.2057
##
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
hpem.df <- as.data.frame(hpem$emmeans)
hpem.df
## fungicide crithidia response SE df asymp.LCL asymp.UCL
## FALSE FALSE 2.320664 0.4672846 Inf 1.563925 3.443568
## TRUE FALSE 3.584248 0.6400565 Inf 2.525776 5.086291
## FALSE TRUE 1.648959 0.3759594 Inf 1.054721 2.577995
## TRUE TRUE 2.546805 0.5126404 Inf 1.716561 3.778610
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
hpem.df$treatment <- c(1, 2, 4, 3)
hpem.df$treatment <-as.factor(hpem.df$treatment)
hp_sum <- brood %>%
group_by(treatment) %>%
summarize(m = mean(honey_pots),
sd = sd(honey_pots),
l = length(honey_pots)) %>%
mutate(se = sqrt(sd/l))
hp_sum
## # A tibble: 4 × 5
## treatment m sd l se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 3.33 2.5 9 0.527
## 2 2 4.11 3.02 9 0.579
## 3 3 2.56 2.24 9 0.499
## 4 4 1.5 1.60 8 0.448
hpcld <- cld(object = hpem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
hpcld
## fungicide crithidia response SE df asymp.LCL asymp.UCL .group
## FALSE TRUE 1.65 0.376 Inf 0.934 2.91 a
## FALSE FALSE 2.32 0.467 Inf 1.405 3.83 a
## TRUE TRUE 2.55 0.513 Inf 1.543 4.20 a
## TRUE FALSE 3.58 0.640 Inf 2.297 5.59 a
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
plot(brood$treatment, brood$honey_pots)
ggplot(data = hpem.df, aes(x = treatment, y = response, fill = treatment)) +
geom_col_pattern(
aes(pattern_density = treatment),
pattern = "stripe", # Set a common pattern type, but differentiate density
pattern_density = c(0, 0, 0, 0.4), # Add density for the third column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
coord_cartesian(ylim = c(0, 5.5)) +
geom_errorbar(aes(ymin = response - SE, ymax = response + SE), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Honeypots") +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 4,
label = "P = 0.05",
size = 7
) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_density_manual(values = c(0, 0, 0, 0.4)) +
scale_x_discrete(labels = custom_labels) +
theme(legend.position = "none") +
geom_segment(x = 2, xend = 3, y = 4.6, yend = 4.6,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 4.4, yend = 4.8,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 4.4, yend = 4.8,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 4, y = 5.2, yend = 5.2,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 5, yend = 5.4,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 5, yend = 5.4,
lineend = "round", linejoin = "round") +
geom_text(x = 2.5, y = 4.5, label = "b", size = 6, vjust = -0.5) +
geom_text(x = 2.5, y = 5.2, label = "a", size = 6, vjust = -0.5)
plot(brood$treatment, brood$drones)
dronecount.mod <- glm(total_drones ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
summary(dronecount.mod) #overdisp
##
## Call:
## glm(formula = total_drones ~ fungicide + crithidia + workers_alive +
## block + duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.94047 -0.89637 -0.09473 0.24751 2.06268
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.954e+00 1.470e+00 2.009 0.0446 *
## fungicideTRUE 1.074e-02 1.403e-01 0.077 0.9390
## crithidiaTRUE 8.278e-02 1.492e-01 0.555 0.5790
## workers_alive 1.287e-01 1.013e-01 1.270 0.2040
## block4 -2.771e-01 3.848e-01 -0.720 0.4714
## block6 -1.855e+01 2.852e+03 -0.007 0.9948
## block7 4.513e-03 3.590e-01 0.013 0.9900
## block8 -3.604e-01 3.738e-01 -0.964 0.3350
## block9 2.126e-01 3.281e-01 0.648 0.5169
## block10 9.119e-01 3.700e-01 2.465 0.0137 *
## block11 2.401e-01 3.955e-01 0.607 0.5437
## block12 -3.373e-02 4.159e-01 -0.081 0.9354
## duration -7.398e-02 2.991e-02 -2.473 0.0134 *
## avg_pollen 3.271e+00 6.989e-01 4.680 2.87e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 300.409 on 34 degrees of freedom
## Residual deviance: 41.778 on 21 degrees of freedom
## AIC: 166.28
##
## Number of Fisher Scoring iterations: 16
dronecount.mod.nb <- glm.nb(total_drones ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen, data = brood)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(total_drones ~ fungicide + crithidia + workers_alive + :
## alternation limit reached
qqnorm(resid(dronecount.mod));qqline(resid(dronecount.mod))
qqnorm(resid(dronecount.mod.nb));qqline(resid(dronecount.mod.nb))
anova(dronecount.mod, dronecount.mod.nb, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: total_drones ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Model 2: total_drones ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 21 41.778
## 2 21 41.772 0 0.0061185
AIC(dronecount.mod, dronecount.mod.nb)
## df AIC
## dronecount.mod 14 166.2838
## dronecount.mod.nb 15 168.2846
dronecount.mod.int <- glm.nb(total_drones ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen, data = brood)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(total_drones ~ fungicide + crithidia + workers_alive + :
## alternation limit reached
drop1(dronecount.mod.int, test = "Chisq")
## Single term deletions
##
## Model:
## total_drones ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 41.772 166.28
## fungicide 1 41.777 164.29 0.006 0.9391
## crithidia 1 42.079 164.59 0.307 0.5794
## workers_alive 1 43.411 165.92 1.639 0.2004
## block 8 79.000 187.51 37.228 1.045e-05 ***
## duration 1 47.542 170.06 5.770 0.0163 *
## avg_pollen 1 65.217 187.73 23.445 1.285e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dc1 <- update(dronecount.mod.int, .~. -workers_alive)
## 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(dc1, test = "Chisq")
## Single term deletions
##
## Model:
## total_drones ~ fungicide + crithidia + block + duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.412 165.92
## fungicide 1 43.475 163.99 0.063 0.801766
## crithidia 1 43.687 164.20 0.275 0.599952
## block 8 83.402 189.91 39.991 3.216e-06 ***
## duration 1 51.319 171.83 7.907 0.004925 **
## avg_pollen 1 93.953 214.47 50.541 1.167e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dronecount.mod.int)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_drones
## LR Chisq Df Pr(>Chisq)
## fungicide 0.006 1 0.9391
## crithidia 0.307 1 0.5794
## workers_alive 1.639 1 0.2004
## block 37.228 8 1.045e-05 ***
## duration 5.770 1 0.0163 *
## avg_pollen 23.445 1 1.285e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(brood$treatment, brood$total_drones)
qqnorm(resid(dc1));qqline(resid(dc1))
summary(dc1)
##
## Call:
## glm.nb(formula = total_drones ~ fungicide + crithidia + block +
## duration + avg_pollen, data = brood, init.theta = 43374.54429,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8968 -1.0530 -0.0842 0.3321 2.1666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.64096 1.33609 2.725 0.00643 **
## fungicideTRUE -0.03409 0.13573 -0.251 0.80169
## crithidiaTRUE 0.07817 0.14893 0.525 0.59968
## block4 -0.42159 0.37217 -1.133 0.25730
## block6 -19.09387 3649.61836 -0.005 0.99583
## block7 -0.11553 0.35020 -0.330 0.74148
## block8 -0.46003 0.37022 -1.243 0.21402
## block9 0.23826 0.32755 0.727 0.46697
## block10 0.81299 0.35742 2.275 0.02293 *
## block11 0.21506 0.39491 0.545 0.58604
## block12 -0.30541 0.36665 -0.833 0.40485
## duration -0.08161 0.02849 -2.864 0.00418 **
## avg_pollen 3.81456 0.56395 6.764 1.34e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(43374.54) family taken to be 1)
##
## Null deviance: 300.364 on 34 degrees of freedom
## Residual deviance: 43.412 on 22 degrees of freedom
## AIC: 167.92
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 43375
## Std. Err.: 765731
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -139.924
Anova(dc1)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_drones
## LR Chisq Df Pr(>Chisq)
## fungicide 0.063 1 0.801766
## crithidia 0.275 1 0.599952
## block 39.991 8 3.216e-06 ***
## duration 7.907 1 0.004925 **
## avg_pollen 50.541 1 1.167e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drones_sum <- brood %>%
group_by(treatment) %>%
summarise(mb = mean(total_drones),
nb = length(total_drones),
sdb = sd(total_drones)) %>%
mutate(seb = (sdb/sqrt(nb)))
drones_sum
## # A tibble: 4 × 5
## treatment mb nb sdb seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 11.3 9 7.84 2.61
## 2 2 7.67 9 7.37 2.46
## 3 3 5.44 9 8.16 2.72
## 4 4 5.62 8 6.52 2.31
drones_sum$plot <- drones_sum$mb + drones_sum$seb
ggplot(data = drones_sum, aes(x = treatment, y = mb, 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, 17)) +
geom_errorbar(aes(ymin = mb - seb, ymax = mb + seb), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Average Adults Males") +
annotate(
geom = "text",
x = c(1, 2, 3, 4),
y = drones_sum$plot + 1, # Adjust the y-position as needed
label = c("a", "a", "a", "a"),
size = 8
) +
theme_classic(base_size = 20) +
annotate(
geom = "text",
x = 1,
y = 16,
label = "P > 0.5",
size = 7
) +
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")
plmod <- glm(cbind(live_larvae, dead_larvae) ~ fungicide + crithidia + avg_pollen + block + duration + workers_alive, data = brood, family = binomial("logit"))
Anova(plmod)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0127 1 0.9103281
## crithidia 0.7930 1 0.3731960
## avg_pollen 1.3008 1 0.2540613
## block 30.5370 8 0.0001698 ***
## duration 8.0601 1 0.0045251 **
## workers_alive 0.8537 1 0.3555162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(plmod, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_larvae, dead_larvae) ~ fungicide + crithidia + avg_pollen +
## block + duration + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.417 106.90
## fungicide 1 42.430 104.91 0.0127 0.9103281
## crithidia 1 43.210 105.69 0.7930 0.3731960
## avg_pollen 1 43.718 106.20 1.3008 0.2540613
## block 8 72.954 121.44 30.5370 0.0001698 ***
## duration 1 50.477 112.96 8.0601 0.0045251 **
## workers_alive 1 43.271 105.75 0.8537 0.3555162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plmod1 <- update(plmod, .~. -workers_alive)
drop1(plmod1, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_larvae, dead_larvae) ~ fungicide + crithidia + avg_pollen +
## block + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.271 105.75
## fungicide 1 43.273 103.75 0.0021 0.9632944
## crithidia 1 43.748 104.23 0.4772 0.4896903
## avg_pollen 1 43.908 104.39 0.6368 0.4248650
## block 8 72.970 119.45 29.6994 0.0002389 ***
## duration 1 56.543 117.02 13.2719 0.0002694 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plmod2 <- update(plmod1, .~. -avg_pollen)
drop1(plmod2, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_larvae, dead_larvae) ~ fungicide + crithidia + block +
## duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.908 104.39
## fungicide 1 43.914 102.39 0.006 0.9391168
## crithidia 1 44.403 102.88 0.495 0.4815381
## block 8 77.421 121.90 33.513 4.976e-05 ***
## duration 1 58.001 116.48 14.094 0.0001739 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(plmod2)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.006 1 0.9391168
## crithidia 0.495 1 0.4815381
## block 33.513 8 4.976e-05 ***
## duration 14.094 1 0.0001739 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(plmod2)
##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ fungicide + crithidia +
## block + duration, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2069 -0.5505 0.0000 0.8991 2.5295
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 17.37268 4.54637 3.821 0.000133 ***
## fungicideTRUE -0.03042 0.39827 -0.076 0.939123
## crithidiaTRUE 0.27725 0.39595 0.700 0.483790
## block4 -1.90856 0.99419 -1.920 0.054894 .
## block6 -1.24483 1.02595 -1.213 0.225000
## block7 -2.48520 0.96780 -2.568 0.010232 *
## block8 -1.63314 0.99628 -1.639 0.101164
## block9 -2.40994 0.94892 -2.540 0.011096 *
## block10 1.72722 0.80604 2.143 0.032126 *
## block11 -1.76714 1.36413 -1.295 0.195171
## block12 -0.11250 0.77348 -0.145 0.884353
## duration -0.32205 0.09431 -3.415 0.000639 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 79.465 on 26 degrees of freedom
## Residual deviance: 43.908 on 15 degrees of freedom
## AIC: 104.39
##
## Number of Fisher Scoring iterations: 5
plot(plmod2)
## Warning: not plotting observations with leverage one:
## 2
ppmod <- glm(cbind(live_pupae, dead_pupae) ~ fungicide + crithidia + avg_pollen + block + duration + workers_alive, data = brood, family = binomial("logit"))
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Anova(ppmod)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0000 1 0.999961
## crithidia 0.0000 1 0.999980
## avg_pollen 0.0000 1 0.999990
## block 7.8251 8 0.450741
## duration 7.3586 1 0.006674 **
## workers_alive 0.0000 1 0.999913
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(ppmod, test = "Chisq")
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ fungicide + crithidia + avg_pollen +
## block + duration + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 0.0000 34.810
## fungicide 1 0.0000 32.810 0.0000 0.999961
## crithidia 1 0.0000 32.810 0.0000 0.999980
## avg_pollen 1 0.0000 32.810 0.0000 0.999990
## block 8 7.8251 26.635 7.8251 0.450741
## duration 1 7.3586 40.168 7.3586 0.006674 **
## workers_alive 1 0.0000 32.810 0.0000 0.999913
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ppmod1 <- update(ppmod, .~. -block)
drop1(ppmod1, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ fungicide + crithidia + avg_pollen +
## duration + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 7.8251 26.635
## fungicide 1 7.9246 24.734 0.0996 0.7524
## crithidia 1 8.1487 24.958 0.3236 0.5694
## avg_pollen 1 9.8300 26.640 2.0050 0.1568
## duration 1 9.1617 25.971 1.3366 0.2476
## workers_alive 1 13.0866 29.896 5.2615 0.0218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ppmod2 <- update(ppmod1, .~. -duration)
drop1(ppmod2, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ fungicide + crithidia + avg_pollen +
## workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 9.1617 25.971
## fungicide 1 9.3080 24.118 0.1464 0.70204
## crithidia 1 9.1617 23.971 0.0001 0.99348
## avg_pollen 1 10.0306 24.840 0.8689 0.35125
## workers_alive 1 15.4089 30.218 6.2473 0.01244 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ppmod3 <- update(ppmod2, .~. -avg_pollen)
drop1(ppmod3, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ fungicide + crithidia + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 10.031 24.840
## fungicide 1 10.111 22.921 0.0804 0.77680
## crithidia 1 10.227 23.037 0.1967 0.65742
## workers_alive 1 16.167 28.977 6.1367 0.01324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(ppmod3)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0804 1 0.77680
## crithidia 0.1967 1 0.65742
## workers_alive 6.1367 1 0.01324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ppmod3)
##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ fungicide + crithidia +
## workers_alive, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9211 0.0000 0.2516 0.3382 0.8204
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03532 2.30205 0.015 0.9878
## fungicideTRUE -0.38188 1.35530 -0.282 0.7781
## crithidiaTRUE 0.62650 1.46413 0.428 0.6687
## workers_alive 0.95585 0.41792 2.287 0.0222 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 18.342 on 23 degrees of freedom
## Residual deviance: 10.031 on 20 degrees of freedom
## AIC: 24.84
##
## Number of Fisher Scoring iterations: 6
brood$live.lp <- brood$live_larvae + brood$live_pupae
brood$dead.lp <- brood$dead_larvae + brood$dead_pupae
lp.mod <- glm(cbind(live.lp, dead.lp) ~ fungicide + crithidia + block + duration, data = brood, family = binomial("logit"))
drop1(lp.mod, test = "Chisq")
## Single term deletions
##
## Model:
## cbind(live.lp, dead.lp) ~ fungicide + crithidia + block + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 44.974 109.37
## fungicide 1 44.979 107.37 0.0056 0.9402700
## crithidia 1 45.383 107.78 0.4091 0.5224156
## block 8 71.399 119.79 26.4247 0.0008882 ***
## duration 1 59.367 121.76 14.3931 0.0001483 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lp.mod)
##
## Call:
## glm(formula = cbind(live.lp, dead.lp) ~ fungicide + crithidia +
## block + duration, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5478 -0.3930 0.0000 0.6606 2.3412
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 17.28245 4.30768 4.012 6.02e-05 ***
## fungicideTRUE -0.02839 0.37895 -0.075 0.94028
## crithidiaTRUE 0.23684 0.37211 0.636 0.52446
## block4 -2.01994 0.95821 -2.108 0.03503 *
## block6 -1.41021 0.96714 -1.458 0.14481
## block7 -2.42735 0.92378 -2.628 0.00860 **
## block8 -1.50318 0.96952 -1.550 0.12104
## block9 -2.29364 0.88581 -2.589 0.00962 **
## block10 1.03795 0.70848 1.465 0.14291
## block11 -1.64692 1.31903 -1.249 0.21182
## block12 -0.56740 0.72634 -0.781 0.43470
## duration -0.31120 0.08941 -3.481 0.00050 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 76.134 on 26 degrees of freedom
## Residual deviance: 44.974 on 15 degrees of freedom
## AIC: 109.37
##
## Number of Fisher Scoring iterations: 5
Anova(lp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live.lp, dead.lp)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0056 1 0.9402700
## crithidia 0.4091 1 0.5224156
## block 26.4247 8 0.0008882 ***
## duration 14.3931 1 0.0001483 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drones$fungicide <- as.factor(drones$fungicide)
drones$crithidia <- as.factor(drones$crithidia)
plot(drones$treatment, drones$dry_weight)
plot(drones_rf$treatment, drones_rf$dry_weight)
plot(brood$treatment, brood$drones)
shapiro.test(drones$dry_weight)
##
## Shapiro-Wilk normality test
##
## data: drones$dry_weight
## W = 0.99166, p-value = 0.1135
hist(drones$dry_weight)
range(drones$dry_weight)
## [1] 0.0166 0.0541
dry <- lmer(dry_weight ~ fungicide + crithidia + workers_alive + block + emerge + (1|colony), data = drones)
drop1(dry, test = "Chisq")
## Single term deletions
##
## Model:
## dry_weight ~ fungicide + crithidia + workers_alive + block +
## emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -1999.9
## fungicide 1 -2001.7 0.1493 0.69923
## crithidia 1 -1997.3 4.5582 0.03276 *
## workers_alive 1 -2000.1 1.7920 0.18069
## block 7 -1997.3 16.6047 0.02013 *
## emerge 1 -1998.7 3.1875 0.07420 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dry1 <- update(dry, .~. -workers_alive)
drop1(dry1, test= "Chisq")
## Single term deletions
##
## Model:
## dry_weight ~ fungicide + crithidia + block + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -2000.1
## fungicide 1 -2001.2 0.9199 0.33749
## crithidia 1 -1998.5 3.5794 0.05850 .
## block 7 -1998.1 16.0008 0.02511 *
## emerge 1 -1998.1 3.9614 0.04656 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dry2 <- update(dry1, .~. -block)
drop1(dry2, test = "Chisq")
## Single term deletions
##
## Model:
## dry_weight ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -1998.1
## fungicide 1 -1999.9 0.1553 0.69350
## crithidia 1 -1999.3 0.7209 0.39585
## emerge 1 -1993.8 6.2922 0.01213 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dry2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## fungicide 0.1985 1 0.65596
## crithidia 0.6241 1 0.42951
## emerge 6.0119 1 0.01421 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dry2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: dry_weight ~ fungicide + crithidia + emerge + (1 | colony)
## Data: drones
##
## REML criterion at convergence: -1958.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.86701 -0.59248 -0.03622 0.60958 2.58281
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 6.923e-06 0.002631
## Residual 4.271e-05 0.006535
## Number of obs: 281, groups: colony, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0522562 0.0059986 8.711
## fungicide1 -0.0006294 0.0014127 -0.446
## crithidia1 0.0011136 0.0014096 0.790
## emerge -0.0003927 0.0001602 -2.452
##
## Correlation of Fixed Effects:
## (Intr) fngcd1 crthd1
## fungicide1 0.073
## crithidia1 -0.072 -0.002
## emerge -0.983 -0.177 -0.023
sum_dry <- drones %>%
group_by(treatment) %>%
summarise(m = mean(dry_weight),
sd = sd(dry_weight),
n = length(dry_weight)) %>%
mutate(se = sd/sqrt(n))
sum_dry
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.0373 0.00752 102 0.000745
## 2 2 0.0376 0.00653 69 0.000786
## 3 3 0.0378 0.00706 49 0.00101
## 4 4 0.0393 0.00665 61 0.000852
drones.na <- na.omit(drones)
drones.na$alive <- as.logical(drones.na$`alive?`)
shapiro.test(drones$radial_cell)
##
## Shapiro-Wilk normality test
##
## data: drones$radial_cell
## W = 0.97715, p-value = 0.0001916
hist(drones$radial_cell)
descdist(drones.na$radial_cell, discrete = FALSE)
## summary statistics
## ------
## min: 2073.526 max: 3083.439
## median: 2710.923
## mean: 2694.575
## estimated sd: 171.6365
## estimated skewness: -0.5726179
## estimated kurtosis: 4.019655
range(drones$radial_cell)
## [1] NA NA
drones.na$square <- drones.na$radial_cell^3
shapiro.test(drones.na$square)
##
## Shapiro-Wilk normality test
##
## data: drones.na$square
## W = 0.99215, p-value = 0.1516
hist(drones.na$square)
drones.na$log <- log(drones.na$radial_cell)
shapiro.test(drones.na$log)
##
## Shapiro-Wilk normality test
##
## data: drones.na$log
## W = 0.95963, p-value = 5.927e-07
hist(drones.na$square)
rad_mod <- lmer(square ~ fungicide + crithidia + workers_alive + block + mean.pollen + emerge + alive + (1|colony), data = drones.na)
drop1(rad_mod, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + block + mean.pollen +
## emerge + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12934
## fungicide 1 12939 7.2355 0.007148 **
## crithidia 1 12933 1.3953 0.237507
## workers_alive 1 12934 1.7286 0.188584
## block 7 12932 11.7617 0.108673
## mean.pollen 1 12935 3.3474 0.067311 .
## emerge 1 12936 4.4149 0.035627 *
## alive 1 12933 1.4090 0.235217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm1 <- update(rad_mod, .~. -block)
drop1(rm1, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + mean.pollen +
## emerge + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12932
## fungicide 1 12932 2.9510 0.085825 .
## crithidia 1 12930 0.7704 0.380101
## workers_alive 1 12930 0.2499 0.617119
## mean.pollen 1 12930 0.2225 0.637150
## emerge 1 12936 6.7496 0.009377 **
## alive 1 12930 0.7498 0.386527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm2 <- update(rm1, .~. -workers_alive)
drop1(rm2, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + mean.pollen + emerge + alive +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12930
## fungicide 1 12931 3.5558 0.05934 .
## crithidia 1 12928 0.6686 0.41353
## mean.pollen 1 12928 0.6137 0.43340
## emerge 1 12935 7.0223 0.00805 **
## alive 1 12929 0.7897 0.37419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm3 <- update(rm2, .~. -mean.pollen)
drop1(rm3, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + emerge + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12928
## fungicide 1 12930 3.6326 0.056657 .
## crithidia 1 12927 0.5522 0.457427
## emerge 1 12935 8.1007 0.004425 **
## alive 1 12927 0.8091 0.368391
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm4 <- update(rm3, .~. -alive)
drop1(rm4, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12927
## fungicide 1 12929 3.8119 0.050889 .
## crithidia 1 12926 0.6516 0.419537
## emerge 1 12933 7.9923 0.004698 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(rm4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: square ~ fungicide + crithidia + emerge + (1 | colony)
## Data: drones.na
##
## REML criterion at convergence: 12752
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9415 -0.5357 0.0455 0.5915 2.9574
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 8.620e+17 9.284e+08
## Residual 1.189e+19 3.448e+09
## Number of obs: 276, groups: colony, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.858e+10 3.044e+09 9.388
## fungicide1 -1.143e+09 6.022e+08 -1.898
## crithidia1 4.190e+08 5.965e+08 0.703
## emerge -2.274e+08 8.179e+07 -2.781
##
## Correlation of Fixed Effects:
## (Intr) fngcd1 crthd1
## fungicide1 0.123
## crithidia1 -0.055 -0.009
## emerge -0.989 -0.206 -0.022
Anova(rm4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: square
## Chisq Df Pr(>Chisq)
## fungicide 3.6022 1 0.057704 .
## crithidia 0.4935 1 0.482353
## emerge 7.7324 1 0.005424 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rad_mod <- lmer(square ~ fungicide + crithidia + workers_alive + block + mean.pollen + days_active + alive + (1|colony), data = drones.na)
drop1(rad_mod, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + block + mean.pollen +
## days_active + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12938
## fungicide 1 12946 9.4396 0.002124 **
## crithidia 1 12938 2.0127 0.155989
## workers_alive 1 12939 2.5934 0.107306
## block 7 12938 13.9266 0.052503 .
## mean.pollen 1 12940 3.3790 0.066032 .
## days_active 1 12936 0.0170 0.896404
## alive 1 12938 1.2506 0.263436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm1 <- update(rad_mod, .~. -days_active)
drop1(rm1, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + block + mean.pollen +
## alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12936
## fungicide 1 12944 9.8595 0.00169 **
## crithidia 1 12936 2.0132 0.15593
## workers_alive 1 12937 2.6202 0.10551
## block 7 12936 14.0964 0.04949 *
## mean.pollen 1 12938 3.4592 0.06290 .
## alive 1 12936 1.2552 0.26255
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm2 <- update(rm1, .~. -block)
drop1(rm2, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + mean.pollen +
## alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12936
## fungicide 1 12939 4.4893 0.03411 *
## crithidia 1 12935 0.7453 0.38796
## workers_alive 1 12935 0.5227 0.46970
## mean.pollen 1 12935 0.6837 0.40833
## alive 1 12935 0.6207 0.43080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm3 <- update(rm2, .~. -workers_alive)
drop1(rm3, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + mean.pollen + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12935
## fungicide 1 12938 5.5520 0.01846 *
## crithidia 1 12934 0.5944 0.44074
## mean.pollen 1 12935 1.6920 0.19333
## alive 1 12934 0.6684 0.41362
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm5 <- update(rm3, .~. -alive)
drop1(rm5, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + mean.pollen + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12934
## fungicide 1 12937 5.7439 0.01655 *
## crithidia 1 12932 0.6878 0.40693
## mean.pollen 1 12933 1.7244 0.18913
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm6 <- update(rm5, .~. -mean.pollen)
summary(rm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: square ~ fungicide + crithidia + (1 | colony)
## Data: drones.na
##
## REML criterion at convergence: 12797.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85719 -0.55413 0.02787 0.61341 2.71516
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 1.111e+18 1.054e+09
## Residual 1.208e+19 3.475e+09
## Number of obs: 276, groups: colony, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.022e+10 4.863e+08 41.589
## fungicide1 -1.526e+09 6.296e+08 -2.424
## crithidia1 3.586e+08 6.373e+08 0.563
##
## Correlation of Fixed Effects:
## (Intr) fngcd1
## fungicide1 -0.559
## crithidia1 -0.518 -0.011
Anova(rm6)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: square
## Chisq Df Pr(>Chisq)
## fungicide 5.8751 1 0.01536 *
## crithidia 0.3165 1 0.57371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rm4, rm6, test = "Chisq")
## Data: drones.na
## Models:
## rm6: square ~ fungicide + crithidia + (1 | colony)
## rm4: square ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rm6 5 12933 12951 -6461.6 12923
## rm4 6 12927 12949 -6457.6 12915 7.9923 1 0.004698 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(rm4, rm6)
## df AIC
## rm4 6 12764.01
## rm6 5 12807.83
qqnorm(resid(rad_mod));qqline(resid(rad_mod))
qqnorm(resid(rm4));qqline(resid(rm4))
qqnorm(resid(rm6));qqline(resid(rm6))
Anova(rm3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: square
## Chisq Df Pr(>Chisq)
## fungicide 5.0629 1 0.02444 *
## crithidia 0.4090 1 0.52249
## mean.pollen 1.5012 1 0.22049
## alive 0.7139 1 0.39815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rem <- emmeans(rm6, pairwise ~ fungicide, type = "response")
rem
## $emmeans
## fungicide emmean SE df lower.CL upper.CL
## 0 2.04e+10 4.22e+08 17.1 1.95e+10 2.13e+10
## 1 1.89e+10 4.82e+08 21.1 1.79e+10 1.99e+10
##
## Results are averaged over the levels of: crithidia
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## fungicide0 - fungicide1 1.53e+09 6.34e+08 19.1 2.405 0.0265
##
## Results are averaged over the levels of: crithidia
## Degrees-of-freedom method: kenward-roger
re <- setDT(as.data.frame(rem$emmeans))
cont_radial <- setDT(as.data.frame(rem$contrasts))
rad.cld <- cld(object =rem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
rad.cld
## fungicide emmean SE df lower.CL upper.CL .group
## 1 1.89e+10 4.82e+08 21.1 1.77e+10 2.00e+10 a
## 0 2.04e+10 4.22e+08 17.1 1.94e+10 2.14e+10 b
##
## Results are averaged over the levels of: crithidia
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
sum_radial <- drones_rf %>%
group_by(treatment) %>%
summarise(m = mean(radial_cell),
sd = sd(radial_cell),
n = length(radial_cell)) %>%
mutate(se = sd/sqrt(n))
sum_radial
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 2699. 182. 99 18.3
## 2 2 2664. 174. 68 21.1
## 3 3 2660. 173. 49 24.7
## 4 4 2750. 136. 60 17.5
sum_radial$plot <- sum_radial$m + sum_radial$se
ggplot(sum_radial, aes(x = treatment, y = m, fill = treatment)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_viridis_d() +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Radial Cell Length (um)") +
theme_classic(base_size = 20) +
coord_cartesian(ylim=c(2620, 2800)) +
annotate(geom = "text",
x = 1, y = 3 ,
label = "P > 0.05",
size = 8) +
theme(legend.position = "none",
axis.text = element_text(size = 20), # Set axis label font size
axis.title = element_text(size = 20)) + # Set axis title font size
theme(text = element_text(size = 20)) +
annotate(geom = "text",
label = "P = 0.06",
x = 1, y = 2762,
size = 7)
drones_rf$fungicide <-as.factor(drones_rf$fungicide)
drones_rf$crithidia <- as.factor(drones_rf$crithidia)
shapiro.test(drones_rf$relative_fat_original)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$relative_fat_original
## W = 0.97273, p-value = 4.049e-05
hist(drones_rf$relative_fat_original)
plot(drones_rf$treatment, drones_rf$relative_fat_original)
range(drones_rf$relative_fat_original)
## [1] 1.78e-07 3.40e-06
drones_rf$log_ref <- log(drones_rf$relative_fat_original)
shapiro.test(drones_rf$log_ref)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$log_ref
## W = 0.94339, p-value = 8.074e-09
drones_rf$squarerf <- sqrt(drones_rf$relative_fat_original)
shapiro.test(drones_rf$squarerf)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$squarerf
## W = 0.99031, p-value = 0.06387
rf_mod <- lmer(squarerf ~ fungicide + crithidia + block + mean.pollen + workers_alive + emerge + (1|colony), data = drones_rf)
drop1(rf_mod, test = "Chisq")
## Single term deletions
##
## Model:
## squarerf ~ fungicide + crithidia + block + mean.pollen + workers_alive +
## emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3926.8
## fungicide 1 -3927.2 1.5763 0.2092949
## crithidia 1 -3926.6 2.1577 0.1418577
## block 7 -3922.1 18.6722 0.0092788 **
## mean.pollen 1 -3928.1 0.6181 0.4317364
## workers_alive 1 -3928.7 0.0348 0.8519805
## emerge 1 -3916.1 12.6298 0.0003796 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf1 <- update(rf_mod, .~. -mean.pollen)
drop1(rf1, test = "Chisq")
## Single term deletions
##
## Model:
## squarerf ~ fungicide + crithidia + block + workers_alive + emerge +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3928.1
## fungicide 1 -3928.8 1.3699 0.2418298
## crithidia 1 -3927.9 2.2779 0.1312301
## block 7 -3922.2 19.9284 0.0057265 **
## workers_alive 1 -3929.8 0.3342 0.5632243
## emerge 1 -3917.5 12.6776 0.0003701 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf2 <- update(rf1, .~. -workers_alive)
drop1(rf2, test = "Chisq")
## Single term deletions
##
## Model:
## squarerf ~ fungicide + crithidia + block + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3929.8
## fungicide 1 -3930.7 1.0627 0.3025989
## crithidia 1 -3929.1 2.6716 0.1021510
## block 7 -3924.2 19.6016 0.0064976 **
## emerge 1 -3919.5 12.3463 0.0004419 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf3 <- update(rf2, .~. -block)
drop1(rf3, test = "Chisq")
## Single term deletions
##
## Model:
## squarerf ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3924.2
## fungicide 1 -3926.2 0.0054 0.94158
## crithidia 1 -3922.0 4.2162 0.04004 *
## emerge 1 -3908.2 17.9955 2.214e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf4 <- update(rf3, .~. -crithidia)
drop1(rf4, test = "Chisq")
## Single term deletions
##
## Model:
## squarerf ~ fungicide + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -3922.0
## fungicide 1 -3924.0 0.0137 0.9068
## emerge 1 -3906.6 17.3510 3.107e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rf_mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: squarerf
## Chisq Df Pr(>Chisq)
## fungicide 0.8451 1 0.3579549
## crithidia 1.7576 1 0.1849263
## block 11.6312 7 0.1133564
## mean.pollen 0.0127 1 0.9101582
## workers_alive 0.0298 1 0.8628689
## emerge 14.1663 1 0.0001673 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rf4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: squarerf
## Chisq Df Pr(>Chisq)
## fungicide 0.0213 1 0.8839
## emerge 17.9437 1 2.275e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rf3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: squarerf
## Chisq Df Pr(>Chisq)
## fungicide 0.0149 1 0.90287
## crithidia 3.9457 1 0.04699 *
## emerge 18.5688 1 1.639e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(rf3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: squarerf ~ fungicide + crithidia + emerge + (1 | colony)
## Data: drones_rf
##
## REML criterion at convergence: -3857.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8428 -0.5043 0.0419 0.5569 3.9576
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 8.144e-09 9.024e-05
## Residual 3.419e-08 1.849e-04
## Number of obs: 276, groups: colony, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.789e-03 1.769e-04 10.113
## fungicide1 -5.541e-06 4.541e-05 -0.122
## crithidia1 9.031e-05 4.546e-05 1.986
## emerge -2.031e-05 4.714e-06 -4.309
##
## Correlation of Fixed Effects:
## (Intr) fngcd1 crthd1
## fungicide1 0.046
## crithidia1 -0.070 0.002
## emerge -0.980 -0.162 -0.035
anova(rf3, rf4, test = "Chisq")
## Data: drones_rf
## Models:
## rf4: squarerf ~ fungicide + emerge + (1 | colony)
## rf3: squarerf ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rf4 5 -3922.0 -3903.9 1966.0 -3932.0
## rf3 6 -3924.2 -3902.5 1968.1 -3936.2 4.2162 1 0.04004 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(rf3, rf4)
## df AIC
## rf3 6 -3845.295
## rf4 5 -3861.588
qqnorm(resid(rf3));qqline(resid(rf3))
qqnorm(resid(rf4));qqline(resid(rf4))
rf_em <- emmeans(rf3, pairwise ~ crithidia*fungicide, type = "response")
rf_em
## $emmeans
## crithidia fungicide emmean SE df lower.CL upper.CL
## 0 0 0.00103 3.55e-05 17.8 0.000955 0.0011
## 1 0 0.00112 4.07e-05 19.3 0.001035 0.0012
## 0 1 0.00102 3.84e-05 22.1 0.000944 0.0011
## 1 1 0.00111 4.33e-05 22.1 0.001024 0.0012
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## crithidia0 fungicide0 - crithidia1 fungicide0 -9.03e-05 4.56e-05 19.8 -1.979
## crithidia0 fungicide0 - crithidia0 fungicide1 5.54e-06 4.56e-05 20.9 0.122
## crithidia0 fungicide0 - crithidia1 fungicide1 -8.48e-05 6.46e-05 20.1 -1.312
## crithidia1 fungicide0 - crithidia0 fungicide1 9.59e-05 6.44e-05 20.6 1.488
## crithidia1 fungicide0 - crithidia1 fungicide1 5.54e-06 4.56e-05 20.9 0.122
## crithidia0 fungicide1 - crithidia1 fungicide1 -9.03e-05 4.56e-05 19.8 -1.979
## p.value
## 0.2291
## 0.9993
## 0.5660
## 0.4621
## 0.9993
## 0.2291
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
rf_e <- setDT(as.data.frame(rf_em$emmeans))
rf_ce <- setDT(as.data.frame(rf_em$contrasts))
rf_em
## $emmeans
## crithidia fungicide emmean SE df lower.CL upper.CL
## 0 0 0.00103 3.55e-05 17.8 0.000955 0.0011
## 1 0 0.00112 4.07e-05 19.3 0.001035 0.0012
## 0 1 0.00102 3.84e-05 22.1 0.000944 0.0011
## 1 1 0.00111 4.33e-05 22.1 0.001024 0.0012
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio
## crithidia0 fungicide0 - crithidia1 fungicide0 -9.03e-05 4.56e-05 19.8 -1.979
## crithidia0 fungicide0 - crithidia0 fungicide1 5.54e-06 4.56e-05 20.9 0.122
## crithidia0 fungicide0 - crithidia1 fungicide1 -8.48e-05 6.46e-05 20.1 -1.312
## crithidia1 fungicide0 - crithidia0 fungicide1 9.59e-05 6.44e-05 20.6 1.488
## crithidia1 fungicide0 - crithidia1 fungicide1 5.54e-06 4.56e-05 20.9 0.122
## crithidia0 fungicide1 - crithidia1 fungicide1 -9.03e-05 4.56e-05 19.8 -1.979
## p.value
## 0.2291
## 0.9993
## 0.5660
## 0.4621
## 0.9993
## 0.2291
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
rf_e
## crithidia fungicide emmean SE df lower.CL
## 1: 0 0 0.001029495 3.550520e-05 17.79935 0.0009548407
## 2: 1 0 0.001119806 4.067936e-05 19.26722 0.0010347432
## 3: 0 1 0.001023954 3.837426e-05 22.04877 0.0009443804
## 4: 1 1 0.001114265 4.334502e-05 22.07399 0.0010243905
## upper.CL
## 1: 0.001104149
## 2: 0.001204869
## 3: 0.001103527
## 4: 0.001204140
rf_ce
## contrast estimate SE
## 1: crithidia0 fungicide0 - crithidia1 fungicide0 -9.031152e-05 4.563038e-05
## 2: crithidia0 fungicide0 - crithidia0 fungicide1 5.541136e-06 4.558328e-05
## 3: crithidia0 fungicide0 - crithidia1 fungicide1 -8.477039e-05 6.459089e-05
## 4: crithidia1 fungicide0 - crithidia0 fungicide1 9.585266e-05 6.440459e-05
## 5: crithidia1 fungicide0 - crithidia1 fungicide1 5.541136e-06 4.558328e-05
## 6: crithidia0 fungicide1 - crithidia1 fungicide1 -9.031152e-05 4.563038e-05
## df t.ratio p.value
## 1: 19.82692 -1.9791972 0.2290557
## 2: 20.88428 0.1215607 0.9993362
## 3: 20.10602 -1.3124202 0.5659652
## 4: 20.58874 1.4882893 0.4621385
## 5: 20.88428 0.1215607 0.9993362
## 6: 19.82692 -1.9791972 0.2290557
rf_e$treatment <- c(1, 4, 2, 3)
rf_e$treatment <- as.factor(rf_e$treatment)
rf_sum <- drones_rf %>%
group_by(treatment) %>%
summarise(m = mean(relative_fat_original),
sd = sd(relative_fat_original),
n = length(relative_fat_original)) %>%
mutate(se = sd/sqrt(n))
rf_sum
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.00000115 0.000000438 99 0.0000000440
## 2 2 0.00000111 0.000000408 68 0.0000000495
## 3 3 0.00000129 0.000000508 49 0.0000000726
## 4 4 0.00000128 0.000000435 60 0.0000000561
rf_sum$plot <- rf_sum$m + rf_sum$se
ggplot(rf_e, aes(x = treatment, y = emmean, fill = treatment)) +
geom_bar(stat = "identity", color = "black") +
geom_col_pattern(
aes(pattern = treatment),
pattern_density = c(0, 0, 0, 0.4), # Add density for the fourth column
pattern_spacing = 0.03,
position = position_dodge(0.9)
) +
scale_fill_viridis_d() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Square Root (Relative Fat (g/mm))") +
theme_classic(base_size = 20) +
coord_cartesian(ylim=c(0.0009, 0.00121)) +
annotate(geom = "text",
x = 1, y = 0.00119,
label = "P = 0.05",
size = 8) +
theme(legend.position = "none",
axis.text = element_text(size = 16), # Set axis label font size
axis.title = element_text(size = 16)) + # Set axis title font size
theme(text = element_text(size = 16)) +
scale_x_discrete(labels = custom_labels) +
scale_fill_manual(values = c("lightgreen", "lightblue", "lightblue", "grey")) +
scale_pattern_manual(values = c("none", "none", "stripe", "none")) +
geom_segment(x = 1, xend = 2, y = 0.0011, yend = 0.0011,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 0.00109, yend = 0.00111,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 0.00109, yend = 0.00111,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 4, y = 0.0012, yend = 0.0012,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 0.00119, yend = 0.00121,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 0.00119, yend = 0.00121,
lineend = "round", linejoin = "round") +
geom_text(x = 1.5, y = 0.00111, label = "a", size = 6, vjust = -0.5) +
geom_text(x = 3.5, y = 0.001201, label = "b", size = 6, vjust = -0.5) +
theme(legend.position = "none")
drones$fungicide <- as.logical(drones$fungicide)
drones$crithidia <- as.logical(drones$crithidia)
em.mod <- glm(emerge ~ fungicide + crithidia + dry_weight + live_weight + workers_alive + mean.pollen, data = drones.na, family = "poisson")
summary(em.mod)
##
## Call:
## glm(formula = emerge ~ fungicide + crithidia + dry_weight + live_weight +
## workers_alive + mean.pollen, family = "poisson", data = drones.na)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.98926 -0.33225 0.00052 0.28013 1.76734
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.739837 0.083374 44.856 <2e-16 ***
## fungicide1 0.033579 0.020613 1.629 0.103
## crithidia1 -0.003914 0.020394 -0.192 0.848
## dry_weight -1.691845 1.430164 -1.183 0.237
## live_weight 0.015777 0.023522 0.671 0.502
## workers_alive -0.008052 0.014627 -0.551 0.582
## mean.pollen -0.073468 0.063707 -1.153 0.249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 58.303 on 275 degrees of freedom
## Residual deviance: 48.891 on 269 degrees of freedom
## AIC: 1570.1
##
## Number of Fisher Scoring iterations: 3
em.mod <- glm.nb(emerge ~ fungicide + crithidia + dry_weight + live_weight + workers_alive + mean.pollen, data = drones.na)
## 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
summary(em.mod)
##
## Call:
## glm.nb(formula = emerge ~ fungicide + crithidia + dry_weight +
## live_weight + workers_alive + mean.pollen, data = drones.na,
## init.theta = 3809154.841, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.98925 -0.33225 0.00052 0.28013 1.76733
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.739837 0.083375 44.856 <2e-16 ***
## fungicide1 0.033579 0.020613 1.629 0.103
## crithidia1 -0.003914 0.020394 -0.192 0.848
## dry_weight -1.691845 1.430174 -1.183 0.237
## live_weight 0.015777 0.023523 0.671 0.502
## workers_alive -0.008052 0.014627 -0.551 0.582
## mean.pollen -0.073468 0.063708 -1.153 0.249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3809155) family taken to be 1)
##
## Null deviance: 58.303 on 275 degrees of freedom
## Residual deviance: 48.890 on 269 degrees of freedom
## AIC: 1572.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3809155
## Std. Err.: 44034635
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -1556.135
drop1(em.mod, test = "Chisq")
## Single term deletions
##
## Model:
## emerge ~ fungicide + crithidia + dry_weight + live_weight + workers_alive +
## mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 48.890 1570.1
## fungicide 1 51.540 1570.8 2.64962 0.1036
## crithidia 1 48.927 1568.2 0.03685 0.8478
## dry_weight 1 50.288 1569.5 1.39821 0.2370
## live_weight 1 49.341 1568.6 0.45110 0.5018
## workers_alive 1 49.193 1568.4 0.30289 0.5821
## mean.pollen 1 50.217 1569.5 1.32744 0.2493
em1 <- update(em.mod, .~. -workers_alive)
## 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(em1, test = "Chisq")
## Single term deletions
##
## Model:
## emerge ~ fungicide + crithidia + dry_weight + live_weight + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.193 1568.4
## fungicide 1 52.544 1569.8 3.3507 0.06718 .
## crithidia 1 49.214 1566.5 0.0209 0.88504
## dry_weight 1 50.643 1567.9 1.4496 0.22859
## live_weight 1 49.705 1567.0 0.5125 0.47407
## mean.pollen 1 51.782 1569.0 2.5893 0.10759
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em2 <- update(em1, .~. -live_weight)
## 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(em2, test = "Chisq")
## Single term deletions
##
## Model:
## emerge ~ fungicide + crithidia + dry_weight + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 49.705 1567.0
## fungicide 1 53.137 1568.4 3.4317 0.06396 .
## crithidia 1 49.720 1565.0 0.0147 0.90334
## dry_weight 1 51.140 1566.4 1.4348 0.23098
## mean.pollen 1 52.707 1568.0 3.0014 0.08319 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em3 <- update(em2, .~. -dry_weight)
## 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(em3, test = "Chisq")
## Single term deletions
##
## Model:
## emerge ~ fungicide + crithidia + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.140 1566.4
## fungicide 1 54.822 1568.1 3.6820 0.05500 .
## crithidia 1 51.192 1564.4 0.0521 0.81952
## mean.pollen 1 54.631 1567.9 3.4905 0.06172 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em4 <- update(em3, .~. -mean.pollen)
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## Warning in glm.nb(formula = emerge ~ fungicide + crithidia, data = drones.na, :
## alternation limit reached
Anova(em4)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## fungicide 3.6718 1 0.05534 .
## crithidia 0.0054 1 0.94153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(em4)
##
## Call:
## glm.nb(formula = emerge ~ fungicide + crithidia, data = drones.na,
## init.theta = 3366696.841, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.86749 -0.30028 -0.02929 0.28290 1.97226
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.605679 0.015123 238.421 <2e-16 ***
## fungicide1 0.038133 0.019883 1.918 0.0551 .
## crithidia1 -0.001478 0.020155 -0.073 0.9415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3366782) family taken to be 1)
##
## Null deviance: 58.302 on 275 degrees of freedom
## Residual deviance: 54.631 on 273 degrees of freedom
## AIC: 1569.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3366697
## Std. Err.: 37068676
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -1561.876
emem <- emmeans(em4, pairwise ~ fungicide, type = "response")
emem
## $emmeans
## fungicide response SE df asymp.LCL asymp.UCL
## 0 36.8 0.489 Inf 35.8 37.8
## 1 38.2 0.575 Inf 37.1 39.4
##
## Results are averaged over the levels of: crithidia
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## fungicide0 / fungicide1 0.963 0.0191 Inf 1 -1.918 0.0551
##
## Results are averaged over the levels of: crithidia
## Tests are performed on the log scale
emcld <- cld(object = emem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
emcld
## fungicide response SE df asymp.LCL asymp.UCL .group
## 0 36.8 0.489 Inf 35.7 37.9 a
## 1 38.2 0.575 Inf 36.9 39.5 a
##
## Results are averaged over the levels of: crithidia
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
## Intervals are back-transformed from the log scale
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
pairs(emem)
## contrast ratio SE df null z.ratio p.value
## fungicide0 / fungicide1 0.963 0.0191 Inf 1 -1.918 0.0551
##
## Results are averaged over the levels of: crithidia
## Tests are performed on the log scale
em_sum <- drones %>%
group_by(treatment) %>%
summarise(m = mean(emerge),
sd = sd(emerge),
n = length(emerge)) %>%
mutate(se = sd/sqrt(n))
em_sum
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 37.0 2.41 102 0.238
## 2 2 38.1 2.54 69 0.306
## 3 3 38.4 4.08 49 0.583
## 4 4 36.5 2.23 61 0.286
em_sum$plot <- em_sum$m + em_sum$se
ggplot(em_sum, aes(x = treatment, y = m, fill = treatment)) +
geom_bar(stat = "identity", color = "black") +
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)
) +
scale_fill_viridis_d() +
geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Days") +
theme_classic(base_size = 20) +
coord_cartesian(ylim=c(35, 41)) +
annotate(geom = "text",
x = 1, y = 39,
label = "P = 0.05",
size = 8) +
theme(legend.position = "none",
axis.text = element_text(size = 20), # Set axis label font size
axis.title = element_text(size = 20)) + # Set axis title font size
theme(text = element_text(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 = 2, xend = 3, y = 39.6, yend = 39.6,
lineend = "round", linejoin = "round") +
geom_segment(x = 2, xend = 2, y = 39.4, yend = 39.8,
lineend = "round", linejoin = "round") +
geom_segment(x = 3, xend = 3, y = 39.4, yend = 39.8,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 4, y = 40.1, yend = 40.1,
lineend = "round", linejoin = "round") +
geom_segment(x = 1, xend = 1, y = 39.9, yend = 40.3,
lineend = "round", linejoin = "round") +
geom_segment(x = 4, xend = 4, y = 39.9, yend = 40.3,
lineend = "round", linejoin = "round") +
geom_text(x = 2.5, y = 39.65, label = "b", size = 6, vjust = -0.5) +
geom_text(x = 2.5, y = 40.15, label = "a", size = 6, vjust = -0.5)
p <- ggplot(qpcr, aes(x = days_since_innoculation, y = spores, color = colony)) +
geom_point() +
geom_line(aes(group = bee_id), alpha = 0.5) +
labs(title = "Spores per Bee Over Time",
x = "Time",
y = "Number of Spores",
color = "Colony") +
facet_wrap(~bee_id)
interactive_plot <- ggplotly(p)
interactive_plot
library(htmlwidgets)
saveWidget(interactive_plot, file = "interactive_plot.html")
unique(qpcr$colony)
## [1] T3.1 T3.10 T3.11 T3.12 T3.4 T3.6 T3.7 T3.8 T3.9 T4.1 T4.11 T4.12
## [13] T4.4 T4.6 T4.7 T4.8 T4.9
## 17 Levels: T3.1 T3.10 T3.11 T3.12 T3.4 T3.6 T3.7 T3.8 T3.9 T4.1 T4.11 ... T4.9
subset_t3.1 <- subset(qpcr, colony == "T3.1")
subset_t3.10 <- subset(qpcr, colony == "T3.10")
subset_t3.11 <- subset(qpcr, colony == "T3.11")
subset_t3.12 <- subset(qpcr, colony == "T3.12")
subset_t3.4 <- subset(qpcr, colony == "T3.4")
subset_t3.6 <- subset(qpcr, colony == "T3.6")
subset_t3.7 <- subset(qpcr, colony == "T3.7")
subset_t3.8 <- subset(qpcr, colony == "T3.8")
subset_t3.9 <- subset(qpcr, colony == "T3.9")
subset_t4.1 <- subset(qpcr, colony == "T4.1")
subset_t4.11 <- subset(qpcr, colony == "T4.11")
subset_t4.12 <- subset(qpcr, colony == "T4.12")
subset_t4.4 <- subset(qpcr, colony == "T4.4")
subset_t4.6 <- subset(qpcr, colony == "T4.6")
subset_t4.7 <- subset(qpcr, colony == "T4.7")
subset_t4.8 <- subset(qpcr, colony == "T4.8")
subset_t4.9 <- subset(qpcr, colony == "T4.9")
ggplot(subset_t3.1, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.10, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.11, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.12, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.4, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.6, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.7, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.8, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t3.9, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.1, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.11, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.12, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.4, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.6, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.7, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.8, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)
ggplot(subset_t4.9, aes(x = days_since_innoculation, y = spores, color = inoculate)) +
geom_point() +
coord_cartesian(ylim=c(0, 370)) +
geom_line(aes(group = bee_id), alpha = 0.5) +
geom_text(aes(label = spores), vjust = -1, size = 5) + # Adding data labels
labs(title = "Crithidia spores per bee over time",
x = "Days Since Inoculation",
y = "Number of Spores") +
facet_wrap(~bee_id) +
theme(legend.position = "none") +
theme_gray(base_size = 15)