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$whole_dif <- as.double(pollen$whole_dif)
pollen <- na.omit(pollen)
range(pollen$difference)
## [1] -0.98780 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
## Average pollen consumption per colony
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.5268 -0.2476 -0.1363 0.1874 1.0576
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.460679 0.021143 21.789 < 2e-16 ***
## treatment2 0.047174 0.030208 1.562 0.118630
## treatment3 0.100480 0.029931 3.357 0.000812 ***
## treatment4 0.053390 0.029931 1.784 0.074703 .
## treatment5 -0.001879 0.030272 -0.062 0.950524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3356 on 1231 degrees of freedom
## Multiple R-squared: 0.01281, Adjusted R-squared: 0.009604
## F-statistic: 3.994 on 4 and 1231 DF, p-value: 0.003177
wd.emm <- emmeans(wd.1, "treatment")
summary(wd.emm)
## treatment emmean SE df lower.CL upper.CL
## 1 0.461 0.0211 1231 0.419 0.502
## 2 0.508 0.0216 1231 0.466 0.550
## 3 0.561 0.0212 1231 0.520 0.603
## 4 0.514 0.0212 1231 0.473 0.556
## 5 0.459 0.0217 1231 0.416 0.501
##
## 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.1R1 0.8049667
## 4 1.1R2 0.5213458
## 5 1.2R1 0.4704294
## 6 1.2R2 0.3374200
## 7 1.3R1 0.3910053
## 8 1.3R2 0.4512891
## 9 1.4R2 0.6063016
## 10 1.5R2 0.7079516
## 11 1.7R2 0.7400381
## 12 1.9R2 0.2240081
## 13 2.11R2 0.4178270
## 14 2.12R2 0.4035568
## 15 2.1R1 0.7282895
## 16 2.1R2 0.6101589
## 17 2.2R1 0.4663045
## 18 2.2R2 0.5109234
## 19 2.3R1 0.4052000
## 20 2.3R2 0.3280036
## 21 2.4R2 0.3881394
## 22 2.5R2 0.7194915
## 23 2.7R2 0.5299685
## 24 2.9R2 0.5857152
## 25 3.11R2 0.4216410
## 26 3.12R2 0.3390993
## 27 3.1R1 0.8014682
## 28 3.1R2 0.3711948
## 29 3.2R1 0.8020500
## 30 3.2R2 0.3461010
## 31 3.3R1 0.5873429
## 32 3.3R2 0.8465806
## 33 3.4R2 0.4120433
## 34 3.5R2 0.8906211
## 35 3.7R2 0.6266680
## 36 3.9R2 0.4409511
## 37 4.11R2 0.3456011
## 38 4.12R2 0.2496295
## 39 4.1R1 0.8837867
## 40 4.1R2 0.7074755
## 41 4.2R1 0.3319238
## 42 4.2R2 0.3871742
## 43 4.3R1 0.3944500
## 44 4.3R2 0.5800074
## 45 4.4R2 0.8356247
## 46 4.5R2 0.2335967
## 47 4.7R2 0.6237400
## 48 4.9R2 0.5322950
## 49 5.11R2 0.4113960
## 50 5.12R2 0.2077741
## 51 5.1R1 0.6799737
## 52 5.1R2 0.3782286
## 53 5.2R1 0.7140056
## 54 5.2R2 0.4912468
## 55 5.3R1 0.2857654
## 56 5.3R2 0.2128778
## 57 5.4R2 0.4563045
## 58 5.5R2 0.7095479
## 59 5.7R2 0.6113189
## 60 5.9R2 0.4188290
# add queenright original colony column
qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)
merge <- merge(wd.summary, qro, by.x = "colony")
Input emerge and drone count data
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")))
trunc <- merge(merge, trunc.csv, by.x = "colony")
mortality <- read_csv("surviving workers per colony.csv")
mortality$colony <- as.factor(mortality$colony)
trunc <- merge(trunc, mortality, by.x = "colony")
trunc$qro <- as.factor(trunc$qro)
trunc <- trunc[-c(16)]
Get rid of colonies shut down
trunc$count[trunc$count < 0] <- NA
trunc <- na.omit(trunc)
range(trunc$count)
## [1] 1 27
hist(trunc$emerge)
mod1 <- glm(emerge ~ treatment + whole.mean + alive + round + qro, data = trunc, family = "poisson")
summary(mod1)
##
## Call:
## glm(formula = emerge ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = trunc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0123 -0.4334 -0.0705 0.1221 2.8545
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.742720 0.260236 14.382 <2e-16 ***
## treatment2 0.008746 0.075684 0.116 0.9080
## treatment3 0.048595 0.074081 0.656 0.5118
## treatment4 0.029443 0.075419 0.390 0.6962
## treatment5 -0.001636 0.076655 -0.021 0.9830
## whole.mean -0.310729 0.153178 -2.029 0.0425 *
## alive 0.020952 0.038937 0.538 0.5905
## round2 -0.117260 0.166319 -0.705 0.4808
## qroB3 -0.047805 0.094183 -0.508 0.6118
## qroB4 0.010764 0.108758 0.099 0.9212
## qroB5 0.046902 0.076473 0.613 0.5397
## qroK1 -0.006226 0.170374 -0.037 0.9709
## qroK2/K1 -0.119402 0.238927 -0.500 0.6173
## qroK3 -0.090795 0.237102 -0.383 0.7018
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 35.029 on 53 degrees of freedom
## Residual deviance: 26.561 on 40 degrees of freedom
## AIC: 348.43
##
## Number of Fisher Scoring iterations: 4
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 0.6940 4 0.95207
## whole.mean 4.1179 1 0.04243 *
## alive 0.2917 1 0.58914
## round 0
## qro 1.2918 6 0.97211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod2 <- glm.nb(emerge ~ treatment + whole.mean + alive + round + qro, data = trunc)
summary(mod2)
##
## Call:
## glm.nb(formula = emerge ~ treatment + whole.mean + alive + round +
## qro, data = trunc, init.theta = 1245956.089, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0122 -0.4334 -0.0705 0.1221 2.8544
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.742720 0.260240 14.382 <2e-16 ***
## treatment2 0.008746 0.075685 0.116 0.9080
## treatment3 0.048595 0.074083 0.656 0.5119
## treatment4 0.029443 0.075420 0.390 0.6963
## treatment5 -0.001636 0.076656 -0.021 0.9830
## whole.mean -0.310729 0.153180 -2.029 0.0425 *
## alive 0.020952 0.038937 0.538 0.5905
## round2 -0.117260 0.166322 -0.705 0.4808
## qroB3 -0.047805 0.094184 -0.508 0.6118
## qroB4 0.010764 0.108760 0.099 0.9212
## qroB5 0.046902 0.076474 0.613 0.5397
## qroK1 -0.006226 0.170377 -0.037 0.9709
## qroK2/K1 -0.119402 0.238930 -0.500 0.6173
## qroK3 -0.090795 0.237106 -0.383 0.7018
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1245956) family taken to be 1)
##
## Null deviance: 35.028 on 53 degrees of freedom
## Residual deviance: 26.560 on 40 degrees of freedom
## AIC: 350.43
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1245956
## Std. Err.: 24925247
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -320.432
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 0.6940 4 0.95207
## whole.mean 4.1179 1 0.04243 *
## alive 0.2917 1 0.58914
## round 0
## qro 1.2918 6 0.97211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3 <- glm(emerge ~ treatment + whole.mean + alive + round + qro, data = trunc, family = "quasipoisson")
summary(mod3)
##
## Call:
## glm(formula = emerge ~ treatment + whole.mean + alive + round +
## qro, family = "quasipoisson", data = trunc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0123 -0.4334 -0.0705 0.1221 2.8545
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.742720 0.219263 17.070 <2e-16 ***
## treatment2 0.008746 0.063768 0.137 0.8916
## treatment3 0.048595 0.062418 0.779 0.4408
## treatment4 0.029443 0.063544 0.463 0.6456
## treatment5 -0.001636 0.064586 -0.025 0.9799
## whole.mean -0.310729 0.129060 -2.408 0.0208 *
## alive 0.020952 0.032806 0.639 0.5267
## round2 -0.117260 0.140133 -0.837 0.4077
## qroB3 -0.047805 0.079354 -0.602 0.5503
## qroB4 0.010764 0.091635 0.117 0.9071
## qroB5 0.046902 0.064433 0.728 0.4709
## qroK1 -0.006226 0.143550 -0.043 0.9656
## qroK2/K1 -0.119402 0.201309 -0.593 0.5564
## qroK3 -0.090795 0.199771 -0.454 0.6519
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.7098969)
##
## Null deviance: 35.029 on 53 degrees of freedom
## Residual deviance: 26.561 on 40 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 0.6940 4 0.95207
## whole.mean 4.1179 1 0.04243 *
## alive 0.2917 1 0.58914
## round 0
## qro 1.2918 6 0.97211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod4 <- glm(emerge ~ treatment, data = trunc, family = "poisson")
summary(mod4)
##
## Call:
## glm(formula = emerge ~ treatment, family = "poisson", data = trunc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2047 -0.5727 -0.1825 0.2989 2.7418
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.594569 0.052414 68.580 <2e-16 ***
## treatment2 0.001497 0.072395 0.021 0.983
## treatment3 0.036417 0.070390 0.517 0.605
## treatment4 0.018803 0.072099 0.261 0.794
## treatment5 0.005479 0.074024 0.074 0.941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 35.029 on 53 degrees of freedom
## Residual deviance: 34.633 on 49 degrees of freedom
## AIC: 338.5
##
## Number of Fisher Scoring iterations: 4
Anova(mod4)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 0.39575 4 0.9828
trunc.2 <- subset(trunc, trunc$round != 1)
# Models no round 1
mod1 <- glm(emerge ~ treatment + whole.mean + alive + qro, data = trunc.2, family = "poisson")
summary(mod1)
##
## Call:
## glm(formula = emerge ~ treatment + whole.mean + alive + qro,
## family = "poisson", data = trunc.2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.82873 -0.43465 -0.05986 0.23330 1.48750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.76549 0.20951 17.973 <2e-16 ***
## treatment2 -0.02846 0.08823 -0.323 0.7470
## treatment3 0.01450 0.08506 0.170 0.8646
## treatment4 -0.07553 0.08888 -0.850 0.3954
## treatment5 -0.03243 0.09165 -0.354 0.7235
## whole.mean -0.48475 0.19215 -2.523 0.0116 *
## alive 0.01780 0.04177 0.426 0.6700
## qroB3 -0.03500 0.09482 -0.369 0.7121
## qroB4 0.04105 0.11659 0.352 0.7248
## qroB5 0.04293 0.07849 0.547 0.5844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 23.538 on 39 degrees of freedom
## Residual deviance: 12.768 on 30 degrees of freedom
## AIC: 249.75
##
## Number of Fisher Scoring iterations: 4
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 1.3532 4 0.85228
## whole.mean 6.4239 1 0.01126 *
## alive 0.1826 1 0.66918
## qro 0.5453 3 0.90882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod2 <- glm.nb(emerge ~ treatment + whole.mean + alive + qro, data = trunc.2)
summary(mod2)
##
## Call:
## glm.nb(formula = emerge ~ treatment + whole.mean + alive + qro,
## data = trunc.2, init.theta = 2119057.252, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.82873 -0.43464 -0.05986 0.23330 1.48749
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.76549 0.20952 17.972 <2e-16 ***
## treatment2 -0.02846 0.08823 -0.323 0.7471
## treatment3 0.01450 0.08506 0.170 0.8646
## treatment4 -0.07553 0.08888 -0.850 0.3954
## treatment5 -0.03243 0.09165 -0.354 0.7235
## whole.mean -0.48475 0.19215 -2.523 0.0116 *
## alive 0.01780 0.04177 0.426 0.6701
## qroB3 -0.03500 0.09482 -0.369 0.7121
## qroB4 0.04105 0.11659 0.352 0.7248
## qroB5 0.04293 0.07849 0.547 0.5844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2119057) family taken to be 1)
##
## Null deviance: 23.538 on 39 degrees of freedom
## Residual deviance: 12.767 on 30 degrees of freedom
## AIC: 251.75
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2119057
## Std. Err.: 54674070
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -229.752
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 1.3532 4 0.85228
## whole.mean 6.4239 1 0.01126 *
## alive 0.1826 1 0.66918
## qro 0.5453 3 0.90882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3 <- glm(emerge ~ treatment + whole.mean + alive + qro, data = trunc.2, family = "quasipoisson")
summary(mod3)
##
## Call:
## glm(formula = emerge ~ treatment + whole.mean + alive + qro,
## family = "quasipoisson", data = trunc.2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.82873 -0.43465 -0.05986 0.23330 1.48750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.76549 0.13828 27.230 < 2e-16 ***
## treatment2 -0.02846 0.05823 -0.489 0.628630
## treatment3 0.01450 0.05614 0.258 0.797924
## treatment4 -0.07553 0.05866 -1.288 0.207755
## treatment5 -0.03243 0.06049 -0.536 0.595853
## whole.mean -0.48475 0.12682 -3.822 0.000621 ***
## alive 0.01780 0.02757 0.646 0.523477
## qroB3 -0.03500 0.06258 -0.559 0.580176
## qroB4 0.04105 0.07695 0.533 0.597682
## qroB5 0.04293 0.05181 0.829 0.413867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.4356321)
##
## Null deviance: 23.538 on 39 degrees of freedom
## Residual deviance: 12.768 on 30 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
Anova(mod1)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 1.3532 4 0.85228
## whole.mean 6.4239 1 0.01126 *
## alive 0.1826 1 0.66918
## qro 0.5453 3 0.90882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod4 <- glm(emerge ~ treatment, data = trunc.2, family = "poisson")
summary(mod4)
##
## Call:
## glm(formula = emerge ~ treatment, family = "poisson", data = trunc.2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2101 -0.5203 -0.1296 0.2110 2.2253
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.60317 0.06238 57.763 <2e-16 ***
## treatment2 -0.02313 0.08588 -0.269 0.788
## treatment3 0.02266 0.08276 0.274 0.784
## treatment4 -0.08419 0.08715 -0.966 0.334
## treatment5 0.01112 0.08519 0.131 0.896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 23.538 on 39 degrees of freedom
## Residual deviance: 21.468 on 35 degrees of freedom
## AIC: 248.45
##
## Number of Fisher Scoring iterations: 4
Anova(mod4)
## Analysis of Deviance Table (Type II tests)
##
## Response: emerge
## LR Chisq Df Pr(>Chisq)
## treatment 2.0699 4 0.7229
emsum <- trunc %>%
group_by(treatment) %>%
summarise( met = mean(emerge),
sdet = sd(emerge),
net = length(emerge)) %>%
mutate( set = sdet/ sqrt(net))
emsum
## # A tibble: 5 × 5
## treatment met sdet net set
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 36.4 5.64 10 1.78
## 2 2 36.5 2.02 11 0.608
## 3 3 37.8 4.63 12 1.34
## 4 4 37.1 7.12 11 2.15
## 5 5 36.6 5.70 10 1.80
emsum.2 <- trunc.2 %>%
group_by(treatment) %>%
summarise( met = mean(emerge),
sdet = sd(emerge),
net = length(emerge)) %>%
mutate( set = sdet/ sqrt(net))
emsum.2
## # A tibble: 5 × 5
## treatment met sdet net set
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 36.7 6.82 7 2.58
## 2 2 35.9 1.89 8 0.666
## 3 3 37.6 5.08 9 1.69
## 4 4 33.8 2.55 8 0.901
## 5 5 37.1 6.29 8 2.22
emp <- ggplot(data = emsum, aes(x = treatment, y = met, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(33, 40)) +
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 (round 1 and 2") +
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)
emp.2 <- ggplot(data = emsum.2, aes(x = treatment, y = met, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(33, 40)) +
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 (round 2 only)") +
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)
emp
emp.2
# Drone Counts
hist(trunc$count)
range(trunc$count)
## [1] 1 27
Models with both rounds
d1 <- glm(count ~ treatment + whole.mean + alive + round + qro, data = trunc, family = "poisson")
summary(d1)
##
## Call:
## glm(formula = count ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = trunc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4799 -1.0154 -0.3271 0.8876 2.9649
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.89661 0.53366 1.680 0.0929 .
## treatment2 0.11924 0.15345 0.777 0.4372
## treatment3 -0.32248 0.15751 -2.047 0.0406 *
## treatment4 0.03111 0.15131 0.206 0.8371
## treatment5 0.31289 0.15739 1.988 0.0468 *
## whole.mean 2.50248 0.32335 7.739 1e-14 ***
## alive -0.05044 0.06502 -0.776 0.4379
## round2 0.28740 0.42747 0.672 0.5014
## qroB3 0.29768 0.15757 1.889 0.0589 .
## qroB4 -0.13880 0.17514 -0.793 0.4281
## qroB5 0.33766 0.13346 2.530 0.0114 *
## qroK1 -1.22585 0.45692 -2.683 0.0073 **
## qroK2/K1 -1.08452 0.83082 -1.305 0.1918
## qroK3 -1.66264 1.09259 -1.522 0.1281
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 291.565 on 53 degrees of freedom
## Residual deviance: 91.832 on 40 degrees of freedom
## AIC: 320.16
##
## Number of Fisher Scoring iterations: 5
d2 <- glm.nb(count ~ treatment + whole.mean + alive + round + qro, data = trunc)
summary(d2)
##
## Call:
## glm.nb(formula = count ~ treatment + whole.mean + alive + round +
## qro, data = trunc, init.theta = 13.4380213, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1127 -0.7363 -0.2593 0.6246 2.1625
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7876079 0.7034563 1.120 0.2629
## treatment2 0.1448808 0.2058901 0.704 0.4816
## treatment3 -0.3097078 0.2081349 -1.488 0.1367
## treatment4 -0.0003659 0.2058872 -0.002 0.9986
## treatment5 0.3541784 0.2102854 1.684 0.0921 .
## whole.mean 2.7078231 0.4396612 6.159 7.33e-10 ***
## alive -0.0553122 0.0940265 -0.588 0.5564
## round2 0.2829419 0.5176550 0.547 0.5847
## qroB3 0.3371566 0.2206565 1.528 0.1265
## qroB4 -0.1647450 0.2517763 -0.654 0.5129
## qroB5 0.3921952 0.1901719 2.062 0.0392 *
## qroK1 -1.2464954 0.5469443 -2.279 0.0227 *
## qroK2/K1 -1.0599902 0.9260832 -1.145 0.2524
## qroK3 -1.5787816 1.1684994 -1.351 0.1767
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(13.438) family taken to be 1)
##
## Null deviance: 182.988 on 53 degrees of freedom
## Residual deviance: 58.386 on 40 degrees of freedom
## AIC: 315.07
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 13.44
## Std. Err.: 7.16
##
## 2 x log-likelihood: -285.068
Anova(d2)
## Analysis of Deviance Table (Type II tests)
##
## Response: count
## LR Chisq Df Pr(>Chisq)
## treatment 11.511 4 0.02138 *
## whole.mean 38.189 1 6.421e-10 ***
## alive 0.320 1 0.57174
## round 0
## qro 13.009 6 0.04289 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d2em <- emmeans(d2, "treatment")
pairs(d2em)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.144881 0.206 Inf -0.704 0.9557
## treatment1 - treatment3 0.309708 0.208 Inf 1.488 0.5703
## treatment1 - treatment4 0.000366 0.206 Inf 0.002 1.0000
## treatment1 - treatment5 -0.354178 0.210 Inf -1.684 0.4437
## treatment2 - treatment3 0.454589 0.200 Inf 2.271 0.1544
## treatment2 - treatment4 0.145247 0.200 Inf 0.725 0.9509
## treatment2 - treatment5 -0.209298 0.193 Inf -1.083 0.8157
## treatment3 - treatment4 -0.309342 0.200 Inf -1.543 0.5343
## treatment3 - treatment5 -0.663886 0.203 Inf -3.267 0.0096
## treatment4 - treatment5 -0.354544 0.204 Inf -1.740 0.4095
##
## Results are averaged over the levels of: qro, round
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
Models with only round 2
d3 <- glm(count ~ treatment + whole.mean + alive + qro, data = trunc.2, family = "poisson")
summary(d3)
##
## Call:
## glm(formula = count ~ treatment + whole.mean + alive + qro, family = "poisson",
## data = trunc.2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4301 -1.1053 -0.4060 0.9594 2.9916
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.13922 0.32392 3.517 0.000437 ***
## treatment2 0.01134 0.15962 0.071 0.943348
## treatment3 -0.39057 0.16260 -2.402 0.016300 *
## treatment4 0.01959 0.15446 0.127 0.899097
## treatment5 0.20977 0.16672 1.258 0.208305
## whole.mean 2.46157 0.33607 7.325 2.4e-13 ***
## alive -0.02429 0.06679 -0.364 0.716095
## qroB3 0.28191 0.15766 1.788 0.073764 .
## qroB4 -0.10252 0.17806 -0.576 0.564755
## qroB5 0.35418 0.13412 2.641 0.008270 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 166.848 on 39 degrees of freedom
## Residual deviance: 73.299 on 30 degrees of freedom
## AIC: 255.3
##
## Number of Fisher Scoring iterations: 5
Anova(d3)
## Analysis of Deviance Table (Type II tests)
##
## Response: count
## LR Chisq Df Pr(>Chisq)
## treatment 15.327 4 0.004069 **
## whole.mean 54.254 1 1.762e-13 ***
## alive 0.131 1 0.717040
## qro 12.314 3 0.006383 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d4 <- glm.nb(count ~ treatment + whole.mean + alive + qro, data = trunc.2)
summary(d4)
##
## Call:
## glm.nb(formula = count ~ treatment + whole.mean + alive + qro,
## data = trunc.2, 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(d4)
## 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
AIC(d3, d4)
## df AIC
## d3 10 255.3033
## d4 11 251.5959
d4em <- emmeans(d4, "treatment")
pairs(d4em)
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.00185 0.213 Inf 0.009 1.0000
## treatment1 - treatment3 0.40510 0.214 Inf 1.890 0.3225
## treatment1 - treatment4 0.00740 0.209 Inf 0.035 1.0000
## treatment1 - treatment5 -0.21493 0.223 Inf -0.962 0.8723
## treatment2 - treatment3 0.40324 0.207 Inf 1.948 0.2919
## treatment2 - treatment4 0.00554 0.205 Inf 0.027 1.0000
## treatment2 - treatment5 -0.21678 0.203 Inf -1.066 0.8240
## treatment3 - treatment4 -0.39770 0.205 Inf -1.939 0.2965
## treatment3 - treatment5 -0.62003 0.214 Inf -2.903 0.0304
## treatment4 - treatment5 -0.22232 0.213 Inf -1.045 0.8346
##
## 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))
d2sum <- trunc.2 %>%
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 9.1 7.19 10 2.27
## 2 2 9.09 6.27 11 1.89
## 3 3 7.58 6.16 12 1.78
## 4 4 10.3 9.24 11 2.79
## 5 5 11 5.91 10 1.87
d2sum
## # 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
cp1 <- ggplot(data = d1sum, aes(x = treatment, y = md, 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 = md - sed,
ymax = md + sed),
position = position_dodge(0.9), width = 0.4) +
labs(y = "Mean Count",) +
ggtitle("Average Count of Drones per Treatment (round 1 and 2") +
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)
cp2 <- ggplot(data = d2sum, aes(x = treatment, y = md, fill = treatment)) +
geom_col(position = "dodge", col = "black")+
coord_cartesian(ylim=c(0,17)) +
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 (round 2 only)") +
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)
cp1
cp2