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)
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("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"),
round = col_factor(levels = c("1",
"2", "3")), trial = col_skip()))
qpcr$colony <- as.factor(qpcr$colony)
qpcr$bee_id <- as.factor(qpcr$bee_id)
qpcr$spores <- as.double(qpcr$spores)
## 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)
pollen$count <- as.factor(pollen$`pollen ball id`)
#drones <- na.omit(drones)
#brood <- na.omit(brood)
shapiro.test(pollen$whole_dif)
##
## Shapiro-Wilk normality test
##
## data: pollen$whole_dif
## W = 0.77483, p-value < 2.2e-16
hist(pollen$whole_dif)
range(pollen$whole_dif)
## [1] 0.03316 1.39545
pollen$box <- bcPower(pollen$whole_dif, -4, gamma=1)
shapiro.test(pollen$box)
##
## Shapiro-Wilk normality test
##
## data: pollen$box
## W = 0.94829, p-value = 1.356e-15
hist(pollen$box)
pollen$log <- log(pollen$whole_dif)
shapiro.test(pollen$log)
##
## Shapiro-Wilk normality test
##
## data: pollen$log
## W = 0.93499, 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.62065, 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.86233, p-value < 2.2e-16
hist(pollen$root)
descdist(pollen$whole_dif, discrete = FALSE)
## summary statistics
## ------
## min: 0.03316 max: 1.39545
## median: 0.27686
## mean: 0.4042769
## estimated sd: 0.3016975
## estimated skewness: 1.532669
## estimated kurtosis: 4.289505
ggplot(pollen, aes(x = box, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.01, 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)")
pollen$count <- as.numeric(pollen$count)
pollen$fungicide <- as.logical(pollen$fungicide)
pollen$crithidia <- as.logical(pollen$crithidia)
pollen$id <- as.factor(pollen$`pollen ball id`)
pol.mod <- lmer(whole_dif ~ fungicide*crithidia + id + block + count + workers_alive + (1|colony), data = pollen)
drop1(pol.mod, test = "Chisq")
## Single term deletions
##
## Model:
## whole_dif ~ fungicide * crithidia + id + block + count + workers_alive +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -460.06
## id 23 -286.07 219.985 < 2.2e-16 ***
## block 8 -435.85 40.209 2.929e-06 ***
## count 0 -460.06 0.000
## workers_alive 1 -404.13 57.931 2.714e-14 ***
## fungicide:crithidia 1 -462.02 0.035 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pm1 <- update(pol.mod, .~. -count)
drop1(pm1, test = "Chisq")
## Single term deletions
##
## Model:
## whole_dif ~ fungicide + crithidia + id + block + workers_alive +
## (1 | colony) + fungicide:crithidia
## npar AIC LRT Pr(Chi)
## <none> -460.06
## id 24 -86.35 421.71 < 2.2e-16 ***
## block 8 -435.85 40.21 2.929e-06 ***
## workers_alive 1 -404.13 57.93 2.714e-14 ***
## fungicide:crithidia 1 -462.02 0.04 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pol.mod1 <- lmer(whole_dif ~ fungicide + crithidia + block + id + workers_alive + (1|colony), data = pollen)
drop1(pol.mod1, test = "Chisq")
## Single term deletions
##
## Model:
## whole_dif ~ fungicide + crithidia + block + id + workers_alive +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -462.02
## fungicide 1 -462.21 1.82 0.17744
## crithidia 1 -460.01 4.02 0.04506 *
## block 8 -437.84 40.19 2.957e-06 ***
## id 24 -88.33 421.70 < 2.2e-16 ***
## workers_alive 1 -406.11 57.91 2.738e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(pol.mod1));qqline(resid(pol.mod1))
#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.464 0.0319 27.9 0.399 0.530
## TRUE 0.388 0.0320 28.2 0.323 0.454
##
## Results are averaged over the levels of: fungicide, block, id
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## FALSE - TRUE 0.0758 0.044 25.4 1.722 0.0973
##
## Results are averaged over the levels of: fungicide, block, id
## Degrees-of-freedom method: kenward-roger
library(glmmTMB)
pol.glmm1 <- glmmTMB(whole_dif ~ fungicide*crithidia + workers_alive + block + id + (1|colony), data = pollen)
summary(pol.glmm1)
## Family: gaussian ( identity )
## Formula:
## whole_dif ~ fungicide * crithidia + workers_alive + block + id +
## (1 | colony)
## Data: pollen
##
## AIC BIC logLik deviance df.resid
## -460.1 -279.6 269.0 -538.1 716
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.01066 0.1032
## Residual 0.02577 0.1605
## Number of obs: 755, groups: colony, 36
##
## Dispersion estimate for gaussian family (sigma^2): 0.0258
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.181489 0.090823 -1.998 0.045687 *
## fungicideTRUE -0.042846 0.051573 -0.831 0.406094
## crithidiaTRUE -0.068435 0.051804 -1.321 0.186489
## workers_alive 0.083882 0.010786 7.777 7.42e-15 ***
## block4 0.433802 0.077369 5.607 2.06e-08 ***
## block6 -0.117723 0.076878 -1.531 0.125694
## block7 0.084247 0.077590 1.086 0.277572
## block8 0.164267 0.077256 2.126 0.033482 *
## block9 0.036582 0.077197 0.474 0.635589
## block10 0.266068 0.077299 3.442 0.000577 ***
## block11 -0.027382 0.077413 -0.354 0.723553
## block12 0.099556 0.077150 1.290 0.196904
## id3 -0.038660 0.037864 -1.021 0.307244
## id4 -0.036851 0.038449 -0.958 0.337851
## id5 -0.040024 0.038150 -1.049 0.294131
## id6 -0.051089 0.038163 -1.339 0.180658
## id7 -0.007163 0.037940 -0.189 0.850253
## id8 0.042733 0.037974 1.125 0.260448
## id9 0.092730 0.038025 2.439 0.014742 *
## id10 0.202909 0.038054 5.332 9.71e-08 ***
## id11 0.346610 0.038193 9.075 < 2e-16 ***
## id12 0.389216 0.038040 10.232 < 2e-16 ***
## id13 0.359250 0.038369 9.363 < 2e-16 ***
## id14 0.309928 0.038582 8.033 9.51e-16 ***
## id15 0.293200 0.038582 7.599 2.97e-14 ***
## id16 0.293770 0.038641 7.603 2.90e-14 ***
## id17 0.292566 0.038765 7.547 4.45e-14 ***
## id18 0.314074 0.039040 8.045 8.63e-16 ***
## id19 0.306038 0.039190 7.809 5.76e-15 ***
## id20 0.288617 0.039517 7.304 2.80e-13 ***
## id21 0.247867 0.043494 5.699 1.21e-08 ***
## id22 0.263248 0.052692 4.996 5.85e-07 ***
## id23 0.283535 0.057731 4.911 9.05e-07 ***
## id24 0.298766 0.062538 4.777 1.78e-06 ***
## id25 0.312301 0.065417 4.774 1.81e-06 ***
## id26 0.161324 0.167332 0.964 0.334997
## fungicideTRUE:crithidiaTRUE -0.013680 0.072800 -0.188 0.850941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(pol.glmm1, test = "Chisq")
## Single term deletions
##
## Model:
## whole_dif ~ fungicide * crithidia + workers_alive + block + id +
## (1 | colony)
## Df AIC LRT Pr(>Chi)
## <none> -460.06
## workers_alive 1 -404.13 57.93 2.714e-14 ***
## block 8 -435.85 40.21 2.929e-06 ***
## id 24 -86.35 421.71 < 2.2e-16 ***
## fungicide:crithidia 1 -462.02 0.04 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pol.glmm <- glmmTMB(whole_dif ~ fungicide + crithidia + workers_alive + block + id + (1|colony), data = pollen)
summary(pol.glmm)
## Family: gaussian ( identity )
## Formula:
## whole_dif ~ fungicide + crithidia + workers_alive + block + id +
## (1 | colony)
## Data: pollen
##
## AIC BIC logLik deviance df.resid
## -462.0 -286.2 269.0 -538.0 717
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.01067 0.1033
## Residual 0.02577 0.1605
## Number of obs: 755, groups: colony, 36
##
## Dispersion estimate for gaussian family (sigma^2): 0.0258
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.177568 0.088424 -2.008 0.044628 *
## fungicideTRUE -0.049708 0.036436 -1.364 0.172488
## crithidiaTRUE -0.075321 0.036644 -2.055 0.039834 *
## workers_alive 0.083791 0.010774 7.777 7.44e-15 ***
## block4 0.433844 0.077408 5.605 2.09e-08 ***
## block6 -0.117707 0.076918 -1.530 0.125943
## block7 0.084195 0.077630 1.085 0.278111
## block8 0.164291 0.077296 2.125 0.033547 *
## block9 0.036601 0.077237 0.474 0.635592
## block10 0.266073 0.077339 3.440 0.000581 ***
## block11 -0.027412 0.077453 -0.354 0.723404
## block12 0.099555 0.077190 1.290 0.197140
## id3 -0.038664 0.037864 -1.021 0.307189
## id4 -0.036878 0.038449 -0.959 0.337494
## id5 -0.040040 0.038150 -1.050 0.293939
## id6 -0.051128 0.038162 -1.340 0.180322
## id7 -0.007217 0.037939 -0.190 0.849135
## id8 0.042690 0.037973 1.124 0.260924
## id9 0.092682 0.038024 2.437 0.014791 *
## id10 0.202858 0.038053 5.331 9.77e-08 ***
## id11 0.346550 0.038192 9.074 < 2e-16 ***
## id12 0.389160 0.038039 10.231 < 2e-16 ***
## id13 0.359180 0.038367 9.362 < 2e-16 ***
## id14 0.309847 0.038579 8.031 9.64e-16 ***
## id15 0.293119 0.038579 7.598 3.01e-14 ***
## id16 0.293687 0.038638 7.601 2.94e-14 ***
## id17 0.292479 0.038762 7.546 4.51e-14 ***
## id18 0.313975 0.039036 8.043 8.75e-16 ***
## id19 0.305936 0.039186 7.807 5.85e-15 ***
## id20 0.288504 0.039512 7.302 2.84e-13 ***
## id21 0.247715 0.043487 5.696 1.22e-08 ***
## id22 0.263106 0.052687 4.994 5.92e-07 ***
## id23 0.283290 0.057716 4.908 9.19e-07 ***
## id24 0.298489 0.062521 4.774 1.80e-06 ***
## id25 0.312040 0.065402 4.771 1.83e-06 ***
## id26 0.160790 0.167306 0.961 0.336524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(pol.glmm, test = "Chisq")
## Single term deletions
##
## Model:
## whole_dif ~ fungicide + crithidia + workers_alive + block + id +
## (1 | colony)
## Df AIC LRT Pr(>Chi)
## <none> -462.02
## fungicide 1 -462.21 1.82 0.17744
## crithidia 1 -460.01 4.02 0.04506 *
## workers_alive 1 -406.11 57.91 2.738e-14 ***
## block 8 -437.84 40.19 2.957e-06 ***
## id 24 -88.33 421.70 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(pol.glmm));qqline(resid(pol.glmm))
pe <- emmeans(pol.glmm, pairwise ~ crithidia, type = "response")
pe
## $emmeans
## crithidia emmean SE df lower.CL upper.CL
## FALSE 0.464 0.0268 717 0.411 0.517
## TRUE 0.389 0.0269 717 0.336 0.441
##
## Results are averaged over the levels of: fungicide, block, id
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## FALSE - TRUE 0.0753 0.0366 717 2.055 0.0402
##
## Results are averaged over the levels of: fungicide, block, id
AIC(pol.glmm, pol.glmm1)
## df AIC
## pol.glmm 38 -462.0243
## pol.glmm1 39 -460.0596
Anova(pol.glmm)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: whole_dif
## Chisq Df Pr(>Chisq)
## fungicide 1.8612 1 0.17249
## crithidia 4.2249 1 0.03983 *
## workers_alive 60.4788 1 7.437e-15 ***
## block 73.4485 8 1.008e-12 ***
## id 566.6711 24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(data = pollen, aes(x = count, 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.376 0.280 188 0.0204
pollen_box_sum
## # A tibble: 4 × 5
## treatment mean sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.178 0.0430 184 0.00317
## 2 2 0.167 0.0428 188 0.00312
## 3 3 0.152 0.0380 195 0.00272
## 4 4 0.160 0.0418 188 0.00305
pollen_sum$plot <- pollen_sum$mean + pollen_sum$se
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)
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, 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.6516 1 0.41956
## crithidia 2.9046 1 0.08832 .
## mean.pollen 18.2550 1 1.932e-05 ***
## block 19.1268 8 0.01420 *
## days_active 1.4265 1 0.23234
## fungicide:crithidia 1.3554 1 0.24434
## ---
## 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 24.125 91.566
## mean.pollen 1 42.380 107.821 18.2550 1.932e-05 ***
## block 8 43.252 94.692 19.1268 0.0142 *
## days_active 1 25.551 90.992 1.4265 0.2323
## fungicide:crithidia 1 25.480 90.921 1.3554 0.2443
## ---
## 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, 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.6516 1 0.41956
## crithidia 2.9046 1 0.08832 .
## mean.pollen 18.4715 1 1.725e-05 ***
## block 19.6722 8 0.01165 *
## days_active 0.8626 1 0.35301
## ---
## 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 25.480 90.921
## fungicide 1 26.132 89.573 0.6516 0.41956
## crithidia 1 28.385 91.826 2.9046 0.08832 .
## mean.pollen 1 43.952 107.392 18.4715 1.725e-05 ***
## block 8 45.152 94.593 19.6722 0.01165 *
## days_active 1 26.343 89.784 0.8626 0.35301
## ---
## 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 26.343 89.784
## fungicide 1 27.282 88.723 0.9395 0.332394
## crithidia 1 30.172 91.613 3.8291 0.050370 .
## mean.pollen 1 52.833 114.274 26.4903 2.649e-07 ***
## block 8 47.521 94.962 21.1782 0.006689 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(cbw2)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.9395 1 0.332394
## crithidia 3.8291 1 0.050370 .
## mean.pollen 26.4903 1 2.649e-07 ***
## block 21.1782 8 0.006689 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(cbw2));qqline(resid(cbw2))
Anova(cbw2)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(workers_alive, workers_dead)
## LR Chisq Df Pr(>Chisq)
## fungicide 0.9395 1 0.332394
## crithidia 3.8291 1 0.050370 .
## mean.pollen 26.4903 1 2.649e-07 ***
## block 21.1782 8 0.006689 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(cbw2)
##
## 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.9252 -0.5366 0.1348 0.6577 1.5162
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1264 1.0823 -1.965 0.0495 *
## fungicideTRUE 0.4218 0.4386 0.962 0.3362
## crithidiaTRUE -0.8727 0.4483 -1.947 0.0516 .
## mean.pollen 10.4521 2.5014 4.179 2.93e-05 ***
## block4 14.8345 3600.7582 0.004 0.9967
## block6 0.3860 0.7888 0.489 0.6246
## block7 -1.1692 0.7926 -1.475 0.1402
## block8 0.6117 0.8994 0.680 0.4965
## block9 0.9506 0.9531 0.997 0.3186
## block10 -2.3981 0.9979 -2.403 0.0163 *
## block11 -0.1471 0.7723 -0.191 0.8489
## block12 -1.9234 0.9288 -2.071 0.0384 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 104.128 on 35 degrees of freedom
## Residual deviance: 26.343 on 24 degrees of freedom
## AIC: 89.784
##
## 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.89 1.96
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.89 1.96 9 0.655
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 = 1,
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")
workers$fungicide <- as.logical(workers$fungicide)
workers$crithidia <- as.logical(workers$crithidia)
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> 1420.9
## avg_pollen 1 1420.5 1.6542 0.1984
## inoculate 1 1418.9 0.0009 0.9755
## block 8 1416.8 11.8861 0.1564
## fungicide:crithidia 1 1419.6 0.7004 0.4027
dayswrk <- glmer.nb(days_alive ~ fungicide + crithidia + block + inoculate + (1|colony), data = workers)
drop1(dayswrk, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + inoculate + (1 |
## colony)
## npar AIC LRT Pr(Chi)
## <none> 1419.2
## fungicide 1 1417.3 0.1000 0.7518
## crithidia 1 1419.6 2.3437 0.1258
## block 8 1413.5 10.2617 0.2471
## inoculate 1 1417.2 0.0012 0.9719
dayswrk1 <- update(dayswrk, .~. -inoculate)
drop1(dayswrk1, test = "Chisq")
## Single term deletions
##
## Model:
## days_alive ~ fungicide + crithidia + block + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 1417.2
## fungicide 1 1415.3 0.1013 0.7503
## crithidia 1 1417.6 2.3458 0.1256
## block 8 1411.5 10.2611 0.2472
Anova(dayswrk1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: days_alive
## Chisq Df Pr(>Chisq)
## fungicide 0.1013 1 0.7503
## crithidia 2.3459 1 0.1256
## block 10.2518 8 0.2478
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.251 218.69
## fungicide 1 11.254 216.70 0.00314 0.9553
## crithidia 1 11.870 217.31 0.61882 0.4315
## mean.pollen 1 12.311 217.75 1.05948 0.3033
## days_first_ov 1 12.169 217.61 0.91814 0.3380
durmod <- glm.nb(days_active ~ fungicide + crithidia + mean.pollen, 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 12.230 223.43
## fungicide 1 12.233 221.43 0.0027 0.95862
## crithidia 1 12.726 221.92 0.4959 0.48131
## mean.pollen 1 16.626 225.82 4.3955 0.03603 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(durmod)
##
## Call:
## glm.nb(formula = days_active ~ fungicide + crithidia + mean.pollen,
## data = duration, init.theta = 2545370.79, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2314 -0.2677 0.1415 0.2977 1.1130
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.895011 0.073638 52.894 <2e-16 ***
## fungicideTRUE -0.002607 0.050247 -0.052 0.9586
## crithidiaTRUE 0.036099 0.051271 0.704 0.4814
## mean.pollen -0.250560 0.120236 -2.084 0.0372 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2545371) family taken to be 1)
##
## Null deviance: 18.303 on 35 degrees of freedom
## Residual deviance: 12.230 on 32 degrees of freedom
## AIC: 225.43
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2545371
## Std. Err.: 66958888
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -215.427
Anova(durmod)
## Analysis of Deviance Table (Type II tests)
##
## Response: days_active
## LR Chisq Df Pr(>Chisq)
## fungicide 0.0027 1 0.95862
## crithidia 0.4959 1 0.48131
## mean.pollen 4.3955 1 0.03603 *
## ---
## 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
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
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> 17.299 193.75
## avg.pol 1 20.216 194.67 2.9167 0.08767 .
## workers_alive 1 18.139 192.59 0.8399 0.35944
## block 8 26.334 186.79 9.0353 0.33933
## fungicide:crithidia 1 17.727 192.18 0.4283 0.51283
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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> 17.727 192.18
## fungicide 1 17.835 190.29 0.1076 0.7429
## crithidia 1 19.227 191.68 1.4995 0.2208
## avg.pol 1 20.392 192.84 2.6647 0.1026
## workers_alive 1 18.797 191.25 1.0696 0.3010
## block 8 26.773 185.23 9.0457 0.3385
ov1 <- update(ov, .~. -block)
## 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(ov1, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide + crithidia + avg.pol + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 26.772 185.23
## fungicide 1 26.878 183.33 0.1053 0.74551
## crithidia 1 28.305 184.76 1.5327 0.21571
## avg.pol 1 32.704 189.16 5.9317 0.01487 *
## workers_alive 1 28.931 185.38 2.1591 0.14173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ov2 <- update(ov1, .~. -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(ov2, test = "Chisq")
## Single term deletions
##
## Model:
## days_first_ov ~ fungicide + crithidia + avg.pol
## Df Deviance AIC LRT Pr(>Chi)
## <none> 28.931 185.39
## fungicide 1 29.045 183.50 0.1146 0.7350
## crithidia 1 29.519 183.97 0.5882 0.4431
## avg.pol 1 49.820 204.27 20.8889 4.867e-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.1146 1 0.7350
## crithidia 0.5882 1 0.4431
## avg.pol 20.8889 1 4.867e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ov2)
##
## Call:
## glm.nb(formula = days_first_ov ~ fungicide + crithidia + avg.pol,
## data = duration, init.theta = 151930.6281, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.71248 -0.65163 -0.07036 0.53669 1.82330
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.95946 0.14611 20.255 < 2e-16 ***
## fungicideTRUE -0.03434 0.10142 -0.339 0.735
## crithidiaTRUE -0.07821 0.10197 -0.767 0.443
## avg.pol -1.14088 0.25754 -4.430 9.43e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(151930.6) family taken to be 1)
##
## Null deviance: 50.159 on 34 degrees of freedom
## Residual deviance: 28.931 on 31 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 187.39
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 151931
## Std. Err.: 5542697
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -177.386
plot(duration$treatment, duration$days_first_ov)
brood$fungicide <- as.logical(brood$fungicide)
brood$crithidia <- as.logical(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 9.473 1 0.002085 **
## crithidia 0.483 1 0.487155
## block 87.471 8 1.515e-15 ***
## workers_alive 35.058 1 3.200e-09 ***
## duration 0.253 1 0.615167
## avg_pollen 19.007 1 1.302e-05 ***
## fungicide:crithidia 0.286 1 0.592710
## ---
## 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> 49.836 249.43
## block 8 137.307 320.90 87.471 1.515e-15 ***
## workers_alive 1 84.894 282.48 35.058 3.200e-09 ***
## duration 1 50.089 247.68 0.253 0.6152
## avg_pollen 1 68.843 266.43 19.007 1.302e-05 ***
## fungicide:crithidia 1 50.122 247.71 0.286 0.5927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
brood.mod <- glm.nb(brood_cells ~ fungicide + crithidia + block + workers_alive + duration, data = brood)
Anova(brood.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 3.316 1 0.06862 .
## crithidia 0.007 1 0.93175
## block 102.672 8 < 2.2e-16 ***
## workers_alive 64.517 1 9.572e-16 ***
## duration 0.751 1 0.38607
## ---
## 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 45.437 259.54
## fungicide 1 48.753 260.86 3.316 0.06862 .
## crithidia 1 45.445 257.55 0.007 0.93175
## block 8 148.110 346.22 102.672 < 2.2e-16 ***
## workers_alive 1 109.954 322.06 64.517 9.572e-16 ***
## duration 1 46.189 258.29 0.751 0.38607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bm1 <- update(brood.mod, .~. -duration)
drop1(bm1, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + crithidia + block + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 46.20 258.29
## fungicide 1 49.27 259.36 3.070 0.07973 .
## crithidia 1 46.21 256.30 0.010 0.91928
## block 8 148.91 345.00 102.709 < 2e-16 ***
## workers_alive 1 139.03 349.13 92.831 < 2e-16 ***
## ---
## 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, data = brood, init.theta = 10.1805701, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.61598 -1.10509 0.04918 0.40418 2.05004
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.80904 0.38129 2.122 0.033852 *
## fungicideTRUE 0.26020 0.14524 1.792 0.073211 .
## crithidiaTRUE -0.01541 0.14808 -0.104 0.917136
## block4 -0.23235 0.26721 -0.870 0.384542
## block6 -1.78611 0.38966 -4.584 4.57e-06 ***
## block7 -0.63042 0.29537 -2.134 0.032816 *
## block8 -1.01617 0.28261 -3.596 0.000324 ***
## block9 -1.04473 0.27873 -3.748 0.000178 ***
## block10 0.50706 0.25340 2.001 0.045390 *
## block11 -1.70311 0.36066 -4.722 2.33e-06 ***
## block12 0.12012 0.29090 0.413 0.679659
## workers_alive 0.68405 0.07325 9.339 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(10.1806) family taken to be 1)
##
## Null deviance: 305.7 on 35 degrees of freedom
## Residual deviance: 46.2 on 24 degrees of freedom
## AIC: 260.29
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 10.18
## Std. Err.: 4.15
##
## 2 x log-likelihood: -234.294
bm2 <- update(bm1, .~. -crithidia)
drop1(bm2, test = "Chisq")
## Single term deletions
##
## Model:
## brood_cells ~ fungicide + block + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 46.193 256.30
## fungicide 1 49.298 257.41 3.105 0.07806 .
## block 8 149.523 343.63 103.330 < 2e-16 ***
## workers_alive 1 153.451 361.56 107.258 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(bm1)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## fungicide 3.070 1 0.07973 .
## crithidia 0.010 1 0.91928
## block 102.709 8 < 2e-16 ***
## workers_alive 92.831 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(resid(bm1));qqline(resid(bm1))
broodem <- emmeans(bm1, 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
## TRUE 13.7 1.58 Inf 10.6 17.7 a
## FALSE 13.9 1.54 Inf 10.9 17.8 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.
brood_sum <- brood %>%
group_by(treatment) %>%
summarise(mb = mean(brood_cells),
nb = length(brood_cells),
sdb = sd(brood_cells)) %>%
mutate(seb = (sdb/sqrt(nb)))
brood_sum
## # A tibble: 4 × 5
## treatment mb nb sdb seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 30.7 9 23.1 7.71
## 2 2 29.8 9 21.1 7.04
## 3 3 21.7 9 24.9 8.30
## 4 4 23.9 9 26.3 8.77
brood_sum$plot <- brood_sum$mb + brood_sum$seb
plot(brood$treatment, brood$brood_cells)
ggplot(data = brood_sum, aes(x = treatment, y = mb, 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, 42)) +
geom_errorbar(aes(ymin = mb - seb, ymax = mb + seb), 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 = 41,
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")
#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.967 153.63
## workers_alive 1 47.737 164.40 12.7701 0.0003522 ***
## block 8 46.453 149.12 11.4865 0.1756262
## duration 1 39.330 156.00 4.3633 0.0367209 *
## avg_pollen 1 58.179 174.85 23.2126 1.45e-06 ***
## fungicide:crithidia 1 35.089 151.76 0.1227 0.7260822
## ---
## 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.2574 1 0.611898
## crithidia 3.1498 1 0.075934 .
## workers_alive 12.7658 1 0.000353 ***
## block 11.4846 8 0.175721
## duration 4.3615 1 0.036760 *
## avg_pollen 23.2025 1 1.458e-06 ***
## fungicide:crithidia 0.1222 1 0.726688
## ---
## 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.959 153.64
## workers_alive 1 47.725 164.40 12.7658 0.000353 ***
## block 8 46.444 149.12 11.4846 0.175721
## duration 1 39.321 156.00 4.3615 0.036760 *
## avg_pollen 1 58.162 174.84 23.2025 1.458e-06 ***
## fungicide:crithidia 1 35.081 151.76 0.1222 0.726688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
livepup.mod <- glm(live_pupae ~ fungicide + crithidia ++ workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
livepup.mod.nb <- 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 + :
## alternation limit reached
summary(livepup.mod)
##
## Call:
## glm(formula = live_pupae ~ fungicide + crithidia + +workers_alive +
## block + duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0545 -0.8359 -0.1715 0.4672 2.2654
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.56936 2.71700 -2.418 0.015612 *
## fungicideTRUE -0.08776 0.17284 -0.508 0.611621
## crithidiaTRUE -0.30903 0.17470 -1.769 0.076907 .
## workers_alive 0.52057 0.15220 3.420 0.000625 ***
## block4 -0.65692 0.43189 -1.521 0.128253
## block6 -2.18079 1.05286 -2.071 0.038330 *
## block7 -0.12475 0.38798 -0.322 0.747807
## block8 -0.61961 0.42220 -1.468 0.142222
## block9 -0.34465 0.38437 -0.897 0.369898
## block10 -0.81683 0.44468 -1.837 0.066226 .
## block11 -0.80235 0.52146 -1.539 0.123887
## block12 -0.49818 0.47864 -1.041 0.297956
## duration 0.10408 0.05153 2.020 0.043413 *
## avg_pollen 4.11880 0.90588 4.547 5.45e-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.752 on 35 degrees of freedom
## Residual deviance: 35.089 on 22 degrees of freedom
## AIC: 151.76
##
## 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 14 151.7561
## livepup.mod.nb 15 153.7595
anova(livepup.mod, livepup.mod.nb, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: live_pupae ~ fungicide + crithidia + +workers_alive + block +
## duration + avg_pollen
## Model 2: live_pupae ~ fungicide + crithidia + +workers_alive + block +
## duration + avg_pollen
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 22 35.089
## 2 22 35.081 0 0.0084061
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> 35.089 151.76
## fungicide 1 35.347 150.01 0.2579 0.6115536
## crithidia 1 38.243 152.91 3.1538 0.0757492 .
## workers_alive 1 47.746 162.41 12.6561 0.0003743 ***
## block 8 46.667 147.33 11.5771 0.1710932
## duration 1 39.519 154.19 4.4296 0.0353208 *
## avg_pollen 1 58.189 172.85 23.0991 1.539e-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 + duration +
## avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 46.667 147.33
## fungicide 1 47.277 145.94 0.610 0.4347
## crithidia 1 48.451 147.12 1.784 0.1816
## workers_alive 1 67.688 166.35 21.021 4.542e-06 ***
## duration 1 51.610 150.28 4.943 0.0262 *
## avg_pollen 1 104.321 202.99 57.655 3.124e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lp1)
##
## Call:
## glm(formula = live_pupae ~ fungicide + crithidia + workers_alive +
## duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2883 -0.9998 -0.4256 0.4991 3.2325
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.82576 1.52462 -3.165 0.00155 **
## fungicideTRUE -0.11609 0.14903 -0.779 0.43599
## crithidiaTRUE -0.20712 0.15674 -1.321 0.18634
## workers_alive 0.51273 0.12063 4.251 2.13e-05 ***
## duration 0.06014 0.02677 2.247 0.02467 *
## avg_pollen 3.44614 0.48583 7.093 1.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 252.752 on 35 degrees of freedom
## Residual deviance: 46.667 on 30 degrees of freedom
## AIC: 147.33
##
## Number of Fisher Scoring iterations: 5
Anova(lp1)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.610 1 0.4347
## crithidia 1.784 1 0.1816
## workers_alive 21.021 1 4.542e-06 ***
## duration 4.943 1 0.0262 *
## avg_pollen 57.655 1 3.124e-14 ***
## ---
## 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 + block +
## duration + avg_pollen
## Model 2: live_pupae ~ fungicide + crithidia + workers_alive + duration +
## avg_pollen
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 22 35.089
## 2 30 46.667 -8 -11.577 0.1711
AIC(livepup.mod, lp1)
## df AIC
## livepup.mod 14 151.7561
## lp1 6 147.3332
be <- emmeans(lp1, "crithidia")
pairs(be)
## contrast estimate SE df z.ratio p.value
## FALSE - TRUE 0.207 0.157 Inf 1.321 0.1863
##
## 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")
broodcld <- cld(object = broodem,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
broodcld
## crithidia rate SE df asymp.LCL asymp.UCL .group
## TRUE 2.27 0.370 Inf 1.57 3.27 a
## FALSE 2.79 0.429 Inf 1.98 3.93 a
##
## Results are averaged over the levels of: fungicide
## 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.
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 4 9 5.17 1.72
livepup_sum$plot <- livepup_sum$mb + livepup_sum$seb
plot(brood$treatment, brood$live_pupae)
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, 0.4), # 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 = 1,
y = 12.5,
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")
plot(brood$treatment, brood$live_larvae)
livelar.mod <- glm(live_larvae ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
summary(livelar.mod) #overdisp
##
## Call:
## glm(formula = live_larvae ~ fungicide + crithidia + workers_alive +
## block + duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6068 -1.6238 -0.5069 1.0864 2.5547
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.90936 1.11950 1.706 0.088094 .
## fungicideTRUE 0.19053 0.09850 1.934 0.053080 .
## crithidiaTRUE 0.18432 0.09419 1.957 0.050357 .
## workers_alive 0.34362 0.07345 4.678 2.89e-06 ***
## block4 -0.35446 0.23717 -1.495 0.135039
## block6 -1.82877 0.60694 -3.013 0.002586 **
## block7 -0.55917 0.23774 -2.352 0.018674 *
## block8 -0.90571 0.24315 -3.725 0.000195 ***
## block9 -0.66470 0.22504 -2.954 0.003140 **
## block10 0.75755 0.22272 3.401 0.000671 ***
## block11 -1.31312 0.38511 -3.410 0.000650 ***
## block12 0.32055 0.23934 1.339 0.180472
## duration -0.03840 0.02207 -1.740 0.081935 .
## avg_pollen 2.24259 0.48644 4.610 4.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 690.91 on 35 degrees of freedom
## Residual deviance: 91.83 on 22 degrees of freedom
## AIC: 237.57
##
## Number of Fisher Scoring iterations: 5
livelar.mod.int <- glm(live_larvae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
summary(livelar.mod.int) #overdisp
##
## Call:
## glm(formula = live_larvae ~ fungicide * crithidia + workers_alive +
## block + duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6955 -1.6402 -0.5423 1.1193 2.2554
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.17431 1.15107 1.889 0.058898 .
## fungicideTRUE 0.12449 0.12521 0.994 0.320099
## crithidiaTRUE 0.11112 0.12715 0.874 0.382157
## workers_alive 0.34537 0.07376 4.682 2.83e-06 ***
## block4 -0.36849 0.23880 -1.543 0.122817
## block6 -1.83346 0.60674 -3.022 0.002513 **
## block7 -0.56421 0.23834 -2.367 0.017924 *
## block8 -0.91356 0.24198 -3.775 0.000160 ***
## block9 -0.70003 0.22921 -3.054 0.002257 **
## block10 0.76821 0.22112 3.474 0.000512 ***
## block11 -1.35595 0.38821 -3.493 0.000478 ***
## block12 0.29780 0.24254 1.228 0.219513
## duration -0.04277 0.02246 -1.904 0.056896 .
## avg_pollen 2.16957 0.49298 4.401 1.08e-05 ***
## fungicideTRUE:crithidiaTRUE 0.15563 0.18144 0.858 0.391032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 690.912 on 35 degrees of freedom
## Residual deviance: 91.094 on 21 degrees of freedom
## AIC: 238.84
##
## Number of Fisher Scoring iterations: 5
livelar.mod.nb.int <- glm.nb(live_larvae ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen, 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, data = brood, init.theta = 9.07466479,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9298 -1.2359 -0.3805 0.5865 2.1165
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.36901 1.86239 0.735 0.462291
## fungicideTRUE 0.12018 0.22743 0.528 0.597224
## crithidiaTRUE -0.02186 0.24449 -0.089 0.928765
## workers_alive 0.43366 0.12208 3.552 0.000382 ***
## block4 -0.65537 0.43381 -1.511 0.130857
## block6 -1.75205 0.68248 -2.567 0.010252 *
## block7 -0.65454 0.39712 -1.648 0.099312 .
## block8 -1.20924 0.39409 -3.068 0.002152 **
## block9 -0.75635 0.35561 -2.127 0.033425 *
## block10 0.70621 0.37204 1.898 0.057666 .
## block11 -1.48747 0.51072 -2.912 0.003586 **
## block12 0.37671 0.39256 0.960 0.337239
## duration -0.03602 0.03572 -1.009 0.313178
## avg_pollen 2.70545 0.87460 3.093 0.001979 **
## fungicideTRUE:crithidiaTRUE 0.27237 0.34618 0.787 0.431409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(9.0747) family taken to be 1)
##
## Null deviance: 311.724 on 35 degrees of freedom
## Residual deviance: 51.376 on 21 degrees of freedom
## AIC: 229.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 9.07
## Std. Err.: 5.05
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -197.70
livelar.mod.nb <- glm.nb(live_larvae ~ fungicide + crithidia + workers_alive + block + duration + avg_pollen, 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, data = brood, init.theta = 8.987768981,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8127 -1.1505 -0.4229 0.5430 2.1928
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.98943 1.81342 0.546 0.58533
## fungicideTRUE 0.23197 0.17923 1.294 0.19557
## crithidiaTRUE 0.10121 0.18221 0.555 0.57858
## workers_alive 0.42852 0.12171 3.521 0.00043 ***
## block4 -0.60665 0.43116 -1.407 0.15942
## block6 -1.77937 0.68529 -2.597 0.00942 **
## block7 -0.63088 0.39774 -1.586 0.11270
## block8 -1.16564 0.39233 -2.971 0.00297 **
## block9 -0.71906 0.35055 -2.051 0.04024 *
## block10 0.70056 0.37496 1.868 0.06171 .
## block11 -1.42670 0.50319 -2.835 0.00458 **
## block12 0.38682 0.39093 0.989 0.32243
## duration -0.02891 0.03499 -0.826 0.40872
## avg_pollen 2.73864 0.87600 3.126 0.00177 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(8.9878) family taken to be 1)
##
## Null deviance: 310.310 on 35 degrees of freedom
## Residual deviance: 51.831 on 22 degrees of freedom
## AIC: 228.31
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 8.99
## Std. Err.: 5.08
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -198.31
drop1(livelar.mod.nb.int, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide * crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.376 227.70
## workers_alive 1 63.882 238.21 12.506 0.0004056 ***
## block 8 107.289 267.61 55.913 2.933e-09 ***
## duration 1 52.267 226.59 0.891 0.3452583
## avg_pollen 1 60.587 234.91 9.211 0.0024050 **
## fungicide:crithidia 1 51.987 226.31 0.611 0.4344570
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(livelar.mod.nb, test = "Chisq")
## Single term deletions
##
## Model:
## live_larvae ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.831 226.31
## fungicide 1 53.501 225.98 1.670 0.1961955
## crithidia 1 52.126 224.60 0.295 0.5868366
## workers_alive 1 64.093 236.57 12.262 0.0004623 ***
## block 8 107.087 265.57 55.256 3.936e-09 ***
## duration 1 52.440 224.92 0.609 0.4352187
## avg_pollen 1 61.252 233.73 9.421 0.0021450 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ll1 <- update(livelar.mod.nb, .~. -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> 51.138 224.90
## fungicide 1 52.618 224.38 1.480 0.223743
## crithidia 1 51.325 223.09 0.188 0.664965
## workers_alive 1 66.612 238.38 15.474 8.363e-05 ***
## block 8 107.238 265.00 56.100 2.697e-09 ***
## avg_pollen 1 60.958 232.72 9.820 0.001726 **
## ---
## 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 = 8.31087449,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8906 -1.0503 -0.3603 0.6087 2.1456
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.48898 0.50908 -0.961 0.336787
## fungicideTRUE 0.22321 0.18244 1.224 0.221140
## crithidiaTRUE 0.08251 0.18540 0.445 0.656290
## workers_alive 0.46356 0.11923 3.888 0.000101 ***
## block4 -0.53894 0.43472 -1.240 0.215077
## block6 -1.86957 0.68175 -2.742 0.006101 **
## block7 -0.52319 0.38118 -1.373 0.169885
## block8 -1.11085 0.39494 -2.813 0.004913 **
## block9 -0.63959 0.35012 -1.827 0.067733 .
## block10 0.62654 0.36247 1.729 0.083891 .
## block11 -1.41876 0.50792 -2.793 0.005218 **
## block12 0.44021 0.39404 1.117 0.263929
## avg_pollen 2.84629 0.88950 3.200 0.001375 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(8.3109) family taken to be 1)
##
## Null deviance: 298.876 on 35 degrees of freedom
## Residual deviance: 51.138 on 23 degrees of freedom
## AIC: 226.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 8.31
## Std. Err.: 4.54
##
## 2 x log-likelihood: -198.903
Anova(ll1)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 1.480 1 0.223743
## crithidia 0.188 1 0.664965
## workers_alive 15.474 1 8.363e-05 ***
## block 56.100 8 2.697e-09 ***
## avg_pollen 9.820 1 0.001726 **
## ---
## 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.31 0.956 Inf 4.49 8.85 a
## TRUE 6.85 1.057 Inf 4.85 9.67 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 14.9 9 18.7 6.23
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 <- glm.nb(dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
drop1(dl.mod, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 34.069 118.39
## fungicide 1 34.491 116.81 0.42200 0.5159
## crithidia 1 34.310 116.63 0.24113 0.6234
## avg_pollen 1 35.877 118.20 1.80740 0.1788
## workers_alive 1 34.709 117.03 0.63955 0.4239
dl1 <- update(dl.mod, .~. -workers_alive)
drop1(dl1, test = "Chisq")
## Single term deletions
##
## Model:
## dead_larvae ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 33.976 117.02
## fungicide 1 34.462 115.51 0.4862 0.48563
## crithidia 1 34.007 115.06 0.0316 0.85896
## avg_pollen 1 39.815 120.86 5.8390 0.01567 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(brood$treatment, brood$dead_larvae)
Anova(dl1)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.4862 1 0.48563
## crithidia 0.0316 1 0.85896
## avg_pollen 5.8390 1 0.01567 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dl1)
##
## Call:
## glm.nb(formula = dead_larvae ~ fungicide + crithidia + avg_pollen,
## data = brood, init.theta = 0.694531903, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5917 -1.0642 -0.8710 0.4214 1.5860
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2997 0.7873 -1.651 0.0988 .
## fungicideTRUE 0.3838 0.5182 0.741 0.4589
## crithidiaTRUE 0.0926 0.5277 0.175 0.8607
## avg_pollen 2.8820 1.1724 2.458 0.0140 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.6945) family taken to be 1)
##
## Null deviance: 39.938 on 35 degrees of freedom
## Residual deviance: 33.976 on 32 degrees of freedom
## AIC: 119.02
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.695
## Std. Err.: 0.331
##
## 2 x log-likelihood: -109.023
dp.mod <- glm.nb(dead_pupae ~ fungicide + crithidia + avg_pollen + workers_alive, 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 13.256 33.256
## fungicide 1 13.489 31.490 0.2336 0.62883
## crithidia 1 13.410 31.411 0.1541 0.69468
## avg_pollen 1 19.183 37.183 5.9270 0.01491 *
## workers_alive 1 17.549 35.550 4.2932 0.03826 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.2336 1 0.62883
## crithidia 0.1541 1 0.69468
## avg_pollen 5.9270 1 0.01491 *
## workers_alive 4.2932 1 0.03826 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dp.mod)
##
## Call:
## glm.nb(formula = dead_pupae ~ fungicide + crithidia + avg_pollen +
## workers_alive, data = brood, init.theta = 7846.438957, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1378 -0.3987 -0.2465 -0.1208 1.3019
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2686 1.6092 -2.031 0.0422 *
## fungicideTRUE 0.4393 0.9185 0.478 0.6324
## crithidiaTRUE 0.4173 1.0597 0.394 0.6938
## avg_pollen 8.6462 4.4149 1.958 0.0502 .
## workers_alive -1.0658 0.5938 -1.795 0.0726 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7846.439) family taken to be 1)
##
## Null deviance: 19.740 on 35 degrees of freedom
## Residual deviance: 13.256 on 31 degrees of freedom
## AIC: 35.256
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7846
## Std. Err.: 363749
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -23.256
plot(brood$treatment, brood$dead_pupae)
Anova(dp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.2336 1 0.62883
## crithidia 0.1541 1 0.69468
## avg_pollen 5.9270 1 0.01491 *
## workers_alive 4.2932 1 0.03826 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tl.mod <- glm.nb(total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
drop1(tl.mod, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.843 247.76
## fungicide 1 47.281 249.20 3.4379 0.06372 .
## crithidia 1 44.368 246.28 0.5250 0.46874
## avg_pollen 1 67.532 269.45 23.6883 1.133e-06 ***
## workers_alive 1 46.292 248.21 2.4490 0.11760
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tl.mod1 <- update(tl.mod, .~. -workers_alive)
drop1(tl.mod1, test = "Chisq")
## Single term deletions
##
## Model:
## total_larvae ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.813 248.13
## fungicide 1 46.281 248.60 2.468 0.1162
## crithidia 1 43.817 246.14 0.004 0.9516
## avg_pollen 1 92.798 295.12 48.985 2.579e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(tl.mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_larvae
## LR Chisq Df Pr(>Chisq)
## fungicide 2.468 1 0.1162
## crithidia 0.004 1 0.9516
## avg_pollen 48.985 1 2.579e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(tl.mod1)
##
## Call:
## glm.nb(formula = total_larvae ~ fungicide + crithidia + avg_pollen,
## data = brood, init.theta = 1.388643983, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8612 -1.3451 -0.4270 0.4856 1.4306
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.52648 0.48115 -1.094 0.2739
## fungicideTRUE 0.52458 0.31805 1.649 0.0991 .
## crithidiaTRUE 0.01889 0.32256 0.059 0.9533
## avg_pollen 6.06181 0.73581 8.238 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.3886) family taken to be 1)
##
## Null deviance: 94.065 on 35 degrees of freedom
## Residual deviance: 43.813 on 32 degrees of freedom
## AIC: 250.13
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.389
## Std. Err.: 0.448
##
## 2 x log-likelihood: -240.133
plot(brood$treatment, brood$total_larvae)
tp.mod <- glm.nb(total_pupae ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
## 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
## Df Deviance AIC LRT Pr(>Chi)
## <none> 48.107 155.06
## fungicide 1 48.150 153.11 0.043 0.836032
## crithidia 1 48.503 153.46 0.396 0.528996
## avg_pollen 1 92.212 197.17 44.105 3.112e-11 ***
## workers_alive 1 58.022 162.98 9.915 0.001639 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(tp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_pupae
## LR Chisq Df Pr(>Chisq)
## fungicide 0.043 1 0.836032
## crithidia 0.396 1 0.528996
## avg_pollen 44.105 1 3.112e-11 ***
## workers_alive 9.915 1 0.001639 **
## ---
## 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, data = brood, init.theta = 22.77068534, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1745 -1.1437 -0.2733 0.1067 3.1363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44348 0.45544 -3.169 0.00153 **
## fungicideTRUE -0.03668 0.17559 -0.209 0.83454
## crithidiaTRUE -0.11590 0.18143 -0.639 0.52294
## avg_pollen 3.31346 0.51336 6.454 1.09e-10 ***
## workers_alive 0.34404 0.11402 3.017 0.00255 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(22.7707) family taken to be 1)
##
## Null deviance: 206.184 on 35 degrees of freedom
## Residual deviance: 48.107 on 31 degrees of freedom
## AIC: 157.06
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 22.8
## Std. Err.: 27.9
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -145.064
plot(brood$treatment, brood$total_pupae)
egg.mod <- glm.nb(eggs ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
drop1(egg.mod, test = "Chisq")
## Single term deletions
##
## Model:
## eggs ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.916 281.18
## fungicide 1 45.286 280.55 1.3697 0.24186
## crithidia 1 43.935 279.20 0.0191 0.89005
## avg_pollen 1 46.478 281.74 2.5623 0.10944
## workers_alive 1 48.379 283.64 4.4628 0.03464 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em1 <- update(egg.mod, .~. -avg_pollen)
drop1(em1, test = "Chisq")
## Single term deletions
##
## Model:
## eggs ~ fungicide + crithidia + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 43.367 281.62
## fungicide 1 43.939 280.19 0.5719 0.4495
## crithidia 1 43.482 279.73 0.1144 0.7352
## workers_alive 1 58.565 294.82 15.1974 9.684e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(em1)
## Analysis of Deviance Table (Type II tests)
##
## Response: eggs
## LR Chisq Df Pr(>Chisq)
## fungicide 0.5719 1 0.4495
## crithidia 0.1144 1 0.7352
## workers_alive 15.1974 1 9.684e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(em1)
##
## Call:
## glm.nb(formula = eggs ~ fungicide + crithidia + workers_alive,
## data = brood, init.theta = 0.9210361052, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5071 -1.0676 -0.1225 0.3527 1.6103
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7006 0.6135 1.142 0.253
## fungicideTRUE 0.2787 0.3643 0.765 0.444
## crithidiaTRUE -0.1241 0.3840 -0.323 0.747
## workers_alive 0.5793 0.1265 4.579 4.67e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.921) family taken to be 1)
##
## Null deviance: 59.183 on 35 degrees of freedom
## Residual deviance: 43.367 on 32 degrees of freedom
## AIC: 283.62
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.921
## Std. Err.: 0.251
##
## 2 x log-likelihood: -273.620
plot(brood$treatment, brood$eggs)
hp.mod <- glm.nb(honey_pots ~ fungicide + crithidia + avg_pollen + workers_alive, data = brood)
drop1(hp.mod, test = "Chisq")
## Single term deletions
##
## Model:
## honey_pots ~ fungicide + crithidia + avg_pollen + workers_alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 35.669 141.06
## fungicide 1 39.320 142.71 3.6505 0.056053 .
## crithidia 1 36.596 139.99 0.9269 0.335659
## avg_pollen 1 43.770 147.16 8.1004 0.004425 **
## workers_alive 1 37.407 140.80 1.7372 0.187489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hp1 <- update(hp.mod, .~. -workers_alive)
drop1(hp1, test = "Chisq")
## Single term deletions
##
## Model:
## honey_pots ~ fungicide + crithidia + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 36.875 140.79
## fungicide 1 40.087 142.00 3.2115 0.07312 .
## crithidia 1 38.392 140.31 1.5173 0.21803
## avg_pollen 1 59.717 161.63 22.8423 1.758e-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.57633 32
## 2 fungicide + crithidia + avg_pollen + workers_alive 19.67711 31
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -132.7911
## 2 -131.0589 1 vs 2 1 1.732166 0.1881346
Anova(hp.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pots
## LR Chisq Df Pr(>Chisq)
## fungicide 3.6505 1 0.056053 .
## crithidia 0.9269 1 0.335659
## avg_pollen 8.1004 1 0.004425 **
## workers_alive 1.7372 1 0.187489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(hp.mod, hp1)
## df AIC
## hp.mod 6 143.0589
## hp1 5 142.7911
Anova(hp1)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pots
## LR Chisq Df Pr(>Chisq)
## fungicide 3.2115 1 0.07312 .
## crithidia 1.5173 1 0.21803
## avg_pollen 22.8423 1 1.758e-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.5763296, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0615 -0.7252 -0.3248 0.4896 2.2243
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.08803 0.33089 -0.266 0.7902
## fungicideTRUE 0.38631 0.21588 1.789 0.0735 .
## crithidiaTRUE -0.27473 0.22410 -1.226 0.2202
## avg_pollen 2.28998 0.46852 4.888 1.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(17.5763) family taken to be 1)
##
## Null deviance: 66.981 on 35 degrees of freedom
## Residual deviance: 36.875 on 32 degrees of freedom
## AIC: 142.79
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 17.6
## Std. Err.: 25.2
##
## 2 x log-likelihood: -132.791
plot(brood$treatment, brood$honey_pots)
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) #not overdisp -> Start with this one
##
## 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.9548 -0.8831 -0.1200 0.2351 2.0003
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.569e+00 1.411e+00 1.821 0.0686 .
## fungicideTRUE 3.633e-02 1.374e-01 0.264 0.7915
## crithidiaTRUE 2.629e-02 1.351e-01 0.195 0.8457
## workers_alive 1.418e-01 1.005e-01 1.411 0.1583
## block4 -2.358e-01 3.810e-01 -0.619 0.5360
## block6 -1.760e+01 1.728e+03 -0.010 0.9919
## block7 5.448e-02 3.538e-01 0.154 0.8776
## block8 -3.568e-01 3.738e-01 -0.955 0.3398
## block9 2.143e-01 3.289e-01 0.652 0.5146
## block10 8.056e-01 3.505e-01 2.298 0.0215 *
## block11 2.199e-01 3.944e-01 0.558 0.5771
## block12 -9.284e-04 4.101e-01 -0.002 0.9982
## duration -6.552e-02 2.833e-02 -2.312 0.0208 *
## avg_pollen 3.210e+00 6.942e-01 4.624 3.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 308.730 on 35 degrees of freedom
## Residual deviance: 42.559 on 22 degrees of freedom
## AIC: 171.75
##
## Number of Fisher Scoring iterations: 15
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 22 42.559
## 2 22 42.554 0 0.0050672
AIC(dronecount.mod, dronecount.mod.nb)
## df AIC
## dronecount.mod 14 171.7458
## dronecount.mod.nb 15 173.7469
dronecount.mod.int <- glm(total_drones ~ fungicide*crithidia + workers_alive + block + duration + avg_pollen, data = brood, family = "poisson")
summary(dronecount.mod.int)
##
## 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.9851 -0.9658 -0.1502 0.4410 1.8628
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.89260 1.45933 1.982 0.0475 *
## fungicideTRUE -0.05110 0.17463 -0.293 0.7698
## crithidiaTRUE -0.06693 0.17700 -0.378 0.7053
## workers_alive 0.13297 0.10201 1.304 0.1924
## block4 -0.29041 0.39015 -0.744 0.4567
## block6 -18.60951 2848.99391 -0.007 0.9948
## block7 0.03311 0.35621 0.093 0.9259
## block8 -0.38187 0.37383 -1.022 0.3070
## block9 0.17892 0.33247 0.538 0.5905
## block10 0.79568 0.34810 2.286 0.0223 *
## block11 0.17341 0.39785 0.436 0.6629
## block12 -0.07470 0.42517 -0.176 0.8605
## duration -0.07070 0.02889 -2.447 0.0144 *
## avg_pollen 3.22393 0.69718 4.624 3.76e-06 ***
## fungicideTRUE:crithidiaTRUE 0.21745 0.26553 0.819 0.4128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 308.730 on 35 degrees of freedom
## Residual deviance: 41.887 on 21 degrees of freedom
## AIC: 173.07
##
## Number of Fisher Scoring iterations: 16
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.887 173.07
## workers_alive 1 43.607 172.79 1.720 0.18969
## block 8 81.181 196.37 39.295 4.334e-06 ***
## duration 1 47.481 176.67 5.594 0.01802 *
## avg_pollen 1 64.714 193.90 22.827 1.772e-06 ***
## fungicide:crithidia 1 42.559 171.75 0.672 0.41236
## ---
## 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.070 1 0.79135
## crithidia 0.038 1 0.84573
## workers_alive 1.720 1 0.18969
## block 39.295 8 4.334e-06 ***
## duration 5.594 1 0.01802 *
## avg_pollen 22.827 1 1.772e-06 ***
## fungicide:crithidia 0.672 1 0.41236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(brood$treatment, brood$total_drones)
drop1(dronecount.mod, test = "Chisq")
## Single term deletions
##
## Model:
## total_drones ~ fungicide + crithidia + workers_alive + block +
## duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 42.559 171.75
## fungicide 1 42.629 169.82 0.070 0.79135
## crithidia 1 42.597 169.78 0.038 0.84573
## workers_alive 1 44.584 171.77 2.025 0.15471
## block 8 81.316 194.50 38.757 5.453e-06 ***
## duration 1 47.573 174.76 5.014 0.02515 *
## avg_pollen 1 65.399 192.59 22.840 1.760e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dc1 <- update(dronecount.mod, .~. -workers_alive)
drop1(dc1, test = "Chisq")
## Single term deletions
##
## Model:
## total_drones ~ fungicide + crithidia + block + duration + avg_pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 44.584 171.77
## fungicide 1 44.586 169.77 0.002 0.968341
## crithidia 1 44.587 169.77 0.003 0.953344
## block 8 85.588 196.78 41.004 2.081e-06 ***
## duration 1 51.406 176.59 6.822 0.009002 **
## avg_pollen 1 94.536 219.72 49.952 1.576e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(dc1)
##
## Call:
## glm(formula = total_drones ~ fungicide + crithidia + block +
## duration + avg_pollen, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.91124 -0.96690 -0.04182 0.30982 2.00093
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.254e+00 1.288e+00 2.526 0.01154 *
## fungicideTRUE -5.302e-03 1.336e-01 -0.040 0.96834
## crithidiaTRUE 7.863e-03 1.344e-01 0.059 0.95334
## block4 -3.821e-01 3.708e-01 -1.031 0.30277
## block6 -1.767e+01 1.733e+03 -0.010 0.99187
## block7 -6.462e-02 3.471e-01 -0.186 0.85231
## block8 -4.637e-01 3.712e-01 -1.249 0.21158
## block9 2.513e-01 3.280e-01 0.766 0.44356
## block10 6.737e-01 3.344e-01 2.015 0.04394 *
## block11 1.909e-01 3.937e-01 0.485 0.62772
## block12 -2.884e-01 3.654e-01 -0.789 0.42993
## duration -7.220e-02 2.697e-02 -2.677 0.00742 **
## avg_pollen 3.800e+00 5.659e-01 6.714 1.89e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 308.730 on 35 degrees of freedom
## Residual deviance: 44.584 on 23 degrees of freedom
## AIC: 171.77
##
## Number of Fisher Scoring iterations: 15
Anova(dc1)
## Analysis of Deviance Table (Type II tests)
##
## Response: total_drones
## LR Chisq Df Pr(>Chisq)
## fungicide 0.002 1 0.968341
## crithidia 0.003 1 0.953344
## block 41.004 8 2.081e-06 ***
## duration 6.822 1 0.009002 **
## avg_pollen 49.952 1 1.576e-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 6.89 9 7.18 2.39
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, 0.4), # 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 Live Pupae") +
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", "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")
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.2335 1 0.628962
## crithidia 0.0073 1 0.931992
## avg_pollen 1.3257 1 0.249572
## block 26.0519 8 0.001029 **
## duration 5.9041 1 0.015106 *
## workers_alive 1.1226 1 0.289357
## ---
## 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> 46.944 114.62
## fungicide 1 47.178 112.85 0.2335 0.628962
## crithidia 1 46.952 112.62 0.0073 0.931992
## avg_pollen 1 48.270 113.94 1.3257 0.249572
## block 8 72.996 124.67 26.0519 0.001029 **
## duration 1 52.848 118.52 5.9041 0.015106 *
## workers_alive 1 48.067 113.74 1.1226 0.289357
## ---
## 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> 48.067 113.74
## fungicide 1 48.310 111.98 0.2435 0.621688
## crithidia 1 48.119 111.79 0.0525 0.818706
## avg_pollen 1 48.596 112.27 0.5290 0.467034
## block 8 73.013 122.69 24.9462 0.001588 **
## duration 1 58.565 122.24 10.4982 0.001195 **
## ---
## 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> 48.596 112.27
## fungicide 1 48.829 110.50 0.2330 0.6292730
## crithidia 1 48.616 110.29 0.0199 0.8877947
## block 8 77.731 125.40 29.1347 0.0003003 ***
## duration 1 59.635 121.31 11.0388 0.0008922 ***
## ---
## 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.2330 1 0.6292730
## crithidia 0.0199 1 0.8877947
## block 29.1347 8 0.0003003 ***
## duration 11.0388 1 0.0008922 ***
## ---
## 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.3237 -0.5519 0.0000 0.7894 2.5069
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.60178 3.97877 3.670 0.000243 ***
## fungicideTRUE 0.17912 0.37191 0.482 0.630077
## crithidiaTRUE -0.04958 0.35142 -0.141 0.887792
## block4 -1.42228 0.88102 -1.614 0.106452
## block6 -1.56753 1.01330 -1.547 0.121875
## block7 -1.97885 0.86417 -2.290 0.022028 *
## block8 -1.37953 0.92912 -1.485 0.137603
## block9 -2.07110 0.86256 -2.401 0.016345 *
## block10 1.04134 0.65380 1.593 0.111216
## block11 -1.44815 1.30648 -1.108 0.267675
## block12 -0.08539 0.76238 -0.112 0.910824
## duration -0.26367 0.08352 -3.157 0.001595 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 79.536 on 27 degrees of freedom
## Residual deviance: 48.596 on 16 degrees of freedom
## AIC: 112.27
##
## Number of Fisher Scoring iterations: 5
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: 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.47 1 0.4948
## crithidia 0.00 1 1.0000
## avg_pollen 0.00 1 1.0000
## block 7.89 8 0.4447
## duration 334.66 1 <2e-16 ***
## workers_alive 5.85 1 0.0156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(ppmod, test = "Chisq")
## 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.00 36.63
## fungicide 1 0.47 35.10 0.47 0.4948
## crithidia 1 0.00 34.63 0.00 1.0000
## avg_pollen 1 0.00 34.63 0.00 1.0000
## block 8 7.89 28.52 7.89 0.4447
## duration 1 334.66 369.29 334.66 <2e-16 ***
## workers_alive 1 5.85 40.48 5.85 0.0156 *
## ---
## 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.8857 28.518
## fungicide 1 8.5302 27.163 0.6445 0.42207
## crithidia 1 9.3781 28.011 1.4925 0.22184
## avg_pollen 1 11.8154 30.448 3.9298 0.04744 *
## duration 1 10.9993 29.632 3.1137 0.07764 .
## workers_alive 1 14.4788 33.112 6.5931 0.01024 *
## ---
## 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> 10.999 29.632
## fungicide 1 11.021 27.653 0.0211 0.884443
## crithidia 1 11.486 28.119 0.4864 0.485529
## avg_pollen 1 12.481 29.114 1.4817 0.223507
## workers_alive 1 18.067 34.699 7.0672 0.007851 **
## ---
## 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> 12.481 29.114
## fungicide 1 12.559 27.192 0.0779 0.78011
## crithidia 1 12.525 27.158 0.0444 0.83313
## workers_alive 1 18.536 33.169 6.0550 0.01387 *
## ---
## 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.0779 1 0.78011
## crithidia 0.0444 1 0.83313
## workers_alive 6.0550 1 0.01387 *
## ---
## 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.7274 0.0000 0.2950 0.4739 0.9934
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.06189 2.02471 0.031 0.976
## fungicideTRUE 0.32433 1.17111 0.277 0.782
## crithidiaTRUE -0.23659 1.11658 -0.212 0.832
## workers_alive 0.86841 0.36701 2.366 0.018 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20.674 on 24 degrees of freedom
## Residual deviance: 12.481 on 21 degrees of freedom
## AIC: 29.114
##
## 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> 49.967 117.75
## fungicide 1 50.195 115.98 0.2287 0.6324816
## crithidia 1 50.106 115.89 0.1398 0.7085134
## block 8 71.401 123.19 21.4348 0.0060778 **
## duration 1 60.986 126.78 11.0197 0.0009015 ***
## ---
## 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.7214 -0.2108 0.0000 0.8444 2.4278
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.4401 3.7762 3.824 0.000131 ***
## fungicideTRUE 0.1698 0.3558 0.477 0.633193
## crithidiaTRUE -0.1214 0.3245 -0.374 0.708384
## block4 -1.5238 0.8558 -1.781 0.074988 .
## block6 -1.7535 0.9540 -1.838 0.066072 .
## block7 -1.9081 0.8343 -2.287 0.022188 *
## block8 -1.2544 0.9113 -1.376 0.168680
## block9 -1.9538 0.8137 -2.401 0.016340 *
## block10 0.4782 0.6109 0.783 0.433707
## block11 -1.3380 1.2701 -1.053 0.292157
## block12 -0.5217 0.7178 -0.727 0.467398
## duration -0.2509 0.0791 -3.172 0.001516 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 76.309 on 27 degrees of freedom
## Residual deviance: 49.967 on 16 degrees of freedom
## AIC: 117.76
##
## 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.2287 1 0.6324816
## crithidia 0.1398 1 0.7085134
## block 21.4348 8 0.0060778 **
## duration 11.0197 1 0.0009015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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
## fungicide -0.0006294 0.0014127 -0.446
## crithidia 0.0011136 0.0014096 0.790
## emerge -0.0003927 0.0001602 -2.452
##
## Correlation of Fixed Effects:
## (Intr) fungcd crithd
## fungicide 0.073
## crithidia -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
shapiro.test(drones$radial_cell)
##
## Shapiro-Wilk normality test
##
## data: drones$radial_cell
## W = 0.97715, p-value = 0.0001916
hist(drones$radial_cell)
drones$log_rad <- log(drones$radial_cell)
shapiro.test(drones$log_rad)
##
## Shapiro-Wilk normality test
##
## data: drones$log_rad
## W = 0.96056, p-value = 6.877e-07
hist(drones$log_rad)
drones$square <- drones$radial_cell^3
shapiro.test(drones$square)
##
## Shapiro-Wilk normality test
##
## data: drones$square
## W = 0.99259, p-value = 0.1797
drones$box_rad <- bcPower(drones$radial_cell, lambda = 3, gamma = 1)
shapiro.test(drones$box_rad)
##
## Shapiro-Wilk normality test
##
## data: drones$box_rad
## W = 0.99259, p-value = 0.1795
hist(drones$box_rad)
plot(drones$treatment, drones$radial_cell)
plot(drones_rf$treatment, drones_rf$radial_cell)
drones_rf$square <-(drones_rf$radial_cell)^3
shapiro.test(drones_rf$square)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$square
## W = 0.99215, p-value = 0.1516
drones_rf$mm <- (drones_rf$radial_cell)/1000
range(drones_rf$mm)
## [1] 2.073526 3.083439
shapiro.test(drones_rf$mm)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$mm
## W = 0.97655, p-value = 0.0001662
rad_mod <- lmer(square ~ fungicide + crithidia + workers_alive + block + mean.pollen + emerge + (1|colony), data = drones_rf)
drop1(rad_mod, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + block + mean.pollen +
## emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12933
## fungicide 1 12938 7.0625 0.007872 **
## crithidia 1 12933 1.6540 0.198421
## workers_alive 1 12933 1.9253 0.165269
## block 7 12930 11.1025 0.134210
## mean.pollen 1 12934 3.1333 0.076707 .
## emerge 1 12936 4.2611 0.038996 *
## ---
## 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 + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 12930
## fungicide 1 12931 3.0785 0.07933 .
## crithidia 1 12929 0.8953 0.34404
## workers_alive 1 12929 0.2898 0.59035
## mean.pollen 1 12929 0.2152 0.64270
## emerge 1 12935 6.6204 0.01008 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm2 <- update(rm1, .~. -mean.pollen)
drop1(rm2, test = "Chisq")
## Single term deletions
##
## Model:
## square ~ fungicide + crithidia + workers_alive + emerge + (1 |
## colony)
## npar AIC LRT Pr(Chi)
## <none> 12929
## fungicide 1 12930 2.9609 0.085300 .
## crithidia 1 12927 0.8763 0.349209
## workers_alive 1 12927 0.7077 0.400224
## emerge 1 12934 7.0805 0.007793 **
## ---
## 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 + 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(rm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: square ~ fungicide + crithidia + emerge + (1 | colony)
## Data: drones_rf
##
## 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
## fungicide -1.143e+09 6.022e+08 -1.898
## crithidia 4.190e+08 5.965e+08 0.703
## emerge -2.274e+08 8.179e+07 -2.781
##
## Correlation of Fixed Effects:
## (Intr) fungcd crithd
## fungicide 0.123
## crithidia -0.055 -0.009
## emerge -0.989 -0.206 -0.022
Anova(rm3)
## 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
qqnorm(resid(rad_mod));qqline(resid(rad_mod))
qqnorm(resid(rm1));qqline(resid(rm1))
qqnorm(resid(rm2));qqline(resid(rm2))
qqnorm(resid(rm3));qqline(resid(rm3))
Anova(rad_mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: square
## Chisq Df Pr(>Chisq)
## fungicide 4.9667 1 0.02584 *
## crithidia 0.3539 1 0.55189
## workers_alive 0.7190 1 0.39646
## block 6.3989 7 0.49402
## mean.pollen 1.0793 1 0.29884
## emerge 4.6697 1 0.03070 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rem <- emmeans(rm3, pairwise ~ fungicide, type = "response")
rem
## $emmeans
## fungicide emmean SE df lower.CL upper.CL
## 0 2.03e+10 3.96e+08 17.1 1.94e+10 2.11e+10
## 1 1.91e+10 4.62e+08 20.5 1.82e+10 2.01e+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.14e+09 6.08e+08 19.4 1.879 0.0754
##
## 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.91e+10 4.62e+08 20.5 1.80e+10 2.03e+10 a
## 0 2.03e+10 3.96e+08 17.1 1.93e+10 2.13e+10 a
##
## 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
sum_radial.mm <- drones_rf %>%
group_by(treatment) %>%
summarise(m = mean(mm),
sd = sd(mm),
n = length(mm)) %>%
mutate(se = sd/sqrt(n))
sum_radial.mm
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 2.70 0.182 99 0.0183
## 2 2 2.66 0.174 68 0.0211
## 3 3 2.66 0.173 49 0.0247
## 4 4 2.75 0.136 60 0.0175
sum_radial.mm$plot <- sum_radial.mm$m + sum_radial.mm$se
ggplot(sum_radial.mm, 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 (mm)") +
theme_classic(base_size = 20) +
coord_cartesian(ylim=c(0,3)) +
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))
shapiro.test(drones_rf$relative_fat)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$relative_fat
## W = 0.97273, p-value = 4.049e-05
hist(drones_rf$relative_fat)
plot(drones_rf$treatment, drones_rf$relative_fat)
plot(drones_rf$treatment, drones_rf$fat_content)
range(drones_rf$relative_fat)
## [1] 0.0178 0.3400
drones_rf$log_ref <- log(drones_rf$relative_fat)
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$box_rf <- bcPower(drones_rf$relative_fat, lambda = -5, gamma = 3)
shapiro.test(drones_rf$box_rf)
##
## Shapiro-Wilk normality test
##
## data: drones_rf$box_rf
## W = 0.98983, p-value = 0.05093
hist(drones_rf$box_rf)
rf_mod <- lmer(box_rf ~ fungicide + crithidia + block + mean.pollen + workers_alive + emerge + (1|colony), data = drones_rf)
drop1(rf_mod, test = "Chisq")
## Single term deletions
##
## Model:
## box_rf ~ fungicide + crithidia + block + mean.pollen + workers_alive +
## emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4745.1
## fungicide 1 -4745.4 1.6685 0.1964544
## crithidia 1 -4744.6 2.4165 0.1200666
## block 7 -4739.0 20.0043 0.0055605 **
## mean.pollen 1 -4745.9 1.1285 0.2880969
## workers_alive 1 -4747.0 0.0185 0.8918694
## emerge 1 -4733.6 13.4465 0.0002455 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf1 <- update(rf_mod, .~. -workers_alive)
drop1(rf1, test = "Chisq")
## Single term deletions
##
## Model:
## box_rf ~ fungicide + crithidia + block + mean.pollen + emerge +
## (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4747.0
## fungicide 1 -4747.3 1.7028 0.1919275
## crithidia 1 -4746.5 2.5201 0.1124013
## block 7 -4740.5 20.5477 0.0045006 **
## mean.pollen 1 -4747.5 1.5334 0.2156066
## emerge 1 -4735.5 13.5735 0.0002294 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf2 <- update(rf1, .~. -mean.pollen)
drop1(rf2, test = "Chisq")
## Single term deletions
##
## Model:
## box_rf ~ fungicide + crithidia + block + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4747.5
## fungicide 1 -4748.5 0.9992 0.3175155
## crithidia 1 -4746.4 3.0840 0.0790674 .
## block 7 -4741.3 20.1800 0.0051937 **
## emerge 1 -4736.4 13.0788 0.0002987 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf3 <- update(rf2, .~. -block)
drop1(rf3, test = "Chisq")
## Single term deletions
##
## Model:
## box_rf ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4741.3
## fungicide 1 -4743.3 0.0000 0.99560
## crithidia 1 -4738.8 4.5172 0.03356 *
## emerge 1 -4724.5 18.7897 1.46e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf4 <- update(rf3, .~. -fungicide)
drop1(rf4, test = "Chisq")
## Single term deletions
##
## Model:
## box_rf ~ crithidia + emerge + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -4743.3
## crithidia 1 -4740.8 4.5197 0.03351 *
## emerge 1 -4726.0 19.3533 1.086e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rf4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: box_rf
## Chisq Df Pr(>Chisq)
## crithidia 4.4038 1 0.03586 *
## emerge 19.9613 1 7.903e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rf3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: box_rf
## Chisq Df Pr(>Chisq)
## fungicide 0.0028 1 0.95772
## crithidia 4.2484 1 0.03929 *
## emerge 19.3367 1 1.096e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(rf3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: box_rf ~ fungicide + crithidia + emerge + (1 | colony)
## Data: drones_rf
##
## REML criterion at convergence: -4662.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4419 -0.5485 -0.0038 0.5572 4.3343
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 4.178e-10 2.044e-05
## Residual 1.771e-09 4.209e-05
## Number of obs: 276, groups: colony, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.995e-01 4.024e-05 4957.603
## fungicide -5.461e-07 1.030e-05 -0.053
## crithidia 2.126e-05 1.031e-05 2.061
## emerge -4.716e-06 1.072e-06 -4.397
##
## Correlation of Fixed Effects:
## (Intr) fungcd crithd
## fungicide 0.046
## crithidia -0.069 0.002
## emerge -0.980 -0.162 -0.035
qqnorm(resid(rf4));qqline(resid(rf4))
qqnorm(resid(rf3));qqline(resid(rf3))
anova(rf3, rf4, test = "Chisq")
## Data: drones_rf
## Models:
## rf4: box_rf ~ crithidia + emerge + (1 | colony)
## rf3: box_rf ~ fungicide + crithidia + emerge + (1 | colony)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rf4 5 -4743.3 -4725.2 2376.7 -4753.3
## rf3 6 -4741.3 -4719.6 2376.7 -4753.3 0 1 0.9956
AIC(rf3, rf4)
## df AIC
## rf3 6 -4650.563
## rf4 5 -4673.716
rf_em <- emmeans(rf3, pairwise ~ crithidia, type = "response")
rf_e <- setDT(as.data.frame(rf_em$emmeans))
rf_ce <- setDT(as.data.frame(rf_em$contrasts))
rf_em
## $emmeans
## crithidia emmean SE df lower.CL upper.CL
## 0 0.19931 6.6023e-06 19.34 0.19930 0.19932
## 1 0.19933 8.0124e-06 20.59 0.19931 0.19935
##
## Results are averaged over the levels of: fungicide
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## crithidia0 - crithidia1 -2.13e-05 1.04e-05 19.8 -2.054 0.0534
##
## Results are averaged over the levels of: fungicide
## Degrees-of-freedom method: kenward-roger
rf_e
## crithidia emmean SE df lower.CL upper.CL
## 1: 0 0.1993103 6.602275e-06 19.34415 0.1992965 0.1993241
## 2: 1 0.1993316 8.012361e-06 20.59147 0.1993149 0.1993483
rf_ce
## contrast estimate SE df t.ratio
## 1: crithidia0 - crithidia1 -2.125763e-05 1.035116e-05 19.81501 -2.053647
## p.value
## 1: 0.05344528
drones_rf$rfgmm <- (drones_rf$relative_fat_original)*1000
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 <- drones_rf %>%
group_by(treatment) %>%
summarise(m = mean(rfgmm),
sd = sd(rfgmm),
n = length(rfgmm)) %>%
mutate(se = sd/sqrt(n))
rf_sum
## # A tibble: 4 × 5
## treatment m sd n se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.00115 0.000438 99 0.0000440
## 2 2 0.00111 0.000408 68 0.0000495
## 3 3 0.00129 0.000508 49 0.0000726
## 4 4 0.00128 0.000435 60 0.0000561
rf_sum$plot <- rf_sum$m + rf_sum$se
ggplot(rf_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, 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 = m - se, ymax = m + se), width = 0.2, position = position_dodge(0.9)) +
labs(x = "Treatment", y = "Relative Fat (g/mm)") +
theme_classic(base_size = 20) +
coord_cartesian(ylim=c(0, 0.0015)) +
annotate(geom = "text",
x = 1, y = 0.05,
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", "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")
drones$fungicide <- as.logical(drones$fungicide)
drones$crithidia <- as.logical(drones$crithidia)
em.mod <- glm.nb(emerge ~ fungicide + crithidia + dry_weight + live_weight + workers_alive + mean.pollen, data = drones)
## 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(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> 50.254 1598.9
## fungicide 1 53.127 1599.8 2.87292 0.09008 .
## crithidia 1 50.358 1597.0 0.10395 0.74714
## dry_weight 1 51.685 1598.3 1.43124 0.23156
## live_weight 1 50.753 1597.4 0.49862 0.48011
## workers_alive 1 50.577 1597.2 0.32243 0.57015
## mean.pollen 1 51.461 1598.1 1.20639 0.27205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em.mod <- glm.nb(emerge ~ fungicide + crithidia + mean.pollen, data = drones)
## 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 + mean.pollen,
## data = drones, init.theta = 3582640.712, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.03374 -0.37053 -0.04538 0.30199 1.82050
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.668832 0.037831 96.980 <2e-16 ***
## fungicideTRUE 0.038467 0.019723 1.950 0.0511 .
## crithidiaTRUE -0.007431 0.020095 -0.370 0.7116
## mean.pollen -0.097250 0.054745 -1.776 0.0757 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3582641) family taken to be 1)
##
## Null deviance: 59.662 on 280 degrees of freedom
## Residual deviance: 52.683 on 277 degrees of freedom
## AIC: 1597.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3582641
## Std. Err.: 40058159
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -1587.309
drop1(em.mod, test = "Chisq")
## Single term deletions
##
## Model:
## emerge ~ fungicide + crithidia + mean.pollen
## Df Deviance AIC LRT Pr(>Chi)
## <none> 52.683 1595.3
## fungicide 1 56.480 1597.1 3.7967 0.05135 .
## crithidia 1 52.820 1593.5 0.1368 0.71148
## mean.pollen 1 55.833 1596.5 3.1506 0.07590 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(em.mod)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## fungicide 3.7967 1 0.05135 .
## crithidia 0.1368 1 0.71148
## mean.pollen 3.1506 1 0.07590 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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, 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 = 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(30, 40)) +
annotate(geom = "text",
x = 1, y = 0.05,
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", "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")
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