Start by imputing pollen data and creating a new data frame with average pollen consumption on a per-colony basis
### Figure out average pollen consumption by treatment
pollen <- read_csv("pollen1.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", "6", "7", "9", "11",
"12")), start_date = col_date(format = "%m/%d/%Y"),
start_time = col_character(), end_date = col_date(format = "%m/%d/%Y"),
end_time = col_character()))
pollen$colony <- as.factor(pollen$colony)
pollen$pid <- as.factor(pollen$pid)
pollen$count <- as.factor(pollen$count)
pollen <- subset(pollen, pollen$round != 1)
pollen <- na.omit(pollen)
range(pollen$difference)
## [1] -0.30646 1.56542
# get rid of negative numbers
pollen$difference[pollen$difference < 0] <- NA
pollen <- na.omit(pollen)
range(pollen$difference)
## [1] 0.002715 1.565420
# add queenright original colony column
qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)
pollen <- merge(pollen, qro, by.x = "colony")
pollen$pid <- as.factor(pollen$pid)
pollen$whole_dif <- as.numeric(pollen$difference)
wd.1 <- lm(difference ~ treatment, data = pollen)
summary(wd.1)
##
## Call:
## lm(formula = difference ~ treatment, data = pollen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4853 -0.2416 -0.1504 0.1576 1.0638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.429541 0.024170 17.772 <2e-16 ***
## treatment2 0.072099 0.034886 2.067 0.0390 *
## treatment3 0.078078 0.034406 2.269 0.0235 *
## treatment4 0.058490 0.034988 1.672 0.0949 .
## treatment5 0.004982 0.035040 0.142 0.8870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3375 on 915 degrees of freedom
## Multiple R-squared: 0.009927, Adjusted R-squared: 0.005599
## F-statistic: 2.294 on 4 and 915 DF, p-value: 0.05773
wd.emm <- emmeans(wd.1, "treatment")
summary(wd.emm)
## treatment emmean SE df lower.CL upper.CL
## 1 0.430 0.0242 915 0.382 0.477
## 2 0.502 0.0252 915 0.452 0.551
## 3 0.508 0.0245 915 0.460 0.556
## 4 0.488 0.0253 915 0.438 0.538
## 5 0.435 0.0254 915 0.385 0.484
##
## Confidence level used: 0.95
wd.summary <- pollen %>%
group_by(colony) %>%
summarize(whole.mean = mean(difference))
as.data.frame(wd.summary) # this is the data frame I will merge with subsequent data frames to contain average pollen consumption per colony as a new column
## colony whole.mean
## 1 1.11R2 0.2829509
## 2 1.12R2 0.1697964
## 3 1.1R2 0.5213458
## 4 1.2R2 0.3374200
## 5 1.3R2 0.4512891
## 6 1.4R2 0.6063016
## 7 1.5R2 0.7079516
## 8 1.7R2 0.7400381
## 9 1.9R2 0.2240081
## 10 2.11R2 0.4178270
## 11 2.12R2 0.4035568
## 12 2.1R2 0.6101589
## 13 2.2R2 0.5109234
## 14 2.3R2 0.3280036
## 15 2.4R2 0.3881394
## 16 2.5R2 0.7194915
## 17 2.7R2 0.5299685
## 18 2.9R2 0.5857152
## 19 3.11R2 0.4216410
## 20 3.12R2 0.3390993
## 21 3.1R2 0.3711948
## 22 3.2R2 0.3461010
## 23 3.3R2 0.8465806
## 24 3.4R2 0.4120433
## 25 3.5R2 0.8906211
## 26 3.7R2 0.6266680
## 27 3.9R2 0.4409511
## 28 4.11R2 0.3456011
## 29 4.12R2 0.2496295
## 30 4.1R2 0.7074755
## 31 4.2R2 0.3871742
## 32 4.3R2 0.5800074
## 33 4.4R2 0.8356247
## 34 4.5R2 0.2335967
## 35 4.7R2 0.6237400
## 36 4.9R2 0.5322950
## 37 5.11R2 0.4113960
## 38 5.12R2 0.2077741
## 39 5.1R2 0.3782286
## 40 5.2R2 0.4912468
## 41 5.3R2 0.2128778
## 42 5.4R2 0.4563045
## 43 5.5R2 0.7095479
## 44 5.7R2 0.6113189
## 45 5.9R2 0.4188290
mortality <- read_csv("surviving workers per colony.csv")
mortality$colony <- as.factor(mortality$colony)
ggplot(pollen, aes(x=difference, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.05,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Pollen Consumption (g)") +
labs(y = "Number of Pollen Balls", x = "Pollen Consumed (g)") +
facet_wrap(vars(treatment))
shapiro.test(pollen$difference)
##
## Shapiro-Wilk normality test
##
## data: pollen$difference
## W = 0.84265, p-value < 2.2e-16
pollen$boxp <- bcPower(pollen$difference, -3, gamma=1)
shapiro.test(pollen$boxp)
##
## Shapiro-Wilk normality test
##
## data: pollen$boxp
## W = 0.9588, p-value = 2.044e-15
ggplot(pollen, aes(x=boxp, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.009,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Pollen Consumption (g) - BoxCox Transformed") +
labs(y = "Number of Pollen Balls", x = "Pollen Consumed (g), BoxCox power transformation")
#### Pollen Model
pbox <- glm(boxp ~ treatment + bees_alive + qro + start_date, data = pollen)
summary(pbox)
##
## Call:
## glm(formula = boxp ~ treatment + bees_alive + qro + start_date,
## data = pollen)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.201511 -0.038005 -0.006821 0.044469 0.124757
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.806e+01 2.841e+00 -9.878 < 2e-16 ***
## treatment2 1.655e-02 6.017e-03 2.750 0.006084 **
## treatment3 1.610e-02 5.941e-03 2.710 0.006865 **
## treatment4 2.167e-02 6.064e-03 3.573 0.000371 ***
## treatment5 -2.725e-03 6.077e-03 -0.448 0.653937
## bees_alive 2.848e-02 3.123e-03 9.121 < 2e-16 ***
## qroB3 2.024e-02 6.818e-03 2.969 0.003065 **
## qroB4 5.463e-02 6.767e-03 8.073 2.17e-15 ***
## qroB5 1.531e-02 4.852e-03 3.155 0.001655 **
## start_date 1.461e-03 1.474e-04 9.910 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.003361313)
##
## Null deviance: 3.6750 on 919 degrees of freedom
## Residual deviance: 3.0588 on 910 degrees of freedom
## AIC: -2617
##
## Number of Fisher Scoring iterations: 2
Anova(pbox)
## Analysis of Deviance Table (Type II tests)
##
## Response: boxp
## LR Chisq Df Pr(>Chisq)
## treatment 25.542 4 3.914e-05 ***
## bees_alive 83.186 1 < 2.2e-16 ***
## qro 68.697 3 8.113e-15 ***
## start_date 98.200 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(pbox, test = "F")
## Single term deletions
##
## Model:
## boxp ~ treatment + bees_alive + qro + start_date
## Df Deviance AIC F value Pr(>F)
## <none> 3.0588 -2617.0
## treatment 4 3.1446 -2599.5 6.3855 4.548e-05 ***
## bees_alive 1 3.3384 -2538.5 83.1862 < 2.2e-16 ***
## qro 3 3.2897 -2556.0 22.8991 2.691e-14 ***
## start_date 1 3.3889 -2524.7 98.1998 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(pbox)
boxemm <- emmeans(pbox, "treatment")
pairs(boxemm)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 -0.016546 0.00602 910 -2.750 0.0478
## treatment1 - treatment3 -0.016098 0.00594 910 -2.710 0.0533
## treatment1 - treatment4 -0.021668 0.00606 910 -3.573 0.0034
## treatment1 - treatment5 0.002725 0.00608 910 0.448 0.9916
## treatment2 - treatment3 0.000448 0.00603 910 0.074 1.0000
## treatment2 - treatment4 -0.005123 0.00625 910 -0.819 0.9247
## treatment2 - treatment5 0.019271 0.00615 910 3.132 0.0154
## treatment3 - treatment4 -0.005570 0.00620 910 -0.899 0.8974
## treatment3 - treatment5 0.018823 0.00607 910 3.103 0.0169
## treatment4 - treatment5 0.024394 0.00634 910 3.849 0.0012
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
# Assumption 1 - No Outliers
plot(pollen$treatment, pollen$difference)
summary(pollen$difference)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002715 0.234403 0.323385 0.472047 0.632838 1.565420
boxplot(pollen$difference)
out <- boxplot.stats(pollen$difference)$out
out_ind <- which(pollen$difference %in% c(out))
out_ind
## [1] 81 122 148 162 163 164 247 264 268 284 326 328 329 330 367 369 370 387 439
## [20] 476 483 514 515 516 517 518 521 613 653 669 673 678 720 875 891
pollen.out <- pollen[-c(81,122,148,162,163,164,247,264,268,284,
326,328,329,330,367,369,370,387,439,476,
483,514,515,516,517,518,521,613,29,653,
669,673,678,720,875,891),]
boxplot(pollen.out$difference)
library(EnvStats)
test <- rosnerTest(pollen$difference,
k = 35
)
test
## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7 R.8
## 3.230400 3.227370 3.175798 3.017167 3.030928 3.036798 2.946971 2.930484
## R.9 R.10 R.11 R.12 R.13 R.14 R.15 R.16
## 2.897848 2.848433 2.818843 2.802060 2.787905 2.794322 2.800783 2.798089
## R.17 R.18 R.19 R.20 R.21 R.22 R.23 R.24
## 2.794140 2.802657 2.802424 2.789264 2.802324 2.816092 2.794116 2.749785
## R.25 R.26 R.27 R.28 R.29 R.30 R.31 R.32
## 2.738855 2.750513 2.688535 2.694709 2.698672 2.685078 2.689350 2.697842
## R.33 R.34 R.35
## 2.708858 2.694204 2.702961
##
## $sample.size
## [1] 920
##
## $parameters
## k
## 35
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8
## 4.019343 4.019073 4.018802 4.018531 4.018260 4.017989 4.017717 4.017445
## lambda.9 lambda.10 lambda.11 lambda.12 lambda.13 lambda.14 lambda.15 lambda.16
## 4.017172 4.016899 4.016626 4.016353 4.016079 4.015805 4.015531 4.015256
## lambda.17 lambda.18 lambda.19 lambda.20 lambda.21 lambda.22 lambda.23 lambda.24
## 4.014981 4.014705 4.014430 4.014153 4.013877 4.013600 4.013323 4.013046
## lambda.25 lambda.26 lambda.27 lambda.28 lambda.29 lambda.30 lambda.31 lambda.32
## 4.012768 4.012490 4.012211 4.011932 4.011653 4.011374 4.011094 4.010814
## lambda.33 lambda.34 lambda.35
## 4.010533 4.010252 4.009971
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 35 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 0.205450 0.211640 0.351300 0.180440 0.150420 0.166610 0.214700 0.239390
## [9] 0.235980 0.278630 0.435890 0.422590 0.408180 0.350460 0.445960 0.603040
## [17] 0.333540 0.278300 0.182380 0.175990 0.165760 0.188270 0.232330 0.201980
## [25] 0.108300 0.175750 0.139230 0.164590 0.210730 0.180380 0.160270 0.167570
## [33] 0.200150 0.154060 0.177430 0.167300 0.161910 0.148860 0.156700 0.158880
## [41] 0.164800 0.175660 0.175490 0.153150 0.288840 0.229220 0.260440 0.112780
## [49] 0.092160 0.229520 0.261760 0.276360 0.432240 1.065710 0.397080 0.585310
## [57] 0.660780 0.698440 0.891430 0.838130 0.760310 0.977340 0.847720 0.213130
## [65] 0.299240 0.053210 0.107150 0.223970 0.203480 0.139000 0.185460 0.191400
## [73] 0.380010 0.285180 0.212490 0.222130 0.293170 0.359630 0.378720 0.565590
## [81] 1.247980 1.071640 0.353490 0.222090 0.274360 0.314590 0.354880 0.399050
## [89] 0.308830 0.250470 0.295540 0.274880 0.287600 0.197020 0.234150 0.274950
## [97] 0.303640 0.253240 0.259870 0.267040 0.910440 0.826240 0.961760 0.787230
## [105] 0.610410 0.709510 0.583580 0.534920 0.371600 0.368680 0.383300 0.232760
## [113] 0.367810 0.057410 0.291320 0.122730 0.074550 0.312200 0.441260 0.655330
## [121] 0.949270 1.372910 1.047990 1.105490 0.728460 0.945370 1.091710 0.841690
## [129] 0.437940 0.388310 0.287980 0.304980 0.338530 0.257850 0.045500 0.271410
## [137] 0.229260 0.308320 0.497180 1.082170 1.157790 1.131620 1.201010 0.744970
## [145] 0.537090 1.170700 1.088990 1.323870 1.043250 0.716590 0.305630 0.250800
## [153] 0.294200 0.210120 0.154150 0.305480 0.411660 0.461440 0.880260 1.090730
## [161] 0.831600 1.238180 1.423773 1.354040 0.826110 1.088440 1.061540 1.173070
## [169] 0.699500 0.259480 0.245510 0.190650 0.103610 0.137490 0.362010 0.213860
## [177] 0.222590 0.185160 0.256280 0.298270 0.240040 0.256830 0.231030 0.267870
## [185] 0.224050 0.216650 0.221410 0.202960 0.230960 0.221710 0.224240 0.198830
## [193] 0.197720 0.207950 0.207050 0.272640 0.263370 0.453940 0.306050 0.214870
## [201] 0.280830 0.381160 0.365720 0.512730 1.123150 1.136430 0.710010 0.467580
## [209] 0.285030 0.221730 0.278700 0.260190 0.274250 0.249860 0.298300 0.111640
## [217] 0.254170 0.316830 0.317650 0.291950 0.278650 0.342470 0.423300 0.636120
## [225] 0.926040 1.145690 0.717420 0.465180 0.268170 0.243600 0.250620 0.258720
## [233] 0.261980 0.256960 0.197700 0.509833 0.239410 0.203620 0.224010 0.305390
## [241] 0.226260 0.312760 0.347800 0.502350 0.632490 0.754300 1.324110 1.086720
## [249] 0.912790 1.031540 0.829260 0.692310 0.687940 0.785480 0.494480 0.098640
## [257] 0.252770 0.137930 0.084460 0.316020 0.247800 0.243700 0.276270 1.332559
## [265] 0.400870 0.580150 1.111260 1.565420 0.997240 0.553460 0.465500 0.402350
## [273] 0.403320 0.409060 0.574930 0.449190 0.440060 0.408280 0.300530 0.218920
## [281] 0.279550 0.298480 0.284660 1.408069 0.294880 0.305410 0.226270 0.237310
## [289] 0.258960 0.266810 0.238150 0.243230 0.239900 0.281490 0.257350 0.333100
## [297] 0.259000 0.234090 0.188650 0.238560 0.176770 0.247480 0.129210 0.540830
## [305] 0.439450 0.568420 1.026180 1.173030 0.439070 0.268420 0.234420 0.220480
## [313] 0.273220 0.200090 0.227680 0.315420 0.238270 0.414200 0.363290 0.456840
## [321] 0.423780 0.504500 0.670940 0.946830 0.887000 1.264060 1.007000 1.344220
## [329] 1.473180 1.387130 0.776940 0.814340 0.481190 0.393020 0.279380 0.255770
## [337] 0.229600 0.250150 0.268840 0.292090 0.294500 0.348960 0.326320 0.464980
## [345] 0.793600 0.910840 1.085900 1.185000 1.020100 0.681440 0.559980 0.413880
## [353] 0.388730 0.549310 0.276080 0.215990 0.209900 0.063690 0.616930 0.265180
## [361] 0.306970 0.292100 0.327270 0.534010 0.916730 0.607560 1.469500 1.012030
## [369] 1.434330 1.287410 0.905680 0.633880 0.337140 0.297019 0.290620 0.254380
## [377] 0.282210 0.163290 0.234950 0.220060 0.267880 0.327400 0.357200 0.461710
## [385] 0.989830 0.512510 1.474180 0.191470 0.220580 0.233730 0.383920 0.579170
## [393] 0.721780 0.420730 0.245500 0.311980 0.073390 0.303950 0.236760 0.263430
## [401] 0.243390 0.275030 0.268970 0.283230 0.383560 0.493530 0.755010 0.548010
## [409] 0.303940 0.246150 0.355970 0.418600 0.652990 0.575680 0.487970 0.338490
## [417] 0.252200 0.229230 0.292010 0.266520 0.262030 0.234350 0.111290 0.244430
## [425] 0.237650 0.199050 0.249640 0.151400 0.208110 0.211430 0.253390 0.265140
## [433] 0.302360 0.334000 0.419870 0.445010 0.553250 0.451860 1.363000 0.511440
## [441] 0.417790 0.348600 0.300080 0.289360 0.360170 0.420450 0.343240 0.224790
## [449] 0.214030 0.299710 0.256620 0.282730 0.265840 0.238810 0.219800 0.448630
## [457] 0.779960 0.453270 0.294490 0.354930 0.445030 0.463820 0.303720 0.361640
## [465] 0.316420 0.354540 0.338690 0.336270 0.358180 0.265500 0.290200 0.565730
## [473] 0.679300 0.846040 1.156570 1.338620 1.136460 0.964570 1.190350 1.190570
## [481] 1.102180 1.003710 1.323910 1.151600 0.177340 0.037450 0.034320 0.277560
## [489] 0.163160 0.116370 0.182830 0.309320 0.332210 0.477610 0.724430 0.748500
## [497] 0.628090 0.765000 0.696030 0.383740 0.339420 0.414130 0.619180 0.772060
## [505] 0.454160 0.242920 0.325640 0.395240 0.329380 0.418310 0.479340 0.708760
## [513] 0.942510 1.247490 1.351750 1.557570 1.533510 1.286970 1.080230 1.190970
## [521] 1.312720 1.104520 0.818360 0.595610 0.252050 0.255770 0.233530 0.230070
## [529] 0.213180 0.261670 0.249040 0.319990 0.440240 0.631850 0.714310 1.168130
## [537] 1.121740 1.115360 1.064480 1.032270 1.000810 0.874490 0.614330 0.516900
## [545] 0.462880 1.013605 0.153880 0.203910 0.199850 0.252770 0.239410 0.276630
## [553] 0.303720 0.413260 0.410040 0.342290 0.436490 0.389690 0.647130 0.744390
## [561] 0.811730 0.836950 0.831780 0.545070 0.339080 0.237930 0.229170 0.331360
## [569] 0.002715 0.032775 0.354110 0.433720 0.498110 0.295700 0.411050 0.402300
## [577] 0.445070 0.374850 0.489200 0.442010 0.379090 0.349750 0.511910 0.200640
## [585] 0.263270 0.234080 0.278850 0.323450 0.285030 0.391030 0.216510 0.197010
## [593] 0.205560 0.239730 0.220710 0.266990 0.213350 0.255390 0.274910 0.230840
## [601] 0.224170 0.221440 0.250280 0.286120 0.238070 0.247790 0.189730 0.196050
## [609] 0.232570 0.428880 0.635590 1.155420 1.259550 1.069710 1.229890 1.022630
## [617] 0.940580 1.009420 1.058660 0.825600 1.048520 0.824450 0.268240 0.239240
## [625] 0.214440 0.240490 0.286090 0.312740 0.008865 0.045465 0.225395 1.156240
## [633] 1.114360 0.741010 0.957060 0.608770 0.237380 0.208750 0.315690 0.195850
## [641] 0.198330 0.169080 0.304060 0.370970 0.191920 0.337910 0.164840 0.224290
## [649] 0.489440 0.663680 0.680785 0.788305 1.294820 1.154350 0.652770 0.456300
## [657] 0.367640 0.471720 0.518040 0.745330 1.142970 0.529800 0.483320 0.505930
## [665] 0.263420 0.447040 0.321010 0.480610 1.336980 0.762315 0.972915 1.043700
## [673] 1.262150 1.016080 1.167380 0.921070 0.821430 1.349450 1.024280 1.167990
## [681] 0.184910 0.222460 0.217310 0.220810 0.219960 0.247700 0.212790 0.220460
## [689] 0.202820 0.234270 0.317510 0.186860 0.190290 0.240950 0.201700 0.208140
## [697] 0.260570 0.243430 0.261300 0.276070 0.264460 0.274580 0.212260 0.238750
## [705] 0.240190 0.260110 0.246450 0.182320 0.278870 0.217690 0.257690 0.178080
## [713] 0.282570 0.365900 0.591090 0.761150 0.978050 0.986870 1.148810 1.239330
## [721] 0.881800 0.561250 0.635050 1.044640 0.636160 0.323810 0.335760 0.241800
## [729] 0.307850 0.325140 0.339830 0.415250 0.386660 0.752270 1.167410 1.071420
## [737] 1.014370 0.323320 0.337960 0.512120 0.625650 0.491670 0.609020 0.298700
## [745] 0.309420 0.240050 0.238780 0.182260 0.272800 0.284250 0.331640 0.446800
## [753] 0.972120 0.869940 0.348410 0.269610 0.217190 0.234980 0.296180 0.403730
## [761] 0.697210 1.016370 0.297480 0.184500 0.199100 0.182390 0.163050 0.200660
## [769] 0.206330 0.202850 0.214660 0.181430 0.204970 0.191360 0.159100 0.135380
## [777] 0.143170 0.374020 0.494840 0.312920 0.224850 0.163550 0.184170 0.183490
## [785] 0.180920 0.161850 0.183650 0.202000 0.193450 0.181240 0.163750 0.210360
## [793] 0.086550 0.212340 0.055810 0.061700 0.141610 0.211220 0.208540 0.236660
## [801] 0.440090 0.671460 0.851170 0.868790 0.462170 0.504330 0.573430 0.617570
## [809] 0.601240 0.326440 0.437570 0.174280 0.233800 0.171470 0.046350 0.186520
## [817] 0.209290 0.250630 0.220900 0.254290 0.418060 0.765460 1.052550 1.092220
## [825] 0.996360 0.877900 1.000510 0.842070 0.565530 0.599940 0.309430 0.325010
## [833] 0.214860 0.278460 0.261310 0.031550 0.195320 0.043660 0.037750 0.425930
## [841] 0.433650 0.208270 0.170350 0.261390 0.264070 0.106610 0.045590 0.113350
## [849] 0.241570 0.362320 0.347940 0.841530 1.215710 0.725750 0.435190 0.439440
## [857] 0.556610 0.766590 0.920220 0.793740 0.338860 0.179260 0.263610 0.281980
## [865] 0.277340 0.224310 0.165260 0.432010 0.336050 0.624870 0.994370 1.107250
## [873] 1.148710 1.049630 1.251650 1.075830 1.228060 0.478750 0.821990 1.136120
## [881] 0.583620 0.318040 0.331790 0.255800 0.247470 0.246880 0.164140 0.597700
## [889] 0.923430 1.101790 1.249190 1.075580 1.183480 0.771790 0.458730 0.583520
## [897] 0.723820 0.548080 0.222510 0.235870 0.162130 0.186220 0.127080 0.241400
## [905] 0.284800 0.276330 0.280410 0.422650 0.536930 0.807410 0.885320 0.837720
## [913] 0.917150 0.544180 0.278840 0.305250 0.359430 0.366890 0.322890 0.416510
##
## $data.name
## [1] "pollen$difference"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 0.4720471 0.3384637 1.565420 268 3.230400 4.019343 FALSE
## 2 1 0.4708573 0.3367177 1.557570 516 3.227370 4.019073 FALSE
## 3 2 0.4696735 0.3349824 1.533510 517 3.175798 4.018802 FALSE
## 4 3 0.4685134 0.3333149 1.474180 387 3.017167 4.018531 FALSE
## 5 4 0.4674155 0.3318339 1.473180 329 3.030928 4.018260 FALSE
## 6 5 0.4663163 0.3303426 1.469500 367 3.036798 4.017989 FALSE
## 7 6 0.4652187 0.3288499 1.434330 369 2.946971 4.017717 FALSE
## 8 7 0.4641573 0.3274598 1.423773 163 2.930484 4.017445 FALSE
## 9 8 0.4631051 0.3260916 1.408069 284 2.897848 4.017172 FALSE
## 10 9 0.4620678 0.3247618 1.387130 330 2.848433 4.016899 FALSE
## 11 10 0.4610512 0.3234869 1.372910 122 2.818843 4.016626 FALSE
## 12 11 0.4600481 0.3222457 1.363000 439 2.802060 4.016353 FALSE
## 13 12 0.4590537 0.3210247 1.354040 164 2.787905 4.016079 FALSE
## 14 13 0.4580669 0.3198211 1.351750 515 2.794322 4.015805 FALSE
## 15 14 0.4570805 0.3186143 1.349450 678 2.800783 4.015531 FALSE
## 16 15 0.4560945 0.3174043 1.344220 328 2.798089 4.015256 FALSE
## 17 16 0.4551120 0.3162003 1.338620 476 2.794140 4.014981 FALSE
## 18 17 0.4541336 0.3150034 1.336980 669 2.802657 4.014705 FALSE
## 19 18 0.4531548 0.3138013 1.332559 264 2.802424 4.014430 FALSE
## 20 19 0.4521788 0.3126026 1.324110 247 2.789264 4.014153 FALSE
## 21 20 0.4512100 0.3114201 1.323910 483 2.802324 4.013877 FALSE
## 22 21 0.4502392 0.3102280 1.323870 148 2.816092 4.013600 FALSE
## 23 22 0.4492664 0.3090257 1.312720 521 2.794116 4.013323 FALSE
## 24 23 0.4483038 0.3078481 1.294820 653 2.749785 4.013046 FALSE
## 25 24 0.4473590 0.3067162 1.287410 370 2.738855 4.012768 FALSE
## 26 25 0.4464204 0.3055974 1.286970 518 2.750513 4.012490 FALSE
## 27 26 0.4454802 0.3044705 1.264060 326 2.688535 4.012211 FALSE
## 28 27 0.4445635 0.3034043 1.262150 673 2.694709 4.011932 FALSE
## 29 28 0.4436470 0.3023350 1.259550 613 2.698672 4.011653 FALSE
## 30 29 0.4427312 0.3012646 1.251650 875 2.685078 4.011374 FALSE
## 31 30 0.4418223 0.3002092 1.249190 891 2.689350 4.011094 FALSE
## 32 31 0.4409142 0.2991524 1.247980 81 2.697842 4.010814 FALSE
## 33 32 0.4400053 0.2980904 1.247490 514 2.708858 4.010533 FALSE
## 34 33 0.4390950 0.2970209 1.239330 720 2.694204 4.010252 FALSE
## 35 34 0.4381918 0.2959674 1.238180 162 2.702961 4.009971 FALSE
##
## attr(,"class")
## [1] "gofOutlier"
test$all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 0.4720471 0.3384637 1.565420 268 3.230400 4.019343 FALSE
## 2 1 0.4708573 0.3367177 1.557570 516 3.227370 4.019073 FALSE
## 3 2 0.4696735 0.3349824 1.533510 517 3.175798 4.018802 FALSE
## 4 3 0.4685134 0.3333149 1.474180 387 3.017167 4.018531 FALSE
## 5 4 0.4674155 0.3318339 1.473180 329 3.030928 4.018260 FALSE
## 6 5 0.4663163 0.3303426 1.469500 367 3.036798 4.017989 FALSE
## 7 6 0.4652187 0.3288499 1.434330 369 2.946971 4.017717 FALSE
## 8 7 0.4641573 0.3274598 1.423773 163 2.930484 4.017445 FALSE
## 9 8 0.4631051 0.3260916 1.408069 284 2.897848 4.017172 FALSE
## 10 9 0.4620678 0.3247618 1.387130 330 2.848433 4.016899 FALSE
## 11 10 0.4610512 0.3234869 1.372910 122 2.818843 4.016626 FALSE
## 12 11 0.4600481 0.3222457 1.363000 439 2.802060 4.016353 FALSE
## 13 12 0.4590537 0.3210247 1.354040 164 2.787905 4.016079 FALSE
## 14 13 0.4580669 0.3198211 1.351750 515 2.794322 4.015805 FALSE
## 15 14 0.4570805 0.3186143 1.349450 678 2.800783 4.015531 FALSE
## 16 15 0.4560945 0.3174043 1.344220 328 2.798089 4.015256 FALSE
## 17 16 0.4551120 0.3162003 1.338620 476 2.794140 4.014981 FALSE
## 18 17 0.4541336 0.3150034 1.336980 669 2.802657 4.014705 FALSE
## 19 18 0.4531548 0.3138013 1.332559 264 2.802424 4.014430 FALSE
## 20 19 0.4521788 0.3126026 1.324110 247 2.789264 4.014153 FALSE
## 21 20 0.4512100 0.3114201 1.323910 483 2.802324 4.013877 FALSE
## 22 21 0.4502392 0.3102280 1.323870 148 2.816092 4.013600 FALSE
## 23 22 0.4492664 0.3090257 1.312720 521 2.794116 4.013323 FALSE
## 24 23 0.4483038 0.3078481 1.294820 653 2.749785 4.013046 FALSE
## 25 24 0.4473590 0.3067162 1.287410 370 2.738855 4.012768 FALSE
## 26 25 0.4464204 0.3055974 1.286970 518 2.750513 4.012490 FALSE
## 27 26 0.4454802 0.3044705 1.264060 326 2.688535 4.012211 FALSE
## 28 27 0.4445635 0.3034043 1.262150 673 2.694709 4.011932 FALSE
## 29 28 0.4436470 0.3023350 1.259550 613 2.698672 4.011653 FALSE
## 30 29 0.4427312 0.3012646 1.251650 875 2.685078 4.011374 FALSE
## 31 30 0.4418223 0.3002092 1.249190 891 2.689350 4.011094 FALSE
## 32 31 0.4409142 0.2991524 1.247980 81 2.697842 4.010814 FALSE
## 33 32 0.4400053 0.2980904 1.247490 514 2.708858 4.010533 FALSE
## 34 33 0.4390950 0.2970209 1.239330 720 2.694204 4.010252 FALSE
## 35 34 0.4381918 0.2959674 1.238180 162 2.702961 4.009971 FALSE
## We can conclude there are no significant outliers
# Assumption 2 - Normality --> achieved with Box-Cox
# Assumption 3 - Sphericity
head(pollen)
## colony round treatment replicate count start_date start_time start_weight
## 1 1.11R2 2 1 11 2 2022-08-22 10:30 1.18198
## 2 1.11R2 2 1 11 3 2022-08-24 12:00 1.14744
## 3 1.11R2 2 1 11 4 2022-08-26 10:30 1.08716
## 4 1.11R2 2 1 11 5 2022-08-28 1:00 1.25522
## 5 1.11R2 2 1 11 6 2022-08-30 11:30 1.14008
## 6 1.11R2 2 1 11 7 2022-09-01 10:30 1.06975
## end_date end_time end_weight whole_dif difference id pid bees_alive
## 1 2022-08-24 13:00 1.08785 0.20545 0.20545 1242 R2.1.11-2 5
## 2 2022-08-26 11:00 1.04712 0.21164 0.21164 1243 R2.1.11-3 5
## 3 2022-08-28 13:40 0.84718 0.35130 0.35130 1244 R2.1.11-4 5
## 4 2022-08-30 12:00 1.14032 0.18044 0.18044 1245 R2.1.11-5 5
## 5 2022-09-01 10:30 1.05520 0.15042 0.15042 1246 R2.1.11-6 5
## 6 2022-09-03 12:20 0.96868 0.16661 0.16661 1247 R2.1.11-7 5
## qro Z1^-3
## 1 B1 0.1430367
## 2 B1 0.1459384
## 3 B1 0.1982433
## 4 B1 0.1306832
## 5 B1 0.1144012
## 6 B1 0.1233902
library(rstatix)
res <- anova_test(data = pollen, dv = difference, wid = id, between = treatment)
res
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 treatment 4 915 2.294 0.058 0.01
# sphericity *barely* not significant
# Pair Wise comparison
pair <- pollen %>%
pairwise_t_test(difference ~ treatment, paired = FALSE, p.adjust.method = "bonferroni")
pair
## # A tibble: 10 × 9
## .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 difference 1 2 195 180 0.039 * 0.39 ns
## 2 difference 1 3 195 190 0.0235 * 0.235 ns
## 3 difference 2 3 180 190 0.865 ns 1 ns
## 4 difference 1 4 195 178 0.0949 ns 0.949 ns
## 5 difference 2 4 180 178 0.703 ns 1 ns
## 6 difference 3 4 190 178 0.578 ns 1 ns
## 7 difference 1 5 195 177 0.887 ns 1 ns
## 8 difference 2 5 180 177 0.0606 ns 0.606 ns
## 9 difference 3 5 190 177 0.0384 * 0.384 ns
## 10 difference 4 5 178 177 0.136 ns 1 ns
## Summary Statistics
sumpol <- pollen %>%
group_by(treatment, count) %>%
get_summary_stats(difference, type = "mean_sd")
data.frame(sumpol)
## treatment count variable n mean sd
## 1 1 2 difference 9 0.275 0.052
## 2 1 3 difference 9 0.234 0.079
## 3 1 4 difference 9 0.233 0.097
## 4 1 5 difference 9 0.139 0.054
## 5 1 6 difference 9 0.164 0.066
## 6 1 7 difference 9 0.250 0.068
## 7 1 8 difference 9 0.278 0.099
## 8 1 9 difference 9 0.330 0.167
## 9 1 10 difference 9 0.486 0.375
## 10 1 11 difference 9 0.671 0.486
## 11 1 12 difference 9 0.615 0.362
## 12 1 13 difference 9 0.665 0.440
## 13 1 14 difference 9 0.620 0.407
## 14 1 15 difference 9 0.596 0.390
## 15 1 16 difference 9 0.647 0.365
## 16 1 17 difference 8 0.712 0.329
## 17 1 18 difference 9 0.605 0.387
## 18 1 19 difference 9 0.669 0.439
## 19 1 20 difference 9 0.507 0.333
## 20 1 21 difference 5 0.255 0.097
## 21 1 22 difference 5 0.236 0.087
## 22 1 23 difference 5 0.214 0.046
## 23 1 24 difference 2 0.269 0.064
## 24 1 25 difference 2 0.277 0.110
## 25 1 26 difference 2 0.298 0.142
## 26 1 27 difference 2 0.258 0.071
## 27 1 28 difference 1 0.207 NA
## 28 1 29 difference 1 0.250 NA
## 29 2 2 difference 9 0.227 0.073
## 30 2 3 difference 9 0.241 0.038
## 31 2 4 difference 9 0.259 0.088
## 32 2 5 difference 9 0.246 0.116
## 33 2 6 difference 9 0.314 0.122
## 34 2 7 difference 9 0.408 0.384
## 35 2 8 difference 9 0.353 0.088
## 36 2 9 difference 9 0.384 0.087
## 37 2 10 difference 9 0.581 0.324
## 38 2 11 difference 9 0.713 0.313
## 39 2 12 difference 9 0.913 0.334
## 40 2 13 difference 9 0.790 0.330
## 41 2 14 difference 9 0.831 0.498
## 42 2 15 difference 9 0.733 0.464
## 43 2 16 difference 8 0.676 0.552
## 44 2 17 difference 9 0.637 0.446
## 45 2 18 difference 9 0.479 0.263
## 46 2 19 difference 8 0.478 0.207
## 47 2 20 difference 8 0.396 0.178
## 48 2 21 difference 7 0.401 0.144
## 49 2 22 difference 3 0.417 0.113
## 50 2 23 difference 1 0.440 NA
## 51 2 24 difference 1 0.408 NA
## 52 3 2 difference 9 0.231 0.086
## 53 3 3 difference 9 0.245 0.090
## 54 3 4 difference 9 0.226 0.105
## 55 3 5 difference 7 0.274 0.033
## 56 3 6 difference 9 0.247 0.078
## 57 3 7 difference 9 0.294 0.140
## 58 3 8 difference 9 0.346 0.200
## 59 3 9 difference 9 0.425 0.269
## 60 3 10 difference 9 0.535 0.385
## 61 3 11 difference 9 0.657 0.399
## 62 3 12 difference 9 0.815 0.381
## 63 3 13 difference 9 0.754 0.392
## 64 3 14 difference 9 0.793 0.471
## 65 3 15 difference 9 0.683 0.382
## 66 3 16 difference 9 0.697 0.356
## 67 3 17 difference 9 0.780 0.422
## 68 3 18 difference 9 0.718 0.368
## 69 3 19 difference 9 0.669 0.265
## 70 3 20 difference 8 0.531 0.140
## 71 3 21 difference 6 0.421 0.177
## 72 3 22 difference 5 0.352 0.125
## 73 3 23 difference 4 0.341 0.098
## 74 3 24 difference 3 0.575 0.385
## 75 3 25 difference 1 0.267 NA
## 76 3 26 difference 1 0.262 NA
## 77 3 27 difference 1 0.234 NA
## 78 3 28 difference 1 0.111 NA
## 79 4 2 difference 9 0.276 0.107
## 80 4 3 difference 9 0.301 0.084
## 81 4 4 difference 8 0.270 0.104
## 82 4 5 difference 8 0.235 0.101
## 83 4 6 difference 7 0.259 0.102
## 84 4 7 difference 9 0.243 0.093
## 85 4 8 difference 9 0.315 0.152
## 86 4 9 difference 9 0.492 0.366
## 87 4 10 difference 9 0.534 0.238
## 88 4 11 difference 9 0.793 0.386
## 89 4 12 difference 9 0.857 0.421
## 90 4 13 difference 9 0.797 0.430
## 91 4 14 difference 8 0.742 0.430
## 92 4 15 difference 9 0.591 0.351
## 93 4 16 difference 9 0.475 0.285
## 94 4 17 difference 9 0.508 0.276
## 95 4 18 difference 9 0.641 0.412
## 96 4 19 difference 9 0.530 0.294
## 97 4 20 difference 8 0.652 0.412
## 98 4 21 difference 4 0.375 0.303
## 99 4 22 difference 2 0.244 0.028
## 100 4 23 difference 2 0.248 0.038
## 101 4 24 difference 1 0.212 NA
## 102 4 25 difference 1 0.239 NA
## 103 4 26 difference 1 0.240 NA
## 104 4 27 difference 1 0.260 NA
## 105 4 28 difference 1 0.246 NA
## 106 5 2 difference 9 0.232 0.060
## 107 5 3 difference 9 0.250 0.054
## 108 5 4 difference 8 0.179 0.084
## 109 5 5 difference 9 0.191 0.069
## 110 5 6 difference 8 0.148 0.072
## 111 5 7 difference 9 0.186 0.127
## 112 5 8 difference 9 0.293 0.149
## 113 5 9 difference 9 0.386 0.243
## 114 5 10 difference 9 0.449 0.351
## 115 5 11 difference 8 0.620 0.419
## 116 5 12 difference 7 0.811 0.259
## 117 5 13 difference 8 0.815 0.385
## 118 5 14 difference 8 0.751 0.374
## 119 5 15 difference 8 0.628 0.362
## 120 5 16 difference 7 0.523 0.295
## 121 5 17 difference 8 0.653 0.317
## 122 5 18 difference 8 0.548 0.184
## 123 5 19 difference 8 0.558 0.256
## 124 5 20 difference 7 0.676 0.331
## 125 5 21 difference 7 0.341 0.125
## 126 5 22 difference 5 0.290 0.109
## 127 5 23 difference 3 0.272 0.126
## 128 5 24 difference 1 0.181 NA
## 129 5 25 difference 1 0.162 NA
## 130 5 26 difference 1 0.184 NA
## 131 5 27 difference 1 0.202 NA
## 132 5 28 difference 1 0.193 NA
## 133 5 29 difference 1 0.181 NA
## effect of treatment on each pollen ball across time
one.way <- pollen %>%
group_by(treatment) %>%
anova_test(dv = difference, wid = pid, between = count) %>%
get_anova_table() %>%
adjust_pvalue(method = "bonferroni")
one.way
## # A tibble: 5 × 9
## treatment Effect DFn DFd F p `p<.05` ges p.adj
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 count 27 167 3.29 1.37e- 6 * 0.347 6.85e- 6
## 2 2 count 22 157 4.13 7.07e- 8 * 0.367 3.54e- 7
## 3 3 count 26 163 3.89 4.81e- 8 * 0.383 2.41e- 7
## 4 4 count 26 151 3.23 3.6 e- 6 * 0.357 1.8 e- 5
## 5 5 count 27 149 4.99 6.95e-11 * 0.475 3.48e-10
pollen %>%
group_by(treatment) %>%
pairwise_t_test(difference ~ count, paired = FALSE,
p.adjust.method = "bonferroni")
## # A tibble: 0 × 10
## # … with 10 variables: treatment <fct>, .y. <chr>, group1 <chr>, group2 <chr>,
## # n1 <int>, n2 <int>, p <dbl>, p.signif <chr>, p.adj <dbl>,
## # p.adj.signif <chr>
pollen%>%
pairwise_t_test(
difference ~ treatment,
paired = FALSE,
p.adjust.method = "bonferroni"
)
## # A tibble: 10 × 9
## .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 difference 1 2 195 180 0.039 * 0.39 ns
## 2 difference 1 3 195 190 0.0235 * 0.235 ns
## 3 difference 2 3 180 190 0.865 ns 1 ns
## 4 difference 1 4 195 178 0.0949 ns 0.949 ns
## 5 difference 2 4 180 178 0.703 ns 1 ns
## 6 difference 3 4 190 178 0.578 ns 1 ns
## 7 difference 1 5 195 177 0.887 ns 1 ns
## 8 difference 2 5 180 177 0.0606 ns 0.606 ns
## 9 difference 3 5 190 177 0.0384 * 0.384 ns
## 10 difference 4 5 178 177 0.136 ns 1 ns
polsum <- pollen %>%
group_by(treatment) %>%
summarise(mp = mean(difference),
sdp = sd(difference),
np = length(difference)) %>%
mutate(sep = sdp/sqrt(np))
polsum
## # A tibble: 5 × 5
## treatment mp sdp np sep
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.430 0.336 195 0.0240
## 2 2 0.502 0.348 180 0.0259
## 3 3 0.508 0.345 190 0.0250
## 4 4 0.488 0.342 178 0.0256
## 5 5 0.435 0.316 177 0.0238
tidy_anova <- broom::tidy(pbox)
knitr::kable(tidy_anova)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -28.0610062 | 2.8408405 | -9.877713 | 0.0000000 |
| treatment2 | 0.0165456 | 0.0060174 | 2.749634 | 0.0060844 |
| treatment3 | 0.0160980 | 0.0059413 | 2.709508 | 0.0068650 |
| treatment4 | 0.0216683 | 0.0060636 | 3.573491 | 0.0003708 |
| treatment5 | -0.0027254 | 0.0060775 | -0.448447 | 0.6539374 |
| bees_alive | 0.0284798 | 0.0031226 | 9.120646 | 0.0000000 |
| qroB3 | 0.0202435 | 0.0068180 | 2.969139 | 0.0030647 |
| qroB4 | 0.0546267 | 0.0067669 | 8.072624 | 0.0000000 |
| qroB5 | 0.0153108 | 0.0048522 | 3.155436 | 0.0016552 |
| start_date | 0.0014607 | 0.0001474 | 9.909582 | 0.0000000 |
anova_summary <- summary(pbox)
tukey_treatment <- HSD.test(pbox,
trt = "treatment",
console = TRUE) # prints the results to console
##
## Study: pbox ~ "treatment"
##
## HSD Test for boxp
##
## Mean Square Error: 0.003361313
##
## treatment, means
##
## boxp std r Min Max
## 1 0.1899525 0.06548217 195 0.041653345 0.3099233
## 2 0.2100114 0.05728009 180 0.056362833 0.3135908
## 3 0.2112073 0.05879453 190 0.032092405 0.3134084
## 4 0.2050092 0.06367876 178 0.002700324 0.3076306
## 5 0.1924815 0.06789212 177 0.029659134 0.3041338
##
## Alpha: 0.05 ; DF Error: 910
## Critical Value of Studentized Range: 3.865405
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## boxp groups
## 3 0.2112073 a
## 2 0.2100114 a
## 4 0.2050092 ab
## 5 0.1924815 b
## 1 0.1899525 b
tidy_tukey <- as.data.frame(tukey_treatment$groups)
tidy_tukey
## boxp groups
## 3 0.2112073 a
## 2 0.2100114 a
## 4 0.2050092 ab
## 5 0.1924815 b
## 1 0.1899525 b
tidier_tukey <- tidy_tukey %>%
rownames_to_column() %>% # converts rownames to columns
rename(treatment = rowname)
p_max <- pollen %>%
group_by(treatment) %>%
summarize(maxp = max(mean(difference)))
p_for_plotting <- full_join(tidier_tukey, p_max,
by = "treatment")
library(ggpmisc)
polp <- ggplot(data = polsum, aes(x = treatment, y = mp, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(0.4,0.55)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = mp - sep,
ymax = mp + sep),
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 = 20) +
annotate(geom = "text",
x = 2, y = 0.55,
label = "ANOVA = 3.914e-05",
size = 8) +
annotate(geom = "text",
x = c(3, 2, 4, 5, 1 ),
y = 0.034 + p_for_plotting$maxp,
label = p_for_plotting$groups,
size = 8)+
theme(legend.position = "none")
polp
drones <- read_csv("drones_rf.csv", col_types = cols(round = col_factor(levels = c("1","2")), treatment = col_factor(levels = c("1","2", "3", "4", "5")), notes = col_skip(), colony_start = col_skip(), colony_count = col_skip(), alive = col_skip(), emerge_date = col_skip()))
#this data set is for drone emerge times and health metrics
drones$order_on_sheet <- as.factor(drones$order_on_sheet)
drones$replicate <- as.factor(drones$replicate)
drones$colony <- as.factor(drones$colony)
drones$id <- as.factor(drones$id)
drones$relative_fat <- as.double(drones$relative_fat)
drones$radial <- as.double(drones$radial)
drones$`alive?` <- as.double(drones$`alive?`)
drf.pollen <- merge(wd.summary, drones, by.x = "colony") # this is a new data frame with average pollen consumption data per colony
drf.pollen$alive <- as.logical(drf.pollen$`alive?`)
drone.sum <- drf.pollen %>%
group_by(colony) %>%
summarise( count.drone = length(id))
drone.sum
## # A tibble: 40 × 2
## colony count.drone
## <fct> <int>
## 1 1.11R2 4
## 2 1.1R2 6
## 3 1.2R2 12
## 4 1.3R2 11
## 5 1.4R2 16
## 6 1.5R2 21
## 7 1.7R2 12
## 8 2.11R2 9
## 9 2.12R2 5
## 10 2.1R2 11
## # … with 30 more rows
drone.sum <- as.data.frame(drone.sum)
qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)
drf.pollen <- merge(drf.pollen, qro, by.x = "colony")
range(drf.pollen$dry_weight)
## [1] 0.0166 0.0826
ggplot(drf.pollen, aes(x=dry_weight, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.0025,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Dry Weight") +
labs(y = "Count", x = "Weight (g) ")
shapiro.test(drf.pollen$dry_weight)
##
## Shapiro-Wilk normality test
##
## data: drf.pollen$dry_weight
## W = 0.98886, p-value = 0.002697
drf.pollen <- na.omit(drf.pollen)
dry.g1 <- lmer(dry_weight~ treatment + whole.mean + workers_alive + alive + emerge_time + qro + (1|colony), data = drf.pollen)
plot(drf.pollen$treatment, drf.pollen$dry_weight)
summary(dry.g1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: dry_weight ~ treatment + whole.mean + workers_alive + alive +
## emerge_time + qro + (1 | colony)
## Data: drf.pollen
##
## REML criterion at convergence: -2472.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66221 -0.60404 -0.01673 0.64927 3.09019
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 3.839e-06 0.001959
## Residual 5.988e-05 0.007738
## Number of obs: 380, groups: colony, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.945e-02 8.827e-03 3.336
## treatment2 -3.700e-03 1.697e-03 -2.180
## treatment3 -6.636e-03 1.817e-03 -3.653
## treatment4 -1.241e-04 1.823e-03 -0.068
## treatment5 -7.731e-04 1.785e-03 -0.433
## whole.mean 3.347e-04 4.310e-03 0.078
## workers_alive -9.352e-04 7.319e-04 -1.278
## aliveTRUE 1.226e-02 5.761e-03 2.128
## emerge_time 1.229e-04 1.341e-04 0.917
## qroB3 -4.817e-04 1.857e-03 -0.259
## qroB4 3.762e-03 1.985e-03 1.895
## qroB5 8.546e-05 1.663e-03 0.051
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn wrkrs_ alTRUE emrg_t
## treatment2 -0.178
## treatment3 -0.070 0.481
## treatment4 -0.100 0.504 0.507
## treatment5 -0.118 0.527 0.487 0.537
## whole.mean -0.262 0.042 -0.129 -0.180 0.114
## workers_alv -0.316 -0.054 -0.168 -0.249 -0.291 -0.047
## aliveTRUE -0.630 0.048 0.098 0.021 -0.006 -0.131 0.017
## emerge_time -0.626 0.077 0.013 0.204 0.160 0.189 -0.093 0.009
## qroB3 -0.097 0.113 0.049 0.080 0.208 0.003 -0.084 0.069 0.050
## qroB4 -0.011 0.058 0.037 0.220 -0.032 -0.512 0.219 0.037 -0.026
## qroB5 -0.160 0.065 -0.085 -0.149 -0.116 0.035 0.623 -0.031 -0.185
## qroB3 qroB4
## treatment2
## treatment3
## treatment4
## treatment5
## whole.mean
## workers_alv
## aliveTRUE
## emerge_time
## qroB3
## qroB4 0.153
## qroB5 0.145 0.279
Anova(dry.g1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 20.7427 4 0.0003561 ***
## whole.mean 0.0060 1 0.9380898
## workers_alive 1.6325 1 0.2013566
## alive 4.5288 1 0.0333296 *
## emerge_time 0.8403 1 0.3593005
## qro 4.0927 3 0.2516228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(resid(dry.g1)) +
abline(h=0, col="red", lwd=2)
## integer(0)
drop1(dry.g1, test = "Chisq")
## Single term deletions
##
## Model:
## dry_weight ~ treatment + whole.mean + workers_alive + alive +
## emerge_time + qro + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -2585.2
## treatment 4 -2571.2 22.0881 0.0001925 ***
## whole.mean 1 -2587.2 0.0000 0.9983807
## workers_alive 1 -2584.8 2.4544 0.1171967
## alive 1 -2582.5 4.7517 0.0292698 *
## emerge_time 1 -2585.9 1.3549 0.2444249
## qro 3 -2584.9 6.2983 0.0979676 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dry2 <- update(dry.g1, .~. -whole.mean)
dry3 <- update(dry2, .~. -emerge_time)
anova(dry.g1, dry2, dry3)
## Data: drf.pollen
## Models:
## dry3: dry_weight ~ treatment + workers_alive + alive + qro + (1 | colony)
## dry2: dry_weight ~ treatment + workers_alive + alive + emerge_time + qro + (1 | colony)
## dry.g1: dry_weight ~ treatment + whole.mean + workers_alive + alive + emerge_time + qro + (1 | colony)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## dry3 12 -2587.8 -2540.6 1305.9 -2611.8
## dry2 13 -2587.2 -2536.0 1306.6 -2613.2 1.4113 1 0.2348
## dry.g1 14 -2585.2 -2530.1 1306.6 -2613.2 0.0000 1 0.9984
Anova(dry3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 20.7668 4 0.0003522 ***
## workers_alive 1.4774 1 0.2241741
## alive 4.5421 1 0.0330709 *
## qro 5.7551 3 0.1241502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dry3
## Linear mixed model fit by REML ['lmerMod']
## Formula: dry_weight ~ treatment + workers_alive + alive + qro + (1 | colony)
## Data: drf.pollen
## REML criterion at convergence: -2496.839
## Random effects:
## Groups Name Std.Dev.
## colony (Intercept) 0.001889
## Residual 0.007733
## Number of obs: 380, groups: colony, 39
## Fixed Effects:
## (Intercept) treatment2 treatment3 treatment4 treatment5
## 0.0343918 -0.0038324 -0.0067044 -0.0005334 -0.0010362
## workers_alive aliveTRUE qroB3 qroB4 qroB5
## -0.0008716 0.0121367 -0.0005733 0.0037085 0.0003828
anova(dry3)
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 4 0.00134261 0.00033565 5.6125
## workers_alive 1 0.00026194 0.00026194 4.3799
## alive 1 0.00030874 0.00030874 5.1625
## qro 3 0.00034418 0.00011473 1.9184
qqnorm(resid(dry.g1));qqline(resid(dry.g1)) #diagnostics look good
qqnorm(resid(dry3));qqline(resid(dry3))
dry.emm <- emmeans(dry3, "treatment")
pairs(dry.emm)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.003832 0.00168 24.6 2.286 0.1831
## treatment1 - treatment3 0.006704 0.00179 26.2 3.746 0.0073
## treatment1 - treatment4 0.000533 0.00172 22.9 0.310 0.9978
## treatment1 - treatment5 0.001036 0.00174 23.0 0.597 0.9742
## treatment2 - treatment3 0.002872 0.00176 28.3 1.635 0.4886
## treatment2 - treatment4 -0.003299 0.00166 23.9 -1.983 0.3042
## treatment2 - treatment5 -0.002796 0.00167 24.2 -1.679 0.4647
## treatment3 - treatment4 -0.006171 0.00176 24.6 -3.515 0.0135
## treatment3 - treatment5 -0.005668 0.00175 24.9 -3.240 0.0255
## treatment4 - treatment5 0.000503 0.00162 19.9 0.310 0.9978
##
## Results are averaged over the levels of: alive, qro
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
plot(dry.emm, comparisons = TRUE)
dry.sum <- drf.pollen %>%
group_by(treatment) %>%
summarise(dry.m = mean(dry_weight),
dry.sd = sd(dry_weight),
dryn = length(dry_weight)) %>%
mutate(sedry = dry.sd/sqrt(dryn))
dry.sum
## # A tibble: 5 × 5
## treatment dry.m dry.sd dryn sedry
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.0445 0.00811 74 0.000942
## 2 2 0.0396 0.00836 75 0.000966
## 3 3 0.0365 0.00844 59 0.00110
## 4 4 0.0417 0.00845 89 0.000895
## 5 5 0.0422 0.00701 83 0.000770
Interactions
#interactions between mean pollen consumed and treatment
m1 <- lm(dry_weight ~ treatment*whole.mean, data = drf.pollen)
summary(m1)
##
## Call:
## lm(formula = dry_weight ~ treatment * whole.mean, data = drf.pollen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0221597 -0.0053496 -0.0004029 0.0053423 0.0265321
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.713e-02 3.566e-03 10.413 <2e-16 ***
## treatment2 -2.545e-03 5.760e-03 -0.442 0.6588
## treatment3 -2.213e-03 4.882e-03 -0.453 0.6506
## treatment4 7.903e-03 5.054e-03 1.564 0.1188
## treatment5 -7.798e-06 5.382e-03 -0.001 0.9988
## whole.mean 1.299e-02 6.104e-03 2.128 0.0340 *
## treatment2:whole.mean -3.633e-03 1.031e-02 -0.352 0.7247
## treatment3:whole.mean -1.043e-02 7.873e-03 -1.325 0.1860
## treatment4:whole.mean -1.839e-02 8.314e-03 -2.212 0.0276 *
## treatment5:whole.mean -3.275e-03 9.740e-03 -0.336 0.7369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.008029 on 370 degrees of freedom
## Multiple R-squared: 0.1088, Adjusted R-squared: 0.08713
## F-statistic: 5.019 on 9 and 370 DF, p-value: 2.138e-06
Anova(m1)
## Anova Table (Type II tests)
##
## Response: dry_weight
## Sum Sq Df F value Pr(>F)
## treatment 0.0024809 4 9.6201 2.071e-07 ***
## whole.mean 0.0001702 1 2.6403 0.1050
## treatment:whole.mean 0.0003854 4 1.4945 0.2032
## Residuals 0.0238542 370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(drf.pollen, aes(x = whole.mean, y = dry_weight, color = treatment)) +
geom_point() +
geom_smooth(method = lm)+
scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "black", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
theme_cowplot()+
labs(title = "Interaction Between Drone Dry Weight and Average Pollen Consumption",
x = "Average Pollen Consumption (g)",
y= "Drone Dry Weight (g)")
dp <- ggplot(data = dry.sum, aes(x = treatment, y = dry.m, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(0.03, 0.046)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = dry.m - sedry,
ymax = dry.m + sedry),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean 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 = 20)
dry3 <- lmer(dry_weight ~ treatment + workers_alive + alive + qro + (1 | colony), data = drf.pollen)
Anova(dry3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 20.7668 4 0.0003522 ***
## workers_alive 1.4774 1 0.2241741
## alive 4.5421 1 0.0330709 *
## qro 5.7551 3 0.1241502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drytidy_anova <- broom::tidy(dry3)
drytidy_anova
## # A tibble: 12 × 6
## effect group term estimate std.error statistic
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 0.0344 0.00672 5.12
## 2 fixed <NA> treatment2 -0.00383 0.00167 -2.30
## 3 fixed <NA> treatment3 -0.00670 0.00177 -3.78
## 4 fixed <NA> treatment4 -0.000533 0.00171 -0.312
## 5 fixed <NA> treatment5 -0.00104 0.00173 -0.599
## 6 fixed <NA> workers_alive -0.000872 0.000717 -1.22
## 7 fixed <NA> aliveTRUE 0.0121 0.00569 2.13
## 8 fixed <NA> qroB3 -0.000573 0.00183 -0.313
## 9 fixed <NA> qroB4 0.00371 0.00167 2.22
## 10 fixed <NA> qroB5 0.000383 0.00161 0.238
## 11 ran_pars colony sd__(Intercept) 0.00189 NA NA
## 12 ran_pars Residual sd__Observation 0.00773 NA NA
dryanova_summary <- summary(dry3)
drytuk.means <- emmeans(object = dry3,
specs = "treatment",
adjust = "Tukey",
type = "response")
drytuk.means
## treatment emmean SE df lower.CL upper.CL
## 1 0.0375 0.00309 224 0.0295 0.0455
## 2 0.0337 0.00301 237 0.0259 0.0415
## 3 0.0308 0.00299 280 0.0231 0.0386
## 4 0.0370 0.00311 211 0.0289 0.0450
## 5 0.0365 0.00310 213 0.0285 0.0445
##
## Results are averaged over the levels of: alive, qro
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
pairs(drytuk.means)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.003832 0.00168 24.6 2.286 0.1831
## treatment1 - treatment3 0.006704 0.00179 26.2 3.746 0.0073
## treatment1 - treatment4 0.000533 0.00172 22.9 0.310 0.9978
## treatment1 - treatment5 0.001036 0.00174 23.0 0.597 0.9742
## treatment2 - treatment3 0.002872 0.00176 28.3 1.635 0.4886
## treatment2 - treatment4 -0.003299 0.00166 23.9 -1.983 0.3042
## treatment2 - treatment5 -0.002796 0.00167 24.2 -1.679 0.4647
## treatment3 - treatment4 -0.006171 0.00176 24.6 -3.515 0.0135
## treatment3 - treatment5 -0.005668 0.00175 24.9 -3.240 0.0255
## treatment4 - treatment5 0.000503 0.00162 19.9 0.310 0.9978
##
## Results are averaged over the levels of: alive, qro
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
dry.cld.model <- cld(object = drytuk.means,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
dry.cld.model
## treatment emmean SE df lower.CL upper.CL .group
## 3 0.0308 0.00299 280 0.0231 0.0386 a
## 2 0.0337 0.00301 237 0.0259 0.0415 ab
## 5 0.0365 0.00310 213 0.0285 0.0445 b
## 4 0.0370 0.00311 211 0.0289 0.0450 b
## 1 0.0375 0.00309 224 0.0295 0.0455 b
##
## Results are averaged over the levels of: alive, qro
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
drytuk.treatment <- as.data.frame(dry.cld.model)
drytuk.treatment
## treatment emmean SE df lower.CL upper.CL .group
## 3 3 0.03080918 0.002993740 280.1516 0.02306649 0.03855187 a
## 2 2 0.03368122 0.003006538 237.0356 0.02589577 0.04146666 ab
## 5 5 0.03647746 0.003095295 213.4428 0.02845502 0.04449990 b
## 4 4 0.03698026 0.003109611 211.1920 0.02891994 0.04504057 b
## 1 1 0.03751361 0.003090340 224.1624 0.02950745 0.04551978 b
dry_max <- drf.pollen %>%
group_by(treatment) %>%
summarize(maxdry = max(mean(dry_weight)))
dry_for_plotting <- full_join(drytuk.treatment, dry_max,
by="treatment")
Anova(dry3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: dry_weight
## Chisq Df Pr(>Chisq)
## treatment 20.7668 4 0.0003522 ***
## workers_alive 1.4774 1 0.2241741
## alive 4.5421 1 0.0330709 *
## qro 5.7551 3 0.1241502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dp <- dp +
annotate(geom = "text",
x = 2.5, y = 0.045,
label = "ANOVA = 0.0004",
size = 8) +
annotate(geom = "text",
x = c(3, 2, 5, 4, 1),
y = dry_for_plotting$maxdry + 0.002,
label = c("a", "ab", "b", "b", "b"),
size = 8) +
theme(legend.position = "none")
dp
ggplot(drf.pollen, aes(x=treatment, y = dry_weight, fill = treatment)) +
geom_boxplot(position = "identity",col=I("black")) +
ggdist::geom_dots(side="both", color = "black") +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
ggtitle("Drone Dry Weight vs. Treatment") +
labs(y = "Count", x = "Treatment") +
theme_cowplot()+
annotate(geom = "text",
x = 2.5, y = 0.07,
label = "ANOVA = 0.0003522 ***",
size = 5) +
annotate(geom = "text",
x = c(3, 2, 5, 4, 1),
y = dry_for_plotting$maxdry + 0.021,
label = c("a", "ab", "b", "b", "b"),
size = 5) +
theme(legend.position = "none")
shapiro.test(drf.pollen$radial)
##
## Shapiro-Wilk normality test
##
## data: drf.pollen$radial
## W = 0.98323, p-value = 0.0002138
drf.pollen <- na.omit(drf.pollen)
ggplot(drf.pollen, aes(x=radial, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.025,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Dry Radial Cell Length") +
labs(y = "Count", x = "Length (mm) ")
range(drf.pollen$radial)
## [1] 1.6924 3.1542
rad.g1 <- lmer(radial~ treatment + whole.mean + workers_alive + alive + emerge_time + qro + (1|colony), data = drf.pollen)
summary(rad.g1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## radial ~ treatment + whole.mean + workers_alive + alive + emerge_time +
## qro + (1 | colony)
## Data: drf.pollen
##
## REML criterion at convergence: -118.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0124 -0.5616 0.0632 0.5988 3.6513
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.004594 0.06778
## Residual 0.035141 0.18746
## Number of obs: 380, groups: colony, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.119510 0.232552 9.114
## treatment2 -0.027847 0.048779 -0.571
## treatment3 -0.109300 0.051806 -2.110
## treatment4 0.026254 0.052329 0.502
## treatment5 -0.015887 0.051475 -0.309
## whole.mean 0.014518 0.121793 0.119
## workers_alive -0.006245 0.021246 -0.294
## aliveTRUE 0.284881 0.142897 1.994
## emerge_time 0.002043 0.003546 0.576
## qroB3 0.029568 0.052820 0.560
## qroB4 0.032359 0.058199 0.556
## qroB5 0.060243 0.047888 1.258
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn wrkrs_ alTRUE emrg_t
## treatment2 -0.177
## treatment3 -0.084 0.491
## treatment4 -0.096 0.507 0.509
## treatment5 -0.107 0.527 0.494 0.530
## whole.mean -0.270 0.040 -0.112 -0.174 0.104
## workers_alv -0.361 -0.060 -0.164 -0.233 -0.284 -0.041
## aliveTRUE -0.595 0.043 0.106 0.022 -0.006 -0.133 0.015
## emerge_time -0.632 0.074 0.024 0.189 0.143 0.170 -0.083 0.019
## qroB3 -0.079 0.102 0.044 0.076 0.197 -0.021 -0.084 0.060 0.044
## qroB4 -0.009 0.044 0.026 0.203 -0.043 -0.512 0.221 0.034 -0.028
## qroB5 -0.194 0.065 -0.090 -0.137 -0.121 0.035 0.620 -0.038 -0.170
## qroB3 qroB4
## treatment2
## treatment3
## treatment4
## treatment5
## whole.mean
## workers_alv
## aliveTRUE
## emerge_time
## qroB3
## qroB4 0.157
## qroB5 0.142 0.277
Anova(rad.g1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 7.8561 4 0.09699 .
## whole.mean 0.0142 1 0.90512
## workers_alive 0.0864 1 0.76880
## alive 3.9745 1 0.04619 *
## emerge_time 0.3321 1 0.56444
## qro 1.7600 3 0.62367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(rad.g1, test = "Chisq")
## Single term deletions
##
## Model:
## radial ~ treatment + whole.mean + workers_alive + alive + emerge_time +
## qro + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> -150.52
## treatment 4 -148.30 10.2223 0.03684 *
## whole.mean 1 -152.51 0.0123 0.91187
## workers_alive 1 -152.36 0.1607 0.68851
## alive 1 -147.86 4.6548 0.03097 *
## emerge_time 1 -151.95 0.5689 0.45069
## qro 3 -154.11 2.4054 0.49262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rad.g1.1 <- update(rad.g1, .~.-alive)
rad.g1.2 <- update(rad.g1.1, .~.-qro)
Anova(rad.g1.1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 8.6395 4 0.07077 .
## whole.mean 0.1488 1 0.69973
## workers_alive 0.0955 1 0.75734
## emerge_time 0.2585 1 0.61113
## qro 1.7620 3 0.62324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rad.g1.2) # selected model
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: radial
## Chisq Df Pr(>Chisq)
## treatment 9.6886 4 0.04601 *
## whole.mean 0.2288 1 0.63244
## workers_alive 2.1658 1 0.14111
## emerge_time 0.5882 1 0.44310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rad.g1.2
## Linear mixed model fit by REML ['lmerMod']
## Formula: radial ~ treatment + whole.mean + workers_alive + emerge_time +
## (1 | colony)
## Data: drf.pollen
## REML criterion at convergence: -126.9243
## Random effects:
## Groups Name Std.Dev.
## colony (Intercept) 0.06686
## Residual 0.18805
## Number of obs: 380, groups: colony, 39
## Fixed Effects:
## (Intercept) treatment2 treatment3 treatment4 treatment5
## 2.460619 -0.037873 -0.115320 0.030813 -0.009825
## whole.mean workers_alive emerge_time
## 0.048139 -0.023609 0.002669
AIC(rad.g1, rad.g1.1, rad.g1.2)
## df AIC
## rad.g1 14 -90.43206
## rad.g1.1 13 -90.57160
## rad.g1.2 10 -106.92426
plot(resid(rad.g1.2)) +
abline(h=0, col="red", lwd=2)
## integer(0)
qqnorm(resid(rad.g1.2));qqline(resid(rad.g1.2))
rad.emm <- emmeans(rad.g1.2, "treatment")
pairs(rad.emm)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.03787 0.0483 26.9 0.784 0.9329
## treatment1 - treatment3 0.11532 0.0512 30.4 2.252 0.1884
## treatment1 - treatment4 -0.03081 0.0500 25.7 -0.617 0.9711
## treatment1 - treatment5 0.00982 0.0496 25.9 0.198 0.9996
## treatment2 - treatment3 0.07745 0.0500 33.1 1.548 0.5399
## treatment2 - treatment4 -0.06869 0.0476 26.8 -1.442 0.6073
## treatment2 - treatment5 -0.02805 0.0470 26.7 -0.596 0.9744
## treatment3 - treatment4 -0.14613 0.0503 29.2 -2.907 0.0499
## treatment3 - treatment5 -0.10550 0.0507 30.6 -2.081 0.2540
## treatment4 - treatment5 0.04064 0.0476 23.8 0.853 0.9108
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
summary(rad.emm)
## treatment emmean SE df lower.CL upper.CL
## 1 2.49 0.0353 26.3 2.42 2.56
## 2 2.45 0.0333 28.4 2.38 2.52
## 3 2.37 0.0370 35.6 2.30 2.45
## 4 2.52 0.0339 23.5 2.45 2.59
## 5 2.48 0.0341 23.4 2.41 2.55
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
plot(rad.emm)
plot(rad.emm, comparisons = TRUE)
rad.sum <- drf.pollen %>%
group_by(treatment) %>%
summarise(rad.m = mean(radial),
sd.rad = sd(radial),
nrad = length(radial)) %>%
mutate(serad = sd.rad / sqrt(nrad))
plot(drf.pollen$treatment, drf.pollen$radial)
rad.sum
## # A tibble: 5 × 5
## treatment rad.m sd.rad nrad serad
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 2.52 0.162 74 0.0189
## 2 2 2.45 0.183 75 0.0212
## 3 3 2.37 0.245 59 0.0319
## 4 4 2.50 0.222 89 0.0236
## 5 5 2.47 0.166 83 0.0182
#interactions between mean pollen consumed and treatment
m2 <- lm(radial ~ treatment*whole.mean, data = drf.pollen)
summary(m2)
##
## Call:
## lm(formula = radial ~ treatment * whole.mean, data = drf.pollen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80870 -0.11789 0.01653 0.12491 0.72588
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.42946 0.08747 27.774 <2e-16 ***
## treatment2 -0.04333 0.14130 -0.307 0.7593
## treatment3 0.06575 0.11976 0.549 0.5833
## treatment4 0.04964 0.12399 0.400 0.6891
## treatment5 0.01138 0.13202 0.086 0.9314
## whole.mean 0.15272 0.14973 1.020 0.3084
## treatment2:whole.mean -0.02840 0.25286 -0.112 0.9106
## treatment3:whole.mean -0.34600 0.19315 -1.791 0.0741 .
## treatment4:whole.mean -0.11141 0.20396 -0.546 0.5852
## treatment5:whole.mean -0.09398 0.23893 -0.393 0.6943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.197 on 370 degrees of freedom
## Multiple R-squared: 0.06554, Adjusted R-squared: 0.04281
## F-statistic: 2.883 on 9 and 370 DF, p-value: 0.002641
Anova(m2)
## Anova Table (Type II tests)
##
## Response: radial
## Sum Sq Df F value Pr(>F)
## treatment 0.8434 4 5.4345 0.0002917 ***
## whole.mean 0.0001 1 0.0021 0.9636696
## treatment:whole.mean 0.1594 4 1.0271 0.3930733
## Residuals 14.3553 370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(drf.pollen, aes(x = whole.mean, y = radial, fill = treatment)) +
geom_point() +
geom_smooth(method = lm)+
theme_cowplot()+
labs(title = "Interaction Between Radial Cell Length and Average Pollen Consumption",
x = "Average Pollen Consumption (g)",
y= "Radial Cell Length (mm)")
radtidyanova <- broom.mixed::tidy(rad.g1.2)
tidy(rad.g1.2)
## # A tibble: 10 × 6
## effect group term estimate std.error statistic
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 2.46 0.179 13.8
## 2 fixed <NA> treatment2 -0.0379 0.0481 -0.787
## 3 fixed <NA> treatment3 -0.115 0.0509 -2.27
## 4 fixed <NA> treatment4 0.0308 0.0497 0.620
## 5 fixed <NA> treatment5 -0.00982 0.0495 -0.198
## 6 fixed <NA> whole.mean 0.0481 0.101 0.478
## 7 fixed <NA> workers_alive -0.0236 0.0160 -1.47
## 8 fixed <NA> emerge_time 0.00267 0.00348 0.767
## 9 ran_pars colony sd__(Intercept) 0.0669 NA NA
## 10 ran_pars Residual sd__Observation 0.188 NA NA
radtidyanova
## # A tibble: 10 × 6
## effect group term estimate std.error statistic
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 2.46 0.179 13.8
## 2 fixed <NA> treatment2 -0.0379 0.0481 -0.787
## 3 fixed <NA> treatment3 -0.115 0.0509 -2.27
## 4 fixed <NA> treatment4 0.0308 0.0497 0.620
## 5 fixed <NA> treatment5 -0.00982 0.0495 -0.198
## 6 fixed <NA> whole.mean 0.0481 0.101 0.478
## 7 fixed <NA> workers_alive -0.0236 0.0160 -1.47
## 8 fixed <NA> emerge_time 0.00267 0.00348 0.767
## 9 ran_pars colony sd__(Intercept) 0.0669 NA NA
## 10 ran_pars Residual sd__Observation 0.188 NA NA
knitr::kable(radtidyanova)
| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 2.4606188 | 0.1786937 | 13.7700337 |
| fixed | NA | treatment2 | -0.0378733 | 0.0481126 | -0.7871803 |
| fixed | NA | treatment3 | -0.1153201 | 0.0509091 | -2.2652164 |
| fixed | NA | treatment4 | 0.0308126 | 0.0497138 | 0.6197997 |
| fixed | NA | treatment5 | -0.0098248 | 0.0494982 | -0.1984879 |
| fixed | NA | whole.mean | 0.0481393 | 0.1006490 | 0.4782893 |
| fixed | NA | workers_alive | -0.0236094 | 0.0160428 | -1.4716550 |
| fixed | NA | emerge_time | 0.0026694 | 0.0034805 | 0.7669725 |
| ran_pars | colony | sd__(Intercept) | 0.0668568 | NA | NA |
| ran_pars | Residual | sd__Observation | 0.1880492 | NA | NA |
radan_summ <- summary(rad.g1.2)
rad.model.means <- emmeans(object = rad.g1.2,
specs = "treatment")
rad.model.cld <- cld(object = rad.model.means,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
radtuk_treatment <- as.data.frame(rad.model.cld)
rad_tidy_tukey <- as.data.frame(radtuk_treatment$.group)
tidier_rad.tukey <- rad_tidy_tukey %>%
rownames_to_column() %>%
rename(treatment = rowname )
tidier_rad.tukey
## treatment radtuk_treatment$.group
## 1 1 a
## 2 2 ab
## 3 3 ab
## 4 4 ab
## 5 5 b
rad_max <- drf.pollen %>%
group_by(treatment) %>%
summarize(maxrad = max(mean(radial)))
rad_for_plotting <- full_join(tidier_rad.tukey, rad_max,
by="treatment")
rcp <- ggplot(data = rad.sum, aes(x = treatment, y = rad.m, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(2, 2.55)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = rad.m - serad,
ymax = rad.m + serad),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Length (mm) ",) +
ggtitle("Average Radial Cell Length of Drones 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)+
annotate(geom = "text",
x = 3, y = 2.55,
label = "ANOVA = 0.047",
size = 10) +
annotate(geom = "text",
x = c(1, 2, 3, 4, 5),
y = 0.05 + rad_for_plotting$maxrad,
label = c("a","a","a","a","a"),
size = 10) +
theme(legend.position = "none")
rcp
shapiro.test(drf.pollen$relative_fat)
##
## Shapiro-Wilk normality test
##
## data: drf.pollen$relative_fat
## W = 0.80183, p-value < 2.2e-16
range(drf.pollen$relative_fat)
## [1] 0.0002 0.0072
drf.pollen$logrf <- log(drf.pollen$relative_fat)
shapiro.test(drf.pollen$logrf)
##
## Shapiro-Wilk normality test
##
## data: drf.pollen$logrf
## W = 0.93346, p-value = 5.464e-12
ggplot(drf.pollen, aes(x=relative_fat, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.00025,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Relative Fat") +
labs(y = "Count", x = "Relative Fat(g) ")
ggplot(drf.pollen, aes(x=logrf, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.1,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Drone Relative Fat (log transformed)") +
labs(y = "Count", x = "log(Relative Fat(g)) ")
rfglmer <- glmer(logrf ~ treatment + whole.mean + emerge_time + alive + (1|colony), data = drf.pollen)
summary(rfglmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logrf ~ treatment + whole.mean + emerge_time + alive + (1 | colony)
## Data: drf.pollen
##
## REML criterion at convergence: 459.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9538 -0.4786 0.0427 0.4874 3.4402
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.009578 0.09786
## Residual 0.178185 0.42212
## Number of obs: 380, groups: colony, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -7.807165 0.443804 -17.591
## treatment2 -0.068899 0.088289 -0.780
## treatment3 -0.302609 0.094252 -3.211
## treatment4 -0.101569 0.088730 -1.145
## treatment5 0.137935 0.088180 1.564
## whole.mean 0.330769 0.190444 1.737
## emerge_time 0.012872 0.007016 1.835
## aliveTRUE 0.750949 0.311116 2.414
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn emrg_t
## treatment2 -0.206
## treatment3 -0.130 0.476
## treatment4 -0.220 0.519 0.487
## treatment5 -0.215 0.527 0.465 0.526
## whole.mean -0.313 0.064 -0.127 -0.052 0.124
## emerge_time -0.691 0.088 -0.005 0.193 0.146 0.260
## aliveTRUE -0.665 0.049 0.098 0.012 -0.013 -0.128 -0.009
Anova(rfglmer)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logrf
## Chisq Df Pr(>Chisq)
## treatment 23.2142 4 0.0001147 ***
## whole.mean 3.0166 1 0.0824168 .
## emerge_time 3.3664 1 0.0665375 .
## alive 5.8261 1 0.0157903 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rfem <- emmeans(rfglmer, "treatment")
pairs(rfem)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.0689 0.0888 27.1 0.776 0.9354
## treatment1 - treatment3 0.3026 0.0950 29.0 3.185 0.0263
## treatment1 - treatment4 0.1016 0.0894 23.2 1.137 0.7858
## treatment1 - treatment5 -0.1379 0.0886 24.8 -1.557 0.5371
## treatment2 - treatment3 0.2337 0.0944 33.0 2.477 0.1207
## treatment2 - treatment4 0.0327 0.0873 25.5 0.374 0.9956
## treatment2 - treatment5 -0.2068 0.0862 26.6 -2.400 0.1466
## treatment3 - treatment4 -0.2010 0.0936 27.6 -2.149 0.2289
## treatment3 - treatment5 -0.4405 0.0951 30.7 -4.632 0.0006
## treatment4 - treatment5 -0.2395 0.0865 22.4 -2.768 0.0748
##
## Results are averaged over the levels of: alive
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates
plot(resid(rfglmer)) +
abline(h=0, col="red", lwd=2) #this looks good
## integer(0)
qqnorm(resid(rfglmer));qqline(resid(rfglmer)) # this does not :(
drop1(rfglmer, test = "Chisq")
## Single term deletions
##
## Model:
## logrf ~ treatment + whole.mean + emerge_time + alive + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 449.86
## treatment 4 462.15 20.2841 0.0004389 ***
## whole.mean 1 451.48 3.6215 0.0570376 .
## emerge_time 1 451.90 4.0430 0.0443544 *
## alive 1 453.78 5.9219 0.0149536 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(rfglmer)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logrf
## Chisq Df Pr(>Chisq)
## treatment 23.2142 4 0.0001147 ***
## whole.mean 3.0166 1 0.0824168 .
## emerge_time 3.3664 1 0.0665375 .
## alive 5.8261 1 0.0157903 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rf.sum <- drf.pollen %>%
group_by(treatment) %>%
summarise(mrf = mean(relative_fat),
sdrf = sd(relative_fat),
nrf = length(relative_fat)) %>%
mutate(serf = sdrf/sqrt(nrf))
rf.sum
## # A tibble: 5 × 5
## treatment mrf sdrf nrf serf
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 0.00199 0.00102 74 0.000119
## 2 2 0.00168 0.000559 75 0.0000645
## 3 3 0.00143 0.000582 59 0.0000758
## 4 4 0.00169 0.000755 89 0.0000800
## 5 5 0.00212 0.00110 83 0.000121
rfmeans <- emmeans(object = rfglmer,
specs = "treatment")
rfcld <- cld(object = rfmeans,
adjust = "Tukey",
Letters = letters,
alpha = 0.05,
console = TRUE)
rfcld
## treatment emmean SE df lower.CL upper.CL .group
## 3 -7.04 0.163 298 -7.46 -6.62 a
## 4 -6.84 0.167 249 -7.27 -6.41 ab
## 2 -6.81 0.164 268 -7.23 -6.38 ab
## 1 -6.74 0.169 241 -7.18 -6.30 b
## 5 -6.60 0.169 245 -7.04 -6.16 b
##
## Results are averaged over the levels of: alive
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
Anova(rfglmer)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: logrf
## Chisq Df Pr(>Chisq)
## treatment 23.2142 4 0.0001147 ***
## whole.mean 3.0166 1 0.0824168 .
## emerge_time 3.3664 1 0.0665375 .
## alive 5.8261 1 0.0157903 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rfp <- ggplot(data = rf.sum, aes(x = treatment, y = mrf, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0.0011, 0.00229)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = mrf - serf,
ymax = mrf + serf),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Relative Fat (g)",) +
ggtitle("Average Relative Fat (g) of Drones 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)+
annotate(geom = "text",
x = 2, y = 0.00227,
label = "ANOVA = 0.0001147",
size = 8) +
annotate(geom = "text",
x = c(3,4,2,1,5),
y = c(0.00155, 0.00182, 0.0018, 0.00218, 0.00229),
label = c("a", "ab", "ab", "b", "b"),
size = 8) +
theme(legend.position = "none")
rfp
trunc.csv <- read_csv("duration.fortrunc.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"),
emerge_date = col_date(format = "%m/%d/%Y"),
end_date = col_date(format = "%m/%d/%Y"))) # this data set is for drone counts
trunc <- merge(wd.summary, trunc.csv, by.x = "colony")
trunc <- merge(mortality, trunc, by.x = "colony")
trunc <- na.omit(trunc)
trunc$qro <- as.factor(trunc$origin)
trunc <- trunc[-c(8)]
ggplot(drf.pollen, aes(x=emerge_time, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 1,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Average Drone Emerge Time") +
labs(y = "Count", x = "Days")
mod2 <- glm.nb(emerge_time ~ treatment + whole.mean + alive + qro, data = drf.pollen)
summary(mod2)
##
## Call:
## glm.nb(formula = emerge_time ~ treatment + whole.mean + alive +
## qro, data = drf.pollen, init.theta = 2339445.883, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.29006 -0.37800 0.01526 0.34498 1.63512
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.766509 0.116696 32.276 <2e-16 ***
## treatment2 -0.018623 0.026434 -0.705 0.4811
## treatment3 0.013055 0.028139 0.464 0.6427
## treatment4 -0.059797 0.027111 -2.206 0.0274 *
## treatment5 -0.045899 0.026301 -1.745 0.0810 .
## whole.mean -0.162105 0.068033 -2.383 0.0172 *
## aliveTRUE 0.006324 0.113573 0.056 0.9556
## qroB3 -0.017261 0.030350 -0.569 0.5695
## qroB4 -0.004327 0.029767 -0.145 0.8844
## qroB5 0.039768 0.019756 2.013 0.0441 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2339446) family taken to be 1)
##
## Null deviance: 145.60 on 379 degrees of freedom
## Residual deviance: 117.27 on 370 degrees of freedom
## AIC: 2230
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2339446
## Std. Err.: 19359728
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -2207.969
Anova(mod2)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge_time
## LR Chisq Df Pr(>Chisq)
## treatment 10.0183 4 0.04012 *
## whole.mean 5.6827 1 0.01713 *
## alive 0.0031 1 0.95555
## qro 5.5743 3 0.13426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(mod2, test = "F")
## Single term deletions
##
## Model:
## emerge_time ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 117.27 2228.0
## treatment 4 127.29 2230.0 7.9022 4.034e-06 ***
## whole.mean 1 122.95 2231.7 17.9294 2.896e-05 ***
## alive 1 117.27 2226.0 0.0098 0.9211935
## qro 3 122.84 2227.5 5.8625 0.0006412 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3 <- glm.nb(emerge_time ~ treatment + whole.mean + alive, data = drf.pollen)
summary(mod3)
##
## Call:
## glm.nb(formula = emerge_time ~ treatment + whole.mean + alive,
## data = drf.pollen, init.theta = 2216286.051, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.32701 -0.38265 -0.01687 0.36308 1.78327
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.77116 0.11476 32.862 < 2e-16 ***
## treatment2 -0.02121 0.02605 -0.814 0.415504
## treatment3 0.01274 0.02793 0.456 0.648445
## treatment4 -0.05735 0.02540 -2.258 0.023928 *
## treatment5 -0.04802 0.02559 -1.876 0.060595 .
## whole.mean -0.19710 0.05543 -3.556 0.000376 ***
## aliveTRUE 0.03143 0.11284 0.279 0.780576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2216286) family taken to be 1)
##
## Null deviance: 145.60 on 379 degrees of freedom
## Residual deviance: 122.84 on 373 degrees of freedom
## AIC: 2229.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2216286
## Std. Err.: 18049023
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -2213.544
Anova(mod3)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge_time
## LR Chisq Df Pr(>Chisq)
## treatment 10.3601 4 0.0347794 *
## whole.mean 12.6349 1 0.0003786 ***
## alive 0.0784 1 0.7795195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod4 <- glm.nb(emerge_time ~ treatment + whole.mean, data = drf.pollen)
AIC(mod2, mod3, mod4)
## df AIC
## mod2 11 2229.969
## mod3 8 2229.544
## mod4 7 2227.622
plot(mod2)
plot(mod3)
emrem <- emmeans(mod2, "treatment")
Anova(mod2)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge_time
## LR Chisq Df Pr(>Chisq)
## treatment 10.0183 4 0.04012 *
## whole.mean 5.6827 1 0.01713 *
## alive 0.0031 1 0.95555
## qro 5.5743 3 0.13426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairs(emrem)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.0186 0.0264 Inf 0.705 0.9555
## treatment1 - treatment3 -0.0131 0.0281 Inf -0.464 0.9905
## treatment1 - treatment4 0.0598 0.0271 Inf 2.206 0.1775
## treatment1 - treatment5 0.0459 0.0263 Inf 1.745 0.4062
## treatment2 - treatment3 -0.0317 0.0282 Inf -1.124 0.7938
## treatment2 - treatment4 0.0412 0.0267 Inf 1.541 0.5358
## treatment2 - treatment5 0.0273 0.0257 Inf 1.062 0.8261
## treatment3 - treatment4 0.0729 0.0277 Inf 2.628 0.0654
## treatment3 - treatment5 0.0590 0.0283 Inf 2.084 0.2270
## treatment4 - treatment5 -0.0139 0.0270 Inf -0.515 0.9859
##
## Results are averaged over the levels of: alive, qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
library(VGAM) #https://search.r-project.org/CRAN/refmans/VGAM/html/genpoisson2.html
gfit2 <- vglm(emerge_time ~ treatment, genpoisson2, data = drf.pollen, trace = TRUE)
## VGLM linear loop 1 : loglikelihood = -1223.8615
## VGLM linear loop 2 : loglikelihood = -1143.7176
## VGLM linear loop 3 : loglikelihood = -1114.3618
## VGLM linear loop 4 : loglikelihood = -1113.1005
## VGLM linear loop 5 : loglikelihood = -1113.1005
coef(gfit2, matrix = TRUE)
## loglink(meanpar) loglink(disppar)
## (Intercept) 3.691915382 -69.87686
## treatment2 -0.015783500 0.00000
## treatment3 -0.001765548 0.00000
## treatment4 -0.067724273 0.00000
## treatment5 -0.039225096 0.00000
summary(gfit2)
##
## Call:
## vglm(formula = emerge_time ~ treatment, family = genpoisson2,
## data = drf.pollen, trace = TRUE)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 3.692e+00 1.835e-02 201.167 < 2e-16 ***
## (Intercept):2 -6.988e+01 3.804e+04 -0.002 0.99853
## treatment2 -1.578e-02 2.597e-02 -0.608 0.54334
## treatment3 -1.766e-03 2.757e-02 -0.064 0.94894
## treatment4 -6.772e-02 2.523e-02 -2.684 0.00727 **
## treatment5 -3.923e-02 2.548e-02 -1.540 0.12366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: loglink(meanpar), loglink(disppar)
##
## Log-likelihood: -1113.101 on 754 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
emsum <- drf.pollen %>%
group_by(treatment) %>%
summarise( met = mean(emerge_time),
sdet = sd(emerge_time),
net = length(emerge_time)) %>%
mutate( set = sdet/ sqrt(net))
emsum
## # A tibble: 5 × 5
## treatment met sdet net set
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 40.1 5.69 74 0.662
## 2 2 39.5 2.77 75 0.320
## 3 3 40.1 4.06 59 0.528
## 4 4 37.5 2.38 89 0.252
## 5 5 38.6 3.60 83 0.395
emtidy_anova <- broom::tidy(mod2)
emtidy_anova
## # A tibble: 10 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.77 0.117 32.3 1.50e-228
## 2 treatment2 -0.0186 0.0264 -0.705 4.81e- 1
## 3 treatment3 0.0131 0.0281 0.464 6.43e- 1
## 4 treatment4 -0.0598 0.0271 -2.21 2.74e- 2
## 5 treatment5 -0.0459 0.0263 -1.75 8.10e- 2
## 6 whole.mean -0.162 0.0680 -2.38 1.72e- 2
## 7 aliveTRUE 0.00632 0.114 0.0557 9.56e- 1
## 8 qroB3 -0.0173 0.0304 -0.569 5.70e- 1
## 9 qroB4 -0.00433 0.0298 -0.145 8.84e- 1
## 10 qroB5 0.0398 0.0198 2.01 4.41e- 2
amanova_summary <- summary(mod2)
emtukey_treatment <- HSD.test(mod2,
trt = "treatment",
console = TRUE)
##
## Study: mod2 ~ "treatment"
##
## HSD Test for emerge_time
##
## Mean Square Error: 0.3169475
##
## treatment, means
##
## emerge_time std r Min Max
## 1 40.12162 5.692949 74 33 54
## 2 39.49333 2.772321 75 34 47
## 3 40.05085 4.057448 59 33 48
## 4 37.49438 2.379475 89 30 42
## 5 38.57831 3.602575 83 30 46
##
## Alpha: 0.05 ; DF Error: 370
## Critical Value of Studentized Range: 3.876752
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## emerge_time groups
## 1 40.12162 a
## 3 40.05085 a
## 2 39.49333 b
## 5 38.57831 c
## 4 37.49438 d
emtidy_tukey <- as.data.frame(emtukey_treatment$groups)
emtidy_tukey
## emerge_time groups
## 1 40.12162 a
## 3 40.05085 a
## 2 39.49333 b
## 5 38.57831 c
## 4 37.49438 d
emtidier <- emtidy_tukey %>%
rownames_to_column() %>%
rename(treatment = rowname)
emmax <- drf.pollen %>%
group_by(treatment) %>%
summarize(maxem = max(mean(emerge_time)))
em_plotting <- full_join(emtidier, emmax,
by = "treatment")
emp <- ggplot(data = emsum, aes(x = treatment, y = met, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(33, 42)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = met - set,
ymax = met + set),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Ermerge Time (days)",) +
ggtitle("Average Days Until Emergence of Drones 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) +
annotate(geom = "text",
x = 2, y = 42,
label = "ANOVA = 0.041",
size = 8) +
annotate(geom = "text",
x = c(1,3,2,5,4),
y = 1 + em_plotting$maxem,
label = c("a", "a", "b", "c", "d"),
size = 8) +
theme(legend.position = "none")
emp
d2 <- glm.nb(count ~ treatment + whole.mean + alive + qro, data = trunc)
summary(d2)
##
## Call:
## glm.nb(formula = count ~ treatment + whole.mean + alive + qro,
## data = trunc, init.theta = 15.17538143, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0732 -0.8041 -0.3309 0.6613 2.2367
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.957457 0.469261 2.040 0.0413 *
## treatment2 -0.001854 0.213043 -0.009 0.9931
## treatment3 -0.405097 0.214304 -1.890 0.0587 .
## treatment4 -0.007395 0.209181 -0.035 0.9718
## treatment5 0.214929 0.223494 0.962 0.3362
## whole.mean 2.690807 0.455540 5.907 3.49e-09 ***
## alive -0.012515 0.094777 -0.132 0.8949
## qroB3 0.300145 0.215362 1.394 0.1634
## qroB4 -0.133115 0.250252 -0.532 0.5948
## qroB5 0.418361 0.186376 2.245 0.0248 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(15.1754) family taken to be 1)
##
## Null deviance: 102.677 on 39 degrees of freedom
## Residual deviance: 46.061 on 30 degrees of freedom
## AIC: 251.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 15.18
## Std. Err.: 8.87
##
## 2 x log-likelihood: -229.596
Anova(d2)
## Analysis of Deviance Table (Type II tests)
##
## Response: count
## LR Chisq Df Pr(>Chisq)
## treatment 9.406 4 0.05171 .
## whole.mean 34.365 1 4.569e-09 ***
## alive 0.016 1 0.89881
## qro 8.538 3 0.03611 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(d2, test = "F")
## Single term deletions
##
## Model:
## count ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 46.061 249.60
## treatment 4 55.467 251.00 1.5316 0.2182
## whole.mean 1 80.426 281.96 22.3822 4.975e-05 ***
## alive 1 46.077 247.61 0.0105 0.9189
## qro 3 54.599 252.13 1.8536 0.1588
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d3 <- glm.nb(count ~ treatment + whole.mean + qro, data = trunc)
summary(d3)
##
## Call:
## glm.nb(formula = count ~ treatment + whole.mean + qro, data = trunc,
## init.theta = 15.09279171, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0705 -0.8168 -0.3276 0.6468 2.2327
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.905720 0.278275 3.255 0.00113 **
## treatment2 -0.008507 0.206268 -0.041 0.96710
## treatment3 -0.411096 0.209316 -1.964 0.04953 *
## treatment4 -0.011237 0.206220 -0.054 0.95654
## treatment5 0.204427 0.205286 0.996 0.31934
## whole.mean 2.683868 0.450197 5.962 2.5e-09 ***
## qroB3 0.297343 0.214824 1.384 0.16632
## qroB4 -0.123793 0.238296 -0.519 0.60342
## qroB5 0.431772 0.155322 2.780 0.00544 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(15.0928) family taken to be 1)
##
## Null deviance: 102.471 on 39 degrees of freedom
## Residual deviance: 45.989 on 31 degrees of freedom
## AIC: 249.61
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 15.09
## Std. Err.: 8.77
##
## 2 x log-likelihood: -229.612
Anova(d3)
## Analysis of Deviance Table (Type II tests)
##
## Response: count
## LR Chisq Df Pr(>Chisq)
## treatment 9.390 4 0.05205 .
## whole.mean 34.811 1 3.634e-09 ***
## qro 9.909 3 0.01935 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d3em <- emmeans(d3, "treatment")
pairs(d3em)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.00851 0.206 Inf 0.041 1.0000
## treatment1 - treatment3 0.41110 0.209 Inf 1.964 0.2838
## treatment1 - treatment4 0.01124 0.206 Inf 0.054 1.0000
## treatment1 - treatment5 -0.20443 0.205 Inf -0.996 0.8574
## treatment2 - treatment3 0.40259 0.207 Inf 1.944 0.2938
## treatment2 - treatment4 0.00273 0.205 Inf 0.013 1.0000
## treatment2 - treatment5 -0.21293 0.201 Inf -1.061 0.8263
## treatment3 - treatment4 -0.39986 0.205 Inf -1.949 0.2912
## treatment3 - treatment5 -0.61552 0.210 Inf -2.937 0.0274
## treatment4 - treatment5 -0.21566 0.207 Inf -1.043 0.8354
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
d1sum <- trunc %>%
group_by(treatment) %>%
summarise(md = mean(count),
sd = sd(count),
nd = length(count)) %>%
mutate(sed = sd/sqrt(nd))
d1sum
## # A tibble: 5 × 5
## treatment md sd nd sed
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 12.4 5.86 7 2.21
## 2 2 11 5.95 8 2.10
## 3 3 8.67 6.76 9 2.25
## 4 4 13.8 8.45 8 2.99
## 5 5 12.2 5.97 8 2.11
counttuk.means <- emmeans(object = d3,
specs = "treatment",
adjust = "Tukey",
type = "response")
counttuk.means
## treatment response SE df asymp.LCL asymp.UCL
## 1 11.3 1.74 Inf 7.62 16.8
## 2 11.2 1.70 Inf 7.59 16.6
## 3 7.5 1.18 Inf 5.01 11.2
## 4 11.2 1.82 Inf 7.36 17.0
## 5 13.9 2.08 Inf 9.44 20.4
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the log scale
count.cld.model <- cld(object = counttuk.means,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
count.cld.model
## treatment response SE df asymp.LCL asymp.UCL .group
## 3 7.5 1.18 Inf 5.01 11.2 a
## 4 11.2 1.82 Inf 7.36 17.0 ab
## 2 11.2 1.70 Inf 7.59 16.6 ab
## 1 11.3 1.74 Inf 7.62 16.8 ab
## 5 13.9 2.08 Inf 9.44 20.4 b
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
count_max <- trunc %>%
group_by(treatment) %>%
summarize(maxcount = max(mean(count)))
count_max <- as.data.frame(count_max)
ggplot(data = d1sum, aes(x = treatment, y = md, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0,18)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = md - sed,
ymax = md + sed),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Count",) +
ggtitle("Average Count of Drones 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) +
annotate(geom = "text",
x = 2.5, y = 17.5,
label = "Anova = 0.05205",
size = 10) +
annotate(geom = "text",
x = c(3, 4, 2, 1, 5),
y = c(12, 17.5, 14, 16, 15.2),
label = c("a", "ab", "ab", "ab", "b"),
size = 10)
w <- read_csv("weights1.csv", col_types = cols(round = col_factor(levels = c("1",
"2")), date = col_date(format = "%m/%d/%Y"),
treatment = col_factor(levels = c("1",
"2", "3", "4", "5"))))
w$number <- as.factor(w$number)
w$replicate <- as.factor(w$replicate)
w$id <- as.factor(w$id)
w$colony <- as.factor(w$colony)
w <- merge(wd.summary, w, by.x = "colony")
w <- merge(w, qro, by.x = "colony")
shapiro.test(w$weight)
##
## Shapiro-Wilk normality test
##
## data: w$weight
## W = 0.93378, p-value = 1.752e-11
range(w$weight)
## [1] 46.6 68.0
w$boxw <- bcPower(w$weight, -2, gamma=1)
shapiro.test(w$boxw)
##
## Shapiro-Wilk normality test
##
## data: w$boxw
## W = 0.97397, p-value = 5.078e-06
ggplot(w, aes(x=weight, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.5,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Colony Weight") +
labs(y = "Count", x = "Weight (g)")
ggplot(w, aes(x=boxw, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.000005,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Colony Weight - BoxCox Transformed") +
labs(y = "Count", x = "Weight (g) - Box-Cox")
ggplot(w, aes(x=treatment, y = boxw, fill = treatment)) +
geom_boxplot(position = "identity", binwidth = 0.000005,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Colony Weight - BoxCox Transformed") +
labs(y = "Count", x = "Weight (g) - Box-Cox")
wmod <- glm(boxw ~ treatment + days*treatment + whole.mean + bees_alive + qro + days, data = w)
summary(wmod)
##
## Call:
## glm(formula = boxw ~ treatment + days * treatment + whole.mean +
## bees_alive + qro + days, data = w)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.661e-05 -9.357e-06 -2.180e-07 9.239e-06 2.979e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.998e-01 6.511e-06 76763.177 < 2e-16 ***
## treatment2 -5.070e-06 3.689e-06 -1.374 0.17021
## treatment3 -8.807e-06 3.663e-06 -2.404 0.01673 *
## treatment4 -4.453e-06 3.677e-06 -1.211 0.22671
## treatment5 1.137e-07 3.715e-06 0.031 0.97559
## days 1.010e-06 9.406e-08 10.735 < 2e-16 ***
## whole.mean 4.816e-05 4.607e-06 10.454 < 2e-16 ***
## bees_alive 7.398e-07 1.249e-06 0.592 0.55404
## qroB3 -1.549e-06 2.468e-06 -0.628 0.53064
## qroB4 -3.649e-06 2.676e-06 -1.364 0.17360
## qroB5 9.394e-07 1.761e-06 0.533 0.59414
## treatment2:days 4.535e-07 1.398e-07 3.243 0.00130 **
## treatment3:days 5.424e-07 1.330e-07 4.079 5.64e-05 ***
## treatment4:days 4.089e-07 1.390e-07 2.941 0.00349 **
## treatment5:days 2.692e-07 1.396e-07 1.929 0.05459 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.726317e-10)
##
## Null deviance: 2.3041e-07 on 355 degrees of freedom
## Residual deviance: 5.8867e-08 on 341 degrees of freedom
## AIC: -6975.9
##
## Number of Fisher Scoring iterations: 2
Anova(wmod)
## Analysis of Deviance Table (Type II tests)
##
## Response: boxw
## LR Chisq Df Pr(>Chisq)
## treatment 8.82 4 0.0658725 .
## days 849.01 1 < 2.2e-16 ***
## whole.mean 109.30 1 < 2.2e-16 ***
## bees_alive 0.35 1 0.5536493
## qro 3.01 3 0.3904017
## treatment:days 19.98 4 0.0005034 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(wmod, test = "F")
## Single term deletions
##
## Model:
## boxw ~ treatment + days * treatment + whole.mean + bees_alive +
## qro + days
## Df Deviance AIC F value Pr(>F)
## <none> 5.8867e-08 -6975.9
## whole.mean 1 7.7735e-08 -6878.9 109.2954 < 2.2e-16 ***
## bees_alive 1 5.8928e-08 -6977.5 0.3508 0.5540418
## qro 3 5.9387e-08 -6978.7 1.0026 0.3917579
## treatment:days 4 6.2317e-08 -6963.6 4.9956 0.0006337 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(wmod)
wem <- emmeans(wmod, "treatment")
pairs(wem)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 -4.59e-06 2.21e-06 341 -2.080 0.2311
## treatment1 - treatment3 -2.75e-06 2.18e-06 341 -1.260 0.7161
## treatment1 - treatment4 -4.26e-06 2.23e-06 341 -1.910 0.3142
## treatment1 - treatment5 -5.85e-06 2.24e-06 341 -2.609 0.0710
## treatment2 - treatment3 1.84e-06 2.18e-06 341 0.846 0.9161
## treatment2 - treatment4 3.34e-07 2.26e-06 341 0.147 0.9999
## treatment2 - treatment5 -1.26e-06 2.27e-06 341 -0.553 0.9816
## treatment3 - treatment4 -1.51e-06 2.23e-06 341 -0.678 0.9611
## treatment3 - treatment5 -3.10e-06 2.25e-06 341 -1.375 0.6439
## treatment4 - treatment5 -1.59e-06 2.34e-06 341 -0.680 0.9607
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
Model includes time*treatment, reference https://rcompanion.org/handbook/I_09.html
Treatment does not influence colony weight change
library(nlme)
#Determine auto correction in residuals
model.a = gls(boxw ~ treatment + days + treatment*days + whole.mean + bees_alive + qro, data = w)
ACF(model.a, form = ~ days | colony)
## lag ACF
## 1 0 1.0000000
## 2 1 0.3306201
## 3 2 0.3692427
## 4 3 0.2809314
## 5 4 0.2588012
## 6 5 0.2215001
## 7 6 0.3561549
## 8 7 0.2098859
## 9 8 0.4522623
## 10 9 -0.2595828
## 11 10 0.8854066
model.b = gls(boxw ~ treatment + days + treatment*days + whole.mean, data = w)
ACF(model.a, form = ~ days | colony)
## lag ACF
## 1 0 1.0000000
## 2 1 0.3306201
## 3 2 0.3692427
## 4 3 0.2809314
## 5 4 0.2588012
## 6 5 0.2215001
## 7 6 0.3561549
## 8 7 0.2098859
## 9 8 0.4522623
## 10 9 -0.2595828
## 11 10 0.8854066
#make a couple mixed effect models, pick based on AIC
model = lme(boxw ~ treatment + days + treatment*days + whole.mean + bees_alive + qro, random = ~1|colony, correlation = corAR1(form = ~ days | colony, value = 0.331), data = w, method = "REML")
Anova(model)
## Analysis of Deviance Table (Type II tests)
##
## Response: boxw
## Chisq Df Pr(>Chisq)
## treatment 2.2200 4 0.6954
## days 478.7479 1 < 2.2e-16 ***
## whole.mean 29.9100 1 4.526e-08 ***
## bees_alive 2.2041 1 0.1376
## qro 0.6238 3 0.8910
## treatment:days 7.3462 4 0.1187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model1 = lme(boxw ~ treatment + days + treatment*days + whole.mean, random = ~1|colony, correlation = corAR1(form = ~ days | colony, value = 0.331), data = w, method = "REML")
Anova(model1)
## Analysis of Deviance Table (Type II tests)
##
## Response: boxw
## Chisq Df Pr(>Chisq)
## treatment 2.5557 4 0.63469
## days 496.1699 1 < 2.2e-16 ***
## whole.mean 36.1713 1 1.807e-09 ***
## treatment:days 8.2387 4 0.08321 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model, model1)
## df AIC
## model 18 -6823.474
## model1 14 -6921.841
#check if we want the random effect
model.fixed = gls(boxw ~ treatment + days + treatment*days + whole.mean, data = w, method = "REML")
anova(model1, model.fixed) # we do
## Model df AIC BIC logLik Test L.Ratio p-value
## model1 1 14 -6921.841 -6868.031 3474.920
## model.fixed 2 12 -6679.369 -6633.246 3351.684 1 vs 2 246.4723 <.0001
# post hoc analysis
library(multcompView)
library(lsmeans)
marginal = lsmeans(model1,
~treatment:days)
library(multcomp)
cld(marginal,
alpha = 0.05,
Letters = letters,
adjust = "tukey")
## treatment days lsmean SE df lower.CL upper.CL .group
## 1 21.3062 0.499826 3.29698e-06 39 0.499817 0.499835 a
## 3 21.3062 0.499829 3.31250e-06 39 0.499820 0.499838 a
## 2 21.3062 0.499830 3.34898e-06 39 0.499821 0.499839 a
## 4 21.3062 0.499830 3.34579e-06 39 0.499821 0.499839 a
## 5 21.3062 0.499833 3.38702e-06 39 0.499824 0.499842 a
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
# histogram of residuals
x = residuals(model1)
library(rcompanion)
plotNormalHistogram(x)
plot(model1)
qqnorm(resid(model1));qqline(resid(model1))
ggplot(w) +
aes( x = days, y = weight, color = treatment) +
geom_point(aes(color = treatment)) +
geom_smooth() +
scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ylab("Time (days since colony initiation") +
xlab("Weight (g)")+
ggtitle("Colony Weight Change Over Time (g)") +
theme_classic(base_size = 15)
colsum <- w %>%
group_by(treatment) %>%
summarise( mw = mean(weight),
sdw = sd(weight),
nw = length(weight)) %>%
mutate(sew = sdw/ sqrt(nw))
colsum
## # A tibble: 5 × 5
## treatment mw sdw nw sew
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 52.9 3.99 74 0.464
## 2 2 54.0 4.23 71 0.502
## 3 3 54.3 4.96 75 0.572
## 4 4 53.8 4.76 70 0.569
## 5 5 53.7 4.28 66 0.527
w47 <- subset(w, days<47)
ggplot(w47) +
aes( x = days, y = weight, color = treatment) +
geom_point(aes(color = treatment)) +
geom_smooth() +
scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ylab("Time (days since colony initiation") +
xlab("Weight (g) ")+
ggtitle("Colony Weight (g) Each Week") +
facet_grid(vars(treatment)) +
theme_classic(base_size = 12)
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 <- subset(brood, brood$round != 1)
brood <- merge(wd.summary, brood, by.x = "colony")
brood <- merge(brood, mortality, by.x = "colony")
b.gnb <- glm.nb(brood_cells ~ treatment + whole.mean + alive + qro + duration, data = brood)
drop1(b.gnb, test = "F")
## Single term deletions
##
## Model:
## brood_cells ~ treatment + whole.mean + alive + qro + duration
## Df Deviance AIC F value Pr(>F)
## <none> 54.366 345.60
## treatment 4 57.885 341.12 0.5502 0.70011
## whole.mean 1 151.249 440.48 60.5907 4.666e-09 ***
## alive 1 66.211 355.44 7.4080 0.01017 *
## qro 3 71.787 357.02 3.6319 0.02244 *
## duration 1 56.729 345.96 1.4783 0.23242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b.gnb1 <- glm.nb(brood_cells ~ treatment + whole.mean + alive + qro, data = brood)
AIC(b.gnb, b.gnb1)
## df AIC
## b.gnb 12 347.5991
## b.gnb1 11 347.9357
Anova(b.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## treatment 3.519 4 0.4749610
## whole.mean 96.884 1 < 2.2e-16 ***
## alive 11.845 1 0.0005781 ***
## qro 17.422 3 0.0005787 ***
## duration 2.364 1 0.1241867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(b.gnb, test = "F")
## Single term deletions
##
## Model:
## brood_cells ~ treatment + whole.mean + alive + qro + duration
## Df Deviance AIC F value Pr(>F)
## <none> 54.366 345.60
## treatment 4 57.885 341.12 0.5502 0.70011
## whole.mean 1 151.249 440.48 60.5907 4.666e-09 ***
## alive 1 66.211 355.44 7.4080 0.01017 *
## qro 3 71.787 357.02 3.6319 0.02244 *
## duration 1 56.729 345.96 1.4783 0.23242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(b.gnb)
brood.emm <- emmeans(b.gnb, "treatment")
pairs(brood.emm)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.06328 0.133 Inf -0.475 0.9896
## treatment1 - treatment3 -0.06441 0.130 Inf -0.496 0.9877
## treatment1 - treatment4 0.05004 0.131 Inf 0.382 0.9955
## treatment1 - treatment5 0.13947 0.142 Inf 0.985 0.8623
## treatment2 - treatment3 -0.00112 0.123 Inf -0.009 1.0000
## treatment2 - treatment4 0.11332 0.129 Inf 0.880 0.9045
## treatment2 - treatment5 0.20275 0.130 Inf 1.557 0.5251
## treatment3 - treatment4 0.11444 0.125 Inf 0.915 0.8912
## treatment3 - treatment5 0.20388 0.129 Inf 1.580 0.5102
## treatment4 - treatment5 0.08943 0.138 Inf 0.648 0.9671
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
summary(brood.emm)
## treatment emmean SE df asymp.LCL asymp.UCL
## 1 3.41 0.1010 Inf 3.21 3.61
## 2 3.47 0.0966 Inf 3.28 3.66
## 3 3.47 0.0920 Inf 3.29 3.65
## 4 3.36 0.1029 Inf 3.16 3.56
## 5 3.27 0.1024 Inf 3.07 3.47
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
sum_brood <- brood %>%
group_by(treatment) %>%
summarise( mean.b = mean(brood_cells),
n.b = length(brood_cells),
sd.b = sd(brood_cells)) %>%
mutate(seb = sd.b / sqrt(n.b))
sum_brood
## # A tibble: 5 × 5
## treatment mean.b n.b sd.b seb
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 33.8 9 22.6 7.53
## 2 2 36.9 9 11.2 3.74
## 3 3 45.6 9 26.2 8.73
## 4 4 36.7 9 18.3 6.09
## 5 5 29.6 9 17.5 5.82
ggplot(data = sum_brood, aes(x = treatment, y = mean.b, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(20,60)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = mean.b - seb,
ymax = mean.b + seb),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Brood Cellls",) +
ggtitle("Average of All 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)
hp.gnb <- glm.nb(honey_pot ~ treatment + whole.mean + alive + qro + duration, data = brood)
summary(hp.gnb)
##
## Call:
## glm.nb(formula = honey_pot ~ treatment + whole.mean + alive +
## qro + duration, data = brood, init.theta = 18.15285757, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4708 -1.0404 0.1104 0.4711 2.0842
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.001796 0.559486 0.003 0.9974
## treatment2 0.457373 0.245041 1.867 0.0620 .
## treatment3 0.021826 0.253828 0.086 0.9315
## treatment4 0.370384 0.240900 1.537 0.1242
## treatment5 0.231539 0.258995 0.894 0.3713
## whole.mean 1.208730 0.526496 2.296 0.0217 *
## alive 0.182876 0.076220 2.399 0.0164 *
## qroB3 -0.161323 0.256370 -0.629 0.5292
## qroB4 0.089812 0.278004 0.323 0.7466
## qroB5 0.204025 0.196376 1.039 0.2988
## duration 0.003274 0.011574 0.283 0.7773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(18.1529) family taken to be 1)
##
## Null deviance: 88.334 on 44 degrees of freedom
## Residual deviance: 54.400 on 34 degrees of freedom
## AIC: 243.39
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 18.2
## Std. Err.: 17.1
##
## 2 x log-likelihood: -219.385
plot(hp.gnb)
Anova(hp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pot
## LR Chisq Df Pr(>Chisq)
## treatment 6.3720 4 0.17304
## whole.mean 5.0936 1 0.02401 *
## alive 6.1060 1 0.01347 *
## qro 1.9423 3 0.58447
## duration 0.0791 1 0.77854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(hp.gnb, test = "F")
## Single term deletions
##
## Model:
## honey_pot ~ treatment + whole.mean + alive + qro + duration
## Df Deviance AIC F value Pr(>F)
## <none> 54.400 241.38
## treatment 4 60.772 239.76 0.9956 0.42331
## whole.mean 1 59.494 244.48 3.1835 0.08331 .
## alive 1 60.506 245.49 3.8162 0.05903 .
## qro 3 56.343 237.33 0.4046 0.75060
## duration 1 54.479 239.46 0.0494 0.82540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hp.gnb1 <- glm.nb(honey_pot ~ treatment + whole.mean + alive + duration, data = brood)
hp.gnb2 <- glm.nb(honey_pot ~ treatment + whole.mean + alive + qro , data = brood)
hp.gnb3 <- glm.nb(honey_pot ~ treatment + whole.mean + alive, data = brood)
anova(hp.gnb, hp.gnb1, hp.gnb2, hp.gnb3)
## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pot
## Model theta Resid. df
## 1 treatment + whole.mean + alive 17.17565 38
## 2 treatment + whole.mean + alive + duration 17.52117 37
## 3 treatment + whole.mean + alive + qro 17.58049 35
## 4 treatment + whole.mean + alive + qro + duration 18.15286 34
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -221.3473
## 2 -221.3263 1 vs 2 1 0.02101159 0.8847474
## 3 -219.4632 2 vs 3 2 1.86305773 0.3939510
## 4 -219.3853 3 vs 4 1 0.07794046 0.7801081
AIC(hp.gnb, hp.gnb1, hp.gnb2, hp.gnb3)
## df AIC
## hp.gnb 12 243.3853
## hp.gnb1 9 239.3263
## hp.gnb2 11 241.4632
## hp.gnb3 8 237.3473
Anova(hp.gnb3)
## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pot
## LR Chisq Df Pr(>Chisq)
## treatment 6.9398 4 0.139103
## whole.mean 9.7702 1 0.001774 **
## alive 5.5759 1 0.018209 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hp.emm <- emmeans(hp.gnb3, "treatment")
summary(hp.emm)
## treatment emmean SE df asymp.LCL asymp.UCL
## 1 1.50 0.184 Inf 1.14 1.86
## 2 1.99 0.146 Inf 1.70 2.27
## 3 1.54 0.169 Inf 1.21 1.87
## 4 1.87 0.155 Inf 1.57 2.17
## 5 1.78 0.162 Inf 1.47 2.10
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
pairs(hp.emm)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.4889 0.235 Inf -2.079 0.2292
## treatment1 - treatment3 -0.0422 0.250 Inf -0.169 0.9998
## treatment1 - treatment4 -0.3742 0.239 Inf -1.566 0.5195
## treatment1 - treatment5 -0.2879 0.247 Inf -1.164 0.7719
## treatment2 - treatment3 0.4467 0.219 Inf 2.035 0.2491
## treatment2 - treatment4 0.1147 0.210 Inf 0.546 0.9825
## treatment2 - treatment5 0.2010 0.216 Inf 0.931 0.8849
## treatment3 - treatment4 -0.3319 0.224 Inf -1.480 0.5758
## treatment3 - treatment5 -0.2457 0.231 Inf -1.061 0.8263
## treatment4 - treatment5 0.0862 0.224 Inf 0.384 0.9954
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
sum_hp <- brood %>%
group_by(treatment) %>%
summarise( mean.h = mean(honey_pot),
n.h = length(honey_pot),
sd.h = sd(honey_pot)) %>%
mutate(sehp = sd.h / sqrt(n.h))
sum_hp
## # A tibble: 5 × 5
## treatment mean.h n.h sd.h sehp
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 4.22 9 3.27 1.09
## 2 2 7.89 9 3.59 1.20
## 3 3 5.56 9 2.96 0.988
## 4 4 7 9 3.61 1.20
## 5 5 6.22 9 4.35 1.45
ggplot(data = sum_hp, aes(x = treatment, y = mean.h, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(2,10)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = mean.h - sehp,
ymax = mean.h + sehp),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Honey Pots",) +
ggtitle("Average of Honey Pots 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)
ll.gnb <- glm.nb(live_larvae ~ treatment + whole.mean + alive + qro + duration, data = brood)
summary(ll.gnb)
##
## Call:
## glm.nb(formula = live_larvae ~ treatment + whole.mean + alive +
## qro + duration, data = brood, init.theta = 1.845431183, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5413 -1.2555 -0.1094 0.3057 1.8630
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.186964 0.904280 -0.207 0.83620
## treatment2 -0.573712 0.401076 -1.430 0.15259
## treatment3 -0.071946 0.390249 -0.184 0.85373
## treatment4 -0.372603 0.392000 -0.951 0.34185
## treatment5 -0.238070 0.416208 -0.572 0.56732
## whole.mean 5.269608 0.880623 5.984 2.18e-09 ***
## alive 0.174995 0.115932 1.509 0.13118
## qroB3 -0.657878 0.446918 -1.472 0.14101
## qroB4 -0.413773 0.480121 -0.862 0.38879
## qroB5 1.000189 0.326062 3.067 0.00216 **
## duration -0.005932 0.019009 -0.312 0.75500
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8454) family taken to be 1)
##
## Null deviance: 111.962 on 44 degrees of freedom
## Residual deviance: 56.018 on 34 degrees of freedom
## AIC: 352.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.845
## Std. Err.: 0.506
##
## 2 x log-likelihood: -328.499
plot(ll.gnb)
Anova(ll.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 2.8019 4 0.591497
## whole.mean 30.0575 1 4.194e-08 ***
## alive 2.4880 1 0.114717
## qro 15.0416 3 0.001781 **
## duration 0.0996 1 0.752266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(ll.gnb, test = "F")
## Single term deletions
##
## Model:
## live_larvae ~ treatment + whole.mean + alive + qro + duration
## Df Deviance AIC F value Pr(>F)
## <none> 56.018 350.50
## treatment 4 58.820 345.30 0.4252 0.7893787
## whole.mean 1 86.076 378.56 18.2432 0.0001479 ***
## alive 1 58.506 350.99 1.5101 0.2275653
## qro 3 71.060 359.54 3.0431 0.0420193 *
## duration 1 56.118 348.60 0.0605 0.8072275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ll.sum <- brood %>%
group_by(treatment) %>%
summarise(ll.m= mean(live_larvae),
ll.n = length(live_larvae),
sd.ll = sd(live_larvae)) %>%
mutate(sell = sd.ll/ sqrt(ll.n))
ll.sum
## # A tibble: 5 × 5
## treatment ll.m ll.n sd.ll sell
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 26.7 9 26.6 8.88
## 2 2 13.8 9 9.72 3.24
## 3 3 31.2 9 23.1 7.70
## 4 4 17.6 9 14.3 4.75
## 5 5 16.4 9 13.4 4.47
ggplot(brood, aes(x=treatment, y = live_larvae, fill = treatment)) +
geom_boxplot(position = "identity",col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Live Larvae vs. Treatment") +
labs(y = "Count", x = "Treatment") +
theme_cowplot()
### Dead Larvae
dl.gnb <- glm.nb(dead_larvae ~ treatment + whole.mean + alive + qro + duration, data = brood)
summary(dl.gnb)
##
## Call:
## glm.nb(formula = dead_larvae ~ treatment + whole.mean + alive +
## qro + duration, data = brood, init.theta = 0.9840035759,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2046 -1.1468 -0.4028 0.4132 1.7539
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.939940 1.379407 -1.406 0.1596
## treatment2 0.225236 0.632830 0.356 0.7219
## treatment3 0.858769 0.597933 1.436 0.1509
## treatment4 0.636937 0.600881 1.060 0.2891
## treatment5 -0.486894 0.687222 -0.708 0.4786
## whole.mean 1.333415 1.322945 1.008 0.3135
## alive 0.420083 0.202293 2.077 0.0378 *
## qroB3 -0.003499 0.664646 -0.005 0.9958
## qroB4 1.282109 0.712712 1.799 0.0720 .
## qroB5 1.163773 0.505870 2.301 0.0214 *
## duration -0.005819 0.028251 -0.206 0.8368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.984) family taken to be 1)
##
## Null deviance: 78.573 on 44 degrees of freedom
## Residual deviance: 47.165 on 34 degrees of freedom
## AIC: 213.27
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.984
## Std. Err.: 0.321
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -189.269
plot(dl.gnb)
Anova(dl.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 5.4735 4 0.24207
## whole.mean 0.9967 1 0.31810
## alive 3.7838 1 0.05175 .
## qro 5.4077 3 0.14426
## duration 0.0414 1 0.83877
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(dl.gnb, test = "F")
## Single term deletions
##
## Model:
## dead_larvae ~ treatment + whole.mean + alive + qro + duration
## Df Deviance AIC F value Pr(>F)
## <none> 47.165 211.27
## treatment 4 52.639 208.74 0.9864 0.4281
## whole.mean 1 48.162 210.27 0.7185 0.4026
## alive 1 50.949 213.05 2.7276 0.1078
## qro 3 52.573 210.68 1.2994 0.2905
## duration 1 47.207 209.31 0.0298 0.8639
dl.sum <- brood %>%
group_by(treatment) %>%
summarise(dl.m= mean(dead_larvae),
dl.n = length(dead_larvae),
sd.dl = sd(dead_larvae)) %>%
mutate(sedl = sd.dl/ sqrt(dl.n))
dl.sum
## # A tibble: 5 × 5
## treatment dl.m dl.n sd.dl sedl
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 1.78 9 1.99 0.662
## 2 2 3.22 9 5.26 1.75
## 3 3 5.78 9 6.22 2.07
## 4 4 6.67 9 14.8 4.94
## 5 5 1.56 9 1.88 0.626
ggplot(brood, aes(x=treatment, y = dead_larvae, fill = treatment)) +
geom_boxplot(position = "identity",col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Dead Larvae vs. Treatment") +
labs(y = "Count", x = "Treatment") +
theme_cowplot()
brood$all_larvae <- brood$dead_larvae + brood$live_larvae
al.gnb <- glm.nb(all_larvae ~ treatment + whole.mean + alive + qro, data = brood)
summary(al.gnb)
##
## Call:
## glm.nb(formula = all_larvae ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 4.184750533, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3956 -0.9302 -0.1428 0.4582 2.1567
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.29312 0.38332 0.765 0.444468
## treatment2 -0.49453 0.26707 -1.852 0.064066 .
## treatment3 0.04672 0.26509 0.176 0.860098
## treatment4 -0.32028 0.26978 -1.187 0.235153
## treatment5 -0.31376 0.28092 -1.117 0.264044
## whole.mean 4.18091 0.58926 7.095 1.29e-12 ***
## alive 0.16181 0.07994 2.024 0.042940 *
## qroB3 -0.56940 0.30241 -1.883 0.059716 .
## qroB4 -0.07923 0.32694 -0.242 0.808526
## qroB5 0.84286 0.22470 3.751 0.000176 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1848) family taken to be 1)
##
## Null deviance: 157.42 on 44 degrees of freedom
## Residual deviance: 52.03 on 35 degrees of freedom
## AIC: 348.91
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.18
## Std. Err.: 1.20
##
## 2 x log-likelihood: -326.908
Anova(al.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: all_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 6.272 4 0.17976
## whole.mean 47.400 1 5.789e-12 ***
## alive 4.367 1 0.03664 *
## qro 22.395 3 5.397e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(al.gnb, test = "F")
## Single term deletions
##
## Model:
## all_larvae ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 52.030 346.91
## treatment 4 58.302 345.18 1.0547 0.393376
## whole.mean 1 99.430 392.31 31.8850 2.257e-06 ***
## alive 1 56.397 349.28 2.9376 0.095384 .
## qro 3 74.426 363.30 5.0217 0.005331 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
al.sum <- brood %>%
group_by(treatment) %>%
summarise(al.m= mean(all_larvae),
al.n = length(all_larvae),
sd.al = sd(all_larvae))
al.sum
## # A tibble: 5 × 4
## treatment al.m al.n sd.al
## <fct> <dbl> <int> <dbl>
## 1 1 28.4 9 27.0
## 2 2 17 9 12.5
## 3 3 37 9 26.3
## 4 4 24.2 9 23.8
## 5 5 18 9 13.6
larvae.mortality <- glm(cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + qro, data = brood, family = binomial("logit"))
summary(larvae.mortality)
##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ treatment + whole.mean +
## qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.603 -1.646 0.000 1.690 4.054
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5754 0.4250 8.412 < 2e-16 ***
## treatment2 -1.3544 0.3359 -4.032 5.52e-05 ***
## treatment3 -1.0149 0.3083 -3.292 0.000994 ***
## treatment4 -1.7935 0.3205 -5.596 2.19e-08 ***
## treatment5 -0.4505 0.3853 -1.169 0.242319
## whole.mean -0.7257 0.6306 -1.151 0.249779
## qroB3 -0.4609 0.3392 -1.359 0.174302
## qroB4 -0.5556 0.3188 -1.743 0.081415 .
## qroB5 -0.6609 0.2122 -3.114 0.001843 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 329.89 on 42 degrees of freedom
## Residual deviance: 263.65 on 34 degrees of freedom
## AIC: 352.64
##
## Number of Fisher Scoring iterations: 5
Anova(larvae.mortality)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## treatment 43.321 4 8.874e-09 ***
## whole.mean 1.347 1 0.24585
## qro 10.093 3 0.01779 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(larvae.mortality, test = "F")
## Single term deletions
##
## Model:
## cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + qro
## Df Deviance AIC F value Pr(>F)
## <none> 263.65 352.64
## treatment 4 306.97 387.96 1.3967 0.2560
## whole.mean 1 265.00 351.99 0.1737 0.6795
## qro 3 273.74 356.74 0.4339 0.7301
plot(larvae.mortality)
lm.emm <- emmeans(larvae.mortality, "treatment")
pairs(lm.emm)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 1.354 0.336 Inf 4.032 0.0005
## treatment1 - treatment3 1.015 0.308 Inf 3.292 0.0088
## treatment1 - treatment4 1.793 0.320 Inf 5.596 <.0001
## treatment1 - treatment5 0.450 0.385 Inf 1.169 0.7690
## treatment2 - treatment3 -0.339 0.270 Inf -1.258 0.7170
## treatment2 - treatment4 0.439 0.295 Inf 1.488 0.5705
## treatment2 - treatment5 -0.904 0.351 Inf -2.574 0.0753
## treatment3 - treatment4 0.779 0.244 Inf 3.187 0.0125
## treatment3 - treatment5 -0.564 0.335 Inf -1.684 0.4436
## treatment4 - treatment5 -1.343 0.356 Inf -3.777 0.0015
##
## Results are averaged over the levels of: qro
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
summary(lm.emm)
## treatment emmean SE df asymp.LCL asymp.UCL
## 1 2.81 0.281 Inf 2.256 3.36
## 2 1.45 0.230 Inf 1.002 1.90
## 3 1.79 0.193 Inf 1.414 2.17
## 4 1.01 0.255 Inf 0.514 1.51
## 5 2.36 0.294 Inf 1.780 2.93
##
## Results are averaged over the levels of: qro
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
emmeans(larvae.mortality, pairwise ~ treatment, type = "response")
## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.943 0.0151 Inf 0.905 0.966
## 2 0.810 0.0353 Inf 0.732 0.870
## 3 0.857 0.0236 Inf 0.804 0.898
## 4 0.734 0.0498 Inf 0.626 0.820
## 5 0.913 0.0232 Inf 0.856 0.949
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.874 1.3014 Inf 1 4.032 0.0005
## treatment1 / treatment3 2.759 0.8506 Inf 1 3.292 0.0088
## treatment1 / treatment4 6.010 1.9263 Inf 1 5.596 <.0001
## treatment1 / treatment5 1.569 0.6045 Inf 1 1.169 0.7690
## treatment2 / treatment3 0.712 0.1922 Inf 1 -1.258 0.7170
## treatment2 / treatment4 1.551 0.4579 Inf 1 1.488 0.5705
## treatment2 / treatment5 0.405 0.1422 Inf 1 -2.574 0.0753
## treatment3 / treatment4 2.178 0.5320 Inf 1 3.187 0.0125
## treatment3 / treatment5 0.569 0.1906 Inf 1 -1.684 0.4436
## treatment4 / treatment5 0.261 0.0928 Inf 1 -3.777 0.0015
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emm1 <- emmeans(larvae.mortality, ~ treatment, type = "response")
emm1
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.943 0.0151 Inf 0.905 0.966
## 2 0.810 0.0353 Inf 0.732 0.870
## 3 0.857 0.0236 Inf 0.804 0.898
## 4 0.734 0.0498 Inf 0.626 0.820
## 5 0.913 0.0232 Inf 0.856 0.949
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
pairs(emm1)
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.874 1.3014 Inf 1 4.032 0.0005
## treatment1 / treatment3 2.759 0.8506 Inf 1 3.292 0.0088
## treatment1 / treatment4 6.010 1.9263 Inf 1 5.596 <.0001
## treatment1 / treatment5 1.569 0.6045 Inf 1 1.169 0.7690
## treatment2 / treatment3 0.712 0.1922 Inf 1 -1.258 0.7170
## treatment2 / treatment4 1.551 0.4579 Inf 1 1.488 0.5705
## treatment2 / treatment5 0.405 0.1422 Inf 1 -2.574 0.0753
## treatment3 / treatment4 2.178 0.5320 Inf 1 3.187 0.0125
## treatment3 / treatment5 0.569 0.1906 Inf 1 -1.684 0.4436
## treatment4 / treatment5 0.261 0.0928 Inf 1 -3.777 0.0015
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
larvcld <- cld(object = emm1,
adjust = "TUkey",
alpha = 0.05,
Letters = letters)
larvcld
## treatment prob SE df asymp.LCL asymp.UCL .group
## 4 0.734 0.0498 Inf 0.589 0.841 a
## 2 0.810 0.0353 Inf 0.703 0.885 ab
## 3 0.857 0.0236 Inf 0.785 0.908 b
## 5 0.913 0.0232 Inf 0.832 0.957 bc
## 1 0.943 0.0151 Inf 0.889 0.972 c
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
emmdf <- as.data.frame(emm1)
p <- ggplot(data = emmdf, aes(x=treatment, y=prob, fill=treatment)) +
geom_col(position = "dodge", color = "black") +
coord_cartesian(ylim = c(0.5,1.06))
p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
labs(y = "Probability of Survival",) +
ggtitle("Probability of larvae being alive upon colony dissection") +
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)+
annotate(geom = "text",
x = 1.5, y = 1.05,
label = "Anova = 8.874e-09 ***",
size = 5) +
annotate(geom = "text",
x = c(4, 2, 3, 5, 1),
y = c(0.85, 0.89, 0.92, 0.97, 0.99),
label = c("a", "ab", "b", "bc", "c"),
size = 5)
lp.gnb <- glm.nb(live_pupae ~ treatment + whole.mean + alive + qro, data = brood)
summary(lp.gnb)
##
## Call:
## glm.nb(formula = live_pupae ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 1.349988185, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6035 -1.2236 -0.1603 0.4227 1.4934
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.63174 0.68918 -0.917 0.359321
## treatment2 0.09101 0.48104 0.189 0.849935
## treatment3 -0.35995 0.50305 -0.716 0.474276
## treatment4 -0.46695 0.50718 -0.921 0.357210
## treatment5 -0.39231 0.52658 -0.745 0.456264
## whole.mean 4.00378 1.10653 3.618 0.000297 ***
## alive 0.03326 0.14411 0.231 0.817447
## qroB3 -0.51721 0.55444 -0.933 0.350896
## qroB4 -0.39441 0.60398 -0.653 0.513739
## qroB5 0.31255 0.42525 0.735 0.462362
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.35) family taken to be 1)
##
## Null deviance: 71.716 on 44 degrees of freedom
## Residual deviance: 51.808 on 35 degrees of freedom
## AIC: 231.53
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.350
## Std. Err.: 0.470
##
## 2 x log-likelihood: -209.526
Anova(lp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: live_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 1.9477 4 0.7453869
## whole.mean 11.2243 1 0.0008073 ***
## alive 0.0517 1 0.8201473
## qro 2.1706 3 0.5377697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(lp.gnb, test = "F")
## Single term deletions
##
## Model:
## live_pupae ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 51.808 229.53
## treatment 4 53.756 223.47 0.3289 0.856621
## whole.mean 1 63.032 238.75 7.5828 0.009281 **
## alive 1 51.860 227.58 0.0349 0.852841
## qro 3 53.979 225.70 0.4888 0.692302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lp.sum <- brood %>%
group_by(treatment) %>%
summarise(lp.m= mean(live_pupae),
sd.lp = sd(live_pupae))
lp.sum
## # A tibble: 5 × 3
## treatment lp.m sd.lp
## <fct> <dbl> <dbl>
## 1 1 4.78 4.60
## 2 2 5.56 5.00
## 3 3 4.11 5.75
## 4 4 3 2.87
## 5 5 2.89 2.98
dp.gnb <- glm.nb(dead_pupae ~ treatment + whole.mean + alive + qro, data = brood) # stick with this model overall
summary(dp.gnb)
##
## Call:
## glm.nb(formula = dead_pupae ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 2.150708229, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2088 -1.0582 -0.3397 0.4293 2.5982
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.51927 0.73174 -2.076 0.037871 *
## treatment2 0.98992 0.44291 2.235 0.025415 *
## treatment3 1.10607 0.44341 2.494 0.012616 *
## treatment4 -0.01406 0.48434 -0.029 0.976840
## treatment5 0.51263 0.47812 1.072 0.283647
## whole.mean 3.36553 0.96271 3.496 0.000472 ***
## alive 0.19137 0.14866 1.287 0.197994
## qroB3 -0.56450 0.47315 -1.193 0.232837
## qroB4 -0.31027 0.52545 -0.590 0.554860
## qroB5 -0.44380 0.39058 -1.136 0.255849
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.1507) family taken to be 1)
##
## Null deviance: 99.516 on 44 degrees of freedom
## Residual deviance: 50.381 on 35 degrees of freedom
## AIC: 234.06
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.151
## Std. Err.: 0.758
##
## 2 x log-likelihood: -212.065
Anova(dp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 12.0782 4 0.0167790 *
## whole.mean 12.5620 1 0.0003937 ***
## alive 1.5575 1 0.2120246
## qro 1.9988 3 0.5726506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(dp.gnb, test = "F")
## Single term deletions
##
## Model:
## dead_pupae ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 50.381 232.06
## treatment 4 62.459 236.14 2.0977 0.102004
## whole.mean 1 62.943 242.63 8.7269 0.005574 **
## alive 1 51.939 231.62 1.0820 0.305376
## qro 3 52.380 228.06 0.4629 0.710015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
em.dp <- emmeans(dp.gnb, "treatment")
summary(em.dp)
## treatment emmean SE df asymp.LCL asymp.UCL
## 1 0.555 0.372 Inf -0.173 1.28
## 2 1.545 0.300 Inf 0.958 2.13
## 3 1.661 0.292 Inf 1.088 2.23
## 4 0.541 0.378 Inf -0.199 1.28
## 5 1.068 0.328 Inf 0.424 1.71
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
pairs(em.dp)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.9899 0.443 Inf -2.235 0.1668
## treatment1 - treatment3 -1.1061 0.443 Inf -2.494 0.0918
## treatment1 - treatment4 0.0141 0.484 Inf 0.029 1.0000
## treatment1 - treatment5 -0.5126 0.478 Inf -1.072 0.8209
## treatment2 - treatment3 -0.1162 0.377 Inf -0.308 0.9980
## treatment2 - treatment4 1.0040 0.431 Inf 2.329 0.1356
## treatment2 - treatment5 0.4773 0.411 Inf 1.161 0.7734
## treatment3 - treatment4 1.1201 0.426 Inf 2.632 0.0646
## treatment3 - treatment5 0.5934 0.406 Inf 1.463 0.5866
## treatment4 - treatment5 -0.5267 0.470 Inf -1.121 0.7954
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
dp.sum <- brood %>%
group_by(treatment) %>%
summarise(dp.m= mean(dead_pupae),
sd.dp = sd(dead_pupae),
n.dp = length(dead_pupae)) %>%
mutate(sedp = sd.dp/ sqrt(n.dp))
dp.sum
## # A tibble: 5 × 5
## treatment dp.m sd.dp n.dp sedp
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 2 2.06 9 0.687
## 2 2 7.89 12.1 9 4.04
## 3 3 8.89 7.27 9 2.42
## 4 4 2.89 3.02 9 1.01
## 5 5 3.89 4.28 9 1.43
ggplot(data = dp.sum, aes(x = treatment, y = dp.m, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0,14)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = dp.m - sedp,
ymax = dp.m + sedp),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Dead Pupae",) +
ggtitle("Average of Dead Pupae 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)
brood$all_pupae <- brood$dead_pupae + brood$live_pupae
ap.gnb <- glm.nb(all_pupae ~ treatment + whole.mean + alive + qro, data = brood)
summary(ap.gnb)
##
## Call:
## glm.nb(formula = all_pupae ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 3.199279779, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4442 -0.8742 -0.1522 0.4555 2.2193
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.11626 0.47735 -0.244 0.808
## treatment2 0.54919 0.32393 1.695 0.090 .
## treatment3 0.37121 0.33224 1.117 0.264
## treatment4 -0.32605 0.34898 -0.934 0.350
## treatment5 0.03050 0.35262 0.086 0.931
## whole.mean 3.78765 0.73353 5.164 2.42e-07 ***
## alive 0.05939 0.09840 0.604 0.546
## qroB3 -0.50493 0.36204 -1.395 0.163
## qroB4 -0.44186 0.40153 -1.100 0.271
## qroB5 -0.09061 0.28764 -0.315 0.753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.1993) family taken to be 1)
##
## Null deviance: 104.241 on 44 degrees of freedom
## Residual deviance: 51.328 on 35 degrees of freedom
## AIC: 281.34
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.20
## Std. Err.: 1.04
##
## 2 x log-likelihood: -259.344
Anova(ap.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: all_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 8.9012 4 0.06362 .
## whole.mean 27.4182 1 1.639e-07 ***
## alive 0.3515 1 0.55328
## qro 2.6716 3 0.44508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(ap.gnb, test = "F")
## Single term deletions
##
## Model:
## all_pupae ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 51.328 279.34
## treatment 4 60.229 280.25 1.5174 0.2183889
## whole.mean 1 78.747 304.76 18.6961 0.0001211 ***
## alive 1 51.680 277.70 0.2397 0.6274993
## qro 3 54.000 276.02 0.6072 0.6147503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ap.sum <- brood %>%
group_by(treatment) %>%
summarise(al.m= mean(all_pupae),
al.n = length(all_pupae),
sd.al = sd(all_pupae)) %>%
mutate(sep = sd.al/sqrt(al.n))
ap.sum
## # A tibble: 5 × 5
## treatment al.m al.n sd.al sep
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 6.78 9 6.22 2.07
## 2 2 13.4 9 14.1 4.68
## 3 3 13 9 9.5 3.17
## 4 4 5.89 9 5.16 1.72
## 5 5 6.78 9 5.89 1.96
ggplot(data = ap.sum, aes(x = treatment, y = al.m, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0,20)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = al.m - sep,
ymax = al.m + sep),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean of all Pupae",) +
ggtitle("Average of All (dead and alive) Pupae 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)
pupae.mortality <- glm(cbind(live_pupae, dead_pupae) ~ treatment + whole.mean + qro, data = brood, family = binomial("logit"))
summary(pupae.mortality)
##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ treatment + whole.mean +
## qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.507 -1.388 0.374 1.075 3.360
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1605 0.5381 2.157 0.031025 *
## treatment2 -1.3652 0.3537 -3.859 0.000114 ***
## treatment3 -1.7369 0.3605 -4.818 1.45e-06 ***
## treatment4 -0.8751 0.4075 -2.148 0.031750 *
## treatment5 -1.6576 0.4124 -4.019 5.84e-05 ***
## whole.mean -0.6907 0.8159 -0.847 0.397251
## qroB3 -0.5694 0.4330 -1.315 0.188433
## qroB4 0.9133 0.3583 2.549 0.010805 *
## qroB5 0.8723 0.3042 2.868 0.004135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.65 on 40 degrees of freedom
## Residual deviance: 127.41 on 32 degrees of freedom
## AIC: 215.17
##
## Number of Fisher Scoring iterations: 4
Anova(pupae.mortality)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## treatment 29.8204 4 5.324e-06 ***
## whole.mean 0.7204 1 0.396014
## qro 18.8401 3 0.000295 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(pupae.mortality, test = "F")
## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ treatment + whole.mean + qro
## Df Deviance AIC F value Pr(>F)
## <none> 127.41 215.18
## treatment 4 157.23 237.00 1.8724 0.1395
## whole.mean 1 128.13 213.90 0.1809 0.6734
## qro 3 146.25 228.01 1.5773 0.2140
plot(pupae.mortality)
lm.emm <- emmeans(pupae.mortality, "treatment")
pairs(lm.emm)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 1.3652 0.354 Inf 3.859 0.0011
## treatment1 - treatment3 1.7369 0.361 Inf 4.818 <.0001
## treatment1 - treatment4 0.8751 0.407 Inf 2.148 0.2001
## treatment1 - treatment5 1.6576 0.412 Inf 4.019 0.0006
## treatment2 - treatment3 0.3717 0.289 Inf 1.286 0.6999
## treatment2 - treatment4 -0.4900 0.346 Inf -1.417 0.6165
## treatment2 - treatment5 0.2924 0.340 Inf 0.860 0.9113
## treatment3 - treatment4 -0.8618 0.360 Inf -2.393 0.1172
## treatment3 - treatment5 -0.0794 0.352 Inf -0.225 0.9994
## treatment4 - treatment5 0.7824 0.408 Inf 1.916 0.3085
##
## Results are averaged over the levels of: qro
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
summary(lm.emm)
## treatment emmean SE df asymp.LCL asymp.UCL
## 1 1.133 0.319 Inf 0.507 1.7587
## 2 -0.233 0.243 Inf -0.708 0.2432
## 3 -0.604 0.248 Inf -1.091 -0.1176
## 4 0.258 0.341 Inf -0.411 0.9259
## 5 -0.525 0.293 Inf -1.099 0.0495
##
## Results are averaged over the levels of: qro
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
plot(brood$treatment, cbind(brood$dead_pupae, brood$live_pupae))
emmeans(pupae.mortality, pairwise ~ treatment, type = "response")
## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.756 0.0589 Inf 0.624 0.853
## 2 0.442 0.0599 Inf 0.330 0.561
## 3 0.353 0.0567 Inf 0.251 0.471
## 4 0.564 0.0839 Inf 0.399 0.716
## 5 0.372 0.0685 Inf 0.250 0.512
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.916 1.385 Inf 1 3.859 0.0011
## treatment1 / treatment3 5.680 2.048 Inf 1 4.818 <.0001
## treatment1 / treatment4 2.399 0.978 Inf 1 2.148 0.2001
## treatment1 / treatment5 5.246 2.164 Inf 1 4.019 0.0006
## treatment2 / treatment3 1.450 0.419 Inf 1 1.286 0.6999
## treatment2 / treatment4 0.613 0.212 Inf 1 -1.417 0.6165
## treatment2 / treatment5 1.340 0.455 Inf 1 0.860 0.9113
## treatment3 / treatment4 0.422 0.152 Inf 1 -2.393 0.1172
## treatment3 / treatment5 0.924 0.325 Inf 1 -0.225 0.9994
## treatment4 / treatment5 2.187 0.893 Inf 1 1.916 0.3085
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emm2 <- emmeans(pupae.mortality, ~ treatment, type = "response")
emm2
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.756 0.0589 Inf 0.624 0.853
## 2 0.442 0.0599 Inf 0.330 0.561
## 3 0.353 0.0567 Inf 0.251 0.471
## 4 0.564 0.0839 Inf 0.399 0.716
## 5 0.372 0.0685 Inf 0.250 0.512
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
pairs(emm2)
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.916 1.385 Inf 1 3.859 0.0011
## treatment1 / treatment3 5.680 2.048 Inf 1 4.818 <.0001
## treatment1 / treatment4 2.399 0.978 Inf 1 2.148 0.2001
## treatment1 / treatment5 5.246 2.164 Inf 1 4.019 0.0006
## treatment2 / treatment3 1.450 0.419 Inf 1 1.286 0.6999
## treatment2 / treatment4 0.613 0.212 Inf 1 -1.417 0.6165
## treatment2 / treatment5 1.340 0.455 Inf 1 0.860 0.9113
## treatment3 / treatment4 0.422 0.152 Inf 1 -2.393 0.1172
## treatment3 / treatment5 0.924 0.325 Inf 1 -0.225 0.9994
## treatment4 / treatment5 2.187 0.893 Inf 1 1.916 0.3085
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
cldpup <- cld(object = emm2,
adjust = "Tukey",
alpha = 0.05,
Letters =letters)
cldpup
## treatment prob SE df asymp.LCL asymp.UCL .group
## 3 0.353 0.0567 Inf 0.224 0.508 a
## 5 0.372 0.0685 Inf 0.218 0.557 a
## 2 0.442 0.0599 Inf 0.298 0.597 a
## 4 0.564 0.0839 Inf 0.350 0.756 ab
## 1 0.756 0.0589 Inf 0.577 0.876 b
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
emmdf2 <- as.data.frame(emm2)
p <- ggplot(data = emmdf2, aes(x=treatment, y=prob, fill=treatment)) +
geom_col(position = "dodge", color = "black") +
coord_cartesian(ylim = c(0,1))
p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
labs(y = "Probability of Survival",) +
ggtitle("Probability of pupae being alive upon colony dissection") +
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)+
annotate(geom = "text",
x = 2, y = 1,
label = "Anova = 5.324e-06***",
size = 5) +
annotate(geom = "text",
x = c(3, 5, 2, 4, 1),
y = c(0.53, 0.58, 0.62, 0.76, 0.9),
label = c("a", "a", "a", "ab", "b"),
size = 5)
brood$all_d_lp <- brood$dead_larvae + brood$dead_pupae
dlp.gnb <- glm.nb(all_d_lp ~ treatment + whole.mean + alive + qro, data = brood)
summary(dlp.gnb)
##
## Call:
## glm.nb(formula = all_d_lp ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 2.163558771, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3985 -0.8771 -0.2727 0.4112 2.2170
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.20339 0.63228 -1.903 0.05701 .
## treatment2 0.67384 0.39484 1.707 0.08789 .
## treatment3 0.85259 0.39457 2.161 0.03071 *
## treatment4 0.33275 0.40703 0.818 0.41363
## treatment5 -0.02301 0.43036 -0.053 0.95737
## whole.mean 2.62915 0.85547 3.073 0.00212 **
## alive 0.31627 0.13026 2.428 0.01519 *
## qroB3 -0.43342 0.42969 -1.009 0.31313
## qroB4 0.35285 0.47649 0.741 0.45898
## qroB5 0.35311 0.33836 1.044 0.29667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.1636) family taken to be 1)
##
## Null deviance: 101.542 on 44 degrees of freedom
## Residual deviance: 47.974 on 35 degrees of freedom
## AIC: 277.77
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.164
## Std. Err.: 0.622
##
## 2 x log-likelihood: -255.767
Anova(dlp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: all_d_lp
## LR Chisq Df Pr(>Chisq)
## treatment 7.7590 4 0.100815
## whole.mean 10.0531 1 0.001521 **
## alive 5.4642 1 0.019410 *
## qro 2.8749 3 0.411310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(dlp.gnb, test = "F")
## Single term deletions
##
## Model:
## all_d_lp ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 47.974 275.77
## treatment 4 55.733 275.53 1.4152 0.24930
## whole.mean 1 58.027 283.82 7.3344 0.01040 *
## alive 1 53.438 279.23 3.9865 0.05369 .
## qro 3 50.849 272.64 0.6992 0.55890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dlp.sum <- brood %>%
group_by(treatment) %>%
summarise(al.m= mean(all_d_lp),
al.n = length(all_d_lp),
sd.al = sd(all_d_lp))
dlp.sum
## # A tibble: 5 × 4
## treatment al.m al.n sd.al
## <fct> <dbl> <int> <dbl>
## 1 1 3.78 9 2.99
## 2 2 11.1 9 13.2
## 3 3 14.7 9 12.6
## 4 4 9.56 9 16.1
## 5 5 5.44 9 5.77
brood$all_a_lp <- brood$live_larvae + brood$live_pupae
alp.gnb <- glm.nb(all_a_lp ~ treatment + whole.mean + alive + qro, data = brood)
summary(alp.gnb)
##
## Call:
## glm.nb(formula = all_a_lp ~ treatment + whole.mean + alive +
## qro, data = brood, init.theta = 2.19074725, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7689 -1.0334 -0.1992 0.3133 2.0300
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.08091 0.49253 -0.164 0.86951
## treatment2 -0.43961 0.34943 -1.258 0.20837
## treatment3 -0.12495 0.35289 -0.354 0.72329
## treatment4 -0.40621 0.35671 -1.139 0.25480
## treatment5 -0.31110 0.36802 -0.845 0.39793
## whole.mean 5.07445 0.78225 6.487 8.76e-11 ***
## alive 0.16567 0.10281 1.611 0.10711
## qroB3 -0.61995 0.38980 -1.590 0.11174
## qroB4 -0.42167 0.43836 -0.962 0.33609
## qroB5 0.93593 0.29749 3.146 0.00165 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.1907) family taken to be 1)
##
## Null deviance: 116.061 on 44 degrees of freedom
## Residual deviance: 55.473 on 35 degrees of freedom
## AIC: 365.13
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.191
## Std. Err.: 0.585
##
## 2 x log-likelihood: -343.126
Anova(alp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: all_a_lp
## LR Chisq Df Pr(>Chisq)
## treatment 2.358 4 0.670151
## whole.mean 34.947 1 3.388e-09 ***
## alive 2.677 1 0.101823
## qro 15.801 3 0.001245 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(alp.gnb, test = "F")
## Single term deletions
##
## Model:
## all_a_lp ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 55.473 363.13
## treatment 4 57.831 357.48 0.3720 0.82696
## whole.mean 1 90.420 396.07 22.0492 4.009e-05 ***
## alive 1 58.150 363.80 1.6889 0.20225
## qro 3 71.274 372.93 3.3232 0.03074 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alp.sum <- brood %>%
group_by(treatment) %>%
summarise(al.m= mean(all_a_lp),
al.n = length(all_a_lp),
sd.al = sd(all_a_lp))
alp.sum
## # A tibble: 5 × 4
## treatment al.m al.n sd.al
## <fct> <dbl> <int> <dbl>
## 1 1 31.4 9 29.9
## 2 2 19.3 9 13.9
## 3 3 35.3 9 26.6
## 4 4 20.6 9 16.2
## 5 5 19.3 9 14.3
brood$all_lp <- brood$all_larvae + brood$all_pupae
lp.gnb <- glm.nb(all_lp ~ treatment + whole.mean + alive + qro, data = brood)
summary(lp.gnb)
##
## Call:
## glm.nb(formula = all_lp ~ treatment + whole.mean + alive + qro,
## data = brood, init.theta = 6.342196108, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.54151 -1.04455 -0.09369 0.40579 2.25119
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.74733 0.31286 2.389 0.016910 *
## treatment2 -0.17512 0.21576 -0.812 0.416993
## treatment3 0.09679 0.21687 0.446 0.655394
## treatment4 -0.36059 0.22204 -1.624 0.104383
## treatment5 -0.30238 0.23008 -1.314 0.188758
## whole.mean 4.06629 0.48121 8.450 < 2e-16 ***
## alive 0.14930 0.06503 2.296 0.021685 *
## qroB3 -0.56419 0.24420 -2.310 0.020869 *
## qroB4 -0.17260 0.26721 -0.646 0.518323
## qroB5 0.62756 0.18429 3.405 0.000661 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(6.3422) family taken to be 1)
##
## Null deviance: 188.770 on 44 degrees of freedom
## Residual deviance: 52.539 on 35 degrees of freedom
## AIC: 366.6
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 6.34
## Std. Err.: 1.80
##
## 2 x log-likelihood: -344.602
Anova(lp.gnb)
## Analysis of Deviance Table (Type II tests)
##
## Response: all_lp
## LR Chisq Df Pr(>Chisq)
## treatment 6.596 4 0.15882
## whole.mean 69.633 1 < 2.2e-16 ***
## alive 5.516 1 0.01884 *
## qro 22.331 3 5.566e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(lp.gnb, test = "F")
## Single term deletions
##
## Model:
## all_lp ~ treatment + whole.mean + alive + qro
## Df Deviance AIC F value Pr(>F)
## <none> 52.539 364.60
## treatment 4 59.135 363.20 1.0986 0.372591
## whole.mean 1 122.172 432.23 46.3872 6.712e-08 ***
## alive 1 58.055 368.12 3.6746 0.063442 .
## qro 3 74.870 380.93 4.9588 0.005673 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(brood$treatment, brood$all_lp)
lp.sum <- brood %>%
group_by(treatment) %>%
summarise(lp.m= mean(all_lp),
lp.n = length(all_lp),
lp.al = sd(all_lp))
lp.sum
## # A tibble: 5 × 4
## treatment lp.m lp.n lp.al
## <fct> <dbl> <int> <dbl>
## 1 1 35.2 9 30.6
## 2 2 30.4 9 19.4
## 3 3 50 9 34.6
## 4 4 30.1 9 26.8
## 5 5 24.8 9 17.1
em.lp <- emmeans(lp.gnb, "treatment")
pairs(em.lp)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.1751 0.216 Inf 0.812 0.9272
## treatment1 - treatment3 -0.0968 0.217 Inf -0.446 0.9918
## treatment1 - treatment4 0.3606 0.222 Inf 1.624 0.4819
## treatment1 - treatment5 0.3024 0.230 Inf 1.314 0.6823
## treatment2 - treatment3 -0.2719 0.208 Inf -1.306 0.6873
## treatment2 - treatment4 0.1855 0.217 Inf 0.854 0.9134
## treatment2 - treatment5 0.1273 0.221 Inf 0.577 0.9785
## treatment3 - treatment4 0.4574 0.215 Inf 2.126 0.2088
## treatment3 - treatment5 0.3992 0.216 Inf 1.850 0.3446
## treatment4 - treatment5 -0.0582 0.231 Inf -0.252 0.9991
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
pl_bind <- glm(cbind(all_a_lp, all_d_lp) ~ treatment + whole.mean + qro, data = brood, family = binomial("logit"))
summary(pl_bind )
##
## Call:
## glm(formula = cbind(all_a_lp, all_d_lp) ~ treatment + whole.mean +
## qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.653 -1.318 0.000 2.193 4.821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.5891 0.3105 8.337 < 2e-16 ***
## treatment2 -1.6031 0.2240 -7.157 8.22e-13 ***
## treatment3 -1.1485 0.2136 -5.377 7.58e-08 ***
## treatment4 -1.2431 0.2302 -5.401 6.62e-08 ***
## treatment5 -0.9699 0.2474 -3.921 8.82e-05 ***
## whole.mean -0.9760 0.4614 -2.115 0.0344 *
## qroB3 -0.1974 0.2323 -0.850 0.3954
## qroB4 0.3607 0.2132 1.692 0.0906 .
## qroB5 0.1809 0.1510 1.198 0.2310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 411.27 on 42 degrees of freedom
## Residual deviance: 331.16 on 34 degrees of freedom
## AIC: 470.26
##
## Number of Fisher Scoring iterations: 5
Anova(pl_bind )
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(all_a_lp, all_d_lp)
## LR Chisq Df Pr(>Chisq)
## treatment 63.817 4 4.567e-13 ***
## whole.mean 4.553 1 0.03286 *
## qro 6.895 3 0.07531 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(pl_bind, test = "F")
## Single term deletions
##
## Model:
## cbind(all_a_lp, all_d_lp) ~ treatment + whole.mean + qro
## Df Deviance AIC F value Pr(>F)
## <none> 331.16 470.26
## treatment 4 394.98 526.08 1.6380 0.1873
## whole.mean 1 335.72 472.81 0.4674 0.4988
## qro 3 338.06 471.15 0.2360 0.8707
plot(pl_bind )
lm.emm <- emmeans(pl_bind, pairwise ~ treatment, type = "response")
pairs(lm.emm)
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 4.969 1.113 Inf 1 7.157 <.0001
## treatment1 / treatment3 3.154 0.674 Inf 1 5.377 <.0001
## treatment1 / treatment4 3.466 0.798 Inf 1 5.401 <.0001
## treatment1 / treatment5 2.638 0.652 Inf 1 3.921 0.0008
## treatment2 / treatment3 0.635 0.109 Inf 1 -2.644 0.0626
## treatment2 / treatment4 0.698 0.137 Inf 1 -1.832 0.3552
## treatment2 / treatment5 0.531 0.110 Inf 1 -3.048 0.0195
## treatment3 / treatment4 1.099 0.199 Inf 1 0.524 0.9850
## treatment3 / treatment5 0.836 0.172 Inf 1 -0.870 0.9078
## treatment4 / treatment5 0.761 0.174 Inf 1 -1.198 0.7526
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
summary(lm.emm)
## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.901 0.0180 Inf 0.860 0.931
## 2 0.646 0.0345 Inf 0.576 0.711
## 3 0.742 0.0262 Inf 0.688 0.790
## 4 0.724 0.0379 Inf 0.644 0.792
## 5 0.775 0.0305 Inf 0.710 0.829
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 4.969 1.113 Inf 1 7.157 <.0001
## treatment1 / treatment3 3.154 0.674 Inf 1 5.377 <.0001
## treatment1 / treatment4 3.466 0.798 Inf 1 5.401 <.0001
## treatment1 / treatment5 2.638 0.652 Inf 1 3.921 0.0008
## treatment2 / treatment3 0.635 0.109 Inf 1 -2.644 0.0626
## treatment2 / treatment4 0.698 0.137 Inf 1 -1.832 0.3552
## treatment2 / treatment5 0.531 0.110 Inf 1 -3.048 0.0195
## treatment3 / treatment4 1.099 0.199 Inf 1 0.524 0.9850
## treatment3 / treatment5 0.836 0.172 Inf 1 -0.870 0.9078
## treatment4 / treatment5 0.761 0.174 Inf 1 -1.198 0.7526
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
cldlp <- cld(object = lm.emm,
adjust = "Tukey",
alpha = 0.05,
Letters = letters)
cldlp
## treatment prob SE df asymp.LCL asymp.UCL .group
## 2 0.646 0.0345 Inf 0.554 0.729 a
## 4 0.724 0.0379 Inf 0.617 0.810 ab
## 3 0.742 0.0262 Inf 0.670 0.804 ab
## 5 0.775 0.0305 Inf 0.687 0.844 b
## 1 0.901 0.0180 Inf 0.844 0.938 c
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
plot(brood$treatment, cbind(brood$dead_pupae, brood$live_pupae))
emmeans(pl_bind , pairwise ~ treatment, type = "response")
## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.901 0.0180 Inf 0.860 0.931
## 2 0.646 0.0345 Inf 0.576 0.711
## 3 0.742 0.0262 Inf 0.688 0.790
## 4 0.724 0.0379 Inf 0.644 0.792
## 5 0.775 0.0305 Inf 0.710 0.829
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 4.969 1.113 Inf 1 7.157 <.0001
## treatment1 / treatment3 3.154 0.674 Inf 1 5.377 <.0001
## treatment1 / treatment4 3.466 0.798 Inf 1 5.401 <.0001
## treatment1 / treatment5 2.638 0.652 Inf 1 3.921 0.0008
## treatment2 / treatment3 0.635 0.109 Inf 1 -2.644 0.0626
## treatment2 / treatment4 0.698 0.137 Inf 1 -1.832 0.3552
## treatment2 / treatment5 0.531 0.110 Inf 1 -3.048 0.0195
## treatment3 / treatment4 1.099 0.199 Inf 1 0.524 0.9850
## treatment3 / treatment5 0.836 0.172 Inf 1 -0.870 0.9078
## treatment4 / treatment5 0.761 0.174 Inf 1 -1.198 0.7526
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emm3 <- emmeans(pl_bind , ~ treatment, type = "response")
emm3
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.901 0.0180 Inf 0.860 0.931
## 2 0.646 0.0345 Inf 0.576 0.711
## 3 0.742 0.0262 Inf 0.688 0.790
## 4 0.724 0.0379 Inf 0.644 0.792
## 5 0.775 0.0305 Inf 0.710 0.829
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
pairs(emm3)
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 4.969 1.113 Inf 1 7.157 <.0001
## treatment1 / treatment3 3.154 0.674 Inf 1 5.377 <.0001
## treatment1 / treatment4 3.466 0.798 Inf 1 5.401 <.0001
## treatment1 / treatment5 2.638 0.652 Inf 1 3.921 0.0008
## treatment2 / treatment3 0.635 0.109 Inf 1 -2.644 0.0626
## treatment2 / treatment4 0.698 0.137 Inf 1 -1.832 0.3552
## treatment2 / treatment5 0.531 0.110 Inf 1 -3.048 0.0195
## treatment3 / treatment4 1.099 0.199 Inf 1 0.524 0.9850
## treatment3 / treatment5 0.836 0.172 Inf 1 -0.870 0.9078
## treatment4 / treatment5 0.761 0.174 Inf 1 -1.198 0.7526
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emmdf3 <- as.data.frame(emm3)
p <- ggplot(data = emmdf3, aes(x=treatment, y=prob, fill=treatment)) +
geom_col(position = "dodge", color = "black") +
coord_cartesian(ylim = c(0.5,1))
p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
labs(y = "Probability of Survival",) +
ggtitle("Combined Larvae and Pupae probability of being alive upon colony dissection") +
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)+
annotate(geom = "text",
x = 2, y = 1,
label = "Anova = 4.567e-13***",
size = 5) +
annotate(geom = "text",
x = c(2, 4, 3, 5, 1),
y = c(0.73, 0.81, 0.81, 0.85, 0.95),
label = c("a", "ab", "ab", "b", "c"),
size = 5)
egg.zero <- glmmTMB(eggs ~ treatment + whole.mean + qro, family = "nbinom2", brood) # had to take out alive, because it made model not converge -- Anova from glm shows it's not impactful anywas
summary(egg.zero)
## Family: nbinom2 ( log )
## Formula: eggs ~ treatment + whole.mean + qro
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 269.7 287.8 -124.9 249.7 35
##
##
## Dispersion parameter for nbinom2 family (): 0.855
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2479 0.6562 0.378 0.70563
## treatment2 -0.4207 0.5942 -0.708 0.47889
## treatment3 -0.9663 0.5852 -1.651 0.09867 .
## treatment4 -0.7371 0.5806 -1.270 0.20428
## treatment5 -1.0130 0.5952 -1.702 0.08874 .
## whole.mean 3.3460 1.2335 2.713 0.00667 **
## qroB3 0.6962 0.5810 1.198 0.23083
## qroB4 0.9635 0.6375 1.511 0.13073
## qroB5 1.1071 0.4767 2.322 0.02021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(egg.zero)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: eggs
## Chisq Df Pr(>Chisq)
## treatment 4.1574 4 0.385128
## whole.mean 7.3585 1 0.006675 **
## qro 6.7957 3 0.078702 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(egg.zero)
## Single term deletions
##
## Model:
## eggs ~ treatment + whole.mean + qro
## Df AIC
## <none> 269.74
## treatment 4 265.77
## whole.mean 1 274.31
## qro 3 269.88
eggs.sum <- brood %>%
group_by(treatment) %>%
summarise(eggs.m= mean(eggs),
eggs.n = length(eggs),
eggs.al = sd(eggs)) %>%
mutate(se.eg = eggs.al/sqrt(eggs.n))
eggs.sum
## # A tibble: 5 × 5
## treatment eggs.m eggs.n eggs.al se.eg
## <fct> <dbl> <int> <dbl> <dbl>
## 1 1 14.8 9 27.7 9.22
## 2 2 9.11 9 11.7 3.91
## 3 3 5.56 9 6.56 2.19
## 4 4 6.56 9 5.90 1.97
## 5 5 4.33 9 4.39 1.46
em.egg <- emmeans(egg.zero, "treatment")
pairs(em.egg)
## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.4207 0.594 35 0.708 0.9533
## treatment1 - treatment3 0.9663 0.585 35 1.651 0.4761
## treatment1 - treatment4 0.7371 0.581 35 1.269 0.7111
## treatment1 - treatment5 1.0130 0.595 35 1.702 0.4458
## treatment2 - treatment3 0.5456 0.596 35 0.916 0.8888
## treatment2 - treatment4 0.3163 0.594 35 0.532 0.9834
## treatment2 - treatment5 0.5923 0.571 35 1.038 0.8361
## treatment3 - treatment4 -0.2293 0.576 35 -0.398 0.9945
## treatment3 - treatment5 0.0467 0.592 35 0.079 1.0000
## treatment4 - treatment5 0.2760 0.588 35 0.470 0.9896
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
summary(em.egg)
## treatment emmean SE df lower.CL upper.CL
## 1 2.55 0.404 35 1.727 3.37
## 2 2.13 0.419 35 1.275 2.98
## 3 1.58 0.448 35 0.671 2.49
## 4 1.81 0.439 35 0.920 2.70
## 5 1.53 0.450 35 0.621 2.45
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
# Worker Data
workers <- read_csv("workers2.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")),
biobest_origin = col_factor(levels = c("B1",
"B2", "B3", "B4", "B5", "K1", "K2",
"K3")), start_date = col_skip(), mortality_date = col_skip(),
replaced = col_skip(), mortality = col_skip(), `replacement_wet weight` = col_skip(),
ending_wet_weight = col_skip()))
workers$qro <- as.factor(workers$biobest_origin)
workers$dry <- as.double(workers$`dry weight`)
workers$colony <- as.factor(workers$colony)
workers$proportion <- as.double(workers$proportion)
wrfp <- merge(wd.summary, workers, by.x = "colony")
wrpm <- merge(wrfp, mortality, by.x = "colony")
wrpm <- na.omit(wrpm)
alive.sum <- wrpm %>%
group_by(treatment) %>%
summarise(a.m= mean(alive),
a.n = length(alive),
sd.a = sd(alive),
n.a = length(alive)) %>%
mutate(sea = sd.a / sqrt(n.a))
alive.sum
## # A tibble: 5 × 6
## treatment a.m a.n sd.a n.a sea
## <fct> <dbl> <int> <dbl> <int> <dbl>
## 1 1 3.44 45 1.85 45 0.276
## 2 2 4.22 45 1.57 45 0.233
## 3 3 4.67 45 0.674 45 0.101
## 4 4 3.89 44 1.63 44 0.246
## 5 5 4.33 45 1.58 45 0.236
ggplot(data = alive.sum, aes(x = treatment, y = a.m, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(3,5)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = a.m - sea,
ymax = a.m + sea),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Workers Alive at End of Experiment",) +
ggtitle("Average Workers Alive per Treatment at End of Experiment") +
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)
ggplot(wrpm, aes(x=dry, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.005,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Worker Dry Weight") +
labs(y = "Count", x = "Weight (g) ")
range(wrpm$dry)
## [1] 0.01761 0.12047
shapiro.test(wrpm$dry)
##
## Shapiro-Wilk normality test
##
## data: wrpm$dry
## W = 0.906, p-value = 1.153e-10
which(wrpm$dry %in% c(0.12047))
## [1] 111
wrpm <- wrpm[-114,]
shapiro.test(wrpm$dry)
##
## Shapiro-Wilk normality test
##
## data: wrpm$dry
## W = 0.90487, p-value = 1.03e-10
ggplot(wrpm, aes(x=dry, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.005,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Worker Dry Weight (without highest value)") +
labs(y = "Count", x = "Weight (g) ")
wrpm$log.dry <- log(wrpm$dry)
shapiro.test(wrpm$log.dry)
##
## Shapiro-Wilk normality test
##
## data: wrpm$log.dry
## W = 0.98735, p-value = 0.0456
ggplot(wrpm, aes(x=log.dry, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.05,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Worker Dry Weight (Log Scale)") +
labs(y = "Count", x = "log(Weight (g)) ")
dry.lmer <- lmer(log.dry ~ treatment + whole.mean + qro + (1|colony), data = wrpm)
summary(dry.lmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log.dry ~ treatment + whole.mean + qro + (1 | colony)
## Data: wrpm
##
## REML criterion at convergence: 142.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.78063 -0.68432 0.05282 0.59718 3.14300
##
## Random effects:
## Groups Name Variance Std.Dev.
## colony (Intercept) 0.00000 0.0000
## Residual 0.09885 0.3144
## Number of obs: 223, groups: colony, 45
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.65348 0.07340 -49.776
## treatment2 -0.05444 0.06659 -0.818
## treatment3 -0.01085 0.06719 -0.161
## treatment4 -0.02273 0.06695 -0.340
## treatment5 0.04716 0.06631 0.711
## whole.mean 0.92366 0.12781 7.227
## qroB3 0.23525 0.07020 3.351
## qroB4 0.21422 0.07378 2.903
## qroB5 0.24286 0.05282 4.598
##
## Correlation of Fixed Effects:
## (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn qroB3 qroB4
## treatment2 -0.379
## treatment3 -0.358 0.503
## treatment4 -0.378 0.502 0.500
## treatment5 -0.474 0.494 0.489 0.492
## whole.mean -0.726 -0.097 -0.124 -0.093 0.031
## qroB3 -0.126 0.004 0.025 0.000 -0.001 -0.043
## qroB4 0.112 0.034 0.044 0.030 -0.011 -0.357 0.169
## qroB5 -0.151 0.007 0.009 0.002 -0.002 -0.075 0.218 0.231
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Anova(dry.lmer)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log.dry
## Chisq Df Pr(>Chisq)
## treatment 2.4419 4 0.6551
## whole.mean 52.2235 1 4.953e-13 ***
## qro 29.4964 3 1.761e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(dry.lmer, test = "Chisq")
## Single term deletions
##
## Model:
## log.dry ~ treatment + whole.mean + qro + (1 | colony)
## npar AIC LRT Pr(Chi)
## <none> 129.60
## treatment 4 124.13 2.530 0.6392
## whole.mean 1 171.00 43.398 4.467e-11 ***
## qro 3 151.75 28.145 3.387e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(resid(dry.lmer)) +
abline(h=0, col="red", lwd=2)
## integer(0)
qqnorm(resid(dry.lmer));qqline(resid(dry.lmer))
wdsum <- wrpm %>%
group_by(treatment) %>%
summarise(md = mean(dry),
sdd = sd(dry),
ndry = length(dry))
wdsum
## # A tibble: 5 × 4
## treatment md sdd ndry
## <fct> <dbl> <dbl> <int>
## 1 1 0.0474 0.0205 45
## 2 2 0.0456 0.0159 45
## 3 3 0.0492 0.0200 44
## 4 4 0.0487 0.0227 44
## 5 5 0.0478 0.0176 45
plot(wrpm$treatment, wrpm$dry)
mod1 <- glm(cbind(alive, dead) ~ treatment + start_wet_weight + qro, data = wrpm, family = binomial("logit"))
summary(mod1)
##
## Call:
## glm(formula = cbind(alive, dead) ~ treatment + start_wet_weight +
## qro, family = binomial("logit"), data = wrpm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.0986 0.2549 0.6978 1.0911 2.7303
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2092 0.3067 0.682 0.4952
## treatment2 1.1512 0.2534 4.544 5.52e-06 ***
## treatment3 2.0531 0.3203 6.409 1.46e-10 ***
## treatment4 0.6203 0.2368 2.620 0.0088 **
## treatment5 1.3492 0.2628 5.134 2.83e-07 ***
## start_wet_weight 8.8234 1.6533 5.337 9.45e-08 ***
## qroB3 -1.2300 0.2852 -4.312 1.62e-05 ***
## qroB4 -2.2616 0.2807 -8.057 7.82e-16 ***
## qroB5 -2.1178 0.2275 -9.309 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 716.28 on 222 degrees of freedom
## Residual deviance: 540.95 on 214 degrees of freedom
## AIC: 671.07
##
## Number of Fisher Scoring iterations: 5
anova(mod1)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: cbind(alive, dead)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 222 716.28
## treatment 4 52.966 218 663.31
## start_wet_weight 1 3.129 217 660.18
## qro 3 119.233 214 540.95
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(alive, dead)
## LR Chisq Df Pr(>Chisq)
## treatment 62.025 4 1.088e-12 ***
## start_wet_weight 31.860 1 1.657e-08 ***
## qro 119.233 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
drop1(mod1, test = "F")
## Single term deletions
##
## Model:
## cbind(alive, dead) ~ treatment + start_wet_weight + qro
## Df Deviance AIC F value Pr(>F)
## <none> 540.95 671.07
## treatment 4 602.97 725.09 6.1343 0.0001084 ***
## start_wet_weight 1 572.81 700.93 12.6038 0.0004733 ***
## qro 3 660.18 784.30 15.7229 2.816e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod1) # this doesn't look great
drop1(mod1, test = "F")
## Single term deletions
##
## Model:
## cbind(alive, dead) ~ treatment + start_wet_weight + qro
## Df Deviance AIC F value Pr(>F)
## <none> 540.95 671.07
## treatment 4 602.97 725.09 6.1343 0.0001084 ***
## start_wet_weight 1 572.81 700.93 12.6038 0.0004733 ***
## qro 3 660.18 784.30 15.7229 2.816e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mod1, pairwise ~ treatment, type = "response")
## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.591 0.0391 Inf 0.513 0.665
## 2 0.821 0.0292 Inf 0.756 0.871
## 3 0.918 0.0211 Inf 0.866 0.951
## 4 0.729 0.0357 Inf 0.654 0.793
## 5 0.848 0.0273 Inf 0.786 0.894
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 0.316 0.0801 Inf 1 -4.544 0.0001
## treatment1 / treatment3 0.128 0.0411 Inf 1 -6.409 <.0001
## treatment1 / treatment4 0.538 0.1273 Inf 1 -2.620 0.0668
## treatment1 / treatment5 0.259 0.0682 Inf 1 -5.134 <.0001
## treatment2 / treatment3 0.406 0.1375 Inf 1 -2.662 0.0597
## treatment2 / treatment4 1.700 0.4484 Inf 1 2.013 0.2594
## treatment2 / treatment5 0.820 0.2331 Inf 1 -0.697 0.9572
## treatment3 / treatment4 4.190 1.3768 Inf 1 4.361 0.0001
## treatment3 / treatment5 2.022 0.6975 Inf 1 2.040 0.2468
## treatment4 / treatment5 0.482 0.1313 Inf 1 -2.677 0.0573
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emm1 <- emmeans(mod1, ~ treatment, type = "response")
emm1
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.591 0.0391 Inf 0.513 0.665
## 2 0.821 0.0292 Inf 0.756 0.871
## 3 0.918 0.0211 Inf 0.866 0.951
## 4 0.729 0.0357 Inf 0.654 0.793
## 5 0.848 0.0273 Inf 0.786 0.894
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
pairs(emm1)
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 0.316 0.0801 Inf 1 -4.544 0.0001
## treatment1 / treatment3 0.128 0.0411 Inf 1 -6.409 <.0001
## treatment1 / treatment4 0.538 0.1273 Inf 1 -2.620 0.0668
## treatment1 / treatment5 0.259 0.0682 Inf 1 -5.134 <.0001
## treatment2 / treatment3 0.406 0.1375 Inf 1 -2.662 0.0597
## treatment2 / treatment4 1.700 0.4484 Inf 1 2.013 0.2594
## treatment2 / treatment5 0.820 0.2331 Inf 1 -0.697 0.9572
## treatment3 / treatment4 4.190 1.3768 Inf 1 4.361 0.0001
## treatment3 / treatment5 2.022 0.6975 Inf 1 2.040 0.2468
## treatment4 / treatment5 0.482 0.1313 Inf 1 -2.677 0.0573
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
emmdf <- as.data.frame(emm1)
workcld <- cld(object = emm1,
adjust = "Tukey",
alpha = 0.05,
Letters = letters)
workcld
## treatment prob SE df asymp.LCL asymp.UCL .group
## 1 0.591 0.0391 Inf 0.488 0.687 a
## 4 0.729 0.0357 Inf 0.628 0.811 ab
## 2 0.821 0.0292 Inf 0.733 0.884 bc
## 5 0.848 0.0273 Inf 0.764 0.906 bc
## 3 0.918 0.0211 Inf 0.845 0.959 c
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the logit scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
p <- ggplot(data = emmdf, aes(x=treatment, y=prob, fill=treatment)) +
geom_col(position = "dodge", color = "black") +
coord_cartesian(ylim = c(0.5,1))
p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
labs(y = "Probability of Survival",) +
ggtitle("Probability of workers surviving the whole experiment") +
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)+
annotate(geom = "text",
x = 1.5, y = 0.99,
label = "Anova = 1.300e-12***",
size = 5) +
annotate(geom = "text",
x = c(1, 4, 2, 5, 3),
y = c(0.7, 0.82, 0.9, 0.92, 0.97 ),
label = c("a", "ab", "bc", "bc", "c"),
size = 5)
treatment_table <- read.csv("tt.csv")
kable(treatment_table,
align = "ccccc",
caption = "Table 1: Experimental levels of Pristine® fed to each treatment group throughout the study, reported in PPB, milligrams added to pollen-sucrose mixture, and quantities of sucrose solution and pollen needed for each treatment",
booktabs = TRUE, valign = 't')
| Treatment | Pristine..Concentration..PPB. | Pristine..mg..per.treatment | Sucrose.Solution.per.treatment..ml. | Honey.Bee.Collected.pollen.per.treatment..g. |
|---|---|---|---|---|
| 1 (Control) | 0 | 0.000000 | 35 | 90 |
| 2 | 150 | 0.020157 | 35 | 90 |
| 3 | 1,500 | 0.201570 | 35 | 90 |
| 4 | 15,000 | 2.015700 | 35 | 90 |
| 5 | 150,000 | 20.157000 | 35 | 90 |
library(patchwork)
(dp / rfp / emp/ polp)
library(gridExtra)
grid.arrange(dp, rfp, emp, polp, nrow = 2)