Bees are frequently exposed to fungicides in agricultural landscapes, and while these chemicals are generally not considered to be harmful to insect pollinators, the sublethal effects of fungicides are not well understood. We investigated the non-target effects of exposure to field-realistic concentrations of the fungicide, Pristine®, for the common eastern bumble bee (Bombus impatiens). Over two years, We fed 60 queenless microcolonies of bumble bees pollen patties contaminated with Pristine® ranging from 0 ppb-150,000 ppb, and evaluated the effects of chronic fungicide exposure on brood production, colony weight change, pollen consumption, and drone body quality. We found that while microcolonies fed pollen containing 1,500 ppb consumed the most pollen and produced the most brood, microcolonies exposed to 15,000 ppb gained weight at significantly lower rates. Fungicide exposure also negatively impacted brood development time and drone body quality. Our findings indicate that Pristine® has sublethal effects on bumble bee microcolony growth and brood development, highlighting the potential negative impacts of chronic fungicide exposure to bee health. In the face of the global decline of pollinator species, it is important for us to broaden our understanding of how bees react towards exposure of often overlooked stressors, such as fungicides, to account for potential consequences for pollinator health and pollination services.
This data set includes information from experiments conducted in the summer of 2021 and summer of 2022. I fed 60 bumble bee microcolonies (15 in 2021, 45 in 2022) pollen contaminated with 1 of 5 doses of a common fungicide. Bumble bee microcolonies are comprised of 5 female worker bees and no queen. One worker will establish dominance, act as a pseudo-queen, and lay haploid drone eggs. In this way microcolonies act as a proxy system, allowing us to look at the effects of contaminants on bee health and fecundity, while using a large enough replication size not achievable with full queen-right colonies. I measured worker mortality and ending dry weights, time until first drone emergence, drone size (radial cell size in wing, and dry weight), colony weight change over time, and brood production and brood stage mortality. Honey pots are also counted within brood production, and these are the wax containers bees store their nectar in.
Bumble bee brood stages develop in this order: egg > larvae > pupae > drone. Uncapped brood cells are often turned into honey pots following drone emergence.
Brood cells, in the context of this analysis, are the total count of every reproductive unit seen in the dissected colonies. This includes all honey pots, pupal cells, larval cells, and clumps of eggs. The total brood cell count will not equal the totals of larvae + pupae + eggs + honey pots + drones for several reasons. First, an “egg clump” may contain anywhere between 3-6 eggs, but will count as only one “cell”, and a brood cell may also contain multiple larvae at once. Drones are only counted as such, rather than as pupae, if they are inside uncapped cells or already emerged inside the colony. Cells that are still capped are counted as pupae and will be equivalent to a single brood cell.
Let’s begin by looking at total brood cells.
brood <- read_csv("brood.csv", col_types = cols(round = col_factor(levels = c("1", "2")),
treatment = col_factor(levels = c("1","2", "3", "4", "5")),
replicate = col_factor(levels = c("1", "2", "3", "4", "5", "7", "9", "11", "12"))))
brood$colony <- as.factor(brood$colony)
brood$pupae <- brood$live_pupae + brood$dead_pupae
brood$larvae <- brood$dead_larvae + brood$live_larvae
head(brood)
## # A tibble: 6 × 17
## round dose treatment replicate colony brood_c…¹ honey…² eggs dead_…³ live_…⁴
## <fct> <dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 0 1 11 1.11R2 18 1 14 1 0
## 2 2 0 1 12 1.12R2 3 3 0 0 0
## 3 1 0 1 1 1.1R1 63 5 3 8 40
## 4 2 0 1 1 1.1R2 54 4 5 1 49
## 5 1 0 1 2 1.2R1 18 11 5 0 9
## 6 2 0 1 2 1.2R2 29 10 0 6 21
## # … with 7 more variables: dead_pupae <dbl>, live_pupae <dbl>,
## # dead_drones <dbl>, live_drones <dbl>, drones <dbl>, pupae <dbl>,
## # larvae <dbl>, and abbreviated variable names ¹brood_cells, ²honey_pot,
## # ³dead_larvae, ⁴live_larvae
Here I will look at shape of data. The histogram shows our count data has a slight right tail, but not too severe so we will probably use the negative binomial distribution rather than Poisson.
ggplot(brood, aes(x=brood_cells, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Brood Cells Counted")
Now I will test different distributions to handle this count data. AIC value for negative binomial distribution is much smaller, I will stick with that.
fit1.brood<-fitdistr(brood$brood_cells, "Poisson")
fit2.brood<-fitdistr(brood$brood_cells, "negative binomial")
AIC(fit1.brood, fit2.brood)
## df AIC
## fit1.brood 1 1152.0504
## fit2.brood 2 544.0061
In the model below we are testing the total number of brood cells as a function of dose of Pristine, with 0 PPB acting as the intercept (which is good, because it is our control), and holding round as a random variable with individual colony nested within. We are using glmer.nb because it is a mixed effects model, with a negative binomial distribution because it is count data.
nb.brood <- glmer.nb(brood_cells~ treatment + (1|round/colony), data = brood)
summary(nb.brood)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(5127.122) ( log )
## Formula: brood_cells ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 563.5 580.3 -273.8 547.5 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.91613 -0.14733 0.06488 0.13899 0.24831
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 4.987e-01 7.062e-01
## round (Intercept) 6.593e-11 8.120e-06
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2418 0.2148 15.091 <2e-16 ***
## treatment2 0.2373 0.3006 0.790 0.4298
## treatment3 0.5079 0.2996 1.695 0.0901 .
## treatment4 0.1539 0.3015 0.511 0.6097
## treatment5 -0.1094 0.3038 -0.360 0.7189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.714
## treatment3 -0.716 0.511
## treatment4 -0.710 0.507 0.508
## treatment5 -0.702 0.501 0.503 0.500
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(nb.brood)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: brood_cells
## Chisq Df Pr(>Chisq)
## treatment 4.997 4 0.2876
nb.brood.means <- emmeans(nb.brood, "treatment")
pairs(nb.brood.means)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.2373 0.301 Inf -0.790 0.9337
## treatment1 - treatment3 -0.5079 0.300 Inf -1.695 0.4370
## treatment1 - treatment4 -0.1539 0.302 Inf -0.511 0.9864
## treatment1 - treatment5 0.1094 0.304 Inf 0.360 0.9964
## treatment2 - treatment3 -0.2706 0.297 Inf -0.912 0.8925
## treatment2 - treatment4 0.0834 0.299 Inf 0.279 0.9987
## treatment2 - treatment5 0.3467 0.302 Inf 1.149 0.7805
## treatment3 - treatment4 0.3540 0.298 Inf 1.188 0.7586
## treatment3 - treatment5 0.6173 0.301 Inf 2.051 0.2416
## treatment4 - treatment5 0.2633 0.303 Inf 0.870 0.9080
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
Let’s look at the Poisson distribution just for fun.
The AIC is slightly lower, but the standard deviation is much higher. Overall output isn’t much different though.
pois.brood <- glmer(brood_cells ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(pois.brood)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: brood_cells ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 561.5 576.2 -273.8 547.5 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.91625 -0.14711 0.06466 0.13806 0.24579
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.4988 0.7063
## round (Intercept) 0.0000 0.0000
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2418 0.2146 15.104 <2e-16 ***
## treatment2 0.2373 0.3004 0.790 0.430
## treatment3 0.5078 0.2995 1.695 0.090 .
## treatment4 0.1539 0.3014 0.511 0.610
## treatment5 -0.1093 0.3035 -0.360 0.719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.714
## treatment3 -0.716 0.511
## treatment4 -0.710 0.507 0.508
## treatment5 -0.702 0.502 0.503 0.500
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(pois.brood)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: brood_cells
## Chisq Df Pr(>Chisq)
## treatment 5.0048 4 0.2868
Let’s see what happens if we take a look at a boxplot, and check for outliers.
summary(brood$brood_cells)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 23.00 33.50 37.17 50.25 107.00
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = brood_cells, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Brood Cells Count vs. Treatment")
boxplot.stats(brood$brood_cells)$out
## [1] 96 107
out <- boxplot.stats(brood$brood_cells)$out
out_ind <- which(brood$brood_cells %in% c(out))
out_ind
## [1] 33 39
temp<-brood[-c(33,39), ]
ggplot(temp) +
geom_boxplot(aes(x = treatment, y = brood_cells, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Brood Cells Count vs. Treatment, with Outlier Removed")
If we remove the noted outliers and re-run our models, there is not much of an effect.
nb.brood.out <- glmer.nb(brood_cells~ treatment + (1|round/colony), data = temp)
summary(nb.brood.out)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(4835.786) ( log )
## Formula: brood_cells ~ treatment + (1 | round/colony)
## Data: temp
##
## AIC BIC logLik deviance df.resid
## 537.7 554.2 -260.9 521.7 50
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96228 -0.15916 0.06712 0.15545 0.25068
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 4.649e-01 6.818e-01
## round (Intercept) 1.487e-12 1.219e-06
## Number of obs: 58, groups: colony:round, 58; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.24719 0.20824 15.593 <2e-16 ***
## treatment2 0.23260 0.29111 0.799 0.424
## treatment3 0.42725 0.29676 1.440 0.150
## treatment4 0.02757 0.29916 0.092 0.927
## treatment5 -0.10824 0.29434 -0.368 0.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.715
## treatment3 -0.701 0.501
## treatment4 -0.693 0.496 0.486
## treatment5 -0.701 0.502 0.492 0.489
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(nb.brood.out)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: brood_cells
## Chisq Df Pr(>Chisq)
## treatment 4.1042 4 0.3921
pois.brood.out <- glmer(brood_cells ~ treatment + (1|round/colony), data = temp, family = "poisson")
summary(pois.brood.out)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: brood_cells ~ treatment + (1 | round/colony)
## Data: temp
##
## AIC BIC logLik deviance df.resid
## 535.7 550.2 -260.9 521.7 51
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96241 -0.15891 0.06691 0.15430 0.24926
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.465 0.6819
## round (Intercept) 0.000 0.0000
## Number of obs: 58, groups: colony:round, 58; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.24724 0.20785 15.623 <2e-16 ***
## treatment2 0.23255 0.29078 0.800 0.424
## treatment3 0.42721 0.29645 1.441 0.150
## treatment4 0.02756 0.29874 0.092 0.926
## treatment5 -0.10821 0.29379 -0.368 0.713
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.714
## treatment3 -0.700 0.500
## treatment4 -0.693 0.495 0.486
## treatment5 -0.702 0.502 0.492 0.489
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(pois.brood.out)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: brood_cells
## Chisq Df Pr(>Chisq)
## treatment 4.1119 4 0.3911
I will now group the brood count data by treatment in order to plot the averages.
brood.means <- brood %>%
group_by(treatment) %>%
summarize(brd.mean = mean(brood_cells),
brd.sd = sd(brood_cells),
brd.n = length(brood_cells)) %>%
mutate(brd.se = brd.sd/sqrt(brd.n))
as.data.frame(brood.means)
## treatment brd.mean brd.sd brd.n brd.se
## 1 1 33.75000 22.11180 12 6.383128
## 2 2 34.83333 11.98357 12 3.459360
## 3 3 48.83333 26.43632 12 7.631507
## 4 4 38.41667 28.31144 12 8.172810
## 5 5 30.00000 18.18091 12 5.248376
brood.means
## # A tibble: 5 × 5
## treatment brd.mean brd.sd brd.n brd.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 33.8 22.1 12 6.38
## 2 2 34.8 12.0 12 3.46
## 3 3 48.8 26.4 12 7.63
## 4 4 38.4 28.3 12 8.17
## 5 5 30 18.2 12 5.25
ggplot(data = brood.means, aes(x = treatment, y = brd.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(20,60)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = brd.mean - brd.se,
ymax = brd.mean + brd.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Brood Cells",) +
ggtitle("Average Brood Cells per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 12)
Here is a barplot showing average brood cell production per treatment level. We see that Treatment 3 has the highest overall brood production.However, it is not significantly higher according to our mixed effects models.
Look at shape of data.
ggplot(brood, aes(x=eggs, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Eggs Counted")
AIC value for negative binomial distribution is smaller, I will stick with that.
fit1.egg<-fitdistr(brood$eggs, "Poisson")
fit2.egg<-fitdistr(brood$eggs, "negative binomial")
AIC(fit1.egg, fit2.egg)
## df AIC
## fit1.egg 1 897.4400
## fit2.egg 2 367.1369
In the model below we are testing the total number of eggs as a function of dose of pristine, with 0 PPB acting as the intercept (which is good, because it is our control), and holding round as a random variable with individual colony nested within. We are using glmer.nb because it is a mixed effects model, and a negative binomial distribution because it is count data.
nb.egg <- glmer.nb(eggs~ treatment + (1|round/colony), data = brood)
summary(nb.egg)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(775.8463) ( log )
## Formula: eggs ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 388.0 404.8 -186.0 372.0 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.85900 -0.78214 0.04674 0.15731 0.20764
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 2.267 1.505
## round (Intercept) 0.000 0.000
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.34026 0.47945 2.795 0.00518 **
## treatment2 0.02914 0.67009 0.043 0.96532
## treatment3 -0.28786 0.67950 -0.424 0.67183
## treatment4 -0.20255 0.67158 -0.302 0.76296
## treatment5 -0.44438 0.68305 -0.651 0.51531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.698
## treatment3 -0.683 0.494
## treatment4 -0.701 0.499 0.491
## treatment5 -0.677 0.491 0.486 0.488
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(nb.brood)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: brood_cells
## Chisq Df Pr(>Chisq)
## treatment 4.997 4 0.2876
nb.brood.means <- emmeans(nb.brood, "treatment")
pairs(nb.brood.means)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.2373 0.301 Inf -0.790 0.9337
## treatment1 - treatment3 -0.5079 0.300 Inf -1.695 0.4370
## treatment1 - treatment4 -0.1539 0.302 Inf -0.511 0.9864
## treatment1 - treatment5 0.1094 0.304 Inf 0.360 0.9964
## treatment2 - treatment3 -0.2706 0.297 Inf -0.912 0.8925
## treatment2 - treatment4 0.0834 0.299 Inf 0.279 0.9987
## treatment2 - treatment5 0.3467 0.302 Inf 1.149 0.7805
## treatment3 - treatment4 0.3540 0.298 Inf 1.188 0.7586
## treatment3 - treatment5 0.6173 0.301 Inf 2.051 0.2416
## treatment4 - treatment5 0.2633 0.303 Inf 0.870 0.9080
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
Let’s look at the Poisson distribution.
pois.egg <- glmer(eggs ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(pois.egg)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: eggs ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 386.0 400.7 -186.0 372.0 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.85893 -0.78208 0.04659 0.15587 0.20564
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 2.268e+00 1.506e+00
## round (Intercept) 2.742e-10 1.656e-05
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.34019 0.47951 2.795 0.00519 **
## treatment2 0.02911 0.67020 0.043 0.96535
## treatment3 -0.28789 0.67957 -0.424 0.67183
## treatment4 -0.20261 0.67157 -0.302 0.76289
## treatment5 -0.44443 0.68323 -0.650 0.51538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.698
## treatment3 -0.683 0.494
## treatment4 -0.701 0.499 0.491
## treatment5 -0.677 0.491 0.487 0.489
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Let’s see what happens if we take a look at a boxplot, and check for outliers.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 5.000 7.867 13.000 87.000
## [1] 87 36
## [1] 10 21
I removed two outliers above. Let’s see how the models change with these removed.
nb.egg.out <- glmer.nb(eggs~ treatment + (1|round/colony), data = temp)
summary(nb.egg.out)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(758.2865) ( log )
## Formula: eggs ~ treatment + (1 | round/colony)
## Data: temp
##
## AIC BIC logLik deviance df.resid
## 370.5 387.0 -177.2 354.5 50
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.84828 -0.76653 0.03157 0.15271 0.20771
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 2.346e+00 1.532e+00
## round (Intercept) 1.390e-11 3.728e-06
## Number of obs: 58, groups: colony:round, 58; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.33130 0.48731 2.732 0.0063 **
## treatment2 0.02853 0.68098 0.042 0.9666
## treatment3 -0.48397 0.71239 -0.679 0.4969
## treatment4 -0.31919 0.70060 -0.456 0.6487
## treatment5 -0.44937 0.69434 -0.647 0.5175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.698
## treatment3 -0.659 0.479
## treatment4 -0.681 0.486 0.464
## treatment5 -0.676 0.491 0.472 0.476
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
pois.egg.out <- glmer(eggs ~ treatment + (1|round/colony), data = temp, family = "poisson")
summary(pois.egg.out)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: eggs ~ treatment + (1 | round/colony)
## Data: temp
##
## AIC BIC logLik deviance df.resid
## 368.5 382.9 -177.3 354.5 51
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8482 -0.7665 0.0315 0.1521 0.2051
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 2.348e+00 1.532e+00
## round (Intercept) 1.425e-10 1.194e-05
## Number of obs: 58, groups: colony:round, 58; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.33117 0.48739 2.731 0.00631 **
## treatment2 0.02851 0.68110 0.042 0.96661
## treatment3 -0.48392 0.71247 -0.679 0.49701
## treatment4 -0.31922 0.70068 -0.456 0.64869
## treatment5 -0.44935 0.69442 -0.647 0.51757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.698
## treatment3 -0.659 0.479
## treatment4 -0.681 0.486 0.464
## treatment5 -0.676 0.491 0.472 0.476
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
The model outputs didn’t really change.
Let’s look at the averages in a graph.
egg.means <- brood %>%
group_by(treatment) %>%
summarize(egg.mean = mean(eggs),
egg.sd = sd(eggs),
egg.n = length(eggs)) %>%
mutate(egg.se = egg.sd/sqrt(egg.n))
as.data.frame(egg.means)
## treatment egg.mean egg.sd egg.n egg.se
## 1 1 12.000000 24.128444 12 6.965282
## 2 2 8.833333 10.390847 12 2.999579
## 3 3 7.083333 7.609305 12 2.196617
## 4 4 5.833333 5.702206 12 1.646085
## 5 5 5.583333 5.583390 12 1.611786
egg.means
## # A tibble: 5 × 5
## treatment egg.mean egg.sd egg.n egg.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 12 24.1 12 6.97
## 2 2 8.83 10.4 12 3.00
## 3 3 7.08 7.61 12 2.20
## 4 4 5.83 5.70 12 1.65
## 5 5 5.58 5.58 12 1.61
ggplot(data = egg.means, aes(x = treatment, y = egg.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(4,20)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = egg.mean - egg.se,
ymax = egg.mean + egg.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Eggs Counted",) +
ggtitle("Average Eggs Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 12)
ggplot(brood, aes(x= larvae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Larvae")
ggplot(brood, aes(x= live_larvae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Live Larvae")
ggplot(brood, aes(x= dead_larvae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Dead Larvae")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = larvae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Larvae Count vs. Treatment")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = live_larvae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Live Larvae Count vs. Treatment")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = dead_larvae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Dead Larvae Count vs. Treatment")
Let’s look at the poisson distributions. None of them appear to be overdispersed, and all AICs are about 2 points lower than the negative binomial models (below).
ll.p <- glmer(live_larvae ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(ll.p)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: live_larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 523.8 538.4 -254.9 509.8 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.23062 -0.12690 0.05380 0.09937 0.12892
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.716e+00 1.310e+00
## round (Intercept) 1.843e-10 1.357e-05
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.444549 0.400304 6.107 1.02e-09 ***
## treatment2 -0.001154 0.558503 -0.002 0.998
## treatment3 0.568924 0.556113 1.023 0.306
## treatment4 -0.101561 0.560165 -0.181 0.856
## treatment5 -0.227026 0.563122 -0.403 0.687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.713
## treatment3 -0.715 0.511
## treatment4 -0.709 0.507 0.509
## treatment5 -0.701 0.502 0.504 0.501
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
dispersion_glmer(ll.p)
## [1] 1.115205
Anova(ll.p)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: live_larvae
## Chisq Df Pr(>Chisq)
## treatment 2.4466 4 0.6542
dl.p <- glmer(dead_larvae ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(dl.p)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: dead_larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 306.4 321.0 -146.2 292.4 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.95515 -0.70852 -0.04771 0.21096 0.42232
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.626e+00 1.275301
## round (Intercept) 1.598e-06 0.001264
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1271785 0.4651363 0.273 0.7845
## treatment2 0.4835485 0.6296274 0.768 0.4425
## treatment3 1.2648461 0.6105948 2.071 0.0383 *
## treatment4 0.5091232 0.6328365 0.805 0.4211
## treatment5 0.0004361 0.6467050 0.001 0.9995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.711
## treatment3 -0.744 0.542
## treatment4 -0.702 0.524 0.538
## treatment5 -0.687 0.513 0.527 0.511
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00682513 (tol = 0.002, component 1)
dispersion_glmer(dl.p)
## [1] 1.039809
Anova(dl.p)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead_larvae
## Chisq Df Pr(>Chisq)
## treatment 5.9389 4 0.2038
al.p <- glmer(larvae ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(al.p)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 536.4 551.1 -261.2 522.4 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.33186 -0.15569 0.03225 0.12444 0.16265
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.178 1.085
## round (Intercept) 0.000 0.000
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6632356 0.3310252 8.045 8.6e-16 ***
## treatment2 0.0589037 0.4627226 0.127 0.899
## treatment3 0.7341882 0.4594049 1.598 0.110
## treatment4 0.0003391 0.4640456 0.001 0.999
## treatment5 -0.1764449 0.4661457 -0.379 0.705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.713
## treatment3 -0.719 0.513
## treatment4 -0.709 0.506 0.511
## treatment5 -0.703 0.503 0.507 0.501
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
dispersion_glmer(al.p)
## [1] 1.071774
Anova(al.p)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: larvae
## Chisq Df Pr(>Chisq)
## treatment 4.7927 4 0.3092
Let’s look at the negative binomial distributions.
ll <- glmer.nb(live_larvae ~ treatment + (1|round/colony), data = brood)
summary(ll)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1570.91) ( log )
## Formula: live_larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 525.8 542.5 -254.9 509.8 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.22971 -0.12687 0.05413 0.10017 0.13040
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.718e+00 1.311e+00
## round (Intercept) 4.585e-12 2.141e-06
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4434023 0.4040318 6.048 1.47e-09 ***
## treatment2 -0.0004644 0.5617265 -0.001 0.999
## treatment3 0.5697509 0.5591709 1.019 0.308
## treatment4 -0.1010853 0.5635260 -0.179 0.858
## treatment5 -0.2271953 0.5669634 -0.401 0.689
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.714
## treatment3 -0.716 0.513
## treatment4 -0.709 0.508 0.510
## treatment5 -0.699 0.502 0.505 0.501
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(ll)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: live_larvae
## Chisq Df Pr(>Chisq)
## treatment 2.4298 4 0.6572
plot(ll)
dl <- glmer.nb(dead_larvae ~ treatment + (1|round/colony), data = brood)
summary(dl)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(616.6248) ( log )
## Formula: dead_larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 308.4 325.1 -146.2 292.4 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.95503 -0.70833 -0.04769 0.21136 0.42468
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.626e+00 1.275e+00
## round (Intercept) 4.812e-12 2.194e-06
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1276144 0.4653181 0.274 0.7839
## treatment2 0.4834914 0.6298781 0.768 0.4427
## treatment3 1.2651231 0.6108373 2.071 0.0383 *
## treatment4 0.5090631 0.6331008 0.804 0.4214
## treatment5 -0.0005323 0.6469681 -0.001 0.9993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.711
## treatment3 -0.744 0.542
## treatment4 -0.702 0.524 0.538
## treatment5 -0.687 0.513 0.527 0.511
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(dl)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead_larvae
## Chisq Df Pr(>Chisq)
## treatment 5.9407 4 0.2036
al <- glmer.nb(larvae ~ treatment + (1|round/colony), data = brood)
summary(al)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(2583.849) ( log )
## Formula: larvae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 538.4 555.1 -261.2 522.4 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.33142 -0.15582 0.03238 0.12522 0.16504
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.179e+00 1.086e+00
## round (Intercept) 4.915e-14 2.217e-07
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6627975 0.3328142 8.001 1.24e-15 ***
## treatment2 0.0591595 0.4643151 0.127 0.899
## treatment3 0.7345953 0.4608705 1.594 0.111
## treatment4 0.0005223 0.4657910 0.001 0.999
## treatment5 -0.1764631 0.4680737 -0.377 0.706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.714
## treatment3 -0.721 0.515
## treatment4 -0.709 0.507 0.512
## treatment5 -0.703 0.503 0.507 0.502
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(al)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: larvae
## Chisq Df Pr(>Chisq)
## treatment 4.7762 4 0.311
Plot the means.
wide.larv <- gather(brood, life, count, dead_larvae:live_larvae)
wide.larv
## # A tibble: 120 × 17
## round dose treatment replicate colony brood_…¹ honey…² eggs dead_…³ live_…⁴
## <fct> <dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 0 1 11 1.11R2 18 1 14 4 3
## 2 2 0 1 12 1.12R2 3 3 0 0 0
## 3 1 0 1 1 1.1R1 63 5 3 25 4
## 4 2 0 1 1 1.1R2 54 4 5 2 6
## 5 1 0 1 2 1.2R1 18 11 5 4 5
## 6 2 0 1 2 1.2R2 29 10 0 0 0
## 7 1 0 1 3 1.3R1 20 3 3 4 2
## 8 2 0 1 3 1.3R2 32 2 13 3 7
## 9 2 0 1 4 1.4R2 56 4 11 2 11
## 10 2 0 1 5 1.5R2 50 8 87 0 4
## # … with 110 more rows, 7 more variables: dead_drones <dbl>, live_drones <dbl>,
## # drones <dbl>, pupae <dbl>, larvae <dbl>, life <chr>, count <dbl>, and
## # abbreviated variable names ¹brood_cells, ²honey_pot, ³dead_pupae,
## # ⁴live_pupae
wide.sum <- wide.larv %>%
group_by(treatment) %>%
summarize(wide.mean = mean(larvae),
wide.sd = sd(larvae),
wide.n = length(larvae)) %>%
mutate(wide.se = wide.sd/sqrt(wide.n))
wide.sum
## # A tibble: 5 × 5
## treatment wide.mean wide.sd wide.n wide.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 26.9 24.5 24 5.01
## 2 2 20.8 15.9 24 3.24
## 3 3 38.8 24.8 24 5.07
## 4 4 24.7 24.8 24 5.07
## 5 5 19.3 14.4 24 2.93
larv.labs <- c("Live Larvae", "Dead Larvae")
ggplot(data = wide.sum, aes(x = treatment, y = wide.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = wide.mean - wide.se,
ymax = wide.mean + wide.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Larvae Counted",) +
ggtitle("Average Larvae Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 15)
larvae.count <- brood %>%
group_by(treatment) %>%
summarize(larve.c = sum(larvae),
larve.a = mean(larvae),
sd.l = sd(larvae))
larvae.count %>%
kbl(caption = "Total Larvae Produced per Treatment", col.names = c("Treatment", "Total Larvae", "Mean Larvae", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Total Larvae | Mean Larvae | SD |
|---|---|---|---|
| 1 | 323 | 26.92 | 25.08 |
| 2 | 249 | 20.75 | 16.22 |
| 3 | 466 | 38.83 | 25.39 |
| 4 | 296 | 24.67 | 25.41 |
| 5 | 232 | 19.33 | 14.67 |
wide.c.sum <- wide.larv %>%
group_by(treatment, life) %>%
summarize(widec.mean = mean(count),
widec.sd = sd(count),
widec.n = length(count)) %>%
mutate(widec.se = widec.sd/sqrt(widec.n))
larv.labs <- c("Live Larvae", "Dead Larvae")
names(larv.labs) <- c("dead_larvae", "live_larvae")
ggplot(data = wide.c.sum, aes(x = treatment, y = widec.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = widec.mean - widec.se,
ymax = widec.mean + widec.se),
position = position_dodge(0.9), width = 0.6) +
labs(y = "Mean Larvae Counted",) +
ggtitle("Average Larvae Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 20)+
facet_wrap(vars(life),
labeller = labeller(life = larv.labs))
ggplot(brood, aes(x= pupae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Pupae")
ggplot(brood, aes(x= live_pupae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 1, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Live Pupae")
ggplot(brood, aes(x= dead_pupae, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 3, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Dead Pupae")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = pupae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Pupal Cell Count vs. Treatment")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = live_pupae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Live Pupal Cell Count vs. Treatment")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = dead_pupae, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Total Dead Pupal Count vs. Treatment")
Here are the negative binomial distribtuions.
lp <- glmer.nb(live_pupae ~ treatment + (1|round/colony), data = brood)
summary(lp)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(608.8931) ( log )
## Formula: live_pupae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 321.2 337.9 -152.6 305.2 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.02805 -0.90586 0.02918 0.34484 0.49080
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.007 1.004
## round (Intercept) 0.000 0.000
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.12185 0.34282 3.272 0.00107 **
## treatment2 -0.08934 0.48151 -0.186 0.85280
## treatment3 -0.04177 0.47929 -0.087 0.93054
## treatment4 -0.36466 0.48828 -0.747 0.45518
## treatment5 -0.49152 0.49300 -0.997 0.31876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.682
## treatment3 -0.691 0.494
## treatment4 -0.675 0.485 0.487
## treatment5 -0.667 0.481 0.482 0.474
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(lp)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: live_pupae
## Chisq Df Pr(>Chisq)
## treatment 1.542 4 0.8192
dp <- glmer.nb(dead_pupae ~ treatment + (1|round/colony), data = brood)
summary(dp)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(955.2041) ( log )
## Formula: dead_pupae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 364.1 380.9 -174.1 348.1 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.05812 -0.48855 0.01005 0.17369 0.35897
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.1480 1.0715
## round (Intercept) 0.1438 0.3792
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.93488 0.47046 1.987 0.0469 *
## treatment2 0.61386 0.51070 1.202 0.2294
## treatment3 1.26875 0.49928 2.541 0.0110 *
## treatment4 0.08845 0.52593 0.168 0.8664
## treatment5 0.40124 0.51499 0.779 0.4359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.570
## treatment3 -0.594 0.544
## treatment4 -0.548 0.516 0.525
## treatment5 -0.569 0.527 0.540 0.511
Anova(dp)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead_pupae
## Chisq Df Pr(>Chisq)
## treatment 8.5133 4 0.07449 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dpem <- emmeans(dp, "treatment")
pairs(dpem)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.6139 0.511 Inf -1.202 0.7503
## treatment1 - treatment3 -1.2687 0.499 Inf -2.541 0.0817
## treatment1 - treatment4 -0.0885 0.526 Inf -0.168 0.9998
## treatment1 - treatment5 -0.4012 0.515 Inf -0.779 0.9367
## treatment2 - treatment3 -0.6549 0.482 Inf -1.358 0.6545
## treatment2 - treatment4 0.5254 0.510 Inf 1.030 0.8417
## treatment2 - treatment5 0.2126 0.499 Inf 0.426 0.9931
## treatment3 - treatment4 1.1803 0.500 Inf 2.359 0.1268
## treatment3 - treatment5 0.8675 0.487 Inf 1.783 0.3835
## treatment4 - treatment5 -0.3128 0.515 Inf -0.607 0.9740
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
ap <- glmer.nb(pupae ~ treatment + (1|round/colony), data = brood)
summary(ap)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1730.515) ( log )
## Formula: pupae ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 424.1 440.9 -204.1 408.1 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2385 -0.4436 0.0697 0.1893 0.3532
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.85905 0.9268
## round (Intercept) 0.02145 0.1465
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.80423 0.33839 5.332 9.72e-08 ***
## treatment2 0.32065 0.41716 0.769 0.4421
## treatment3 0.74446 0.41178 1.808 0.0706 .
## treatment4 -0.11093 0.42426 -0.261 0.7937
## treatment5 -0.03273 0.42284 -0.077 0.9383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.625
## treatment3 -0.647 0.521
## treatment4 -0.618 0.504 0.511
## treatment5 -0.618 0.505 0.512 0.497
Anova(ap)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: pupae
## Chisq Df Pr(>Chisq)
## treatment 5.9703 4 0.2014
wide.pup <- gather(brood, life, count, dead_pupae:live_pupae)
wide.pup
## # A tibble: 120 × 17
## round dose treatment replicate colony brood_…¹ honey…² eggs dead_…³ live_…⁴
## <fct> <dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 0 1 11 1.11R2 18 1 14 1 0
## 2 2 0 1 12 1.12R2 3 3 0 0 0
## 3 1 0 1 1 1.1R1 63 5 3 8 40
## 4 2 0 1 1 1.1R2 54 4 5 1 49
## 5 1 0 1 2 1.2R1 18 11 5 0 9
## 6 2 0 1 2 1.2R2 29 10 0 6 21
## 7 1 0 1 3 1.3R1 20 3 3 1 9
## 8 2 0 1 3 1.3R2 32 2 13 0 17
## 9 2 0 1 4 1.4R2 56 4 11 2 77
## 10 2 0 1 5 1.5R2 50 8 87 1 46
## # … with 110 more rows, 7 more variables: dead_drones <dbl>, live_drones <dbl>,
## # drones <dbl>, pupae <dbl>, larvae <dbl>, life <chr>, count <dbl>, and
## # abbreviated variable names ¹brood_cells, ²honey_pot, ³dead_larvae,
## # ⁴live_larvae
widep.sum <- wide.pup %>%
group_by(treatment) %>%
summarize(widep.mean = mean(pupae),
widep.sd = sd(pupae),
widep.n = length(pupae)) %>%
mutate(widep.se = widep.sd/sqrt(widep.n))
pup.labs <- c("Live Pupae", "Dead Pupae")
names(pup.labs) <- c("dead_pupae", "live_pupae")
ggplot(data = widep.sum, aes(x = treatment, y = widep.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = widep.mean - widep.se,
ymax = widep.mean + widep.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Pupae Counted",) +
ggtitle("Average Pupae Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 15)
p.count <- brood %>%
group_by(treatment) %>%
summarize(p.c = sum(pupae),
p.a = mean(pupae),
p.l = sd(pupae))
p.count %>%
kbl(caption = "Total Pupae Produced per Treatment", col.names = c("Treatment", "Total Pupae", "Mean Pupae", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Total Pupae | Mean Pupae | SD |
|---|---|---|---|
| 1 | 105 | 8.75 | 8.32 |
| 2 | 146 | 12.17 | 12.63 |
| 3 | 184 | 15.33 | 9.70 |
| 4 | 111 | 9.25 | 13.90 |
| 5 | 102 | 8.50 | 7.83 |
widepc.sum <- wide.pup %>%
group_by(treatment, life) %>%
summarize(widepc.mean = mean(count),
widepc.sd = sd(count),
widepc.n = length(count)) %>%
mutate(widepc.se = widepc.sd/sqrt(widepc.n))
ggplot(data = widepc.sum, aes(x = treatment, y = widepc.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = widepc.mean - widepc.se,
ymax = widepc.mean + widepc.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Pupae Counted",) +
ggtitle("Average Pupae Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 20)+
facet_wrap(vars(life),
labeller = labeller(life = pup.labs))
brood$dead_all <- rowSums(cbind(brood$dead_larvae, brood$dead_pupae, brood$dead_drones))
brood
## # A tibble: 60 × 18
## round dose treatment replicate colony brood_…¹ honey…² eggs dead_…³ live_…⁴
## <fct> <dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 0 1 11 1.11R2 18 1 14 1 0
## 2 2 0 1 12 1.12R2 3 3 0 0 0
## 3 1 0 1 1 1.1R1 63 5 3 8 40
## 4 2 0 1 1 1.1R2 54 4 5 1 49
## 5 1 0 1 2 1.2R1 18 11 5 0 9
## 6 2 0 1 2 1.2R2 29 10 0 6 21
## 7 1 0 1 3 1.3R1 20 3 3 1 9
## 8 2 0 1 3 1.3R2 32 2 13 0 17
## 9 2 0 1 4 1.4R2 56 4 11 2 77
## 10 2 0 1 5 1.5R2 50 8 87 1 46
## # … with 50 more rows, 8 more variables: dead_pupae <dbl>, live_pupae <dbl>,
## # dead_drones <dbl>, live_drones <dbl>, drones <dbl>, pupae <dbl>,
## # larvae <dbl>, dead_all <dbl>, and abbreviated variable names ¹brood_cells,
## # ²honey_pot, ³dead_larvae, ⁴live_larvae
brood %>%
ggplot(aes(x=dead_all, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Dead Brood (all life stages)") +
labs(y = "Count", x = "Dead Brood")
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = dead_all, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("All Dead Brood")
fitdb1 <- fitdistr(brood$dead_all, "normal")
fitdb2 <- fitdistr(brood$dead_all, "Poisson")
fitdb3 <- fitdistr(brood$dead_all, "negative binomial")
AIC(fitdb1, fitdb2, fitdb3)
## df AIC
## fitdb1 2 485.4221
## fitdb2 1 990.5991
## fitdb3 2 418.8662
With this model I get the warning “iteration limit reached” - I will try glmmTMB instead.
db <- glmer.nb(dead_all ~ treatment + (1|round/colony), data = brood)
summary(db)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1780.753) ( log )
## Formula: dead_all ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 421.7 438.5 -202.9 405.7 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.21047 -0.32867 -0.04794 0.15554 0.30562
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.07640 1.0375
## round (Intercept) 0.08034 0.2834
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4486 0.4069 3.560 0.00037 ***
## treatment2 0.6748 0.4687 1.439 0.15001
## treatment3 1.1707 0.4629 2.529 0.01145 *
## treatment4 0.3523 0.4752 0.741 0.45849
## treatment5 0.1724 0.4783 0.361 0.71843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.601
## treatment3 -0.611 0.531
## treatment4 -0.588 0.517 0.523
## treatment5 -0.586 0.513 0.519 0.507
Anova(db)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead_all
## Chisq Df Pr(>Chisq)
## treatment 8.1283 4 0.08699 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The results for this model are very similar, except I now receive no warning message.
library(glmmTMB)
dbtb <- glmmTMB(dead_all ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(dbtb)
## Family: poisson ( log )
## Formula: dead_all ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 419.7 434.4 -202.9 405.7 53
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.07674 1.0377
## round (Intercept) 0.08038 0.2835
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4486 0.4066 3.563 0.000367 ***
## treatment2 0.6749 0.4686 1.440 0.149783
## treatment3 1.1709 0.4627 2.530 0.011392 *
## treatment4 0.3524 0.4750 0.742 0.458095
## treatment5 0.1725 0.4781 0.361 0.718274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(dbtb)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead_all
## Chisq Df Pr(>Chisq)
## treatment 8.1355 4 0.08674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Also, the AIC value for the second model is lower, so I will stick with that.
AIC(db, dbtb)
## df AIC
## db 8 421.7087
## dbtb 7 419.7091
dbm <- emmeans(dbtb, "treatment")
pairs(dbm)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 -0.675 0.469 53 -1.440 0.6049
## treatment1 - treatment3 -1.171 0.463 53 -2.530 0.0992
## treatment1 - treatment4 -0.352 0.475 53 -0.742 0.9455
## treatment1 - treatment5 -0.172 0.478 53 -0.361 0.9963
## treatment2 - treatment3 -0.496 0.451 53 -1.099 0.8064
## treatment2 - treatment4 0.322 0.464 53 0.695 0.9567
## treatment2 - treatment5 0.502 0.467 53 1.075 0.8183
## treatment3 - treatment4 0.819 0.458 53 1.786 0.3923
## treatment3 - treatment5 0.998 0.462 53 2.162 0.2100
## treatment4 - treatment5 0.180 0.473 53 0.380 0.9954
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
db.means <- brood %>%
group_by(treatment) %>%
summarize(db.mean = mean(dead_all),
db.sd = sd(dead_all),
db.n = length(dead_all)) %>%
mutate(db.se = db.sd/sqrt(db.n))
as.data.frame(db.means)
## treatment db.mean db.sd db.n db.se
## 1 1 6.333333 8.783594 12 2.535605
## 2 2 11.916667 11.973139 12 3.456348
## 3 3 17.250000 12.410589 12 3.582629
## 4 4 12.750000 20.937678 12 6.044187
## 5 5 8.000000 9.045340 12 2.611165
ggplot(data = db.means, aes(x = treatment, y = db.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = db.mean - db.se,
ymax = db.mean + db.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Dead Brood",) +
ggtitle("Average Dead Brood of all Life Stages Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 15)
It looks like there is significantly more dead brood of all life stages in treatment 3.
ggplot(brood) +
geom_boxplot(aes(x = treatment, y = honey_pot, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Count of Honey Pots")
fithp1 <- fitdistr(brood$honey_pot, "poisson")
fithp2 <- fitdistr(brood$honey_pot, "negative binomial")
AIC(fithp1, fithp2)
## df AIC
## fithp1 1 353.1115
## fithp2 2 329.8190
Warning: iteration limit reached
hpnb <- glmer.nb(honey_pot ~ treatment + (1|round/colony), data = brood)
summary(hpnb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1120.364) ( log )
## Formula: honey_pot ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 338.4 355.2 -161.2 322.4 52
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79128 -0.48144 0.08576 0.46688 1.23300
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.686e-01 4.106e-01
## round (Intercept) 3.011e-12 1.735e-06
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4641 0.1836 7.974 1.53e-15 ***
## treatment2 0.5448 0.2412 2.259 0.0239 *
## treatment3 0.3771 0.2448 1.540 0.1235
## treatment4 0.3351 0.2460 1.362 0.1731
## treatment5 0.2456 0.2482 0.990 0.3223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.751
## treatment3 -0.732 0.556
## treatment4 -0.733 0.555 0.545
## treatment5 -0.719 0.547 0.539 0.537
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
hptb <- glmmTMB(honey_pot ~ treatment + (1|round/colony), data = brood, family = "poisson")
summary(hptb)
## Family: poisson ( log )
## Formula: honey_pot ~ treatment + (1 | round/colony)
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 336.4 351.1 -161.2 322.4 53
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.696e-01 4.118e-01
## round (Intercept) 1.770e-10 1.331e-05
## Number of obs: 60, groups: colony:round, 60; round, 2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4638 0.1836 7.972 1.57e-15 ***
## treatment2 0.5449 0.2413 2.259 0.0239 *
## treatment3 0.3771 0.2448 1.540 0.1235
## treatment4 0.3352 0.2460 1.362 0.1731
## treatment5 0.2456 0.2482 0.990 0.3224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(hptb)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: honey_pot
## Chisq Df Pr(>Chisq)
## treatment 5.4238 4 0.2465
AIC(hpnb, hptb)
## df AIC
## hpnb 8 338.3956
## hptb 7 336.3973
Anova(hptb)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: honey_pot
## Chisq Df Pr(>Chisq)
## treatment 5.4238 4 0.2465
hpem <- emmeans(hptb, "treatment")
pairs(hpem)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 -0.5449 0.241 53 -2.259 0.1747
## treatment1 - treatment3 -0.3771 0.245 53 -1.540 0.5415
## treatment1 - treatment4 -0.3352 0.246 53 -1.362 0.6540
## treatment1 - treatment5 -0.2456 0.248 53 -0.990 0.8588
## treatment2 - treatment3 0.1678 0.229 53 0.733 0.9479
## treatment2 - treatment4 0.2097 0.230 53 0.912 0.8909
## treatment2 - treatment5 0.2993 0.233 53 1.285 0.7017
## treatment3 - treatment4 0.0419 0.234 53 0.179 0.9998
## treatment3 - treatment5 0.1315 0.237 53 0.556 0.9808
## treatment4 - treatment5 0.0895 0.238 53 0.376 0.9956
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
hp.means <- brood %>%
group_by(treatment) %>%
summarize(hp.mean = mean(honey_pot),
hp.sd = sd(honey_pot),
hp.n = length(honey_pot)) %>%
mutate(hp.se = hp.sd/sqrt(hp.n))
as.data.frame(hp.means)
## treatment hp.mean hp.sd hp.n hp.se
## 1 1 4.750000 3.441062 12 0.9933491
## 2 2 7.916667 3.146667 12 0.9083646
## 3 3 6.916667 4.501683 12 1.2995240
## 4 4 6.500000 3.289100 12 0.9494815
## 5 5 6.083333 4.055486 12 1.1707181
ggplot(data = hp.means, aes(x = treatment, y = hp.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = hp.mean - hp.se,
ymax = hp.mean + hp.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Honey pot Count",) +
ggtitle("Average Honey Pots Counted per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 15)
data_long <- gather(brood, brood_type, count, brood_cells, honey_pot, eggs, larvae, pupae, drones, factor_key=TRUE)
na <- is.na(data_long)
unique(na)
long.means <- data_long %>%
group_by(treatment, brood_type) %>%
summarize(long.mean = mean(count),
long.sd = sd(count),
long.n = length(count)) %>%
mutate(long.se = long.sd/sqrt(long.n))
long.means
## # A tibble: 30 × 6
## # Groups: treatment [5]
## treatment brood_type long.mean long.sd long.n long.se
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 1 brood_cells 33.8 22.1 12 6.38
## 2 1 honey_pot 4.75 3.44 12 0.993
## 3 1 eggs 12 24.1 12 6.97
## 4 1 larvae 26.9 25.1 12 7.24
## 5 1 pupae 8.75 8.32 12 2.40
## 6 1 drones 7.25 7.24 12 2.09
## 7 2 brood_cells 34.8 12.0 12 3.46
## 8 2 honey_pot 7.92 3.15 12 0.908
## 9 2 eggs 8.83 10.4 12 3.00
## 10 2 larvae 20.8 16.2 12 4.68
## # … with 20 more rows
as.data.frame(long.means)
## treatment brood_type long.mean long.sd long.n long.se
## 1 1 brood_cells 33.750000 22.111803 12 6.3831276
## 2 1 honey_pot 4.750000 3.441062 12 0.9933491
## 3 1 eggs 12.000000 24.128444 12 6.9652819
## 4 1 larvae 26.916667 25.076096 12 7.2388455
## 5 1 pupae 8.750000 8.324389 12 2.4030442
## 6 1 drones 7.250000 7.237842 12 2.0893851
## 7 2 brood_cells 34.833333 11.983575 12 3.4593600
## 8 2 honey_pot 7.916667 3.146667 12 0.9083646
## 9 2 eggs 8.833333 10.390847 12 2.9995791
## 10 2 larvae 20.750000 16.220778 12 4.6825352
## 11 2 pupae 12.166667 12.626331 12 3.6449079
## 12 2 drones 7.333333 5.314360 12 1.5341236
## 13 3 brood_cells 48.833333 26.436316 12 7.6315070
## 14 3 honey_pot 6.916667 4.501683 12 1.2995240
## 15 3 eggs 7.083333 7.609305 12 2.1966170
## 16 3 larvae 38.833333 25.387303 12 7.3286831
## 17 3 pupae 15.333333 9.698485 12 2.7997114
## 18 3 drones 6.500000 5.696889 12 1.6445502
## 19 4 brood_cells 38.416667 28.311444 12 8.1728100
## 20 4 honey_pot 6.500000 3.289100 12 0.9494815
## 21 4 eggs 5.833333 5.702206 12 1.6460850
## 22 4 larvae 24.666667 25.406990 12 7.3343663
## 23 4 pupae 9.250000 13.896533 12 4.0115836
## 24 4 drones 9.166667 9.083785 12 2.6222629
## 25 5 brood_cells 30.000000 18.180909 12 5.2483764
## 26 5 honey_pot 6.083333 4.055486 12 1.1707181
## 27 5 eggs 5.583333 5.583390 12 1.6117858
## 28 5 larvae 19.333333 14.674240 12 4.2360883
## 29 5 pupae 8.500000 7.833495 12 2.2613351
## 30 5 drones 8.166667 6.027714 12 1.7400511
long.means
## # A tibble: 30 × 6
## # Groups: treatment [5]
## treatment brood_type long.mean long.sd long.n long.se
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 1 brood_cells 33.8 22.1 12 6.38
## 2 1 honey_pot 4.75 3.44 12 0.993
## 3 1 eggs 12 24.1 12 6.97
## 4 1 larvae 26.9 25.1 12 7.24
## 5 1 pupae 8.75 8.32 12 2.40
## 6 1 drones 7.25 7.24 12 2.09
## 7 2 brood_cells 34.8 12.0 12 3.46
## 8 2 honey_pot 7.92 3.15 12 0.908
## 9 2 eggs 8.83 10.4 12 3.00
## 10 2 larvae 20.8 16.2 12 4.68
## # … with 20 more rows
brood.labs <- c("Total Brood Cells", "Honey Pots", "Eggs", "Total Larvae", "Total Pupae", "Total Drone Count")
names(brood.labs) <- c("brood_cells", "honey_pot", "eggs", "larvae", "pupae", "drones")
ggplot(data = long.means, aes(x = treatment, y = long.mean, fill = treatment)) +
geom_text(aes(y = -1, label = sprintf("%0.1f", round(long.mean, digits = 1))), size = 5)+
geom_col(position = "dodge", col = "black")+
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = long.mean - long.se,
ymax = long.mean + long.se),
position = position_dodge(1), width = 0.6) +
labs(y = "Average Counts",) +
ggtitle("Summary of all Averages of Brood Metrics for all Treatments") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 20)+
facet_wrap(vars(brood_type),
labeller = labeller(brood_type = brood.labs))
Drone emergence marked the end of the experiment. Microcolonies were shut down 1 week after the first drone emerged.
Input the data.
drones <- read_csv("drones.csv", col_types = cols(round = col_factor(levels = c("1",
"2")), treatment = col_factor(levels = c("1",
"2", "3", "4", "5")), replicate = col_factor(levels = c("1",
"2", "3", "4", "5", "7", "9", "11", "12")),
colony_start = col_date(format = "%m/%d/%Y"),
emerge_date = col_date(format = "%m/%d/%Y"),
`alive?` = col_logical()))
drones$colony <- as.factor(drones$colony)
drones$id <- as.factor(drones$id)
drones$radial <- as.numeric(drones$radial)
drones <- subset(drones, select = -c(abdomen_post_ethyl, fat_content, relative_fat, notes))
glimpse(drones)
## Rows: 461
## Columns: 15
## $ round <fct> 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ order_on_sheet <dbl> 31.1, 414.0, 26.0, 28.0, 25.0, 27.0, 4.0, 2.0, 39.0, 40…
## $ treatment <fct> 3, 1, 4, 5, 4, 5, 3, 3, 4, 5, 5, 3, 1, 1, 1, 1, 4, 1, 3…
## $ replicate <fct> 3, 2, 7, 7, 7, 7, 3, 3, 7, 7, 7, 3, 5, 7, 7, 1, 3, 1, 5…
## $ colony_count <dbl> 6, 12, 2, 2, 1, 1, 3, 1, 3, 3, 4, 2, 1, 1, 2, 2, 2, 1, …
## $ colony <fct> 3.3R1, 1.2R2, 4.7R2, 5.7R2, 4.7R2, 5.7R2, 3.3R2, 3.3R2,…
## $ id <fct> 3.3R1.6, 1.2R2.12, 4.7R2.2, 5.7R2.2, 4.7R2.1, 5.7R2.1, …
## $ colony_start <date> 2021-07-24, 2022-08-14, 2022-08-14, 2022-08-14, 2022-0…
## $ emerge_date <date> 2021-09-06, 2022-10-07, 2022-09-13, 2022-09-13, 2022-0…
## $ emerge_time <dbl> 44, 54, 30, 30, 30, 30, 31, 31, 31, 31, 31, 31, 32, 32,…
## $ wet_weight <dbl> 0.08452, 1.16731, 1.02458, 1.20225, 1.17731, 1.12571, 1…
## $ dry_weight <dbl> 0.03446, 0.04294, 0.02455, 0.03071, 0.03613, 0.04358, 0…
## $ `alive?` <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, T…
## $ radial <dbl> 2.556000, 2.532282, 1.986287, 2.169113, 2.398838, 2.640…
## $ abdomen_dry <dbl> 0.01643, 0.02179, 0.01197, 0.01674, 0.01652, 0.02458, 0…
dc = select(drones, -8,-9,-11,-12,-14,-15)
drone.count <- dc %>%
group_by(round, treatment, replicate, colony) %>%
summarize(drone.c = length(colony),
drone.a = mean(colony_count),
sd.dr = sd(colony_count))
ggplot(drone.count, aes(x=drone.c, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 2, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Radial Cell Length(mm)") +
labs(y = "Count", x = "Length(mm)")
drone.t.sum <- drones %>%
group_by(treatment) %>%
summarize(drone.c = length(id),
drone.a = mean(colony_count),
sd.dr = sd(colony_count))
drone.t.sum %>%
kbl(caption = "Total Drones Produced per Treatment", col.names = c("Treatment", "Total Drones", "Mean Drones", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Total Drones | Mean Drones | SD |
|---|---|---|---|
| 1 | 87 | 7.33 | 5.30 |
| 2 | 88 | 7.01 | 12.52 |
| 3 | 78 | 5.83 | 4.47 |
| 4 | 110 | 9.21 | 6.51 |
| 5 | 98 | 6.62 | 4.36 |
The total counts of drones per treatment doe not appear to vary significantly.
dcn <- glmer.nb(drone.c ~ treatment + (1|round) + (1|colony), data = drone.count)
summary(dcn)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1645.602) ( log )
## Formula: drone.c ~ treatment + (1 | round) + (1 | colony)
## Data: drone.count
##
## AIC BIC logLik deviance df.resid
## 342.7 358.6 -163.4 326.7 46
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.50701 -0.60858 -0.00086 0.32065 0.87558
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.3325 0.5767
## round (Intercept) 0.3279 0.5726
## Number of obs: 54, groups: colony, 54; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.62739 0.47027 3.461 0.000539 ***
## treatment2 0.00255 0.30806 0.008 0.993397
## treatment3 -0.27250 0.30748 -0.886 0.375491
## treatment4 0.05116 0.30841 0.166 0.868245
## treatment5 0.17850 0.30998 0.576 0.564717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.351
## treatment3 -0.349 0.529
## treatment4 -0.340 0.525 0.528
## treatment5 -0.345 0.527 0.526 0.522
Anova(dcn)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: drone.c
## Chisq Df Pr(>Chisq)
## treatment 2.4143 4 0.6601
de = select(drones, 1,3,4,5,6,7,10)
de$time <- as.numeric(de$emerge_time)
em.s <- drones %>%
group_by(treatment) %>%
summarize(em.c = length(id),
em.a = mean(emerge_time),
em.s = sd(emerge_time))
em.s
## # A tibble: 5 × 4
## treatment em.c em.a em.s
## <fct> <int> <dbl> <dbl>
## 1 1 87 39.3 5.68
## 2 2 88 39.4 2.67
## 3 3 78 39.9 4.14
## 4 4 110 37.2 3.31
## 5 5 98 38.4 3.77
em.s %>%
kbl(caption = "Total Drones Produced per Treatment", col.names = c("Treatment", "Total Drones", "Mean Emergence Time", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Total Drones | Mean Emergence Time | SD |
|---|---|---|---|
| 1 | 87 | 39.26 | 5.68 |
| 2 | 88 | 39.44 | 2.67 |
| 3 | 78 | 39.87 | 4.14 |
| 4 | 110 | 37.18 | 3.31 |
| 5 | 98 | 38.44 | 3.77 |
em.s2 <- drones %>%
group_by(treatment) %>%
summarize(em.c = length(id),
em.a = mean(emerge_time),
em.s = sd(emerge_time)) %>%
mutate(em.se = em.s/sqrt(em.c))
ggplot(data = em.s2, aes(x = treatment, y = em.a, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(35,42)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = em.a - em.se,
ymax = em.a + em.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Average Drone Emergence Time (days)",) +
ggtitle("Average Time Until First Drone Emergence (days) Per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 12)
Treatment level 4 has the shortest time to drone production on average. However, based on our best-fitting model with a poisson distribtuion this does not appear to be statistically significant.
den <- glm.nb(emerge_time ~ treatment, data = de)
summary(den)
##
## Call:
## glm.nb(formula = emerge_time ~ treatment, data = de, init.theta = 1735462.952,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.46262 -0.39763 -0.04224 0.27498 2.72579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.670317 0.017110 214.515 <2e-16 ***
## treatment2 0.004544 0.024101 0.189 0.8505
## treatment3 0.015352 0.024785 0.619 0.5357
## treatment4 -0.054498 0.023179 -2.351 0.0187 *
## treatment5 -0.021251 0.023627 -0.899 0.3684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1735463) family taken to be 1)
##
## Null deviance: 194.68 on 460 degrees of freedom
## Residual deviance: 183.18 on 456 degrees of freedom
## AIC: 2727.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1735463
## Std. Err.: 12111257
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -2715.674
Anova(den)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge_time
## LR Chisq Df Pr(>Chisq)
## treatment 11.502 4 0.02146 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
e1 <- glmer.nb(emerge_time ~ treatment + (1|round/colony), data = de)
summary(e1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(10059300) ( log )
## Formula: emerge_time ~ treatment + (1 | round/colony)
## Data: de
##
## AIC BIC logLik deviance df.resid
## 2695.9 2729.0 -1340.0 2679.9 453
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.06701 -0.27094 0.03519 0.27076 2.39290
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 3.960e-03 6.293e-02
## round (Intercept) 1.024e-10 1.012e-05
## Number of obs: 461, groups: colony:round, 54; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.668305 0.024390 150.405 <2e-16 ***
## treatment2 -0.003217 0.035266 -0.091 0.927
## treatment3 0.009388 0.036549 0.257 0.797
## treatment4 -0.042101 0.033090 -1.272 0.203
## treatment5 -0.015709 0.034658 -0.453 0.650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.660
## treatment3 -0.654 0.426
## treatment4 -0.641 0.417 0.417
## treatment5 -0.664 0.436 0.432 0.426
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(e1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: emerge_time
## Chisq Df Pr(>Chisq)
## treatment 2.4787 4 0.6485
e2 <- glmer(emerge_time ~ treatment + (1|round/colony), data = de, family = "poisson")
summary(e2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: emerge_time ~ treatment + (1 | round/colony)
## Data: de
##
## AIC BIC logLik deviance df.resid
## 2693.9 2722.9 -1340.0 2679.9 454
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.06703 -0.27092 0.03522 0.27078 2.39292
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 3.961e-03 6.293e-02
## round (Intercept) 5.677e-10 2.383e-05
## Number of obs: 461, groups: colony:round, 54; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.668300 0.028794 127.398 <2e-16 ***
## treatment2 -0.003221 0.039322 -0.082 0.935
## treatment3 0.009403 0.040040 0.235 0.814
## treatment4 -0.042099 0.040208 -1.047 0.295
## treatment5 -0.015714 0.039328 -0.400 0.689
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.730
## treatment3 -0.718 0.525
## treatment4 -0.721 0.524 0.516
## treatment5 -0.733 0.535 0.526 0.530
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(e2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: emerge_time
## Chisq Df Pr(>Chisq)
## treatment 2.0219 4 0.7317
dispersion_glmer(e2)
## [1] 0.4818671
AIC(den, e1, e2)
## df AIC
## den 6 2727.674
## e1 8 2695.930
## e2 7 2693.929
Drone radial cell length is a body size metric used to make inferences about overall body quality.
Look at shape of data
ggplot(drones, aes(x=radial, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.05,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Radial Cell Length(mm)") +
labs(y = "Count", x = "Length(mm)")
Looks relatively normal, but a Shapiro test tells us otherwise.
shapiro.test(drones$radial)
##
## Shapiro-Wilk normality test
##
## data: drones$radial
## W = 0.98448, p-value = 0.0001093
Let’s get rid of NAs and check for outliers.
R found us 10 outliers. When we remove them, the Shapiro test looks much better. Our histogram is also more normal.
d.na <- na.omit(drones)
ggplot(d.na) +
geom_boxplot(aes(x = treatment, y = radial, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Drone Radial Length(mm)")
boxplot.stats(d.na$radial)$out
## [1] 1.931965 1.939590 1.692364 3.154232 1.716938 1.874962 1.865854 3.059268
## [9] 1.886897 1.894846
out.rad <- boxplot.stats(d.na$radial)$out
out_ind.rad <- which(d.na$radial %in% c(1.874962, 3.154232, 1.716938, 1.865854, 1.939590, 1.886897, 1.894846, 1.931965, 3.059268, 1.692364))
temp.rad<-d.na[-c(194, 200, 218, 227, 247, 252, 258, 299, 319, 342), ]
shapiro.test(temp.rad$radial)
##
## Shapiro-Wilk normality test
##
## data: temp.rad$radial
## W = 0.98479, p-value = 0.0001689
ggplot(temp.rad, aes(x=radial, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.05,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Radial Cell Length(mm): Outliers Removed") +
labs(y = "Count", x = "Length(mm)")
ggplot(temp.rad) +
geom_boxplot(aes(x = treatment, y = radial, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Drone Radial Length(mm): Outliers Removed")
If we treat our original data as normal (lmer), in a mixed effects model we see that treatment three continues the trend of having the most impact, but our statistics are not showing any significance.
rad.mod <- lmer(radial ~ treatment + (1|round/colony), data = d.na)
summary(rad.mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: radial ~ treatment + (1 | round/colony)
## Data: d.na
##
## REML criterion at convergence: -151.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8823 -0.5816 0.0513 0.6762 3.6159
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.004700 0.06856
## round (Intercept) 0.003018 0.05493
## Residual 0.036616 0.19135
## Number of obs: 442, groups: colony:round, 51; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.52721 0.05314 47.553
## treatment2 -0.02840 0.04505 -0.630
## treatment3 -0.09692 0.04693 -2.065
## treatment4 -0.04030 0.04591 -0.878
## treatment5 -0.03511 0.04536 -0.774
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.469
## treatment3 -0.455 0.518
## treatment4 -0.434 0.523 0.502
## treatment5 -0.460 0.534 0.513 0.521
Anova(rad.mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 4.5553 4 0.336
mean.rad.normal <- emmeans(rad.mod, "treatment")
pairs((mean.rad.normal))
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.02840 0.0454 38.0 0.626 0.9699
## treatment1 - treatment3 0.09692 0.0473 40.3 2.047 0.2628
## treatment1 - treatment4 0.04030 0.0463 32.4 0.871 0.9054
## treatment1 - treatment5 0.03511 0.0456 34.2 0.770 0.9375
## treatment2 - treatment3 0.06852 0.0454 42.4 1.509 0.5624
## treatment2 - treatment4 0.01190 0.0449 34.3 0.265 0.9989
## treatment2 - treatment5 0.00671 0.0438 35.9 0.153 0.9999
## treatment3 - treatment4 -0.05662 0.0470 36.7 -1.205 0.7485
## treatment3 - treatment5 -0.06181 0.0458 38.3 -1.350 0.6625
## treatment4 - treatment5 -0.00519 0.0451 31.0 -0.115 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
This is also still true if we run a mixed effects model using our outlier-excluded data as well.
rad.mod.out <- lmer(radial ~ treatment + (1|round/colony), data = temp.rad)
summary(rad.mod.out)
## Linear mixed model fit by REML ['lmerMod']
## Formula: radial ~ treatment + (1 | round/colony)
## Data: temp.rad
##
## REML criterion at convergence: -146.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9143 -0.5859 0.0482 0.6694 3.6213
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.005319 0.07293
## round (Intercept) 0.003301 0.05746
## Residual 0.036467 0.19096
## Number of obs: 432, groups: colony:round, 51; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.52517 0.05526 45.698
## treatment2 -0.02378 0.04665 -0.510
## treatment3 -0.10136 0.04863 -2.084
## treatment4 -0.03776 0.04769 -0.792
## treatment5 -0.03524 0.04712 -0.748
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.468
## treatment3 -0.453 0.518
## treatment4 -0.431 0.523 0.501
## treatment5 -0.457 0.533 0.512 0.518
Anova(rad.mod.out)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 4.809 4 0.3075
rad.out.mean <- emmeans(rad.mod.out, "treatment")
pairs(rad.out.mean)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.02378 0.0469 38.1 0.506 0.9862
## treatment1 - treatment3 0.10136 0.0490 41.1 2.068 0.2535
## treatment1 - treatment4 0.03776 0.0480 33.2 0.786 0.9329
## treatment1 - treatment5 0.03524 0.0473 34.7 0.744 0.9444
## treatment2 - treatment3 0.07759 0.0470 42.9 1.651 0.4743
## treatment2 - treatment4 0.01399 0.0466 34.9 0.300 0.9981
## treatment2 - treatment5 0.01146 0.0454 36.0 0.252 0.9991
## treatment3 - treatment4 -0.06360 0.0488 37.9 -1.304 0.6902
## treatment3 - treatment5 -0.06613 0.0475 38.9 -1.391 0.6369
## treatment4 - treatment5 -0.00253 0.0469 31.8 -0.054 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
If we switch to a glmer (rad.mod1) using our original data, this trend holds.
rad.mod1 <- glmer(radial ~ treatment + (1|round/colony), data = d.na, family = "gaussian" (link = "identity"))
summary(rad.mod1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: radial ~ treatment + (1 | round/colony)
## Data: d.na
##
## REML criterion at convergence: -151.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8823 -0.5816 0.0513 0.6762 3.6159
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.004700 0.06856
## round (Intercept) 0.003018 0.05493
## Residual 0.036616 0.19135
## Number of obs: 442, groups: colony:round, 51; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.52721 0.05314 47.553
## treatment2 -0.02840 0.04505 -0.630
## treatment3 -0.09692 0.04693 -2.065
## treatment4 -0.04030 0.04591 -0.878
## treatment5 -0.03511 0.04536 -0.774
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.469
## treatment3 -0.455 0.518
## treatment4 -0.434 0.523 0.502
## treatment5 -0.460 0.534 0.513 0.521
Anova(rad.mod1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 4.5553 4 0.336
mean.rad.abn <- emmeans(rad.mod1, "treatment")
pairs((mean.rad.abn))
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.02840 0.0454 38.0 0.626 0.9699
## treatment1 - treatment3 0.09692 0.0473 40.3 2.047 0.2628
## treatment1 - treatment4 0.04030 0.0463 32.4 0.871 0.9054
## treatment1 - treatment5 0.03511 0.0456 34.2 0.770 0.9375
## treatment2 - treatment3 0.06852 0.0454 42.4 1.509 0.5624
## treatment2 - treatment4 0.01190 0.0449 34.3 0.265 0.9989
## treatment2 - treatment5 0.00671 0.0438 35.9 0.153 0.9999
## treatment3 - treatment4 -0.05662 0.0470 36.7 -1.205 0.7485
## treatment3 - treatment5 -0.06181 0.0458 38.3 -1.350 0.6625
## treatment4 - treatment5 -0.00519 0.0451 31.0 -0.115 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
However, when I run a simple linear model with no mixed effects (rad.mod3), we do see that treatment 3 is correlated with significantly smaller radial cell lengths. An AIC comparison of the model containing mixed effects (rad.mod1) and not containing mixed effects (rad.mod3), shows that removing these better fits the data.
rad.mod3 <- lm(radial ~ treatment, data = d.na)
summary(rad.mod3)
##
## Call:
## lm(formula = radial ~ treatment, data = d.na)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.77354 -0.12232 0.02293 0.12435 0.75783
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.51334 0.02192 114.647 < 2e-16 ***
## treatment2 -0.04283 0.03073 -1.394 0.164128
## treatment3 -0.11694 0.03227 -3.624 0.000324 ***
## treatment4 -0.04744 0.02941 -1.613 0.107492
## treatment5 -0.04257 0.03017 -1.411 0.158884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2009 on 437 degrees of freedom
## Multiple R-squared: 0.02975, Adjusted R-squared: 0.02087
## F-statistic: 3.349 on 4 and 437 DF, p-value: 0.01022
Anova(rad.mod3)
## Anova Table (Type II tests)
##
## Response: radial
## Sum Sq Df F value Pr(>F)
## treatment 0.5409 4 3.3494 0.01022 *
## Residuals 17.6417 437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3.mean <- emmeans(rad.mod3, "treatment")
pairs(mod3.mean)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.04283 0.0307 437 1.394 0.6321
## treatment1 - treatment3 0.11694 0.0323 437 3.624 0.0030
## treatment1 - treatment4 0.04744 0.0294 437 1.613 0.4899
## treatment1 - treatment5 0.04257 0.0302 437 1.411 0.6207
## treatment2 - treatment3 0.07411 0.0320 437 2.315 0.1420
## treatment2 - treatment4 0.00460 0.0291 437 0.158 0.9999
## treatment2 - treatment5 -0.00026 0.0299 437 -0.009 1.0000
## treatment3 - treatment4 -0.06950 0.0307 437 -2.261 0.1598
## treatment3 - treatment5 -0.07437 0.0315 437 -2.363 0.1275
## treatment4 - treatment5 -0.00486 0.0285 437 -0.171 0.9998
##
## P value adjustment: tukey method for comparing a family of 5 estimates
AIC(rad.mod1, rad.mod3)
## df AIC
## rad.mod1 8 -135.3494
## rad.mod3 6 -157.3589
Plot of the averages
rad.means <- d.na %>%
group_by(treatment) %>%
summarize(rad.mean = mean(radial),
rad.sd = sd(radial),
rad.n = length(radial)) %>%
mutate(rad.se = rad.sd/sqrt(rad.n))
ggplot(data = rad.means, aes(x = treatment, y = rad.mean, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(2,2.7)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = rad.mean - rad.se,
ymax = rad.mean + rad.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Average Radial Cell Length (mm)",) +
ggtitle("Drone Radial Cell Length(mm) per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 12)
radial.sum <- d.na %>%
group_by(treatment) %>%
summarize(rad.c = length(id),
rad.a = mean(radial),
rad.sd = sd(radial))
radial.sum %>%
kbl(caption = "Drone Radial Cell Length per Treatment", col.names = c("Treatment", "Drones (N)", "Mean Radial Cell Length (mm)", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Drones (N) | Mean Radial Cell Length (mm) | SD |
|---|---|---|---|
| 1 | 84 | 2.51 | 0.17 |
| 2 | 87 | 2.47 | 0.18 |
| 3 | 72 | 2.40 | 0.24 |
| 4 | 105 | 2.47 | 0.23 |
| 5 | 94 | 2.47 | 0.17 |
The dry weights data column had a couple outliers and a typo, I fixed that in the chunk below.
ggplot(d.na) +
geom_boxplot(aes(x = treatment, y = dry_weight, fill = treatment)) +
coord_cartesian(ylim=c(0, 0.1)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Drone Dry Weight (g)")
ggplot(d.na) +
geom_histogram(aes(x = dry_weight, fill = treatment),
position = "identity", binwidth = 0.025 ,col=I("black")) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Drone Dry Weight (g)")
range(d.na$dry_weight)
## [1] 0.01459 0.45700
which(d.na$dry_weight %in% c(0.457, 0.08262))
## [1] 12 54
temp.dry<- d.na[-c(12, 54),]
range(temp.dry$dry_weight)
## [1] 0.01459 0.06588
ggplot(temp.dry) +
geom_histogram(aes(x = dry_weight, fill = treatment),
position = "identity", binwidth = 0.001 ,col=I("black")) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Drone Dry Weight (g): Outliers Removed") +
scale_x_continuous(breaks = seq(0.01, 0.07, by=0.01))
shapiro.test(temp.dry$dry_weight)
##
## Shapiro-Wilk normality test
##
## data: temp.dry$dry_weight
## W = 0.9941, p-value = 0.08619
We see in a mixed effects model that dry weights for drones are significantly smaller in treatment 3.
dry.mod <- lmer(dry_weight~ treatment + (1|round/colony), data = temp.dry)
summary(dry.mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: dry_weight ~ treatment + (1 | round/colony)
## Data: temp.dry
##
## REML criterion at convergence: -2932.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.98470 -0.60514 -0.03077 0.65998 3.10437
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 6.718e-06 0.002592
## round (Intercept) 1.214e-06 0.001102
## Residual 6.175e-05 0.007858
## Number of obs: 440, groups: colony:round, 51; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0419017 0.0015843 26.449
## treatment2 -0.0027990 0.0017726 -1.579
## treatment3 -0.0055106 0.0018513 -2.977
## treatment4 -0.0013233 0.0018004 -0.735
## treatment5 -0.0007752 0.0017817 -0.435
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.612
## treatment3 -0.590 0.515
## treatment4 -0.581 0.526 0.503
## treatment5 -0.604 0.534 0.512 0.524
Anova(dry.mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 11.167 4 0.02475 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dry.emm <- emmeans(dry.mod, "treatment")
pairs(dry.emm)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.002799 0.00179 37.5 1.564 0.5292
## treatment1 - treatment3 0.005511 0.00188 40.0 2.937 0.0412
## treatment1 - treatment4 0.001323 0.00182 31.4 0.728 0.9483
## treatment1 - treatment5 0.000775 0.00179 33.7 0.432 0.9924
## treatment2 - treatment3 0.002712 0.00180 42.5 1.510 0.5617
## treatment2 - treatment4 -0.001476 0.00177 33.5 -0.831 0.9189
## treatment2 - treatment5 -0.002024 0.00172 35.7 -1.175 0.7654
## treatment3 - treatment4 -0.004187 0.00187 36.0 -2.243 0.1873
## treatment3 - treatment5 -0.004735 0.00181 38.2 -2.619 0.0868
## treatment4 - treatment5 -0.000548 0.00177 30.1 -0.309 0.9979
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
dry.means <- temp.dry %>%
group_by(treatment) %>%
summarize(dry.mean = mean(dry_weight),
dry.sd = sd(dry_weight),
dry.n = length(dry_weight)) %>%
mutate(dry.se = dry.sd/sqrt(dry.n))
ggplot(data = dry.means, aes(x = treatment, y = dry.mean, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(0.03,0.05)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = dry.mean - dry.se,
ymax = dry.mean + dry.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Drone Dry Weight(g)",) +
ggtitle("Average Drone Dry Weight(g) per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 12)
dry.sum <- d.na %>%
group_by(treatment) %>%
summarize(dry.c = length(id),
dry.a = mean(dry_weight),
dry.sd = sd(dry_weight))
dry.sum %>%
kbl(caption = "Drone Dry Weight per Treatment", col.names = c("Treatment", "Drones (N)", "Mean Dry Weight (g)", "SD"), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria", position = "left")
| Treatment | Drones (N) | Mean Dry Weight (g) | SD |
|---|---|---|---|
| 1 | 84 | 0.05 | 0.05 |
| 2 | 87 | 0.04 | 0.01 |
| 3 | 72 | 0.04 | 0.01 |
| 4 | 105 | 0.04 | 0.01 |
| 5 | 94 | 0.04 | 0.01 |
content pending
We were not expecting to see significant mortality in workers. We expected we might see some sub-lethal effects, such as decreased dry weight.
I still need to analyze the last replicate from round 2
workers <- read_csv("workers.analyze.csv",
col_types = cols(round = col_factor(levels = c("1",
"2")), treatment = col_factor(levels = c("1",
"2", "3", "4", "5")), replicate = col_factor(levels = c("1",
"2", "3", "4", "5", "7", "9", "11",
"12")), start_date = col_date(format = "%m/%d/%Y"),
mortality_date = col_date(format = "%m/%d/%Y"),
survived = col_logical(), shut_down_dt = col_datetime(format = "%m/%d/%Y %H:%M")))
workers$q_origin <- as.factor(workers$q_origin)
workers$dry_weight <- as.numeric(workers$`dry weight`)
workers$colony <- as.factor(workers$colony)
workers <- subset(workers, select = -c(replaced, mortality_date, `dry weight`, shut_down_dt, ending_wet_weight, start_wet_weight, start_date))
glimpse(workers)
## Rows: 300
## Columns: 8
## $ round <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2,…
## $ treatment <fct> 1, 3, 3, 1, 2, 4, 1, 3, 1, 3, 1, 4, 5, 5, 5, 1, 2, 3, 1, 4,…
## $ replicate <fct> 9, 2, 2, 9, 3, 11, 11, 9, 9, 1, 3, 5, 1, 2, 11, 3, 4, 9, 2,…
## $ colony <fct> 1.9R2, 3.2R2, 3.2R2, 1.9R2, 2.3R2, 4.11R2, 1.11R2, 3.9R2, 1…
## $ q_origin <fct> B1, B5, B5, B1, B3, B1, B1, B1, B1, B1, B3, B4, B1, B5, B1,…
## $ days_alive <dbl> 55, 44, 44, 55, 42, 40, 45, 41, 54, 48, 45, 4, 44, 45, 41, …
## $ survived <lgl> FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE, TR…
## $ dry_weight <dbl> 0.017610, 0.018760, 0.019900, 0.021080, 0.021800, 0.022340,…
work.mean <- workers %>%
group_by(round, treatment, colony, q_origin, survived) %>%
summarize(w.m = length(survived))
as.data.frame(work.mean)
## round treatment colony q_origin survived w.m
## 1 1 1 1.1R1 K1 TRUE 5
## 2 1 1 1.2R1 K1 TRUE 5
## 3 1 1 1.3R1 K1 TRUE 5
## 4 1 2 2.1R1 K1 TRUE 5
## 5 1 2 2.2R1 K1 TRUE 5
## 6 1 2 2.3R1 K1 TRUE 4
## 7 1 2 2.3R1 K2 TRUE 1
## 8 1 3 3.1R1 K1 TRUE 5
## 9 1 3 3.2R1 K1 TRUE 5
## 10 1 3 3.3R1 K2 TRUE 1
## 11 1 3 3.3R1 K3 TRUE 4
## 12 1 4 4.1R1 K1 TRUE 5
## 13 1 4 4.2R1 K1 TRUE 5
## 14 1 4 4.3R1 K3 TRUE 5
## 15 1 5 5.1R1 K1 FALSE 1
## 16 1 5 5.1R1 K1 TRUE 4
## 17 1 5 5.2R1 K1 TRUE 5
## 18 1 5 5.3R1 K2 FALSE 1
## 19 1 5 5.3R1 K2 TRUE 4
## 20 2 1 1.11R2 B1 TRUE 5
## 21 2 1 1.12R2 B1 FALSE 5
## 22 2 1 1.1R2 B1 TRUE 5
## 23 2 1 1.2R2 B5 FALSE 1
## 24 2 1 1.2R2 B5 TRUE 4
## 25 2 1 1.3R2 B3 TRUE 5
## 26 2 1 1.4R2 B5 FALSE 4
## 27 2 1 1.4R2 B5 TRUE 1
## 28 2 1 1.5R2 B4 FALSE 1
## 29 2 1 1.5R2 B4 TRUE 4
## 30 2 1 1.7R2 B1 TRUE 5
## 31 2 1 1.9R2 B1 FALSE 3
## 32 2 1 1.9R2 B1 TRUE 2
## 33 2 2 2.11R2 B1 TRUE 5
## 34 2 2 2.12R2 B1 TRUE 5
## 35 2 2 2.1R2 B1 TRUE 5
## 36 2 2 2.2R2 B5 FALSE 1
## 37 2 2 2.2R2 B5 TRUE 4
## 38 2 2 2.3R2 B3 TRUE 5
## 39 2 2 2.4R2 B5 FALSE 5
## 40 2 2 2.5R2 B4 FALSE 1
## 41 2 2 2.5R2 B4 TRUE 4
## 42 2 2 2.7R2 B1 TRUE 5
## 43 2 2 2.9R2 B1 TRUE 5
## 44 2 3 3.11R2 B1 TRUE 6
## 45 2 3 3.12R2 B1 TRUE 5
## 46 2 3 3.1R2 B1 TRUE 5
## 47 2 3 3.2R2 B5 TRUE 5
## 48 2 3 3.3R2 B3 TRUE 5
## 49 2 3 3.4R2 B5 FALSE 2
## 50 2 3 3.4R2 B5 TRUE 3
## 51 2 3 3.5R2 B4 FALSE 1
## 52 2 3 3.5R2 B4 TRUE 4
## 53 2 3 3.7R2 B1 TRUE 5
## 54 2 3 3.9R2 B1 TRUE 5
## 55 2 4 4.11R2 B1 TRUE 5
## 56 2 4 4.12R2 B1 FALSE 1
## 57 2 4 4.12R2 B1 TRUE 4
## 58 2 4 4.1R2 B1 TRUE 5
## 59 2 4 4.2R2 B5 FALSE 2
## 60 2 4 4.2R2 B5 TRUE 3
## 61 2 4 4.3R2 B3 TRUE 5
## 62 2 4 4.4R2 B5 TRUE 5
## 63 2 4 4.5R2 B4 FALSE 4
## 64 2 4 4.5R2 B4 TRUE 1
## 65 2 4 4.7R2 B1 TRUE 5
## 66 2 4 4.9R2 B1 TRUE 4
## 67 2 5 5.11R2 B1 TRUE 5
## 68 2 5 5.12R2 B1 TRUE 5
## 69 2 5 5.1R2 B1 TRUE 5
## 70 2 5 5.2R2 B5 FALSE 1
## 71 2 5 5.2R2 B5 TRUE 4
## 72 2 5 5.3R2 B3 FALSE 5
## 73 2 5 5.4R2 B5 TRUE 5
## 74 2 5 5.5R2 B4 TRUE 5
## 75 2 5 5.7R2 B1 TRUE 5
## 76 2 5 5.9R2 B1 TRUE 5
mortality <- spread(work.mean, survived, w.m)
mortality
## # A tibble: 62 × 6
## # Groups: round, treatment, colony, q_origin [62]
## round treatment colony q_origin `FALSE` `TRUE`
## <fct> <fct> <fct> <fct> <int> <int>
## 1 1 1 1.1R1 K1 NA 5
## 2 1 1 1.2R1 K1 NA 5
## 3 1 1 1.3R1 K1 NA 5
## 4 1 2 2.1R1 K1 NA 5
## 5 1 2 2.2R1 K1 NA 5
## 6 1 2 2.3R1 K1 NA 4
## 7 1 2 2.3R1 K2 NA 1
## 8 1 3 3.1R1 K1 NA 5
## 9 1 3 3.2R1 K1 NA 5
## 10 1 3 3.3R1 K2 NA 1
## # … with 52 more rows
mortality[is.na(mortality)] = 0
as.data.frame(mortality)
## round treatment colony q_origin FALSE TRUE
## 1 1 1 1.1R1 K1 0 5
## 2 1 1 1.2R1 K1 0 5
## 3 1 1 1.3R1 K1 0 5
## 4 1 2 2.1R1 K1 0 5
## 5 1 2 2.2R1 K1 0 5
## 6 1 2 2.3R1 K1 0 4
## 7 1 2 2.3R1 K2 0 1
## 8 1 3 3.1R1 K1 0 5
## 9 1 3 3.2R1 K1 0 5
## 10 1 3 3.3R1 K2 0 1
## 11 1 3 3.3R1 K3 0 4
## 12 1 4 4.1R1 K1 0 5
## 13 1 4 4.2R1 K1 0 5
## 14 1 4 4.3R1 K3 0 5
## 15 1 5 5.1R1 K1 1 4
## 16 1 5 5.2R1 K1 0 5
## 17 1 5 5.3R1 K2 1 4
## 18 2 1 1.11R2 B1 0 5
## 19 2 1 1.12R2 B1 5 0
## 20 2 1 1.1R2 B1 0 5
## 21 2 1 1.2R2 B5 1 4
## 22 2 1 1.3R2 B3 0 5
## 23 2 1 1.4R2 B5 4 1
## 24 2 1 1.5R2 B4 1 4
## 25 2 1 1.7R2 B1 0 5
## 26 2 1 1.9R2 B1 3 2
## 27 2 2 2.11R2 B1 0 5
## 28 2 2 2.12R2 B1 0 5
## 29 2 2 2.1R2 B1 0 5
## 30 2 2 2.2R2 B5 1 4
## 31 2 2 2.3R2 B3 0 5
## 32 2 2 2.4R2 B5 5 0
## 33 2 2 2.5R2 B4 1 4
## 34 2 2 2.7R2 B1 0 5
## 35 2 2 2.9R2 B1 0 5
## 36 2 3 3.11R2 B1 0 6
## 37 2 3 3.12R2 B1 0 5
## 38 2 3 3.1R2 B1 0 5
## 39 2 3 3.2R2 B5 0 5
## 40 2 3 3.3R2 B3 0 5
## 41 2 3 3.4R2 B5 2 3
## 42 2 3 3.5R2 B4 1 4
## 43 2 3 3.7R2 B1 0 5
## 44 2 3 3.9R2 B1 0 5
## 45 2 4 4.11R2 B1 0 5
## 46 2 4 4.12R2 B1 1 4
## 47 2 4 4.1R2 B1 0 5
## 48 2 4 4.2R2 B5 2 3
## 49 2 4 4.3R2 B3 0 5
## 50 2 4 4.4R2 B5 0 5
## 51 2 4 4.5R2 B4 4 1
## 52 2 4 4.7R2 B1 0 5
## 53 2 4 4.9R2 B1 0 4
## 54 2 5 5.11R2 B1 0 5
## 55 2 5 5.12R2 B1 0 5
## 56 2 5 5.1R2 B1 0 5
## 57 2 5 5.2R2 B5 1 4
## 58 2 5 5.3R2 B3 5 0
## 59 2 5 5.4R2 B5 0 5
## 60 2 5 5.5R2 B4 0 5
## 61 2 5 5.7R2 B1 0 5
## 62 2 5 5.9R2 B1 0 5
mortality <- (mortality %>%
rename("dead" = "FALSE", "alive" = "TRUE"))
ggplot(mortality, aes(x=dead, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 1 ,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Number of Workers that Died") +
labs(y = "Count", x = "Dead Workers")
ggplot(mortality, aes(x=alive, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 1 ,col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Number of Workers that Lived") +
labs(y = "Count", x = "Surviving Workers")
ggplot(mortality) +
geom_boxplot(aes(x = treatment, y = dead, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Worker Mortality per Treatment")
fitm1<-fitdistr(mortality$dead, "normal")
fitm2<-fitdistr(mortality$dead, "Poisson")
fitm4<-fitdistr(mortality$dead, "negative binomial")
#Use AIC() to compare the fit to each distribution
AIC(fitm1,fitm2,fitm4)
## df AIC
## fitm1 2 214.6261
## fitm2 1 163.9519
## fitm4 2 130.2186
As with most of the other data, treatment 3 seems like it has something going on but when we follow through with the Anova we do not see significant effects.
mort.mod <- glmer.nb(dead ~ treatment + (1|round/colony) + (1|q_origin), data = mortality)
summary(mort.mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(205.4241) ( log )
## Formula: dead ~ treatment + (1 | round/colony) + (1 | q_origin)
## Data: mortality
##
## AIC BIC logLik deviance df.resid
## 137.3 156.4 -59.6 119.3 53
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6370 -0.3888 -0.3121 -0.1358 0.8380
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 1.934e+00 1.391e+00
## q_origin (Intercept) 7.602e-01 8.719e-01
## round (Intercept) 2.074e-12 1.440e-06
## Number of obs: 62, groups: colony:round, 60; q_origin, 7; round, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6009 0.7491 -0.802 0.4225
## treatment2 -1.2566 0.9274 -1.355 0.1754
## treatment3 -1.9606 1.0340 -1.896 0.0579 .
## treatment4 -1.0063 0.8966 -1.122 0.2617
## treatment5 -0.7588 0.8721 -0.870 0.3842
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.431
## treatment3 -0.396 0.380
## treatment4 -0.468 0.418 0.381
## treatment5 -0.513 0.417 0.376 0.423
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(mort.mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dead
## Chisq Df Pr(>Chisq)
## treatment 4.1282 4 0.3889
mort.emm <- emmeans(mort.mod, "treatment")
pairs(mort.emm)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 1.257 0.927 Inf 1.355 0.6565
## treatment1 - treatment3 1.961 1.034 Inf 1.896 0.3193
## treatment1 - treatment4 1.006 0.897 Inf 1.122 0.7947
## treatment1 - treatment5 0.759 0.872 Inf 0.870 0.9079
## treatment2 - treatment3 0.704 1.095 Inf 0.643 0.9680
## treatment2 - treatment4 -0.250 0.985 Inf -0.254 0.9991
## treatment2 - treatment5 -0.498 0.972 Inf -0.512 0.9862
## treatment3 - treatment4 -0.954 1.080 Inf -0.884 0.9030
## treatment3 - treatment5 -1.202 1.073 Inf -1.120 0.7961
## treatment4 - treatment5 -0.248 0.951 Inf -0.260 0.9990
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
mort.means <- mortality %>%
group_by(treatment) %>%
summarize(dead.m = mean(dead),
dead.sd = sd(dead),
dead.n = length(dead)) %>%
mutate(dead.se = dead.sd/sqrt(dead.n))
mort.means
## # A tibble: 5 × 5
## treatment dead.m dead.sd dead.n dead.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 1.17 1.80 12 0.520
## 2 2 0.538 1.39 13 0.386
## 3 3 0.231 0.599 13 0.166
## 4 4 0.583 1.24 12 0.358
## 5 5 0.667 1.44 12 0.414
ggplot(data = mort.means, aes(x = treatment, y = dead.m, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(0,2)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = dead.m - dead.se ,
ymax = dead.m + dead.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Dead Workers",) +
ggtitle("Average Count of Worker Mortality per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 12)
w.na <- na.omit(workers)
shapiro.test(w.na$dry_weight)
##
## Shapiro-Wilk normality test
##
## data: w.na$dry_weight
## W = 0.93702, p-value = 2.384e-09
ggplot(w.na) +
geom_boxplot(aes(x = treatment, y = dry_weight, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Worker Dry Weights(g)")
ggplot(w.na) +
geom_histogram(aes(x = dry_weight, fill = treatment), binwidth = 0.01, col=I("black")) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Worker Dry Weights(g)")
dry.mod <- glmer(dry_weight ~ treatment + (1|round/colony) + (1|q_origin), data = w.na)
summary(dry.mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: dry_weight ~ treatment + (1 | round/colony) + (1 | q_origin)
## Data: w.na
##
## REML criterion at convergence: -1385.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8403 -0.6386 -0.1503 0.5603 3.2674
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 4.575e-05 6.764e-03
## q_origin (Intercept) 5.065e-05 7.117e-03
## round (Intercept) 1.593e-12 1.262e-06
## Residual 2.559e-04 1.600e-02
## Number of obs: 271, groups: colony:round, 55; q_origin, 7; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.050630 0.004231 11.967
## treatment2 -0.003292 0.004210 -0.782
## treatment3 0.002929 0.004234 0.692
## treatment4 0.001112 0.004253 0.261
## treatment5 0.002642 0.004262 0.620
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.505
## treatment3 -0.524 0.501
## treatment4 -0.522 0.498 0.508
## treatment5 -0.520 0.499 0.500 0.495
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(dry.mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 2.8313 4 0.5864
wd.m <- w.na %>%
group_by(treatment) %>%
summarize(wd.mean = mean(dry_weight),
wd.sd = sd(dry_weight),
wd.n = length(dry_weight)) %>%
mutate(wd.se = wd.sd/sqrt(wd.n))
as.data.frame(wd.m)
## treatment wd.mean wd.sd wd.n wd.se
## 1 1 0.04929944 0.01938415 54 0.002637848
## 2 2 0.04597327 0.01490736 55 0.002010108
## 3 3 0.05158964 0.01925182 55 0.002595915
## 4 4 0.04978296 0.02106748 54 0.002866920
## 5 5 0.05173909 0.01747866 53 0.002400879
wd.m
## # A tibble: 5 × 5
## treatment wd.mean wd.sd wd.n wd.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.0493 0.0194 54 0.00264
## 2 2 0.0460 0.0149 55 0.00201
## 3 3 0.0516 0.0193 55 0.00260
## 4 4 0.0498 0.0211 54 0.00287
## 5 5 0.0517 0.0175 53 0.00240
ggplot(data = wd.m, aes(x = treatment, y = wd.mean, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0.04, 0.055)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = wd.mean - wd.se,
ymax = wd.mean + wd.se),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Dry Weight(g)",) +
ggtitle("Average Worker Dry Weight(g) per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 12)
Interestingly, we see the opposite of what we predicted. Treatment 3 significantly impacted worker mortality, while treatment had overall no effect on worker dry weights.
I measured pollen consumption by weighing pollen balls fed to my bees before being placed into the microcolonies, and after 48 hours. I accounted for water loss by placing evaporative control pollen balls inside the blocks containing the bees, and accounting for average water loss per type of evaporative control pollen ball in the difference between the before and after weights of the experimental pollen balls.
Let’s look at the shape of the pollen data in a histogram.
The data doesn’t seem too abdormal. Definitely a tail though.
A shapiro-wilk test shows that it is definitely not normal.
shapiro.test(pollen$difference)
##
## Shapiro-Wilk normality test
##
## data: pollen$difference
## W = 0.87913, p-value < 2.2e-16
I tried transforming the data using the square-root and log functions, however this does not really help to normalize anything. Instead I will use a nonparametric test to look at the data.
pollen$dif.sq <- sqrt((pollen$difference)+1)
shapiro.test(pollen$dif.sq)
##
## Shapiro-Wilk normality test
##
## data: pollen$dif.sq
## W = 0.9028, p-value < 2.2e-16
pollen$dif.lg <- log((pollen$difference)+1)
shapiro.test(pollen$dif.lg)
##
## Shapiro-Wilk normality test
##
## data: pollen$dif.lg
## W = 0.92359, p-value < 2.2e-16
Let’s do a non-parametric analysis using a Kruskal Wallis Test
kruskal.test(difference ~ treatment, data = pollen)
##
## Kruskal-Wallis rank sum test
##
## data: difference by treatment
## Kruskal-Wallis chi-squared = 27.622, df = 4, p-value = 1.488e-05
The Kruskal-Wallis test indicates that pollen consumption does differ significantly with treatment level.
Let’s start digging deeper into what the data looks like.
ggplot(pollen) +
geom_boxplot(aes(x = treatment, y = difference, fill = treatment)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
theme_cowplot() +
ggtitle("Pollen Consumed per Treatment")
bar.p <- ggplot(pollen, aes(x=treatment, y = difference, color = treatment)) +
scale_color_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5"))+
geom_col()+
ggtitle("Pollen Consumed (g) per Treatment")
ggplotly(bar.p)
I ran a mixed effects model with the transformed data. This model holds round as a random effect with colony nested inside. I also have the number of bees alive in the colony as a random effect to account for consumption difference due to mortality over the experiment.
pol.mod <- glmer(difference ~ treatment + (1|round/colony) + (1|bees_alive), data = pollen, family = gaussian(link = "identity"))
summary(pol.mod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: difference ~ treatment + (1 | round/colony) + (1 | bees_alive)
## Data: pollen
##
## REML criterion at convergence: 602.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5547 -0.5978 -0.1411 0.4841 3.6746
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.032264 0.17962
## bees_alive (Intercept) 0.004497 0.06706
## round (Intercept) 0.001797 0.04239
## Residual 0.083873 0.28961
## Number of obs: 1258, groups: colony:round, 60; bees_alive, 5; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.48961 0.07514 6.516
## treatment2 0.03595 0.07784 0.462
## treatment3 0.10014 0.07779 1.287
## treatment4 0.03208 0.07781 0.412
## treatment5 -0.01540 0.07798 -0.197
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.507
## treatment3 -0.506 0.500
## treatment4 -0.523 0.498 0.498
## treatment5 -0.510 0.498 0.498 0.498
Anova(pol.mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: difference
## Chisq Df Pr(>Chisq)
## treatment 2.6066 4 0.6257
Next, I am going to summarize the data by microcolony in order to avoid pseudoreplication issues
pollen$round <- as.double(pollen$round)
pol.col <- pollen %>%
group_by(colony, treatment) %>%
summarize(dif.pc = mean(difference),
bee.pc = mean(bees_alive),
round.pc = mean(round))
pol.col <- as.data.frame(pol.col)
pol.col$round.pc <- as.factor(pol.col$round.pc)
shapiro.test(pol.col$dif.pc) #it's not NOT normal
##
## Shapiro-Wilk normality test
##
## data: pol.col$dif.pc
## W = 0.96288, p-value = 0.06512
pol.col$log.d <- log(pol.col$dif.pc)
shapiro.test(pol.col$log.d) ##better I guess???
##
## Shapiro-Wilk normality test
##
## data: pol.col$log.d
## W = 0.97004, p-value = 0.1466
Model grouped data
pol.mod3 <- lmer(dif.pc ~ treatment + (1|bee.pc), data = pol.col)
summary(pol.mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: dif.pc ~ treatment + (1 | bee.pc)
## Data: pol.col
##
## REML criterion at convergence: -13.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4167 -0.7896 -0.1951 0.7855 1.9854
##
## Random effects:
## Groups Name Variance Std.Dev.
## bee.pc (Intercept) 0.006702 0.08187
## Residual 0.033992 0.18437
## Number of obs: 60, groups: bee.pc, 17
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.44191 0.06414 6.889
## treatment2 0.02381 0.07723 0.308
## treatment3 0.08474 0.07721 1.098
## treatment4 0.03487 0.07723 0.452
## treatment5 -0.01390 0.07742 -0.180
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.580
## treatment3 -0.557 0.507
## treatment4 -0.580 0.505 0.507
## treatment5 -0.604 0.501 0.501 0.501
Anova(pol.mod3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dif.pc
## Chisq Df Pr(>Chisq)
## treatment 1.9469 4 0.7455
I added a log-transformed column to the data and will re-run my models with this.
pol.mod2 <- lmer(log.d ~ treatment + (1|bee.pc), data = pol.col)
summary(pol.mod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log.d ~ treatment + (1 | bee.pc)
## Data: pol.col
##
## REML criterion at convergence: 66.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15074 -0.67473 -0.05548 0.78469 1.62163
##
## Random effects:
## Groups Name Variance Std.Dev.
## bee.pc (Intercept) 0.1170 0.3421
## Residual 0.1269 0.3563
## Number of obs: 60, groups: bee.pc, 17
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.96196 0.15031 -6.400
## treatment2 0.08967 0.15704 0.571
## treatment3 0.16357 0.15584 1.050
## treatment4 0.09426 0.15704 0.600
## treatment5 -0.01043 0.15868 -0.066
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4
## treatment2 -0.502
## treatment3 -0.477 0.519
## treatment4 -0.502 0.513 0.519
## treatment5 -0.528 0.505 0.509 0.505
pol.log.means <- emmeans(pol.mod2, "treatment")
pairs(pol.log.means)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 -0.08967 0.161 51.8 -0.558 0.9804
## treatment1 - treatment3 -0.16357 0.159 50.2 -1.029 0.8406
## treatment1 - treatment4 -0.09426 0.161 51.8 -0.587 0.9765
## treatment1 - treatment5 0.01043 0.163 53.1 0.064 1.0000
## treatment2 - treatment3 -0.07390 0.156 48.5 -0.474 0.9894
## treatment2 - treatment4 -0.00459 0.158 50.2 -0.029 1.0000
## treatment2 - treatment5 0.10009 0.161 51.8 0.623 0.9707
## treatment3 - treatment4 0.06931 0.156 48.5 0.444 0.9917
## treatment3 - treatment5 0.17400 0.159 50.2 1.095 0.8083
## treatment4 - treatment5 0.10468 0.161 51.8 0.652 0.9655
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
Anova(pol.mod2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log.d
## Chisq Df Pr(>Chisq)
## treatment 1.7255 4 0.7861
Last I will make a bar plot of the averages based on the entire untransformed pollen data. We see that at treatment level 3, there is increased pollen consumption. However our models tells us that it is not significant.
dif.means1 <- pollen %>%
group_by(treatment) %>%
summarize(dif.mean1 = mean(difference),
dif.sd1 = sd(difference),
dif.n1 = length(difference)) %>%
mutate(dif.se1 = dif.sd1/sqrt(dif.n1))
dif.means1
## # A tibble: 5 × 5
## treatment dif.mean1 dif.sd1 dif.n1 dif.se1
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.453 0.338 256 0.0211
## 2 2 0.502 0.330 245 0.0211
## 3 3 0.552 0.354 255 0.0222
## 4 4 0.500 0.350 258 0.0218
## 5 5 0.451 0.322 244 0.0206
ggplot(data = dif.means1, aes(x = treatment, y = dif.mean1, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(0.4,0.6)) +
scale_fill_manual(values = c("pink1", "plum", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = dif.mean1 - dif.se1,
ymax = dif.mean1 + dif.se1),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Pollen Consumed (g)",) +
ggtitle("Average Pollen Consumed (g) per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 12)
Overall, bees in treatment 3 consumed the most pollen. However, based on our models we do not have statistical significance when comparing pollen consumption among groups.
weights <- read_csv("weights.csv", col_types = cols(round = col_factor(levels = c("1",
"2")), treatment = col_factor(levels = c("1",
"2", "3", "4", "5")), replicate = col_factor(levels = c("1",
"2", "3", "4", "5", "7", "9", "11", "12"))))
weights$id <- as.factor(weights$id)
weights$colony <- as.factor(weights$colony)
ggplot(weights, aes(x=diff, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.5, col=I("black")) +
scale_fill_manual(values = c("pink1", "mediumpurple1", "peachpuff4", "pink3", "peachpuff2"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Weight Change Since First Colony Weigh-In (g)") +
labs(y = "Count", x = "Weight Change (g)")
shapiro.test(weights$diff)
##
## Shapiro-Wilk normality test
##
## data: weights$diff
## W = 0.92118, p-value = 4.787e-15
Transforming with a log diesn’t help so let’s use non-parametric.
weights$logw <- log((weights$weight)+1)
shapiro.test(weights$logw)
##
## Shapiro-Wilk normality test
##
## data: weights$logw
## W = 0.97828, p-value = 1.632e-06
The Kruskal-Wallis test does not indicate that weight change from the intial colony weigh in varries significantly with treatment.
kruskal.test(diff ~ treatment, data = weights)
##
## Kruskal-Wallis rank sum test
##
## data: diff by treatment
## Kruskal-Wallis chi-squared = 2.2337, df = 4, p-value = 0.6929
waov <- glmer(diff ~ treatment * days + (1|round/colony) + (1|days), data = weights)
summary(waov)
## Linear mixed model fit by REML ['lmerMod']
## Formula: diff ~ treatment * days + (1 | round/colony) + (1 | days)
## Data: weights
##
## REML criterion at convergence: 2202
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.94008 -0.57584 -0.06979 0.49077 3.05170
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 3.320 1.822
## days (Intercept) 1.291 1.136
## round (Intercept) 2.142 1.464
## Residual 4.012 2.003
## Number of obs: 474, groups: colony:round, 60; days, 57; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.03735 1.27463 1.598
## treatment2 -0.58161 0.90223 -0.645
## treatment3 -0.63150 0.89973 -0.702
## treatment4 -0.74002 0.90229 -0.820
## treatment5 -0.38846 0.90721 -0.428
## days 0.18091 0.01702 10.629
## treatment2:days 0.04780 0.01997 2.394
## treatment3:days 0.06525 0.01894 3.445
## treatment4:days 0.05417 0.01884 2.875
## treatment5:days 0.01187 0.01994 0.596
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 days trtm2: trtm3: trtm4:
## treatment2 -0.349
## treatment3 -0.353 0.496
## treatment4 -0.349 0.498 0.500
## treatment5 -0.346 0.495 0.497 0.500
## days -0.325 0.235 0.248 0.241 0.242
## trtmnt2:dys 0.143 -0.456 -0.213 -0.210 -0.213 -0.491
## trtmnt3:dys 0.156 -0.222 -0.455 -0.222 -0.223 -0.550 0.474
## trtmnt4:dys 0.151 -0.223 -0.224 -0.448 -0.227 -0.550 0.471 0.497
## trtmnt5:dys 0.144 -0.211 -0.212 -0.211 -0.456 -0.525 0.453 0.475 0.481
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00932703 (tol = 0.002, component 1)
Anova(waov)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: diff
## Chisq Df Pr(>Chisq)
## treatment 1.6011 4 0.808591
## days 301.7570 1 < 2.2e-16 ***
## treatment:days 17.3062 4 0.001685 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(waov)
What I see in this graph is that in treatment 4 there is one colony that seems to have remained stagnant over time. I will try removing this colony and see what happens.
wplot <- ggplot(weights) +
aes( x = days, y = diff, color = treatment) +
geom_point(aes(color = colony)) +
geom_smooth() +
facet_wrap(vars(treatment)) +
theme(legend.position = "none") +
ylab("Time (days since colony initiation") +
xlab("Weight Change(g)")+
ggtitle("Colony Weight Change Over Time From Initial Weigh-In (g)")
ggplotly(wplot,
tooltip = c("colour"))
w4 <- weights[weights$colony != "T4.5R2", ]
With the removal of this colony I don’t see that much of a difference in the diagnostic plots.
w4a <- glmer(logw ~ treatment * days + (1|round/colony) + (1|number), data = w4)
summary(w4a)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logw ~ treatment * days + (1 | round/colony) + (1 | number)
## Data: w4
##
## REML criterion at convergence: -1543.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8087 -0.6865 -0.0305 0.5402 3.1847
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.0011854 0.03443
## number (Intercept) 0.0002310 0.01520
## round (Intercept) 0.0005326 0.02308
## Residual 0.0012379 0.03518
## Number of obs: 464, groups: colony:round, 59; number, 11; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.9053392 0.0220113 177.424
## treatment2 -0.0083810 0.0166230 -0.504
## treatment3 -0.0140184 0.0165679 -0.846
## treatment4 -0.0220954 0.0169332 -1.305
## treatment5 0.0025536 0.0166275 0.154
## days 0.0033914 0.0003366 10.076
## treatment2:days 0.0007077 0.0003431 2.063
## treatment3:days 0.0009998 0.0003239 3.087
## treatment4:days 0.0015268 0.0003328 4.587
## treatment5:days 0.0002007 0.0003368 0.596
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 days trtm2: trtm3: trtm4:
## treatment2 -0.368
## treatment3 -0.370 0.497
## treatment4 -0.364 0.485 0.486
## treatment5 -0.367 0.495 0.496 0.485
## days -0.394 0.182 0.186 0.187 0.182
## trtmnt2:dys 0.132 -0.434 -0.201 -0.195 -0.200 -0.383
## trtmnt3:dys 0.146 -0.212 -0.432 -0.205 -0.210 -0.428 0.472
## trtmnt4:dys 0.141 -0.204 -0.204 -0.429 -0.203 -0.429 0.454 0.476
## trtmnt5:dys 0.144 -0.202 -0.202 -0.196 -0.431 -0.423 0.449 0.472 0.454
Anova(w4a)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logw
## Chisq Df Pr(>Chisq)
## treatment 0.589 4 0.9643
## days 200.796 1 < 2.2e-16 ***
## treatment:days 26.605 4 2.389e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(w4a)
However, the AIC analysis shows that the removal of this colony is a better fit for the data overall.
AIC(waov, w4a)
## df AIC
## waov 14 2229.952
## w4a 14 -1515.757
w4plot <- ggplot(w4) +
aes( x = days, y = diff, color = treatment) +
geom_point(aes(color = colony)) +
geom_smooth() +
facet_wrap(vars(treatment)) +
theme(legend.position = "none") +
ylab("Time (days since colony initiation") +
xlab("Weight Change(g)")+
ggtitle("Colony Weight Change Over Time From Initial Weigh-In (g)")
ggplotly(w4plot,
tooltip = c("colour"))
The last time we have data for all the treatments is day 42 - I will truncate from here.
ggplot(w4) +
aes( x = days, y = diff, color = treatment) +
geom_point(aes(color = colony)) +
geom_smooth() +
theme(legend.position = "none") +
ylab("Weight Change (g)") +
scale_x_continuous(name ="Time (days)",
breaks = round(seq(min(w4$days), max(w4$days), by = 2)))+
ggtitle("Colony Weight Change Over Time From Initial Weigh-In (g)")
table(weights['number'])
## number
## 1 2 3 4 5 6 7 8 9 10 11
## 60 60 60 60 59 59 57 40 12 5 2
wt <- subset(weights,days < 42)
ggplot(wt) +
aes( x = days, y = diff, color = treatment) +
geom_smooth() +
geom_point(aes(color = colony)) +
theme(legend.position = "none") +
ylab(" Weight (g)") +
scale_x_continuous(name ="Time (days)", limits = c(0,42),
breaks = round(seq(min(w4$days), max(w4$days), by = 7)))+
ggtitle("Colony Weight Change Over Time From Initial Weigh-In (g)") +
facet_wrap(vars(treatment))
Let’s see how this impacts our model.
I’m still not overjoyed with the fit, but the residuals are much more evenly distributed in the residuals vs. fitted plot.
wta <- glmer(logw ~ treatment*days + (1|round/colony) + (1|days), data = wt)
summary(wta)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logw ~ treatment * days + (1 | round/colony) + (1 | days)
## Data: wt
##
## REML criterion at convergence: -1424.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.68447 -0.60921 0.00349 0.52664 2.82571
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony:round (Intercept) 0.0008741 0.02957
## days (Intercept) 0.0002226 0.01492
## round (Intercept) 0.0009030 0.03005
## Residual 0.0009450 0.03074
## Number of obs: 408, groups: colony:round, 60; days, 42; round, 2
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.888e+00 2.439e-02 159.447
## treatment2 -1.421e-03 1.461e-02 -0.097
## treatment3 -9.496e-03 1.460e-02 -0.650
## treatment4 -8.980e-03 1.468e-02 -0.612
## treatment5 3.712e-03 1.470e-02 0.252
## days 4.586e-03 3.421e-04 13.406
## treatment2:days 1.724e-04 3.772e-04 0.457
## treatment3:days 6.801e-04 3.675e-04 1.851
## treatment4:days 5.546e-04 3.745e-04 1.481
## treatment5:days 8.271e-05 3.775e-04 0.219
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 days trtm2: trtm3: trtm4:
## treatment2 -0.301
## treatment3 -0.302 0.500
## treatment4 -0.301 0.502 0.502
## treatment5 -0.297 0.499 0.499 0.503
## days -0.268 0.249 0.252 0.249 0.245
## trtmnt2:dys 0.137 -0.452 -0.223 -0.228 -0.226 -0.544
## trtmnt3:dys 0.142 -0.228 -0.453 -0.232 -0.230 -0.568 0.500
## trtmnt4:dys 0.138 -0.229 -0.229 -0.452 -0.232 -0.554 0.504 0.518
## trtmnt5:dys 0.134 -0.226 -0.225 -0.229 -0.456 -0.544 0.500 0.507 0.516
Anova(wta)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logw
## Chisq Df Pr(>Chisq)
## treatment 0.1806 4 0.9962
## days 397.8616 1 <2e-16 ***
## treatment:days 5.3180 4 0.2562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(wta)
Treatment 3 produced on average the highest numbers of overall brood cells, total larvae, and total pupae. Treatment 3 also had the highest average counts of overall dead brood cells, with a P Value of 0.08. Treatment 3 had the lowest average drone weight (P = 0.04). The lowest worker mortality was seen in treatment 3 ( P = 0.39). Treatment 3 also had the highest pollen consumption (P=0.78). Treatment 4 produced the most drones overall, and had the highest overall average drone production. Treatment 4 also took the longest time for average drone emergence (P = 0.02). Colony weight change did not appear to vary by treatment.
I think what may be going on is this: Workers in treatment 3 eat the post pollen. They are ingesting higher quantities at once of the fungicide, even though it is the intermediate dose for this experiment. However, they are also feeding higher quantities of the pollen to their brood cells. This might account for the increased brood production. However, the large doses of fungicide being given to the brood might explain the increased brood mortality, and the decreased drone size in the brood that does survive.