brood
Emily Runnion
2023-01-26
1 Input Pollen Data
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 <- 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.5654202 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.003177wd.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.95wd.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.4188290brood <- 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 <- merge(wd.summary, brood, by.x = "colony") # this is good because I calculated average pollen consumption two different ways (avg_pollen was calculated a couple months ago manually) and it's the same numbers as the group_by function that created whole.mean
mortality <- surviving_workers_per_colony <- read_csv("surviving workers per colony.csv")
mortality$colony <- as.factor(mortality$colony)
brood <- merge(brood, mortality, by.x = "colony")3 Brood cells
plot(brood$treatment, brood$brood_cells)
range(brood$brood_cells)## [1] 0 107Total brood count not affected by treatment, Anova = 0.343328
b.gn <- glm(brood_cells ~ treatment + whole.mean + alive + round, data = brood, family = "poisson")
summary(b.gn) # I don't see any overdispersion ##
## Call:
## glm(formula = brood_cells ~ treatment + whole.mean + alive +
## round, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7102 -1.3003 0.0875 0.9822 3.2979
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6559168 0.1321023 12.535 < 2e-16 ***
## treatment2 -0.0007538 0.0703605 -0.011 0.99145
## treatment3 0.0655198 0.0663548 0.987 0.32344
## treatment4 -0.0190518 0.0690270 -0.276 0.78254
## treatment5 -0.0755697 0.0733367 -1.030 0.30280
## whole.mean 2.7153647 0.1228109 22.110 < 2e-16 ***
## alive 0.0598922 0.0227659 2.631 0.00852 **
## round2 0.2479260 0.0502176 4.937 7.93e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 840.31 on 59 degrees of freedom
## Residual deviance: 196.37 on 52 degrees of freedom
## AIC: 522.11
##
## Number of Fisher Scoring iterations: 5Anova(b.gn)## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## treatment 4.49 4 0.343328
## whole.mean 506.50 1 < 2.2e-16 ***
## alive 7.24 1 0.007119 **
## round 25.11 1 5.403e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1b.gnb <- glm.nb(brood_cells ~ treatment + whole.mean + alive + round, data = brood)
AIC(b.gn, b.gnb)## df AIC
## b.gn 8 522.1057
## b.gnb 9 470.4001plot(b.gn)
plot(b.gnb)
brood.emm <- emmeans(b.gn, "treatment")
pairs(brood.emm)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.000754 0.0704 Inf 0.011 1.0000
## treatment1 - treatment3 -0.065520 0.0664 Inf -0.987 0.8611
## treatment1 - treatment4 0.019052 0.0690 Inf 0.276 0.9987
## treatment1 - treatment5 0.075570 0.0733 Inf 1.030 0.8413
## treatment2 - treatment3 -0.066274 0.0658 Inf -1.008 0.8520
## treatment2 - treatment4 0.018298 0.0683 Inf 0.268 0.9989
## treatment2 - treatment5 0.074816 0.0720 Inf 1.040 0.8370
## treatment3 - treatment4 0.084572 0.0625 Inf 1.354 0.6574
## treatment3 - treatment5 0.141090 0.0690 Inf 2.045 0.2447
## treatment4 - treatment5 0.056518 0.0713 Inf 0.792 0.9330
##
## Results are averaged over the levels of: round
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatessummary(brood.emm)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 3.41 0.0523 Inf 3.31 3.52
## 2 3.41 0.0517 Inf 3.31 3.51
## 3 3.48 0.0490 Inf 3.38 3.57
## 4 3.39 0.0520 Inf 3.29 3.50
## 5 3.34 0.0548 Inf 3.23 3.44
##
## Results are averaged over the levels of: round
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95# the poisson qq plot looks better, so with no overdispersion I will keep poisson even though nb AIC is better
sum_brood <- brood %>%
group_by(treatment) %>%
summarise( mean.b = mean(brood_cells),
n.b = length(brood_cells),
sd.b = sd(brood_cells))
sum_brood## # A tibble: 5 × 4
## treatment mean.b n.b sd.b
## <fct> <dbl> <int> <dbl>
## 1 1 33.8 12 22.1
## 2 2 34.8 12 12.0
## 3 3 48.8 12 26.4
## 4 4 38.4 12 28.3
## 5 5 30 12 18.23.0.1 honey pots
Treatment does not impact count of honey pots
hp.gn <- glm(honey_pot ~ treatment + whole.mean + alive + round, data = brood, family = "poisson")
summary(hp.gn) # overdispersion not horrible but not great ##
## Call:
## glm(formula = honey_pot ~ treatment + whole.mean + alive + round,
## family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.43646 -1.22754 -0.03545 0.62655 2.73036
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.32242 0.32088 1.005 0.314991
## treatment2 0.42460 0.16877 2.516 0.011873 *
## treatment3 0.16792 0.17519 0.959 0.337806
## treatment4 0.21789 0.17527 1.243 0.213820
## treatment5 0.20453 0.17799 1.149 0.250506
## whole.mean 1.11024 0.29453 3.769 0.000164 ***
## alive 0.15792 0.05624 2.808 0.004986 **
## round2 0.03829 0.11783 0.325 0.745238
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 143.28 on 59 degrees of freedom
## Residual deviance: 100.36 on 52 degrees of freedom
## AIC: 324.19
##
## Number of Fisher Scoring iterations: 5hp.gnb <- glm.nb(honey_pot ~ treatment + whole.mean + alive + round, data = brood)
summary(hp.gnb)##
## Call:
## glm.nb(formula = honey_pot ~ treatment + whole.mean + alive +
## round, data = brood, init.theta = 11.3790581, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.12424 -1.05468 -0.06497 0.53218 2.09725
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.30059 0.37846 0.794 0.42706
## treatment2 0.41417 0.20993 1.973 0.04851 *
## treatment3 0.14166 0.21664 0.654 0.51318
## treatment4 0.22496 0.21511 1.046 0.29564
## treatment5 0.17380 0.21856 0.795 0.42649
## whole.mean 1.14711 0.37278 3.077 0.00209 **
## alive 0.16208 0.06558 2.472 0.01345 *
## round2 0.03275 0.15062 0.217 0.82787
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(11.3791) family taken to be 1)
##
## Null deviance: 97.271 on 59 degrees of freedom
## Residual deviance: 68.368 on 52 degrees of freedom
## AIC: 320.07
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 11.38
## Std. Err.: 6.26
##
## 2 x log-likelihood: -302.07AIC(hp.gn, hp.gnb)## df AIC
## hp.gn 8 324.1918
## hp.gnb 9 320.0697plot(hp.gn)
plot(hp.gnb)
Anova(hp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pot
## LR Chisq Df Pr(>Chisq)
## treatment 4.2935 4 0.367738
## whole.mean 9.2018 1 0.002418 **
## alive 6.4046 1 0.011383 *
## round 0.0472 1 0.828010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(hp.gn)## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pot
## LR Chisq Df Pr(>Chisq)
## treatment 6.8464 4 0.1442327
## whole.mean 14.1647 1 0.0001675 ***
## alive 8.9967 1 0.0027046 **
## round 0.1061 1 0.7446642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$honey_pot)
hp.emm <- emmeans(hp.gnb, "treatment")
summary(hp.emm)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 1.60 0.163 Inf 1.28 1.92
## 2 2.01 0.141 Inf 1.73 2.29
## 3 1.74 0.153 Inf 1.44 2.04
## 4 1.82 0.150 Inf 1.53 2.12
## 5 1.77 0.153 Inf 1.47 2.07
##
## Results are averaged over the levels of: round
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95pairs(hp.emm)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.4142 0.210 Inf -1.973 0.2793
## treatment1 - treatment3 -0.1417 0.217 Inf -0.654 0.9660
## treatment1 - treatment4 -0.2250 0.215 Inf -1.046 0.8339
## treatment1 - treatment5 -0.1738 0.219 Inf -0.795 0.9321
## treatment2 - treatment3 0.2725 0.197 Inf 1.385 0.6372
## treatment2 - treatment4 0.1892 0.197 Inf 0.960 0.8728
## treatment2 - treatment5 0.2404 0.200 Inf 1.204 0.7491
## treatment3 - treatment4 -0.0833 0.202 Inf -0.413 0.9939
## treatment3 - treatment5 -0.0321 0.207 Inf -0.155 0.9999
## treatment4 - treatment5 0.0512 0.207 Inf 0.247 0.9992
##
## Results are averaged over the levels of: round
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatessum_hp <- brood %>%
group_by(treatment) %>%
summarise( mean.h = mean(honey_pot),
n.h = length(honey_pot),
sd.h = sd(honey_pot))
sum_hp## # A tibble: 5 × 4
## treatment mean.h n.h sd.h
## <fct> <dbl> <int> <dbl>
## 1 1 4.75 12 3.44
## 2 2 7.92 12 3.15
## 3 3 6.92 12 4.50
## 4 4 6.5 12 3.29
## 5 5 6.08 12 4.063.0.2 Larvae
3.0.2.1 Live Larvae
ll.gn <- glm(live_larvae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(ll.gn)##
## Call:
## glm(formula = live_larvae ~ treatment + whole.mean + alive +
## round + qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.936 -2.258 -0.243 1.351 5.076
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.192482 0.413847 -0.465 0.641856
## treatment2 -0.323271 0.092384 -3.499 0.000467 ***
## treatment3 0.043415 0.081313 0.534 0.593397
## treatment4 -0.559050 0.096698 -5.781 7.41e-09 ***
## treatment5 -0.256481 0.094480 -2.715 0.006634 **
## whole.mean 3.611775 0.203193 17.775 < 2e-16 ***
## alive -0.005275 0.031149 -0.169 0.865535
## round2 1.450699 0.384357 3.774 0.000160 ***
## qroB3 -0.422326 0.127784 -3.305 0.000950 ***
## qroB4 -0.155574 0.110540 -1.407 0.159311
## qroB5 0.497129 0.089010 5.585 2.34e-08 ***
## qroK1 1.296553 0.385974 3.359 0.000782 ***
## qroK2/K1 0.746314 0.510107 1.463 0.143453
## qroK3 1.550465 0.514478 3.014 0.002581 **
## 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: 1059.03 on 59 degrees of freedom
## Residual deviance: 377.59 on 46 degrees of freedom
## AIC: 649.99
##
## Number of Fisher Scoring iterations: 5ll.gnb <- glm.nb(live_larvae ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(ll.gnb)##
## Call:
## glm.nb(formula = live_larvae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 2.421856617, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7305 -1.1020 -0.0898 0.4673 2.1460
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.84237 0.91615 -2.011 0.044327 *
## treatment2 -0.26268 0.29567 -0.888 0.374299
## treatment3 -0.07010 0.29974 -0.234 0.815083
## treatment4 -0.42945 0.30799 -1.394 0.163200
## treatment5 -0.27526 0.30312 -0.908 0.363831
## whole.mean 4.92329 0.63113 7.801 6.15e-15 ***
## alive 0.19334 0.09767 1.980 0.047750 *
## round2 1.45856 0.78212 1.865 0.062196 .
## qroB3 -0.61381 0.38078 -1.612 0.106970
## qroB4 -0.36291 0.39986 -0.908 0.364092
## qroB5 0.97252 0.28604 3.400 0.000674 ***
## qroK1 1.09358 0.79891 1.369 0.171051
## qroK2/K1 0.79681 0.99129 0.804 0.421510
## qroK3 1.56033 1.08979 1.432 0.152207
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.4219) family taken to be 1)
##
## Null deviance: 169.268 on 59 degrees of freedom
## Residual deviance: 75.829 on 46 degrees of freedom
## AIC: 466.56
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.422
## Std. Err.: 0.594
##
## 2 x log-likelihood: -436.558AIC(ll.gn, ll.gnb) # poisson overdispersed, nb AIC better ## df AIC
## ll.gn 14 649.9867
## ll.gnb 15 466.5582plot(ll.gnb)
Anova(ll.gn)## Analysis of Deviance Table (Type II tests)
##
## Response: live_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 60.38 4 2.415e-12 ***
## whole.mean 359.34 1 < 2.2e-16 ***
## alive 0.03 1 0.8657
## round 0
## qro 85.08 6 3.180e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(ll.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: live_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 2.526 4 0.639994
## whole.mean 51.533 1 7.04e-13 ***
## alive 4.106 1 0.042742 *
## round 0
## qro 21.089 6 0.001769 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$live_larvae)
ll.sum <- brood %>%
group_by(treatment) %>%
summarise(ll.m= mean(live_larvae),
ll.n = length(live_larvae),
sd.ll = sd(live_larvae))
ll.sum## # A tibble: 5 × 4
## treatment ll.m ll.n sd.ll
## <fct> <dbl> <int> <dbl>
## 1 1 24.8 12 24.2
## 2 2 17.1 12 14.2
## 3 3 32.2 12 22.3
## 4 4 18.2 12 16.3
## 5 5 17.2 12 13.63.0.2.2 Dead Larvae
dl.gn <- glm(dead_larvae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(dl.gn)##
## Call:
## glm(formula = dead_larvae ~ treatment + whole.mean + alive +
## round + qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2045 -1.4558 -0.4256 0.4491 5.3734
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.48594 0.65906 -3.772 0.000162 ***
## treatment2 0.51441 0.25588 2.010 0.044390 *
## treatment3 0.66273 0.23753 2.790 0.005270 **
## treatment4 0.71362 0.25191 2.833 0.004614 **
## treatment5 -0.06032 0.28763 -0.210 0.833905
## whole.mean 3.80155 0.49445 7.688 1.49e-14 ***
## alive 0.30726 0.10880 2.824 0.004740 **
## round2 -0.59516 0.41872 -1.421 0.155214
## qroB3 0.06794 0.28822 0.236 0.813660
## qroB4 0.65803 0.28644 2.297 0.021602 *
## qroB5 1.20230 0.20836 5.770 7.92e-09 ***
## qroK1 -0.27982 0.41887 -0.668 0.504120
## qroK2/K1 0.27087 0.61463 0.441 0.659426
## qroK3 -16.56609 1275.75394 -0.013 0.989639
## 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: 444.72 on 59 degrees of freedom
## Residual deviance: 174.90 on 46 degrees of freedom
## AIC: 336.35
##
## Number of Fisher Scoring iterations: 13dl.gnb <- glm.nb(dead_larvae ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(dl.gnb)##
## Call:
## glm.nb(formula = dead_larvae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 1.309674691, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1769 -1.1220 -0.3143 0.2250 1.9979
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.889e+00 1.309e+00 -1.442 0.14921
## treatment2 4.813e-01 4.688e-01 1.027 0.30465
## treatment3 8.993e-01 4.635e-01 1.940 0.05237 .
## treatment4 6.262e-01 4.734e-01 1.323 0.18591
## treatment5 -3.202e-02 4.950e-01 -0.065 0.94842
## whole.mean 2.623e+00 9.358e-01 2.803 0.00506 **
## alive 2.789e-01 1.693e-01 1.647 0.09951 .
## round2 -3.225e-01 1.014e+00 -0.318 0.75036
## qroB3 -2.459e-01 5.953e-01 -0.413 0.67961
## qroB4 6.964e-01 5.877e-01 1.185 0.23609
## qroB5 9.095e-01 4.476e-01 2.032 0.04214 *
## qroK1 6.941e-02 1.039e+00 0.067 0.94673
## qroK2/K1 7.456e-02 1.320e+00 0.056 0.95497
## qroK3 -3.821e+01 6.711e+07 0.000 1.00000
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.3097) family taken to be 1)
##
## Null deviance: 118.618 on 59 degrees of freedom
## Residual deviance: 66.516 on 46 degrees of freedom
## AIC: 292.61
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.310
## Std. Err.: 0.413
##
## 2 x log-likelihood: -262.614AIC(dl.gn, dl.gnb) #poisson overdispersed and nb AIC better ## df AIC
## dl.gn 14 336.3505
## dl.gnb 15 292.6140plot(dl.gnb)
Anova(dl.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: dead_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 5.4864 4 0.240927
## whole.mean 7.3377 1 0.006752 **
## alive 2.4515 1 0.117412
## round 0
## qro 8.0093 6 0.237423
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$dead_larvae)
dl.sum <- brood %>%
group_by(treatment) %>%
summarise(dl.m= mean(dead_larvae),
dl.n = length(dead_larvae),
sd.dl = sd(dead_larvae))
dl.sum## # A tibble: 5 × 4
## treatment dl.m dl.n sd.dl
## <fct> <dbl> <int> <dbl>
## 1 1 2.08 12 2.57
## 2 2 3.67 12 4.58
## 3 3 6.67 12 5.65
## 4 4 6.5 12 13.4
## 5 5 2.17 12 2.593.0.2.3 All larvae
brood$all_larvae <- brood$dead_larvae + brood$live_larvae
al.gn <- glm(all_larvae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(al.gn)##
## Call:
## glm(formula = all_larvae ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.419 -2.036 0.000 1.112 4.841
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.22118 0.31121 0.711 0.47726
## treatment2 -0.22206 0.08600 -2.582 0.00982 **
## treatment3 0.10593 0.07628 1.389 0.16492
## treatment4 -0.35973 0.08775 -4.100 4.14e-05 ***
## treatment5 -0.24482 0.08962 -2.732 0.00630 **
## whole.mean 3.69409 0.18752 19.699 < 2e-16 ***
## alive 0.02845 0.02984 0.953 0.34041
## round2 0.87855 0.27506 3.194 0.00140 **
## qroB3 -0.35523 0.11659 -3.047 0.00231 **
## qroB4 -0.06486 0.10243 -0.633 0.52659
## qroB5 0.60704 0.08136 7.461 8.59e-14 ***
## qroK1 0.77995 0.27672 2.819 0.00482 **
## qroK2/K1 0.52390 0.38490 1.361 0.17348
## qroK3 0.73639 0.43723 1.684 0.09214 .
## 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: 1145.08 on 59 degrees of freedom
## Residual deviance: 278.92 on 46 degrees of freedom
## AIC: 576.07
##
## Number of Fisher Scoring iterations: 5al.gnb <- glm.nb(all_larvae ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(al.gnb)##
## Call:
## glm.nb(formula = all_larvae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 4.670976872, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.37207 -0.99239 -0.04914 0.45894 2.29714
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.80827 0.66747 -1.211 0.2259
## treatment2 -0.17292 0.22116 -0.782 0.4343
## treatment3 0.08391 0.22257 0.377 0.7062
## treatment4 -0.31830 0.23001 -1.384 0.1664
## treatment5 -0.24325 0.22728 -1.070 0.2845
## whole.mean 4.34503 0.46996 9.246 < 2e-16 ***
## alive 0.16228 0.07365 2.203 0.0276 *
## round2 0.94335 0.56234 1.678 0.0934 .
## qroB3 -0.54522 0.28736 -1.897 0.0578 .
## qroB4 -0.16259 0.29440 -0.552 0.5808
## qroB5 0.83787 0.21338 3.927 8.62e-05 ***
## qroK1 0.68938 0.57500 1.199 0.2306
## qroK2/K1 0.51135 0.72386 0.706 0.4799
## qroK3 0.79850 0.81869 0.975 0.3294
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.671) family taken to be 1)
##
## Null deviance: 228.519 on 59 degrees of freedom
## Residual deviance: 71.273 on 46 degrees of freedom
## AIC: 467.49
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.67
## Std. Err.: 1.20
##
## 2 x log-likelihood: -437.493AIC(al.gn, al.gnb)## df AIC
## al.gn 14 576.0719
## al.gnb 15 467.4929Anova(al.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: all_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 4.519 4 0.3402420
## whole.mean 78.684 1 < 2.2e-16 ***
## alive 5.129 1 0.0235351 *
## round 0
## qro 26.731 6 0.0001626 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$all_larvae)
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 26.9 12 25.1
## 2 2 20.8 12 16.2
## 3 3 38.8 12 25.4
## 4 4 24.7 12 25.4
## 5 5 19.3 12 14.73.0.2.4 Larvae Outliers
out <- boxplot.stats(brood$dead_larvae)$out
out## [1] 15 14 15 18 46out_dl <- which(brood$dead_larvae %in% c(out))
out_dl## [1] 22 32 34 39 45brood[out_dl, ]## colony whole.mean round dose treatment replicate brood_cells honey_pot eggs
## 22 2.5R2 0.7194915 2 150 2 5 44 8 14
## 32 3.3R2 0.8465806 2 1500 3 3 64 9 15
## 34 3.5R2 0.8906211 2 1500 3 5 96 9 16
## 39 4.1R1 0.8837867 1 15000 4 1 107 4 10
## 45 4.4R2 0.8356247 2 15000 4 4 53 10 13
## dead_larvae live_larvae dead_pupae live_pupae dead_drones live_drones drones
## 22 15 29 8 12 0 0 11
## 32 14 49 24 1 0 0 10
## 34 15 71 12 18 1 2 15
## 39 18 49 44 7 0 0 1
## 45 46 32 6 0 0 0 27
## avg_pollen qro alive dead all_larvae
## 22 0.7194915 B4 4 1 44
## 32 0.5044652 B3 5 0 63
## 34 0.4120433 B4 4 1 86
## 39 0.8837867 K1 5 0 67
## 45 0.5800074 B5 5 0 78brood_dl_out <- brood[-c(22, 32, 34, 39, 45), ]
dl.gn <- glm(dead_larvae ~ treatment + whole.mean + alive + round + qro, data = brood_dl_out, family = "poisson")
summary(dl.gn)##
## Call:
## glm(formula = dead_larvae ~ treatment + whole.mean + alive +
## round + qro, family = "poisson", data = brood_dl_out)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6531 -1.5465 -0.3092 0.5875 4.0338
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.25682 0.71647 -0.358 0.72000
## treatment2 0.21032 0.28332 0.742 0.45788
## treatment3 0.69170 0.26146 2.646 0.00816 **
## treatment4 -0.06780 0.34038 -0.199 0.84211
## treatment5 0.02320 0.28525 0.081 0.93519
## whole.mean 1.54099 0.63828 2.414 0.01577 *
## alive 0.12119 0.10903 1.112 0.26634
## round2 -0.51983 0.42496 -1.223 0.22124
## qroB3 -0.40444 0.52509 -0.770 0.44116
## qroB4 0.08812 0.46686 0.189 0.85029
## qroB5 0.20399 0.30459 0.670 0.50304
## qroK1 0.05738 0.42693 0.134 0.89309
## qroK2/K1 -0.05138 0.63331 -0.081 0.93534
## qroK3 -16.19174 1275.75396 -0.013 0.98987
## 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: 199.08 on 54 degrees of freedom
## Residual deviance: 132.87 on 41 degrees of freedom
## AIC: 270.32
##
## Number of Fisher Scoring iterations: 13dl.gnb <- glm.nb(dead_larvae ~ treatment + whole.mean + alive + round + qro, data = brood_dl_out)
summary(dl.gnb)##
## Call:
## glm.nb(formula = dead_larvae ~ treatment + whole.mean + alive +
## round + qro, data = brood_dl_out, init.theta = 1.385441703,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0598 -1.2596 -0.2605 0.3580 1.9501
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.347e-01 1.321e+00 -0.329 0.742
## treatment2 2.252e-01 4.800e-01 0.469 0.639
## treatment3 7.164e-01 4.839e-01 1.481 0.139
## treatment4 1.090e-01 5.145e-01 0.212 0.832
## treatment5 -5.237e-02 4.803e-01 -0.109 0.913
## whole.mean 1.047e+00 1.087e+00 0.963 0.335
## alive 2.099e-01 1.619e-01 1.296 0.195
## round2 -5.713e-01 1.004e+00 -0.569 0.569
## qroB3 -4.966e-01 7.348e-01 -0.676 0.499
## qroB4 2.836e-01 7.501e-01 0.378 0.705
## qroB5 4.915e-01 4.766e-01 1.031 0.302
## qroK1 8.401e-02 1.024e+00 0.082 0.935
## qroK2/K1 -2.322e-01 1.302e+00 -0.178 0.858
## qroK3 -3.744e+01 4.633e+07 0.000 1.000
## qroK3/K2 NA NA NA NA
##
## (Dispersion parameter for Negative Binomial(1.3854) family taken to be 1)
##
## Null deviance: 82.352 on 54 degrees of freedom
## Residual deviance: 60.289 on 41 degrees of freedom
## AIC: 246.44
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.385
## Std. Err.: 0.507
##
## 2 x log-likelihood: -216.444Anova(dl.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: dead_larvae
## LR Chisq Df Pr(>Chisq)
## treatment 3.2360 4 0.5191
## whole.mean 0.8669 1 0.3518
## alive 1.4894 1 0.2223
## round 0
## qro 4.5988 6 0.5962plot(brood_dl_out$treatment, brood_dl_out$dead_larvae)
3.0.2.5 cbind Larvae
larvae.mortality <- glm(cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + round + qro, data = brood, family = binomial("logit"))
summary(larvae.mortality)##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ treatment + whole.mean +
## round + qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.5952 -1.3529 0.0449 1.4375 3.9529
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3003 0.6680 1.946 0.051606 .
## treatment2 -0.9340 0.2772 -3.369 0.000754 ***
## treatment3 -0.7955 0.2518 -3.159 0.001581 **
## treatment4 -1.4739 0.2661 -5.538 3.06e-08 ***
## treatment5 -0.4817 0.2982 -1.615 0.106284
## whole.mean -0.8593 0.5873 -1.463 0.143406
## round2 2.0956 0.5694 3.681 0.000233 ***
## qroB3 -0.4126 0.3304 -1.249 0.211707
## qroB4 -0.4822 0.3038 -1.587 0.112462
## qroB5 -0.6207 0.2097 -2.961 0.003069 **
## qroK1 1.7468 0.5744 3.041 0.002356 **
## qroK2/K1 0.5697 0.8054 0.707 0.479327
## qroK3 17.5385 1006.6076 0.017 0.986099
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 358.67 on 56 degrees of freedom
## Residual deviance: 279.19 on 44 degrees of freedom
## AIC: 413.88
##
## Number of Fisher Scoring iterations: 14Anova(larvae.mortality)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## treatment 36.909 4 1.881e-07 ***
## whole.mean 2.196 1 0.1383741
## round 0
## qro 26.545 6 0.0001762 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(larvae.mortality)
lm.emm <- emmeans(larvae.mortality, "treatment")
pairs(lm.emm)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.934 0.277 Inf 3.369 0.0068
## treatment1 - treatment3 0.796 0.252 Inf 3.159 0.0137
## treatment1 - treatment4 1.474 0.266 Inf 5.538 <.0001
## treatment1 - treatment5 0.482 0.298 Inf 1.615 0.4876
## treatment2 - treatment3 -0.138 0.230 Inf -0.601 0.9750
## treatment2 - treatment4 0.540 0.256 Inf 2.108 0.2166
## treatment2 - treatment5 -0.452 0.274 Inf -1.651 0.4644
## treatment3 - treatment4 0.678 0.205 Inf 3.314 0.0082
## treatment3 - treatment5 -0.314 0.259 Inf -1.213 0.7439
## treatment4 - treatment5 -0.992 0.280 Inf -3.537 0.0037
##
## Results are averaged over the levels of: qro, round
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatessummary(lm.emm)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 4.21 126 Inf -242 251
## 2 3.27 126 Inf -243 250
## 3 3.41 126 Inf -243 250
## 4 2.73 126 Inf -244 249
## 5 3.72 126 Inf -243 250
##
## Results are averaged over the levels of: qro, round
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95plot(brood$treatment, cbind(brood$dead_larvae, brood$live_larvae))
emmeans(larvae.mortality, pairwise ~ treatment, type = "response")## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.985 1.82 Inf 2.22e-16 1
## 2 0.963 4.43 Inf 2.22e-16 1
## 3 0.968 3.90 Inf 2.22e-16 1
## 4 0.939 7.22 Inf 2.22e-16 1
## 5 0.976 2.90 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## 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 2.545 0.705 Inf 1 3.369 0.0068
## treatment1 / treatment3 2.216 0.558 Inf 1 3.159 0.0137
## treatment1 / treatment4 4.366 1.162 Inf 1 5.538 <.0001
## treatment1 / treatment5 1.619 0.483 Inf 1 1.615 0.4876
## treatment2 / treatment3 0.871 0.201 Inf 1 -0.601 0.9750
## treatment2 / treatment4 1.716 0.440 Inf 1 2.108 0.2166
## treatment2 / treatment5 0.636 0.174 Inf 1 -1.651 0.4644
## treatment3 / treatment4 1.971 0.403 Inf 1 3.314 0.0082
## treatment3 / treatment5 0.731 0.189 Inf 1 -1.213 0.7439
## treatment4 / treatment5 0.371 0.104 Inf 1 -3.537 0.0037
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemm1 <- emmeans(larvae.mortality, ~ treatment, type = "response")
emm1## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.985 1.82 Inf 2.22e-16 1
## 2 0.963 4.43 Inf 2.22e-16 1
## 3 0.968 3.90 Inf 2.22e-16 1
## 4 0.939 7.22 Inf 2.22e-16 1
## 5 0.976 2.90 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalepairs(emm1)## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 2.545 0.705 Inf 1 3.369 0.0068
## treatment1 / treatment3 2.216 0.558 Inf 1 3.159 0.0137
## treatment1 / treatment4 4.366 1.162 Inf 1 5.538 <.0001
## treatment1 / treatment5 1.619 0.483 Inf 1 1.615 0.4876
## treatment2 / treatment3 0.871 0.201 Inf 1 -0.601 0.9750
## treatment2 / treatment4 1.716 0.440 Inf 1 2.108 0.2166
## treatment2 / treatment5 0.636 0.174 Inf 1 -1.651 0.4644
## treatment3 / treatment4 1.971 0.403 Inf 1 3.314 0.0082
## treatment3 / treatment5 0.731 0.189 Inf 1 -1.213 0.7439
## treatment4 / treatment5 0.371 0.104 Inf 1 -3.537 0.0037
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemmdf <- 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,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 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)
3.0.3 Pupae
3.0.3.1 Live Pupae
lp.gn <- glm(live_pupae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(lp.gn)##
## Call:
## glm(formula = live_pupae ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.8267 -1.6016 -0.0004 0.6527 3.1322
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.50614 0.54531 0.928 0.3533
## treatment2 0.10356 0.19551 0.530 0.5963
## treatment3 -0.14536 0.20411 -0.712 0.4764
## treatment4 -0.52525 0.22745 -2.309 0.0209 *
## treatment5 -0.23907 0.22911 -1.043 0.2967
## whole.mean 3.56458 0.49313 7.229 4.88e-13 ***
## alive -0.13253 0.06857 -1.933 0.0533 .
## round2 -0.14985 0.44207 -0.339 0.7346
## qroB3 -0.82745 0.32882 -2.516 0.0119 *
## qroB4 -0.30731 0.24554 -1.252 0.2107
## qroB5 -0.15620 0.22207 -0.703 0.4818
## qroK1 -0.74911 0.45904 -1.632 0.1027
## qroK2/K1 -17.39200 1462.38375 -0.012 0.9905
## qroK3 0.66201 0.69248 0.956 0.3391
## 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: 256.88 on 59 degrees of freedom
## Residual deviance: 155.29 on 46 degrees of freedom
## AIC: 329.63
##
## Number of Fisher Scoring iterations: 14lp.gnb <- glm.nb(live_pupae ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(lp.gnb)##
## Call:
## glm.nb(formula = live_pupae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 1.968375486, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.79611 -1.07682 -0.06755 0.40300 1.57298
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.320e-01 1.047e+00 -0.126 0.900
## treatment2 -2.715e-02 3.720e-01 -0.073 0.942
## treatment3 -3.257e-01 3.841e-01 -0.848 0.397
## treatment4 -5.645e-01 3.993e-01 -1.414 0.157
## treatment5 -3.916e-01 3.946e-01 -0.992 0.321
## whole.mean 3.830e+00 8.278e-01 4.626 3.72e-06 ***
## alive -1.642e-05 1.204e-01 0.000 1.000
## round2 -2.117e-01 8.772e-01 -0.241 0.809
## qroB3 -5.433e-01 4.923e-01 -1.104 0.270
## qroB4 -3.699e-01 4.954e-01 -0.747 0.455
## qroB5 2.235e-01 3.689e-01 0.606 0.545
## qroK1 -8.135e-01 9.035e-01 -0.900 0.368
## qroK2/K1 -3.803e+01 4.745e+07 0.000 1.000
## qroK3 5.721e-01 1.273e+00 0.449 0.653
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.9684) family taken to be 1)
##
## Null deviance: 106.308 on 59 degrees of freedom
## Residual deviance: 68.473 on 46 degrees of freedom
## AIC: 300.19
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.968
## Std. Err.: 0.655
##
## 2 x log-likelihood: -270.190AIC(lp.gn, lp.gnb)## df AIC
## lp.gn 14 329.6261
## lp.gnb 15 300.1901Anova(lp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: live_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 2.9818 4 0.5609
## whole.mean 20.0647 1 7.487e-06 ***
## alive 0.0000 1 0.9996
## round 0
## qro 9.5471 6 0.1451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$live_pupae)
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.5 4.01
## 2 2 4.58 4.66
## 3 3 4.83 5.61
## 4 4 3.17 2.89
## 5 5 2.83 2.763.0.3.2 Dead Pupae
dp.gn <- glm(dead_pupae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(dp.gn)##
## Call:
## glm(formula = dead_pupae ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0791 -1.5266 -0.6266 0.4029 7.5972
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.51552 0.59407 -2.551 0.010738 *
## treatment2 0.62634 0.17804 3.518 0.000435 ***
## treatment3 0.48512 0.17219 2.817 0.004842 **
## treatment4 0.16071 0.18843 0.853 0.393726
## treatment5 0.40674 0.18732 2.171 0.029909 *
## whole.mean 4.18839 0.41492 10.094 < 2e-16 ***
## alive 0.12996 0.09467 1.373 0.169829
## round2 0.16160 0.37498 0.431 0.666507
## qroB3 -0.37790 0.21913 -1.724 0.084619 .
## qroB4 -0.69025 0.23674 -2.916 0.003550 **
## qroB5 -0.68194 0.23096 -2.953 0.003150 **
## qroK1 0.29344 0.37996 0.772 0.439941
## qroK2/K1 -0.02194 0.55718 -0.039 0.968588
## qroK3 -16.24971 1275.75392 -0.013 0.989837
## 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: 543.06 on 59 degrees of freedom
## Residual deviance: 206.79 on 46 degrees of freedom
## AIC: 408.07
##
## Number of Fisher Scoring iterations: 13dp.gnb <- glm.nb(dead_pupae ~ treatment + whole.mean + alive + round + qro, data = brood) # stick with this model overall
summary(dp.gnb)##
## Call:
## glm.nb(formula = dead_pupae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 2.499875077, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2392 -0.9125 -0.4007 0.2642 3.0123
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.912e+00 1.023e+00 -1.869 0.0616 .
## treatment2 6.805e-01 3.485e-01 1.952 0.0509 .
## treatment3 7.938e-01 3.497e-01 2.270 0.0232 *
## treatment4 1.759e-02 3.763e-01 0.047 0.9627
## treatment5 4.045e-01 3.592e-01 1.126 0.2602
## whole.mean 3.822e+00 7.303e-01 5.233 1.66e-07 ***
## alive 1.905e-01 1.369e-01 1.392 0.1640
## round2 3.614e-01 7.688e-01 0.470 0.6383
## qroB3 -6.170e-01 4.487e-01 -1.375 0.1692
## qroB4 -4.554e-01 4.582e-01 -0.994 0.3203
## qroB5 -4.893e-01 3.684e-01 -1.328 0.1841
## qroK1 6.435e-01 7.872e-01 0.817 0.4137
## qroK2/K1 1.256e-01 1.006e+00 0.125 0.9006
## qroK3 -3.587e+01 2.810e+07 0.000 1.0000
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.4999) family taken to be 1)
##
## Null deviance: 160.803 on 59 degrees of freedom
## Residual deviance: 64.766 on 46 degrees of freedom
## AIC: 329.55
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.500
## Std. Err.: 0.730
##
## 2 x log-likelihood: -299.551AIC(dp.gn, dp.gnb)## df AIC
## dp.gn 14 408.0690
## dp.gnb 15 329.5505Anova(dp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 8.3369 4 0.07999 .
## whole.mean 27.4054 1 1.65e-07 ***
## alive 1.9372 1 0.16397
## round 0
## qro 8.4882 6 0.20447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$dead_pupae)
em.dp <- emmeans(dp.gnb, "treatment")
summary(em.dp)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 -3.56 3512544 Inf -6884463 6884455
## 2 -2.88 3512544 Inf -6884462 6884456
## 3 -2.76 3512544 Inf -6884462 6884456
## 4 -3.54 3512544 Inf -6884463 6884455
## 5 -3.15 3512544 Inf -6884462 6884456
##
## Results are averaged over the levels of: qro, round
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95pairs(em.dp)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.6805 0.349 Inf -1.952 0.2897
## treatment1 - treatment3 -0.7938 0.350 Inf -2.270 0.1546
## treatment1 - treatment4 -0.0176 0.376 Inf -0.047 1.0000
## treatment1 - treatment5 -0.4045 0.359 Inf -1.126 0.7929
## treatment2 - treatment3 -0.1133 0.320 Inf -0.354 0.9966
## treatment2 - treatment4 0.6629 0.351 Inf 1.891 0.3220
## treatment2 - treatment5 0.2760 0.326 Inf 0.846 0.9163
## treatment3 - treatment4 0.7762 0.343 Inf 2.261 0.1578
## treatment3 - treatment5 0.3894 0.334 Inf 1.167 0.7702
## treatment4 - treatment5 -0.3869 0.364 Inf -1.063 0.8257
##
## 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 estimatesdp.gnb1 <- glm.nb(dead_pupae ~ treatment + whole.mean + qro, data = brood)
dp.gnb2 <- glm.nb(dead_pupae ~ treatment + whole.mean + qro + alive, data = brood)
dp.gnb3 <- glm.nb(dead_pupae ~ treatment + whole.mean + qro + round, data = brood)
AIC(dp.gnb, dp.gnb1, dp.gnb2, dp.gnb3)## df AIC
## dp.gnb 15 329.5505
## dp.gnb1 14 329.4804
## dp.gnb2 15 329.5505
## dp.gnb3 14 329.4804summary(dp.gnb1)##
## Call:
## glm.nb(formula = dead_pupae ~ treatment + whole.mean + qro, data = brood,
## init.theta = 2.437314906, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2125 -0.9343 -0.2655 0.3687 2.9948
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.213e-01 4.252e-01 -1.931 0.0534 .
## treatment2 7.324e-01 3.460e-01 2.117 0.0343 *
## treatment3 8.525e-01 3.483e-01 2.447 0.0144 *
## treatment4 4.607e-02 3.742e-01 0.123 0.9020
## treatment5 4.718e-01 3.554e-01 1.327 0.1844
## whole.mean 4.134e+00 6.962e-01 5.938 2.88e-09 ***
## qroB3 -6.883e-01 4.462e-01 -1.542 0.1230
## qroB4 -6.668e-01 4.303e-01 -1.550 0.1212
## qroB5 -7.433e-01 3.336e-01 -2.228 0.0258 *
## qroK1 2.256e-01 3.056e-01 0.738 0.4605
## qroK2/K1 -2.457e-01 6.659e-01 -0.369 0.7121
## qroK3 -3.516e+01 1.704e+07 0.000 1.0000
## qroK3/K2 -3.800e-01 7.762e-01 -0.490 0.6244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.4373) family taken to be 1)
##
## Null deviance: 158.343 on 59 degrees of freedom
## Residual deviance: 65.828 on 47 degrees of freedom
## AIC: 329.48
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.437
## Std. Err.: 0.715
##
## 2 x log-likelihood: -301.480Anova(dp.gnb1)## Analysis of Deviance Table (Type II tests)
##
## Response: dead_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 9.342 4 0.05309 .
## whole.mean 38.143 1 6.575e-10 ***
## qro 14.148 7 0.04862 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(dp.gnb1)
plot(dp.gnb)
plot(dp.gnb3)
em.dp3 <- emmeans(dp.gnb3, "treatment")
summary(em.dp3)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 -3.44 2130465 Inf -4175639 4175632
## 2 -2.70 2130465 Inf -4175638 4175633
## 3 -2.58 2130465 Inf -4175638 4175633
## 4 -3.39 2130465 Inf -4175639 4175632
## 5 -2.96 2130465 Inf -4175638 4175632
##
## Results are averaged over the levels of: qro, round
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95pairs(em.dp3)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.7324 0.346 Inf -2.117 0.2127
## treatment1 - treatment3 -0.8525 0.348 Inf -2.447 0.1029
## treatment1 - treatment4 -0.0461 0.374 Inf -0.123 0.9999
## treatment1 - treatment5 -0.4718 0.355 Inf -1.327 0.6741
## treatment2 - treatment3 -0.1201 0.323 Inf -0.372 0.9959
## treatment2 - treatment4 0.6863 0.350 Inf 1.962 0.2850
## treatment2 - treatment5 0.2606 0.328 Inf 0.795 0.9321
## treatment3 - treatment4 0.8064 0.345 Inf 2.339 0.1327
## treatment3 - treatment5 0.3807 0.337 Inf 1.129 0.7913
## treatment4 - treatment5 -0.4258 0.363 Inf -1.172 0.7676
##
## 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 estimatesdp.sum <- brood %>%
group_by(treatment) %>%
summarise(dp.m= mean(dead_pupae),
sd.dp = sd(dead_pupae))
dp.sum## # A tibble: 5 × 3
## treatment dp.m sd.dp
## <fct> <dbl> <dbl>
## 1 1 4.25 6.81
## 2 2 7.58 10.7
## 3 3 10.5 7.48
## 4 4 6.08 12.2
## 5 5 5.67 6.623.0.3.3 All Pupae
brood$all_pupae <- brood$dead_pupae + brood$live_pupae
ap.gn <- glm(all_pupae ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(ap.gn)##
## Call:
## glm(formula = all_pupae ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3262 -1.3312 -0.3055 0.7122 6.7655
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.31164 0.38079 0.818 0.413127
## treatment2 0.40400 0.13060 3.093 0.001979 **
## treatment3 0.22836 0.12904 1.770 0.076769 .
## treatment4 -0.12852 0.14228 -0.903 0.366373
## treatment5 0.12159 0.14209 0.856 0.392152
## whole.mean 4.01723 0.31612 12.708 < 2e-16 ***
## alive -0.05209 0.05323 -0.978 0.327832
## round2 0.04969 0.28539 0.174 0.861771
## qroB3 -0.53741 0.18074 -2.973 0.002946 **
## qroB4 -0.59343 0.16888 -3.514 0.000441 ***
## qroB5 -0.44678 0.15774 -2.832 0.004619 **
## qroK1 -0.04431 0.29048 -0.153 0.878756
## qroK2/K1 -0.69107 0.50014 -1.382 0.167052
## qroK3 -0.12098 0.58529 -0.207 0.836242
## 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: 555.23 on 59 degrees of freedom
## Residual deviance: 199.27 on 46 degrees of freedom
## AIC: 445.6
##
## Number of Fisher Scoring iterations: 5ap.gnb <- glm.nb(all_pupae ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(ap.gnb)##
## Call:
## glm.nb(formula = all_pupae ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 4.133353061, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4798 -0.7123 -0.1907 0.4315 2.7362
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.17116 0.72721 -0.235 0.814
## treatment2 0.37200 0.25576 1.455 0.146
## treatment3 0.24456 0.26160 0.935 0.350
## treatment4 -0.26340 0.27580 -0.955 0.340
## treatment5 0.02550 0.26685 0.096 0.924
## whole.mean 4.05407 0.55494 7.305 2.76e-13 ***
## alive 0.03691 0.08649 0.427 0.670
## round2 0.10450 0.59509 0.176 0.861
## qroB3 -0.54106 0.33283 -1.626 0.104
## qroB4 -0.52807 0.34172 -1.545 0.122
## qroB5 -0.15484 0.26160 -0.592 0.554
## qroK1 0.02531 0.61006 0.041 0.967
## qroK2/K1 -0.63435 0.82236 -0.771 0.440
## qroK3 0.03720 0.94252 0.039 0.969
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1334) family taken to be 1)
##
## Null deviance: 171.290 on 59 degrees of freedom
## Residual deviance: 66.918 on 46 degrees of freedom
## AIC: 380.85
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.13
## Std. Err.: 1.19
##
## 2 x log-likelihood: -350.848AIC(ap.gn, ap.gnb)## df AIC
## ap.gn 14 445.5981
## ap.gnb 15 380.8478Anova(ap.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: all_pupae
## LR Chisq Df Pr(>Chisq)
## treatment 6.817 4 0.1459
## whole.mean 54.397 1 1.638e-13 ***
## alive 0.181 1 0.6707
## round 0
## qro 5.873 6 0.4376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$all_pupae)
ap.sum <- brood %>%
group_by(treatment) %>%
summarise(al.m= mean(all_pupae),
al.n = length(all_pupae),
sd.al = sd(all_pupae))
ap.sum## # A tibble: 5 × 4
## treatment al.m al.n sd.al
## <fct> <dbl> <int> <dbl>
## 1 1 8.75 12 8.32
## 2 2 12.2 12 12.6
## 3 3 15.3 12 9.70
## 4 4 9.25 12 13.9
## 5 5 8.5 12 7.833.1 cbind pupae
pupae.mortality <- glm(cbind(live_pupae, dead_pupae) ~ treatment + whole.mean + round + qro, data = brood, family = binomial("logit"))
summary(pupae.mortality)##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ treatment + whole.mean +
## round + qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7164 -1.2193 0.1324 1.0376 3.6265
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2572 0.7052 1.783 0.074623 .
## treatment2 -0.8958 0.2839 -3.155 0.001605 **
## treatment3 -1.0185 0.2774 -3.671 0.000241 ***
## treatment4 -0.7371 0.3068 -2.403 0.016274 *
## treatment5 -1.1075 0.3132 -3.536 0.000406 ***
## whole.mean -0.8962 0.6727 -1.332 0.182795
## round2 -0.3895 0.5811 -0.670 0.502617
## qroB3 -0.5404 0.4114 -1.314 0.188950
## qroB4 0.7982 0.3291 2.425 0.015294 *
## qroB5 0.8355 0.2936 2.846 0.004430 **
## qroK1 -1.1677 0.6037 -1.934 0.053092 .
## qroK2/K1 -17.1748 1340.4106 -0.013 0.989777
## qroK3 17.1987 1789.1830 0.010 0.992330
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 244.71 on 55 degrees of freedom
## Residual deviance: 163.56 on 43 degrees of freedom
## AIC: 290.85
##
## Number of Fisher Scoring iterations: 15Anova(pupae.mortality)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## treatment 18.135 4 0.001161 **
## whole.mean 1.768 1 0.183612
## round 0
## qro 35.491 6 3.461e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(pupae.mortality)
lm.emm <- emmeans(pupae.mortality, "treatment")
pairs(lm.emm)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.896 0.284 Inf 3.155 0.0139
## treatment1 - treatment3 1.019 0.277 Inf 3.671 0.0022
## treatment1 - treatment4 0.737 0.307 Inf 2.403 0.1145
## treatment1 - treatment5 1.107 0.313 Inf 3.536 0.0037
## treatment2 - treatment3 0.123 0.259 Inf 0.474 0.9897
## treatment2 - treatment4 -0.159 0.297 Inf -0.534 0.9838
## treatment2 - treatment5 0.212 0.297 Inf 0.712 0.9537
## treatment3 - treatment4 -0.281 0.289 Inf -0.973 0.8676
## treatment3 - treatment5 0.089 0.297 Inf 0.299 0.9983
## treatment4 - treatment5 0.370 0.330 Inf 1.122 0.7952
##
## Results are averaged over the levels of: qro, round
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatessummary(lm.emm)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 0.603 279 Inf -547 548
## 2 -0.293 279 Inf -548 547
## 3 -0.416 279 Inf -548 547
## 4 -0.135 279 Inf -548 548
## 5 -0.505 279 Inf -548 547
##
## Results are averaged over the levels of: qro, round
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95plot(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.646 63.9 Inf 2.22e-16 1
## 2 0.427 68.4 Inf 2.22e-16 1
## 3 0.397 66.9 Inf 2.22e-16 1
## 4 0.466 69.5 Inf 2.22e-16 1
## 5 0.376 65.6 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## 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 2.449 0.695 Inf 1 3.155 0.0139
## treatment1 / treatment3 2.769 0.768 Inf 1 3.671 0.0022
## treatment1 / treatment4 2.090 0.641 Inf 1 2.403 0.1145
## treatment1 / treatment5 3.027 0.948 Inf 1 3.536 0.0037
## treatment2 / treatment3 1.131 0.293 Inf 1 0.474 0.9897
## treatment2 / treatment4 0.853 0.253 Inf 1 -0.534 0.9838
## treatment2 / treatment5 1.236 0.367 Inf 1 0.712 0.9537
## treatment3 / treatment4 0.755 0.218 Inf 1 -0.973 0.8676
## treatment3 / treatment5 1.093 0.325 Inf 1 0.299 0.9983
## treatment4 / treatment5 1.448 0.478 Inf 1 1.122 0.7952
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemm2 <- emmeans(pupae.mortality, ~ treatment, type = "response")
emm2## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.646 63.9 Inf 2.22e-16 1
## 2 0.427 68.4 Inf 2.22e-16 1
## 3 0.397 66.9 Inf 2.22e-16 1
## 4 0.466 69.5 Inf 2.22e-16 1
## 5 0.376 65.6 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalepairs(emm2)## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 2.449 0.695 Inf 1 3.155 0.0139
## treatment1 / treatment3 2.769 0.768 Inf 1 3.671 0.0022
## treatment1 / treatment4 2.090 0.641 Inf 1 2.403 0.1145
## treatment1 / treatment5 3.027 0.948 Inf 1 3.536 0.0037
## treatment2 / treatment3 1.131 0.293 Inf 1 0.474 0.9897
## treatment2 / treatment4 0.853 0.253 Inf 1 -0.534 0.9838
## treatment2 / treatment5 1.236 0.367 Inf 1 0.712 0.9537
## treatment3 / treatment4 0.755 0.218 Inf 1 -0.973 0.8676
## treatment3 / treatment5 1.093 0.325 Inf 1 0.299 0.9983
## treatment4 / treatment5 1.448 0.478 Inf 1 1.122 0.7952
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemmdf2 <- 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)
3.1.1 Dead Larvae + Dead Pupae
brood$all_d_lp <- brood$dead_larvae + brood$dead_pupae
dlp.gn <- glm(all_d_lp ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(dlp.gn)##
## Call:
## glm(formula = all_d_lp ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.8217 -1.7824 -0.7181 1.0765 6.9007
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.36218 0.44121 -3.087 0.00202 **
## treatment2 0.58775 0.14609 4.023 5.74e-05 ***
## treatment3 0.54439 0.13919 3.911 9.18e-05 ***
## treatment4 0.38567 0.14907 2.587 0.00967 **
## treatment5 0.26285 0.15704 1.674 0.09417 .
## whole.mean 4.11687 0.31851 12.926 < 2e-16 ***
## alive 0.22157 0.07079 3.130 0.00175 **
## round2 -0.11289 0.27760 -0.407 0.68426
## qroB3 -0.24441 0.17381 -1.406 0.15965
## qroB4 -0.17750 0.17837 -0.995 0.31968
## qroB5 0.25355 0.14220 1.783 0.07458 .
## qroK1 0.06237 0.28054 0.222 0.82405
## qroK2/K1 0.09320 0.41196 0.226 0.82101
## qroK3 -16.05780 773.78389 -0.021 0.98344
## 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: 778.56 on 59 degrees of freedom
## Residual deviance: 250.10 on 46 degrees of freedom
## AIC: 487.03
##
## Number of Fisher Scoring iterations: 12dlp.gnb <- glm.nb(all_d_lp ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(dlp.gnb)##
## Call:
## glm.nb(formula = all_d_lp ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 2.681171311, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5978 -0.8191 -0.3884 0.4920 2.4503
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.302e+00 9.027e-01 -1.442 0.1492
## treatment2 5.751e-01 3.119e-01 1.844 0.0652 .
## treatment3 7.172e-01 3.135e-01 2.288 0.0222 *
## treatment4 2.883e-01 3.243e-01 0.889 0.3739
## treatment5 1.549e-01 3.242e-01 0.478 0.6329
## whole.mean 3.486e+00 6.423e-01 5.428 5.69e-08 ***
## alive 2.490e-01 1.152e-01 2.162 0.0306 *
## round2 7.228e-02 7.039e-01 0.103 0.9182
## qroB3 -5.338e-01 4.001e-01 -1.334 0.1821
## qroB4 6.187e-03 4.067e-01 0.015 0.9879
## qroB5 1.693e-01 3.107e-01 0.545 0.5857
## qroK1 3.220e-01 7.212e-01 0.446 0.6553
## qroK2/K1 4.583e-02 9.086e-01 0.050 0.9598
## qroK3 -3.865e+01 6.711e+07 0.000 1.0000
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.6812) family taken to be 1)
##
## Null deviance: 169.560 on 59 degrees of freedom
## Residual deviance: 63.518 on 46 degrees of freedom
## AIC: 380.55
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.681
## Std. Err.: 0.693
##
## 2 x log-likelihood: -350.552AIC(dlp.gn, dlp.gnb)## df AIC
## dlp.gn 14 487.0289
## dlp.gnb 15 380.5516Anova(dlp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: all_d_lp
## LR Chisq Df Pr(>Chisq)
## treatment 6.6974 4 0.15277
## whole.mean 28.8833 1 7.687e-08 ***
## alive 4.4189 1 0.03554 *
## round 0
## qro 9.1072 6 0.16764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$all_d_lp)
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 6.33 12 8.78
## 2 2 11.2 12 11.7
## 3 3 17.2 12 12.3
## 4 4 12.6 12 21.0
## 5 5 7.83 12 9.013.1.2 Live Larvae + Live Pupae
brood$all_a_lp <- brood$live_larvae + brood$live_pupae
alp.gn <- glm(all_a_lp ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(alp.gn)##
## Call:
## glm(formula = all_a_lp ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.3761 -2.3123 -0.2949 1.4448 5.1001
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.55954 0.31658 1.767 0.07715 .
## treatment2 -0.24675 0.08325 -2.964 0.00304 **
## treatment3 0.01915 0.07545 0.254 0.79962
## treatment4 -0.55043 0.08885 -6.195 5.81e-10 ***
## treatment5 -0.25304 0.08731 -2.898 0.00375 **
## whole.mean 3.60038 0.18756 19.196 < 2e-16 ***
## alive -0.02568 0.02828 -0.908 0.36384
## round2 0.99892 0.28470 3.509 0.00045 ***
## qroB3 -0.48224 0.11877 -4.060 4.90e-05 ***
## qroB4 -0.18153 0.10071 -1.802 0.07148 .
## qroB5 0.39841 0.08235 4.838 1.31e-06 ***
## qroK1 0.77775 0.28686 2.711 0.00670 **
## qroK2/K1 0.04520 0.43965 0.103 0.91812
## qroK3 1.26407 0.40368 3.131 0.00174 **
## 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: 1124.81 on 59 degrees of freedom
## Residual deviance: 377.25 on 46 degrees of freedom
## AIC: 665.14
##
## Number of Fisher Scoring iterations: 5alp.gnb <- glm.nb(all_a_lp ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(alp.gnb)##
## Call:
## glm.nb(formula = all_a_lp ~ treatment + whole.mean + alive +
## round + qro, data = brood, init.theta = 2.914593471, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.99638 -0.99952 -0.09627 0.37375 2.38290
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.01494 0.80474 -1.261 0.207234
## treatment2 -0.20063 0.26827 -0.748 0.454534
## treatment3 -0.11769 0.27296 -0.431 0.666335
## treatment4 -0.46884 0.28003 -1.674 0.094076 .
## treatment5 -0.32096 0.27577 -1.164 0.244473
## whole.mean 4.81474 0.57426 8.384 < 2e-16 ***
## alive 0.17394 0.08754 1.987 0.046924 *
## round2 0.99902 0.68300 1.463 0.143554
## qroB3 -0.57428 0.34320 -1.673 0.094266 .
## qroB4 -0.38758 0.36335 -1.067 0.286116
## qroB5 0.90092 0.25970 3.469 0.000522 ***
## qroK1 0.57577 0.69913 0.824 0.410199
## qroK2/K1 0.07779 0.89400 0.087 0.930661
## qroK3 1.27987 0.96413 1.327 0.184347
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.9146) family taken to be 1)
##
## Null deviance: 179.48 on 59 degrees of freedom
## Residual deviance: 75.34 on 46 degrees of freedom
## AIC: 482.89
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.915
## Std. Err.: 0.705
##
## 2 x log-likelihood: -452.887AIC(alp.gn, alp.gnb)## df AIC
## alp.gn 14 665.1389
## alp.gnb 15 482.8868Anova(alp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: all_a_lp
## LR Chisq Df Pr(>Chisq)
## treatment 3.347 4 0.501558
## whole.mean 59.909 1 9.934e-15 ***
## alive 4.079 1 0.043429 *
## round 0
## qro 22.380 6 0.001033 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$all_a_lp)
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 29.3 12 26.9
## 2 2 21.7 12 16.3
## 3 3 37 12 25.3
## 4 4 21.3 12 18.5
## 5 5 20 12 14.83.1.3 All Larvae + All Pupae
brood$all_lp <- brood$all_larvae + brood$all_pupae
lp.gn <- glm(all_lp ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(lp.gn)##
## Call:
## glm(formula = all_lp ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.172 -1.709 0.000 1.294 4.334
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.92801 0.23341 3.976 7.01e-05 ***
## treatment2 -0.03244 0.07116 -0.456 0.64848
## treatment3 0.13823 0.06558 2.108 0.03506 *
## treatment4 -0.29525 0.07446 -3.965 7.33e-05 ***
## treatment5 -0.14057 0.07560 -1.859 0.06297 .
## whole.mean 3.77499 0.16112 23.429 < 2e-16 ***
## alive 0.00975 0.02595 0.376 0.70710
## round2 0.55386 0.19664 2.817 0.00485 **
## qroB3 -0.41125 0.09786 -4.202 2.64e-05 ***
## qroB4 -0.21264 0.08727 -2.437 0.01483 *
## qroB5 0.35323 0.07123 4.959 7.09e-07 ***
## qroK1 0.45997 0.19840 2.318 0.02043 *
## qroK2/K1 0.07311 0.29894 0.245 0.80680
## qroK3 0.39439 0.34459 1.145 0.25240
## 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: 1407.73 on 59 degrees of freedom
## Residual deviance: 264.85 on 46 degrees of freedom
## AIC: 588.98
##
## Number of Fisher Scoring iterations: 5lp.gnb <- glm.nb(all_lp ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(lp.gnb)##
## Call:
## glm.nb(formula = all_lp ~ treatment + whole.mean + alive + round +
## qro, data = brood, init.theta = 7.881729479, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.55614 -0.96646 -0.05738 0.47429 2.58588
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.05116 0.50993 0.100 0.92008
## treatment2 -0.03587 0.17223 -0.208 0.83503
## treatment3 0.10346 0.17401 0.595 0.55214
## treatment4 -0.33515 0.18078 -1.854 0.06375 .
## treatment5 -0.21804 0.17789 -1.226 0.22031
## whole.mean 4.26660 0.36816 11.589 < 2e-16 ***
## alive 0.13432 0.05753 2.335 0.01955 *
## round2 0.62263 0.42577 1.462 0.14364
## qroB3 -0.55062 0.22433 -2.455 0.01411 *
## qroB4 -0.25356 0.23004 -1.102 0.27036
## qroB5 0.60185 0.16761 3.591 0.00033 ***
## qroK1 0.39800 0.43583 0.913 0.36114
## qroK2/K1 0.09522 0.55933 0.170 0.86482
## qroK3 0.49435 0.63486 0.779 0.43617
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(7.8817) family taken to be 1)
##
## Null deviance: 308.449 on 59 degrees of freedom
## Residual deviance: 70.397 on 46 degrees of freedom
## AIC: 487.75
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 7.88
## Std. Err.: 2.02
##
## 2 x log-likelihood: -457.748AIC(lp.gn, lp.gnb)## df AIC
## lp.gn 14 588.9827
## lp.gnb 15 487.7479Anova(lp.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: all_lp
## LR Chisq Df Pr(>Chisq)
## treatment 7.977 4 0.09244 .
## whole.mean 128.847 1 < 2.2e-16 ***
## alive 5.717 1 0.01680 *
## round 0
## qro 28.109 6 8.961e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(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.7 12 30.0
## 2 2 32.9 12 21.7
## 3 3 54.2 12 33.6
## 4 4 33.9 12 36.0
## 5 5 27.8 12 20.6em.lp <- emmeans(lp.gnb, "treatment")
pairs(em.lp)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.0359 0.172 Inf 0.208 0.9996
## treatment1 - treatment3 -0.1035 0.174 Inf -0.595 0.9759
## treatment1 - treatment4 0.3352 0.181 Inf 1.854 0.3426
## treatment1 - treatment5 0.2180 0.178 Inf 1.226 0.7364
## treatment2 - treatment3 -0.1393 0.169 Inf -0.822 0.9239
## treatment2 - treatment4 0.2993 0.178 Inf 1.684 0.4440
## treatment2 - treatment5 0.1822 0.172 Inf 1.059 0.8274
## treatment3 - treatment4 0.4386 0.176 Inf 2.498 0.0911
## treatment3 - treatment5 0.3215 0.174 Inf 1.843 0.3489
## treatment4 - treatment5 -0.1171 0.184 Inf -0.637 0.9691
##
## 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##cbind larvae and pupae
pl_bind <- glm(cbind(all_a_lp, all_d_lp) ~ treatment + whole.mean + round + qro, data = brood, family = binomial("logit"))
summary(pl_bind )##
## Call:
## glm(formula = cbind(all_a_lp, all_d_lp) ~ treatment + whole.mean +
## round + qro, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.8316 -1.2499 0.0007 2.0115 4.1015
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.043243 0.463012 2.253 0.024249 *
## treatment2 -0.854208 0.170465 -5.011 5.41e-07 ***
## treatment3 -0.606427 0.158353 -3.830 0.000128 ***
## treatment4 -0.917245 0.172442 -5.319 1.04e-07 ***
## treatment5 -0.630907 0.179162 -3.521 0.000429 ***
## whole.mean -0.987356 0.400339 -2.466 0.013652 *
## round2 1.094157 0.397580 2.752 0.005923 **
## qroB3 -0.185856 0.219186 -0.848 0.396474
## qroB4 0.316627 0.198336 1.596 0.110395
## qroB5 0.265745 0.149028 1.783 0.074556 .
## qroK1 0.809957 0.403763 2.006 0.044854 *
## qroK2/K1 -0.006826 0.606105 -0.011 0.991014
## qroK3 16.616994 598.418649 0.028 0.977847
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 495.64 on 57 degrees of freedom
## Residual deviance: 384.05 on 45 degrees of freedom
## AIC: 582.73
##
## Number of Fisher Scoring iterations: 13Anova(pl_bind )## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(all_a_lp, all_d_lp)
## LR Chisq Df Pr(>Chisq)
## treatment 37.919 4 1.165e-07 ***
## whole.mean 6.209 1 0.0127086 *
## round 0
## qro 26.661 6 0.0001676 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(pl_bind )
lm.emm <- emmeans(pl_bind , "treatment")
pairs(lm.emm)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.8542 0.170 Inf 5.011 <.0001
## treatment1 - treatment3 0.6064 0.158 Inf 3.830 0.0012
## treatment1 - treatment4 0.9172 0.172 Inf 5.319 <.0001
## treatment1 - treatment5 0.6309 0.179 Inf 3.521 0.0039
## treatment2 - treatment3 -0.2478 0.149 Inf -1.668 0.4541
## treatment2 - treatment4 0.0630 0.166 Inf 0.379 0.9956
## treatment2 - treatment5 -0.2233 0.166 Inf -1.343 0.6644
## treatment3 - treatment4 0.3108 0.145 Inf 2.150 0.1990
## treatment3 - treatment5 0.0245 0.161 Inf 0.152 0.9999
## treatment4 - treatment5 -0.2863 0.178 Inf -1.607 0.4927
##
## Results are averaged over the levels of: qro, round
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatessummary(lm.emm)## treatment emmean SE df asymp.LCL asymp.UCL
## 1 3.32 74.8 Inf -143 150
## 2 2.46 74.8 Inf -144 149
## 3 2.71 74.8 Inf -144 149
## 4 2.40 74.8 Inf -144 149
## 5 2.69 74.8 Inf -144 149
##
## Results are averaged over the levels of: qro, round
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95plot(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.965 2.52 Inf 2.22e-16 1
## 2 0.922 5.41 Inf 2.22e-16 1
## 3 0.938 4.37 Inf 2.22e-16 1
## 4 0.917 5.70 Inf 2.22e-16 1
## 5 0.936 4.47 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## 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 2.350 0.401 Inf 1 5.011 <.0001
## treatment1 / treatment3 1.834 0.290 Inf 1 3.830 0.0012
## treatment1 / treatment4 2.502 0.432 Inf 1 5.319 <.0001
## treatment1 / treatment5 1.879 0.337 Inf 1 3.521 0.0039
## treatment2 / treatment3 0.781 0.116 Inf 1 -1.668 0.4541
## treatment2 / treatment4 1.065 0.177 Inf 1 0.379 0.9956
## treatment2 / treatment5 0.800 0.133 Inf 1 -1.343 0.6644
## treatment3 / treatment4 1.365 0.197 Inf 1 2.150 0.1990
## treatment3 / treatment5 1.025 0.165 Inf 1 0.152 0.9999
## treatment4 / treatment5 0.751 0.134 Inf 1 -1.607 0.4927
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemm3 <- emmeans(pl_bind , ~ treatment, type = "response")
emm3## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.965 2.52 Inf 2.22e-16 1
## 2 0.922 5.41 Inf 2.22e-16 1
## 3 0.938 4.37 Inf 2.22e-16 1
## 4 0.917 5.70 Inf 2.22e-16 1
## 5 0.936 4.47 Inf 2.22e-16 1
##
## Results are averaged over the levels of: qro, round
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scalepairs(emm3)## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 2.350 0.401 Inf 1 5.011 <.0001
## treatment1 / treatment3 1.834 0.290 Inf 1 3.830 0.0012
## treatment1 / treatment4 2.502 0.432 Inf 1 5.319 <.0001
## treatment1 / treatment5 1.879 0.337 Inf 1 3.521 0.0039
## treatment2 / treatment3 0.781 0.116 Inf 1 -1.668 0.4541
## treatment2 / treatment4 1.065 0.177 Inf 1 0.379 0.9956
## treatment2 / treatment5 0.800 0.133 Inf 1 -1.343 0.6644
## treatment3 / treatment4 1.365 0.197 Inf 1 2.150 0.1990
## treatment3 / treatment5 1.025 0.165 Inf 1 0.152 0.9999
## treatment4 / treatment5 0.751 0.134 Inf 1 -1.607 0.4927
##
## Results are averaged over the levels of: qro, round
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleemmdf3 <- 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,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)
3.1.4 Eggs
eggs.gn <- glm(eggs ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(eggs.gn)##
## Call:
## glm(formula = eggs ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7087 -1.7080 -0.0778 1.6872 10.6744
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.468e+15 2.017e+14 7.280 3.33e-13 ***
## treatment2 -7.564e-02 1.540e-01 -0.491 0.62339
## treatment3 -2.261e-02 1.976e-01 -0.114 0.90890
## treatment4 -6.673e-01 1.676e-01 -3.982 6.84e-05 ***
## treatment5 -3.445e-01 1.940e-01 -1.776 0.07580 .
## whole.mean 2.118e+00 4.305e-01 4.920 8.67e-07 ***
## alive -1.396e-01 5.273e-02 -2.647 0.00812 **
## round2 -1.468e+15 2.017e+14 -7.280 3.33e-13 ***
## qroB3 9.282e-01 2.228e-01 4.166 3.11e-05 ***
## qroB4 1.305e+00 1.997e-01 6.534 6.41e-11 ***
## qroB5 8.729e-01 2.030e-01 4.300 1.71e-05 ***
## qroK1 -1.468e+15 2.017e+14 -7.280 3.33e-13 ***
## qroK2/K1 -1.468e+15 2.017e+14 -7.280 3.33e-13 ***
## qroK3 -1.468e+15 2.017e+14 -7.280 3.33e-13 ***
## qroK3/K2 -1.468e+15 2.017e+14 -7.280 3.33e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 730.50 on 59 degrees of freedom
## Residual deviance: 432.95 on 45 degrees of freedom
## AIC: 627.89
##
## Number of Fisher Scoring iterations: 25#eggs.gnb <- glm.nb(eggs ~ treatment + whole.mean + alive + round + qro, data = brood)
#summary(eggs.gnb)
#AIC(eggs.gn, eggs.gnb) #### I think we have too many zeros in eggs, I will switch to a zero-inflated negative binomial model
library(glmmTMB)
egg.zero <- glmmTMB(eggs ~ treatment + whole.mean + round + 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 + round + qro
## Data: brood
##
## AIC BIC logLik deviance df.resid
## 355.8 387.2 -162.9 325.8 45
##
##
## Dispersion parameter for nbinom2 family (): 1.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.089e+00 3.416e+05 0.000 0.999983
## treatment2 -1.134e-01 4.612e-01 -0.246 0.805714
## treatment3 -6.458e-01 4.679e-01 -1.380 0.167491
## treatment4 -6.436e-01 4.620e-01 -1.393 0.163570
## treatment5 -6.235e-01 4.626e-01 -1.348 0.177760
## whole.mean 3.478e+00 9.668e-01 3.597 0.000322 ***
## round2 7.072e+00 3.416e+05 0.000 0.999983
## qroB3 7.363e-01 5.092e-01 1.446 0.148187
## qroB4 1.022e+00 5.510e-01 1.855 0.063644 .
## qroB5 1.018e+00 4.124e-01 2.469 0.013567 *
## qroK1 7.433e+00 3.416e+05 0.000 0.999983
## qroK2/K1 5.381e+00 3.416e+05 0.000 0.999987
## qroK3 -1.186e+01 3.417e+05 0.000 0.999972
## qroK3/K2 -1.511e+01 3.430e+05 0.000 0.999965
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(egg.zero)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: eggs
## Chisq Df Pr(>Chisq)
## treatment 3.7776 4 0.4369391
## whole.mean 12.9384 1 0.0003219 ***
## round 0.0000 1 0.9999835
## qro 10.7811 7 0.1484564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(brood$treatment, brood$eggs)
eggs.sum <- brood %>%
group_by(treatment) %>%
summarise(eggs.m= mean(eggs),
eggs.n = length(eggs),
eggs.al = sd(eggs))
eggs.sum## # A tibble: 5 × 4
## treatment eggs.m eggs.n eggs.al
## <fct> <dbl> <int> <dbl>
## 1 1 12 12 24.1
## 2 2 8.83 12 10.4
## 3 3 7.08 12 7.61
## 4 4 5.83 12 5.70
## 5 5 5.58 12 5.58em.egg <- emmeans(egg.zero, "treatment")
pairs(em.egg)## contrast estimate SE df t.ratio p.value
## treatment1 - treatment2 0.11343 0.461 45 0.246 0.9992
## treatment1 - treatment3 0.64580 0.468 45 1.380 0.6431
## treatment1 - treatment4 0.64364 0.462 45 1.393 0.6351
## treatment1 - treatment5 0.62347 0.463 45 1.348 0.6634
## treatment2 - treatment3 0.53237 0.468 45 1.138 0.7855
## treatment2 - treatment4 0.53021 0.473 45 1.120 0.7952
## treatment2 - treatment5 0.51004 0.443 45 1.150 0.7791
## treatment3 - treatment4 -0.00216 0.460 45 -0.005 1.0000
## treatment3 - treatment5 -0.02232 0.456 45 -0.049 1.0000
## treatment4 - treatment5 -0.02016 0.464 45 -0.043 1.0000
##
## 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 estimatessummary(em.egg)## treatment emmean SE df lower.CL upper.CL
## 1 -3.22 4270 45 -8604 8597
## 2 -3.33 4270 45 -8604 8597
## 3 -3.86 4270 45 -8604 8597
## 4 -3.86 4270 45 -8604 8597
## 5 -3.84 4270 45 -8604 8597
##
## Results are averaged over the levels of: qro, round
## Results are given on the log (not the response) scale.
## Confidence level used: 0.953.1.5 Count of drones
drones.gn <- glm(drones ~ treatment + whole.mean + alive + round + qro, data = brood, family = "poisson")
summary(drones.gn)##
## Call:
## glm(formula = drones ~ treatment + whole.mean + alive + round +
## qro, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0090 -1.3704 -0.3482 0.8664 3.3686
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.19678 0.50238 -0.392 0.69528
## treatment2 -0.07222 0.15332 -0.471 0.63762
## treatment3 -0.50177 0.16384 -3.062 0.00220 **
## treatment4 -0.06883 0.15300 -0.450 0.65281
## treatment5 0.07981 0.15583 0.512 0.60853
## whole.mean 2.78858 0.33720 8.270 < 2e-16 ***
## alive 0.17049 0.05571 3.061 0.00221 **
## round2 0.15649 0.43030 0.364 0.71610
## qroB3 0.09282 0.17198 0.540 0.58940
## qroB4 0.04393 0.17703 0.248 0.80400
## qroB5 0.54024 0.13143 4.110 3.95e-05 ***
## qroK1 -1.27937 0.45886 -2.788 0.00530 **
## qroK2/K1 -1.55300 0.82996 -1.871 0.06132 .
## qroK3 -1.68681 1.09433 -1.541 0.12322
## 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: 365.35 on 59 degrees of freedom
## Residual deviance: 122.62 on 46 degrees of freedom
## AIC: 346.89
##
## Number of Fisher Scoring iterations: 6drones.gnb <- glm.nb(drones ~ treatment + whole.mean + alive + round + qro, data = brood)
summary(drones.gnb)##
## Call:
## glm.nb(formula = drones ~ treatment + whole.mean + alive + round +
## qro, data = brood, init.theta = 6.834486502, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3447 -1.0280 -0.1238 0.6151 2.0921
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.95091 0.72245 -1.316 0.188096
## treatment2 -0.04781 0.23611 -0.203 0.839525
## treatment3 -0.54081 0.24672 -2.192 0.028377 *
## treatment4 -0.04672 0.24191 -0.193 0.846841
## treatment5 0.11772 0.24046 0.490 0.624447
## whole.mean 3.36606 0.51865 6.490 8.58e-11 ***
## alive 0.26129 0.08384 3.117 0.001829 **
## round2 0.13867 0.59353 0.234 0.815269
## qroB3 0.11161 0.27056 0.413 0.679958
## qroB4 -0.01581 0.29738 -0.053 0.957589
## qroB5 0.77258 0.21648 3.569 0.000359 ***
## qroK1 -1.34370 0.62211 -2.160 0.030780 *
## qroK2/K1 -1.47693 0.97588 -1.513 0.130172
## qroK3 -1.63656 1.23745 -1.323 0.185994
## qroK3/K2 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(6.8345) family taken to be 1)
##
## Null deviance: 193.657 on 59 degrees of freedom
## Residual deviance: 65.131 on 46 degrees of freedom
## AIC: 331.19
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 6.83
## Std. Err.: 2.57
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -301.187AIC(drones.gn, drones.gnb)## df AIC
## drones.gn 14 346.8916
## drones.gnb 15 331.1866Anova(drones.gnb)## Analysis of Deviance Table (Type II tests)
##
## Response: drones
## LR Chisq Df Pr(>Chisq)
## treatment 8.604 4 0.071796 .
## whole.mean 41.974 1 9.25e-11 ***
## alive 9.869 1 0.001680 **
## round 0
## qro 18.972 6 0.004211 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drone.em <- emmeans(drones.gnb, "treatment")
pairs(drone.em)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.04781 0.236 Inf 0.203 0.9996
## treatment1 - treatment3 0.54081 0.247 Inf 2.192 0.1826
## treatment1 - treatment4 0.04672 0.242 Inf 0.193 0.9997
## treatment1 - treatment5 -0.11772 0.240 Inf -0.490 0.9884
## treatment2 - treatment3 0.49300 0.241 Inf 2.047 0.2435
## treatment2 - treatment4 -0.00109 0.237 Inf -0.005 1.0000
## treatment2 - treatment5 -0.16553 0.231 Inf -0.718 0.9525
## treatment3 - treatment4 -0.49409 0.242 Inf -2.041 0.2463
## treatment3 - treatment5 -0.65853 0.244 Inf -2.702 0.0536
## treatment4 - treatment5 -0.16445 0.242 Inf -0.678 0.9612
##
## 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 estimatesdrones.sum <- brood %>%
group_by(treatment) %>%
summarise(m.d = mean(drones),
sd.d = sd(drones))
drones.sum## # A tibble: 5 × 3
## treatment m.d sd.d
## <fct> <dbl> <dbl>
## 1 1 7.25 7.24
## 2 2 7.33 5.31
## 3 3 6.5 5.70
## 4 4 9.17 9.08
## 5 5 8.17 6.03LS0tDQp0aXRsZTogImJyb29kIg0KYXV0aG9yOiAiRW1pbHkgUnVubmlvbiINCmRhdGU6ICIyMDIzLTAxLTI2Ig0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIHRvYzogdHJ1ZQ0KICAgIHRvY19kZXB0aDogNA0KICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQ0KICAgIHRvY19mbG9hdDogdHJ1ZQ0KICAgIHRoZW1lOiBqb3VybmFsDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KLS0tDQoNCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KHdhcm5pbmcgPSBGQUxTRSwgbWVzc2FnZSA9IEZBTFNFKQ0KYGBgDQoNCmBgYHtyIGxvYWQgbGlicmFyaWVzLCBpbmNsdWRlPUZBTFNFfQ0KbGlicmFyeShyZWFkcikNCmxpYnJhcnkoa2FibGVFeHRyYSkNCmxpYnJhcnkoc3RhdHMpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGNhcikNCmxpYnJhcnkoZW1tZWFucykNCmxpYnJhcnkoTUFTUykNCmxpYnJhcnkobG1lNCkNCmxpYnJhcnkoYmxtZWNvKQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShjb3dwbG90KQ0KbGlicmFyeShiZXN0Tm9ybWFsaXplKQ0KbGlicmFyeShwbG90bHkpDQpsaWJyYXJ5KGFncmljb2xhZSkgDQpsaWJyYXJ5KGdncHVicikNCmxpYnJhcnkoZ2x1ZSkNCmxpYnJhcnkobXVsdGNvbXBWaWV3KQ0KYGBgDQoNCg0KIyBJbnB1dCBQb2xsZW4gRGF0YSANCg0KU3RhcnQgYnkgaW1wdXRpbmcgcG9sbGVuIGRhdGEgYW5kIGNyZWF0aW5nIGEgbmV3IGRhdGEgZnJhbWUgd2l0aCBhdmVyYWdlIHBvbGxlbiBjb25zdW1wdGlvbiBvbiBhIHBlci1jb2xvbnkgYmFzaXMgDQoNCg0KYGBge3J9DQojIyMgRmlndXJlIG91dCBhdmVyYWdlIHBvbGxlbiBjb25zdW1wdGlvbiBieSB0cmVhdG1lbnQgDQoNCg0KcG9sbGVuIDwtIHJlYWRfY3N2KCJwb2xsZW4xLmNzdiIsIGNvbF90eXBlcyA9IGNvbHMocm91bmQgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjEiLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICIyIikpLCB0cmVhdG1lbnQgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjEiLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIjIiLCAiMyIsICI0IiwgIjUiKSksIHJlcGxpY2F0ZSA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMSIsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICIyIiwgIjMiLCAiNCIsICI1IiwgIjYiLCAiNyIsICI5IiwgIjExIiwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIjEyIikpLCBzdGFydF9kYXRlID0gY29sX2RhdGUoZm9ybWF0ID0gIiVtLyVkLyVZIiksIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RhcnRfdGltZSA9IGNvbF9jaGFyYWN0ZXIoKSwgZW5kX2RhdGUgPSBjb2xfZGF0ZShmb3JtYXQgPSAiJW0vJWQvJVkiKSwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBlbmRfdGltZSA9IGNvbF9jaGFyYWN0ZXIoKSkpDQoNCg0KcG9sbGVuJGNvbG9ueSA8LSBhcy5mYWN0b3IocG9sbGVuJGNvbG9ueSkNCnBvbGxlbiRwaWQgPC0gYXMuZmFjdG9yKHBvbGxlbiRwaWQpDQpwb2xsZW4kY291bnQgPC0gYXMuZmFjdG9yKHBvbGxlbiRjb3VudCkNCg0KcG9sbGVuIDwtIG5hLm9taXQocG9sbGVuKQ0KDQpyYW5nZShwb2xsZW4kZGlmZmVyZW5jZSkNCg0KIyBnZXQgcmlkIG9mIG5lZ2F0aXZlIG51bWJlcnMNCnBvbGxlbiRkaWZmZXJlbmNlW3BvbGxlbiRkaWZmZXJlbmNlIDwgMF0gPC0gTkENCnBvbGxlbiA8LSBuYS5vbWl0KHBvbGxlbikNCnJhbmdlKHBvbGxlbiRkaWZmZXJlbmNlKQ0KDQpgYGANCg0KDQoNCg0KIyBBdmVyYWdlIHBvbGxlbiBjb25zdW1wdGlvbiBwZXIgY29sb255DQoNCmBgYHtyfQ0KcG9sbGVuJHdob2xlX2RpZiA8LSBhcy5udW1lcmljKHBvbGxlbiRkaWZmZXJlbmNlKQ0Kd2QuMSA8LSBsbShkaWZmZXJlbmNlIH4gdHJlYXRtZW50LCBkYXRhID0gcG9sbGVuKQ0Kc3VtbWFyeSh3ZC4xKQ0Kd2QuZW1tIDwtIGVtbWVhbnMod2QuMSwgInRyZWF0bWVudCIpDQpzdW1tYXJ5KHdkLmVtbSkNCg0Kd2Quc3VtbWFyeSA8LSBwb2xsZW4gJT4lIA0KICBncm91cF9ieShjb2xvbnkpICU+JQ0KICBzdW1tYXJpemUod2hvbGUubWVhbiA9IG1lYW4oZGlmZmVyZW5jZSkpDQoNCmFzLmRhdGEuZnJhbWUod2Quc3VtbWFyeSkgICMgdGhpcyBpcyB0aGUgZGF0YSBmcmFtZSBJIHdpbGwgbWVyZ2Ugd2l0aCBzdWJzZXF1ZW50IGRhdGEgZnJhbWVzIHRvIGNvbnRhaW4gYXZlcmFnZSBwb2xsZW4gY29uc3VtcHRpb24gcGVyIGNvbG9ueSBhcyBhIG5ldyBjb2x1bW4gIA0KDQpgYGANCg0KDQpgYGB7cn0NCg0KYnJvb2QgPC0gcmVhZF9jc3YoImJyb29kLmNzdiIsIGNvbF90eXBlcyA9IGNvbHMocm91bmQgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjEiLCANCiAgICAiMiIpKSwgdHJlYXRtZW50ID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCIxIiwgDQogICAgIjIiLCAiMyIsICI0IiwgIjUiKSksIHJlcGxpY2F0ZSA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMSIsIA0KICAgICIyIiwgIjMiLCAiNCIsICI1IiwgIjciLCAiOSIsICIxMSIsICIxMiIpKSkpDQoNCmJyb29kJGNvbG9ueSA8LSBhcy5mYWN0b3IoYnJvb2QkY29sb255KQ0KDQpicm9vZCA8LSBtZXJnZSh3ZC5zdW1tYXJ5LCBicm9vZCwgYnkueCA9ICJjb2xvbnkiKSAgICAjIHRoaXMgaXMgZ29vZCBiZWNhdXNlIEkgY2FsY3VsYXRlZCBhdmVyYWdlIHBvbGxlbiBjb25zdW1wdGlvbiB0d28gZGlmZmVyZW50IHdheXMgKGF2Z19wb2xsZW4gd2FzIGNhbGN1bGF0ZWQgYSBjb3VwbGUgbW9udGhzIGFnbyBtYW51YWxseSkgYW5kIGl0J3MgdGhlIHNhbWUgbnVtYmVycyBhcyB0aGUgZ3JvdXBfYnkgZnVuY3Rpb24gdGhhdCBjcmVhdGVkIHdob2xlLm1lYW4NCg0KbW9ydGFsaXR5IDwtIHN1cnZpdmluZ193b3JrZXJzX3Blcl9jb2xvbnkgPC0gcmVhZF9jc3YoInN1cnZpdmluZyB3b3JrZXJzIHBlciBjb2xvbnkuY3N2IikNCg0KbW9ydGFsaXR5JGNvbG9ueSA8LSBhcy5mYWN0b3IobW9ydGFsaXR5JGNvbG9ueSkNCg0KYnJvb2QgPC0gbWVyZ2UoYnJvb2QsIG1vcnRhbGl0eSwgYnkueCA9ICJjb2xvbnkiKQ0KDQoNCmBgYA0KDQoNCiMgQnJvb2QgY2VsbHMgDQoNCmBgYHtyfQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkYnJvb2RfY2VsbHMpDQoNCnJhbmdlKGJyb29kJGJyb29kX2NlbGxzKQ0KDQpgYGANCg0KDQpUb3RhbCBicm9vZCBjb3VudCBub3QgYWZmZWN0ZWQgYnkgdHJlYXRtZW50LCBBbm92YSA9IDAuMzQzMzI4DQoNCmBgYHtyfQ0KDQpiLmduIDwtIGdsbShicm9vZF9jZWxscyB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoYi5nbikgICMgSSBkb24ndCBzZWUgYW55IG92ZXJkaXNwZXJzaW9uIA0KDQpBbm92YShiLmduKQ0KDQpiLmduYiA8LSBnbG0ubmIoYnJvb2RfY2VsbHMgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCwgZGF0YSA9IGJyb29kKQ0KDQpBSUMoYi5nbiwgYi5nbmIpDQoNCnBsb3QoYi5nbikNCnBsb3QoYi5nbmIpDQoNCmJyb29kLmVtbSA8LSBlbW1lYW5zKGIuZ24sICJ0cmVhdG1lbnQiKQ0KcGFpcnMoYnJvb2QuZW1tKQ0Kc3VtbWFyeShicm9vZC5lbW0pDQoNCiMgdGhlIHBvaXNzb24gcXEgcGxvdCBsb29rcyBiZXR0ZXIsIHNvIHdpdGggbm8gb3ZlcmRpc3BlcnNpb24gSSB3aWxsIGtlZXAgcG9pc3NvbiBldmVuIHRob3VnaCBuYiBBSUMgaXMgYmV0dGVyIA0KDQpzdW1fYnJvb2QgPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZSggbWVhbi5iID0gbWVhbihicm9vZF9jZWxscyksIA0KICAgICAgICAgICAgIG4uYiA9IGxlbmd0aChicm9vZF9jZWxscyksIA0KICAgICAgICAgICAgIHNkLmIgPSBzZChicm9vZF9jZWxscykpDQoNCnN1bV9icm9vZA0KDQpgYGANCg0KIyMjIGhvbmV5IHBvdHMNCg0KVHJlYXRtZW50IGRvZXMgbm90IGltcGFjdCBjb3VudCBvZiBob25leSBwb3RzIA0KDQpgYGB7cn0NCg0KaHAuZ24gPC0gZ2xtKGhvbmV5X3BvdCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoaHAuZ24pICMgb3ZlcmRpc3BlcnNpb24gbm90IGhvcnJpYmxlIGJ1dCBub3QgZ3JlYXQgDQpocC5nbmIgPC0gZ2xtLm5iKGhvbmV5X3BvdCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kLCBkYXRhID0gYnJvb2QpDQpzdW1tYXJ5KGhwLmduYikNCkFJQyhocC5nbiwgaHAuZ25iKQ0KDQpwbG90KGhwLmduKQ0KcGxvdChocC5nbmIpDQoNCkFub3ZhKGhwLmduYikNCkFub3ZhKGhwLmduKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkaG9uZXlfcG90KQ0KDQpocC5lbW0gPC0gZW1tZWFucyhocC5nbmIsICJ0cmVhdG1lbnQiKQ0Kc3VtbWFyeShocC5lbW0pDQpwYWlycyhocC5lbW0pDQoNCg0Kc3VtX2hwIDwtIGJyb29kICU+JQ0KICBncm91cF9ieSh0cmVhdG1lbnQpICU+JQ0KICBzdW1tYXJpc2UoIG1lYW4uaCA9IG1lYW4oaG9uZXlfcG90KSwgDQogICAgICAgICAgICAgbi5oID0gbGVuZ3RoKGhvbmV5X3BvdCksIA0KICAgICAgICAgICAgIHNkLmggPSBzZChob25leV9wb3QpKQ0KDQpzdW1faHANCg0KYGBgDQoNCg0KIyMjIExhcnZhZQ0KDQojIyMjIExpdmUgTGFydmFlDQoNCmBgYHtyfQ0KDQpsbC5nbiA8LSBnbG0obGl2ZV9sYXJ2YWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kLCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGxsLmduKQ0KbGwuZ25iIDwtIGdsbS5uYihsaXZlX2xhcnZhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QpDQpzdW1tYXJ5KGxsLmduYikNCkFJQyhsbC5nbiwgbGwuZ25iKSAgICAgICMgcG9pc3NvbiBvdmVyZGlzcGVyc2VkLCBuYiBBSUMgYmV0dGVyIA0KDQpwbG90KGxsLmduYikNCg0KQW5vdmEobGwuZ24pDQpBbm92YShsbC5nbmIpDQoNCnBsb3QoYnJvb2QkdHJlYXRtZW50LCBicm9vZCRsaXZlX2xhcnZhZSkNCg0KbGwuc3VtIDwtIGJyb29kICU+JQ0KICBncm91cF9ieSh0cmVhdG1lbnQpICU+JQ0KICBzdW1tYXJpc2UobGwubT0gbWVhbihsaXZlX2xhcnZhZSksIA0KICAgICAgICAgICAgbGwubiA9IGxlbmd0aChsaXZlX2xhcnZhZSksIA0KICAgICAgICAgICAgc2QubGwgPSBzZChsaXZlX2xhcnZhZSkpDQoNCmxsLnN1bQ0KDQpgYGANCiMjIyMgRGVhZCBMYXJ2YWUNCg0KYGBge3J9DQoNCmRsLmduIDwtIGdsbShkZWFkX2xhcnZhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoZGwuZ24pDQpkbC5nbmIgPC0gZ2xtLm5iKGRlYWRfbGFydmFlIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCkNCnN1bW1hcnkoZGwuZ25iKQ0KQUlDKGRsLmduLCBkbC5nbmIpICAgICAjcG9pc3NvbiBvdmVyZGlzcGVyc2VkIGFuZCBuYiBBSUMgYmV0dGVyIA0KcGxvdChkbC5nbmIpDQoNCkFub3ZhKGRsLmduYikNCg0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGJyb29kJGRlYWRfbGFydmFlKQ0KDQpkbC5zdW0gPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZShkbC5tPSBtZWFuKGRlYWRfbGFydmFlKSwgDQogICAgICAgICAgICBkbC5uID0gbGVuZ3RoKGRlYWRfbGFydmFlKSwgDQogICAgICAgICAgICBzZC5kbCA9IHNkKGRlYWRfbGFydmFlKSkNCg0KZGwuc3VtDQoNCmBgYA0KDQojIyMjIEFsbCBsYXJ2YWUNCg0KYGBge3J9DQoNCmJyb29kJGFsbF9sYXJ2YWUgPC0gYnJvb2QkZGVhZF9sYXJ2YWUgKyBicm9vZCRsaXZlX2xhcnZhZQ0KDQphbC5nbiA8LSBnbG0oYWxsX2xhcnZhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoYWwuZ24pDQphbC5nbmIgPC0gZ2xtLm5iKGFsbF9sYXJ2YWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShhbC5nbmIpDQpBSUMoYWwuZ24sIGFsLmduYikNCg0KQW5vdmEoYWwuZ25iKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkYWxsX2xhcnZhZSkNCg0KDQphbC5zdW0gPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZShhbC5tPSBtZWFuKGFsbF9sYXJ2YWUpLCANCiAgICAgICAgICAgIGFsLm4gPSBsZW5ndGgoYWxsX2xhcnZhZSksIA0KICAgICAgICAgICAgc2QuYWwgPSBzZChhbGxfbGFydmFlKSkNCg0KYWwuc3VtDQpgYGANCg0KIyMjIyBMYXJ2YWUgT3V0bGllcnMNCg0KYGBge3J9DQoNCm91dCA8LSBib3hwbG90LnN0YXRzKGJyb29kJGRlYWRfbGFydmFlKSRvdXQgDQpvdXQNCm91dF9kbCA8LSB3aGljaChicm9vZCRkZWFkX2xhcnZhZSAlaW4lIGMob3V0KSkgDQpvdXRfZGwNCg0KYnJvb2Rbb3V0X2RsLCBdDQoNCmJyb29kX2RsX291dCA8LSBicm9vZFstYygyMiwgMzIsIDM0LCAzOSwgNDUpLCBdIA0KDQoNCmRsLmduIDwtIGdsbShkZWFkX2xhcnZhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2RfZGxfb3V0LCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGRsLmduKQ0KZGwuZ25iIDwtIGdsbS5uYihkZWFkX2xhcnZhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2RfZGxfb3V0KQ0Kc3VtbWFyeShkbC5nbmIpDQpBbm92YShkbC5nbmIpDQoNCnBsb3QoYnJvb2RfZGxfb3V0JHRyZWF0bWVudCwgYnJvb2RfZGxfb3V0JGRlYWRfbGFydmFlKQ0KDQpgYGANCg0KDQojIyMjIGNiaW5kIExhcnZhZQ0KDQpgYGB7cn0NCmxhcnZhZS5tb3J0YWxpdHkgPC0gZ2xtKGNiaW5kKGxpdmVfbGFydmFlLCBkZWFkX2xhcnZhZSkgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQpzdW1tYXJ5KGxhcnZhZS5tb3J0YWxpdHkpDQoNCkFub3ZhKGxhcnZhZS5tb3J0YWxpdHkpDQoNCnBsb3QobGFydmFlLm1vcnRhbGl0eSkNCg0KbG0uZW1tIDwtIGVtbWVhbnMobGFydmFlLm1vcnRhbGl0eSwgInRyZWF0bWVudCIpDQpwYWlycyhsbS5lbW0pDQpzdW1tYXJ5KGxtLmVtbSkNCg0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGNiaW5kKGJyb29kJGRlYWRfbGFydmFlLCBicm9vZCRsaXZlX2xhcnZhZSkpDQoNCmVtbWVhbnMobGFydmFlLm1vcnRhbGl0eSwgcGFpcndpc2UgfiB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KDQplbW0xIDwtIGVtbWVhbnMobGFydmFlLm1vcnRhbGl0eSwgfiB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KZW1tMQ0KcGFpcnMoZW1tMSkNCg0KZW1tZGYgPC0gYXMuZGF0YS5mcmFtZShlbW0xKQ0KDQpwIDwtIGdncGxvdChkYXRhID0gZW1tZGYsIGFlcyh4PXRyZWF0bWVudCwgeT1wcm9iLCBmaWxsPXRyZWF0bWVudCkpICsgDQogIGdlb21fY29sKHBvc2l0aW9uID0gImRvZGdlIiwgY29sb3IgPSAiYmxhY2siKSArDQogIGNvb3JkX2NhcnRlc2lhbih5bGltID0gYygwLDEpKQ0KDQpwICsgZ2VvbV9lcnJvcmJhcihwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSguOSksd2lkdGg9LjI1LCBhZXMoeW1heD1hc3ltcC5VQ0wsIHltaW49YXN5bXAuTENMKSxhbHBoYT0uNikrDQogIGxhYnMoeD0iVHJlYXRtZW50IiwgeT0iUHJvYmFiaWxpdHkgb2YgU3Vydml2YWwiLCBmaWxsID0gInRyZWF0bWVudCIpICsgdGhlbWVfYncoKSArDQogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGMoInBlYWNocHVmZjMiLCAiZGFya3NlYWdyZWVuIiwgImxpZ2h0c2VhZ3JlZW4iLCAiZGFya29saXZlZ3JlZW4zIiwgImRhcmtvbGl2ZWdyZWVuNCIpKSArDQogICBsYWJzKHkgPSAiUHJvYmFiaWxpdHkgb2YgU3Vydml2YWwiLCkgKw0KICBnZ3RpdGxlKCJQcm9iYWJpbGl0eSBvZiBsYXJ2YWUgYmVpbmcgYWxpdmUgdXBvbiBjb2xvbnkgZGlzc2VjdGlvbiIpICsNCiAgc2NhbGVfeF9kaXNjcmV0ZShuYW1lID0gIlRyZWF0bWVudCIsIA0KICAgICAgICAgICAgICAgICAgIGxhYmVscyA9IGMoIjAgUFBCIiwgIjE1MCBQUEIiLCAiMSw1MDAgUFBCIiwgIjE1LDAwMCBQUEIiLCAiMTUwLDAwMCBQUEIiKSkgKw0KICB0aGVtZV9jb3dwbG90KCkgKw0KICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDEyKQ0KYGBgDQoNCiMjIyBQdXBhZQ0KDQojIyMjIExpdmUgUHVwYWUNCg0KYGBge3J9DQpscC5nbiA8LSBnbG0obGl2ZV9wdXBhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkobHAuZ24pDQpscC5nbmIgPC0gZ2xtLm5iKGxpdmVfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShscC5nbmIpDQpBSUMobHAuZ24sIGxwLmduYikNCg0KQW5vdmEobHAuZ25iKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkbGl2ZV9wdXBhZSkNCg0KDQpscC5zdW0gPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZShscC5tPSBtZWFuKGxpdmVfcHVwYWUpLCANCiAgICAgICAgICAgIHNkLmxwID0gc2QobGl2ZV9wdXBhZSkpDQoNCmxwLnN1bQ0KYGBgDQoNCg0KIyMjIyBEZWFkIFB1cGFlDQoNCmBgYHtyfQ0KDQpkcC5nbiA8LSBnbG0oZGVhZF9wdXBhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoZHAuZ24pDQpkcC5nbmIgPC0gZ2xtLm5iKGRlYWRfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKSAgICAjIHN0aWNrIHdpdGggdGhpcyBtb2RlbCBvdmVyYWxsIA0Kc3VtbWFyeShkcC5nbmIpDQpBSUMoZHAuZ24sIGRwLmduYikNCg0KQW5vdmEoZHAuZ25iKQ0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGJyb29kJGRlYWRfcHVwYWUpDQoNCmVtLmRwIDwtIGVtbWVhbnMoZHAuZ25iLCAidHJlYXRtZW50IikNCnN1bW1hcnkoZW0uZHApDQpwYWlycyhlbS5kcCkNCg0KZHAuZ25iMSA8LSBnbG0ubmIoZGVhZF9wdXBhZSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBxcm8sIGRhdGEgPSBicm9vZCkNCmRwLmduYjIgPC0gZ2xtLm5iKGRlYWRfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgcXJvICsgYWxpdmUsIGRhdGEgPSBicm9vZCkNCmRwLmduYjMgPC0gZ2xtLm5iKGRlYWRfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgcXJvICsgcm91bmQsIGRhdGEgPSBicm9vZCkNCkFJQyhkcC5nbmIsIGRwLmduYjEsIGRwLmduYjIsIGRwLmduYjMpDQpzdW1tYXJ5KGRwLmduYjEpDQpBbm92YShkcC5nbmIxKQ0KDQpwbG90KGRwLmduYjEpDQpwbG90KGRwLmduYikNCnBsb3QoZHAuZ25iMykNCg0KDQplbS5kcDMgPC0gZW1tZWFucyhkcC5nbmIzLCAidHJlYXRtZW50IikNCnN1bW1hcnkoZW0uZHAzKQ0KcGFpcnMoZW0uZHAzKQ0KDQoNCg0KZHAuc3VtIDwtIGJyb29kICU+JQ0KICBncm91cF9ieSh0cmVhdG1lbnQpICU+JQ0KICBzdW1tYXJpc2UoZHAubT0gbWVhbihkZWFkX3B1cGFlKSwgDQogICAgICAgICAgICBzZC5kcCA9IHNkKGRlYWRfcHVwYWUpKQ0KDQpkcC5zdW0NCg0KYGBgDQoNCiMjIyMgQWxsIFB1cGFlDQoNCmBgYHtyfQ0KYnJvb2QkYWxsX3B1cGFlIDwtIGJyb29kJGRlYWRfcHVwYWUgKyBicm9vZCRsaXZlX3B1cGFlDQoNCmFwLmduIDwtIGdsbShhbGxfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kLCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGFwLmduKQ0KYXAuZ25iIDwtIGdsbS5uYihhbGxfcHVwYWUgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShhcC5nbmIpDQpBSUMoYXAuZ24sIGFwLmduYikNCg0KQW5vdmEoYXAuZ25iKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkYWxsX3B1cGFlKQ0KDQoNCmFwLnN1bSA8LSBicm9vZCAlPiUNCiAgZ3JvdXBfYnkodHJlYXRtZW50KSAlPiUNCiAgc3VtbWFyaXNlKGFsLm09IG1lYW4oYWxsX3B1cGFlKSwgDQogICAgICAgICAgICBhbC5uID0gbGVuZ3RoKGFsbF9wdXBhZSksIA0KICAgICAgICAgICAgc2QuYWwgPSBzZChhbGxfcHVwYWUpKQ0KDQphcC5zdW0NCg0KDQpgYGANCg0KIyMgY2JpbmQgcHVwYWUNCg0KYGBge3J9DQpwdXBhZS5tb3J0YWxpdHkgPC0gZ2xtKGNiaW5kKGxpdmVfcHVwYWUsIGRlYWRfcHVwYWUpIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9IGJpbm9taWFsKCJsb2dpdCIpKQ0Kc3VtbWFyeShwdXBhZS5tb3J0YWxpdHkpDQoNCkFub3ZhKHB1cGFlLm1vcnRhbGl0eSkNCg0KcGxvdChwdXBhZS5tb3J0YWxpdHkpDQoNCmxtLmVtbSA8LSBlbW1lYW5zKHB1cGFlLm1vcnRhbGl0eSwgInRyZWF0bWVudCIpDQpwYWlycyhsbS5lbW0pDQpzdW1tYXJ5KGxtLmVtbSkNCg0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGNiaW5kKGJyb29kJGRlYWRfcHVwYWUsIGJyb29kJGxpdmVfcHVwYWUpKQ0KDQplbW1lYW5zKHB1cGFlLm1vcnRhbGl0eSwgcGFpcndpc2UgfiB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KDQplbW0yIDwtIGVtbWVhbnMocHVwYWUubW9ydGFsaXR5LCB+IHRyZWF0bWVudCwgdHlwZSA9ICJyZXNwb25zZSIpDQplbW0yDQpwYWlycyhlbW0yKQ0KDQplbW1kZjIgPC0gYXMuZGF0YS5mcmFtZShlbW0yKQ0KDQpwIDwtIGdncGxvdChkYXRhID0gZW1tZGYyLCBhZXMoeD10cmVhdG1lbnQsIHk9cHJvYiwgZmlsbD10cmVhdG1lbnQpKSArIA0KICBnZW9tX2NvbChwb3NpdGlvbiA9ICJkb2RnZSIsIGNvbG9yID0gImJsYWNrIikgKw0KICBjb29yZF9jYXJ0ZXNpYW4oeWxpbSA9IGMoMCwxKSkNCg0KcCArIGdlb21fZXJyb3JiYXIocG9zaXRpb249cG9zaXRpb25fZG9kZ2UoLjkpLHdpZHRoPS4yNSwgYWVzKHltYXg9YXN5bXAuVUNMLCB5bWluPWFzeW1wLkxDTCksYWxwaGE9LjYpKw0KICBsYWJzKHg9IlRyZWF0bWVudCIsIHk9IlByb2JhYmlsaXR5IG9mIFN1cnZpdmFsIiwgZmlsbCA9ICJ0cmVhdG1lbnQiKSArIHRoZW1lX2J3KCkgKw0KICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjKCJwZWFjaHB1ZmYzIiwgImRhcmtzZWFncmVlbiIsICJsaWdodHNlYWdyZWVuIiwgImRhcmtvbGl2ZWdyZWVuMyIsICJkYXJrb2xpdmVncmVlbjQiKSkgKw0KICAgbGFicyh5ID0gIlByb2JhYmlsaXR5IG9mIFN1cnZpdmFsIiwpICsNCiAgZ2d0aXRsZSgiUHJvYmFiaWxpdHkgb2YgcHVwYWUgYmVpbmcgYWxpdmUgdXBvbiBjb2xvbnkgZGlzc2VjdGlvbiIpICsNCiAgc2NhbGVfeF9kaXNjcmV0ZShuYW1lID0gIlRyZWF0bWVudCIsIA0KICAgICAgICAgICAgICAgICAgIGxhYmVscyA9IGMoIjAgUFBCIiwgIjE1MCBQUEIiLCAiMSw1MDAgUFBCIiwgIjE1LDAwMCBQUEIiLCAiMTUwLDAwMCBQUEIiKSkgKw0KICB0aGVtZV9jb3dwbG90KCkgKw0KICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDEyKQ0KYGBgDQoNCg0KDQojIyMgRGVhZCBMYXJ2YWUgKyBEZWFkIFB1cGFlDQoNCmBgYHtyfQ0KYnJvb2QkYWxsX2RfbHAgPC0gYnJvb2QkZGVhZF9sYXJ2YWUgKyBicm9vZCRkZWFkX3B1cGFlDQoNCmRscC5nbiA8LSBnbG0oYWxsX2RfbHAgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kLCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGRscC5nbikNCmRscC5nbmIgPC0gZ2xtLm5iKGFsbF9kX2xwIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCkNCnN1bW1hcnkoZGxwLmduYikNCkFJQyhkbHAuZ24sIGRscC5nbmIpDQoNCkFub3ZhKGRscC5nbmIpDQoNCnBsb3QoYnJvb2QkdHJlYXRtZW50LCBicm9vZCRhbGxfZF9scCkNCg0KDQpkbHAuc3VtIDwtIGJyb29kICU+JQ0KICBncm91cF9ieSh0cmVhdG1lbnQpICU+JQ0KICBzdW1tYXJpc2UoYWwubT0gbWVhbihhbGxfZF9scCksIA0KICAgICAgICAgICAgYWwubiA9IGxlbmd0aChhbGxfZF9scCksIA0KICAgICAgICAgICAgc2QuYWwgPSBzZChhbGxfZF9scCkpDQoNCmRscC5zdW0NCg0KYGBgDQoNCiMjIyBMaXZlIExhcnZhZSArIExpdmUgUHVwYWUNCg0KYGBge3J9DQpicm9vZCRhbGxfYV9scCA8LSBicm9vZCRsaXZlX2xhcnZhZSArIGJyb29kJGxpdmVfcHVwYWUNCg0KYWxwLmduIDwtIGdsbShhbGxfYV9scCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHJvdW5kICsgcXJvLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoYWxwLmduKQ0KYWxwLmduYiA8LSBnbG0ubmIoYWxsX2FfbHAgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShhbHAuZ25iKQ0KQUlDKGFscC5nbiwgYWxwLmduYikNCg0KQW5vdmEoYWxwLmduYikNCg0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGJyb29kJGFsbF9hX2xwKQ0KDQoNCmFscC5zdW0gPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZShhbC5tPSBtZWFuKGFsbF9hX2xwKSwgDQogICAgICAgICAgICBhbC5uID0gbGVuZ3RoKGFsbF9hX2xwKSwgDQogICAgICAgICAgICBzZC5hbCA9IHNkKGFsbF9hX2xwKSkNCg0KYWxwLnN1bQ0KYGBgDQoNCg0KIyMjIEFsbCBMYXJ2YWUgKyBBbGwgUHVwYWUNCg0KYGBge3J9DQoNCmJyb29kJGFsbF9scCA8LSBicm9vZCRhbGxfbGFydmFlICsgYnJvb2QkYWxsX3B1cGFlDQoNCmxwLmduIDwtIGdsbShhbGxfbHAgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kLCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGxwLmduKQ0KbHAuZ25iIDwtIGdsbS5uYihhbGxfbHAgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShscC5nbmIpDQpBSUMobHAuZ24sIGxwLmduYikNCg0KQW5vdmEobHAuZ25iKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkYWxsX2xwKQ0KDQoNCmxwLnN1bSA8LSBicm9vZCAlPiUNCiAgZ3JvdXBfYnkodHJlYXRtZW50KSAlPiUNCiAgc3VtbWFyaXNlKGxwLm09IG1lYW4oYWxsX2xwKSwgDQogICAgICAgICAgICBscC5uID0gbGVuZ3RoKGFsbF9scCksIA0KICAgICAgICAgICAgbHAuYWwgPSBzZChhbGxfbHApKQ0KDQpscC5zdW0NCg0KZW0ubHAgPC0gZW1tZWFucyhscC5nbmIsICJ0cmVhdG1lbnQiKQ0KcGFpcnMoZW0ubHApDQpgYGANCg0KIyNjYmluZCBsYXJ2YWUgYW5kIHB1cGFlDQoNCmBgYHtyfQ0KcGxfYmluZCA8LSBnbG0oY2JpbmQoYWxsX2FfbHAsIGFsbF9kX2xwKSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kLCBmYW1pbHkgPSBiaW5vbWlhbCgibG9naXQiKSkNCnN1bW1hcnkocGxfYmluZCApDQoNCkFub3ZhKHBsX2JpbmQgKQ0KDQpwbG90KHBsX2JpbmQgKQ0KDQpsbS5lbW0gPC0gZW1tZWFucyhwbF9iaW5kICwgInRyZWF0bWVudCIpDQpwYWlycyhsbS5lbW0pDQpzdW1tYXJ5KGxtLmVtbSkNCg0KcGxvdChicm9vZCR0cmVhdG1lbnQsIGNiaW5kKGJyb29kJGRlYWRfcHVwYWUsIGJyb29kJGxpdmVfcHVwYWUpKQ0KDQplbW1lYW5zKHBsX2JpbmQgLCBwYWlyd2lzZSB+IHRyZWF0bWVudCwgdHlwZSA9ICJyZXNwb25zZSIpDQoNCmVtbTMgPC0gZW1tZWFucyhwbF9iaW5kICwgfiB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KZW1tMw0KcGFpcnMoZW1tMykNCg0KZW1tZGYzIDwtIGFzLmRhdGEuZnJhbWUoZW1tMykNCg0KcCA8LSBnZ3Bsb3QoZGF0YSA9IGVtbWRmMywgYWVzKHg9dHJlYXRtZW50LCB5PXByb2IsIGZpbGw9dHJlYXRtZW50KSkgKyANCiAgZ2VvbV9jb2wocG9zaXRpb24gPSAiZG9kZ2UiLCBjb2xvciA9ICJibGFjayIpICsNCiAgY29vcmRfY2FydGVzaWFuKHlsaW0gPSBjKDAsMSkpDQoNCnAgKyBnZW9tX2Vycm9yYmFyKHBvc2l0aW9uPXBvc2l0aW9uX2RvZGdlKC45KSx3aWR0aD0uMjUsIGFlcyh5bWF4PWFzeW1wLlVDTCwgeW1pbj1hc3ltcC5MQ0wpLGFscGhhPS42KSsNCiAgbGFicyh4PSJUcmVhdG1lbnQiLCB5PSJQcm9iYWJpbGl0eSBvZiBTdXJ2aXZhbCIsIGZpbGwgPSAidHJlYXRtZW50IikgKyB0aGVtZV9idygpICsNCiAgc2NhbGVfZmlsbF9tYW51YWwodmFsdWVzID0gYygicGVhY2hwdWZmMyIsICJkYXJrc2VhZ3JlZW4iLCAibGlnaHRzZWFncmVlbiIsICJkYXJrb2xpdmVncmVlbjMiLCAiZGFya29saXZlZ3JlZW40IikpICsNCiAgIGxhYnMoeSA9ICJQcm9iYWJpbGl0eSBvZiBTdXJ2aXZhbCIsKSArDQogIGdndGl0bGUoIkNvbWJpbmVkIExhcnZhZSBhbmQgUHVwYWUgcHJvYmFiaWxpdHkgb2YgYmVpbmcgYWxpdmUgdXBvbiBjb2xvbnkgZGlzc2VjdGlvbiIpICsNCiAgc2NhbGVfeF9kaXNjcmV0ZShuYW1lID0gIlRyZWF0bWVudCIsIA0KICAgICAgICAgICAgICAgICAgIGxhYmVscyA9IGMoIjAgUFBCIiwgIjE1MCBQUEIiLCAiMSw1MDAgUFBCIiwgIjE1LDAwMCBQUEIiLCAiMTUwLDAwMCBQUEIiKSkgKw0KICB0aGVtZV9jb3dwbG90KCkgKw0KICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDEyKQ0KYGBgDQoNCg0KDQoNCiMjIyBFZ2dzDQoNCmBgYHtyfQ0KDQplZ2dzLmduIDwtIGdsbShlZ2dzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gInBvaXNzb24iKQ0Kc3VtbWFyeShlZ2dzLmduKQ0KI2VnZ3MuZ25iIDwtIGdsbS5uYihlZ2dzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCkNCiNzdW1tYXJ5KGVnZ3MuZ25iKQ0KI0FJQyhlZ2dzLmduLCBlZ2dzLmduYikgICMjIyMgSSB0aGluayB3ZSBoYXZlIHRvbyBtYW55IHplcm9zIGluIGVnZ3MsIEkgd2lsbCBzd2l0Y2ggdG8gYSB6ZXJvLWluZmxhdGVkIG5lZ2F0aXZlIGJpbm9taWFsIG1vZGVsDQoNCg0KbGlicmFyeShnbG1tVE1CKQ0KZWdnLnplcm8gPC0gZ2xtbVRNQihlZ2dzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIHJvdW5kICsgcXJvLCBmYW1pbHkgPSAibmJpbm9tMiIsIGJyb29kKSAgIyBoYWQgdG8gdGFrZSBvdXQgYWxpdmUsIGJlY2F1c2UgaXQgbWFkZSBtb2RlbCBub3QgY29udmVyZ2UgLS0gQW5vdmEgZnJvbSBnbG0gc2hvd3MgaXQncyBub3QgaW1wYWN0ZnVsIGFueXdhcyANCnN1bW1hcnkoZWdnLnplcm8pDQpBbm92YShlZ2cuemVybykNCg0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkZWdncykNCg0KDQplZ2dzLnN1bSA8LSBicm9vZCAlPiUNCiAgZ3JvdXBfYnkodHJlYXRtZW50KSAlPiUNCiAgc3VtbWFyaXNlKGVnZ3MubT0gbWVhbihlZ2dzKSwgDQogICAgICAgICAgICBlZ2dzLm4gPSBsZW5ndGgoZWdncyksIA0KICAgICAgICAgICAgZWdncy5hbCA9IHNkKGVnZ3MpKQ0KDQplZ2dzLnN1bQ0KDQplbS5lZ2cgPC0gZW1tZWFucyhlZ2cuemVybywgInRyZWF0bWVudCIpDQpwYWlycyhlbS5lZ2cpDQpzdW1tYXJ5KGVtLmVnZykNCg0KYGBgDQoNCiMjIyBDb3VudCBvZiBkcm9uZXMNCg0KYGBge3J9DQoNCmRyb25lcy5nbiA8LSBnbG0oZHJvbmVzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcm91bmQgKyBxcm8sIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gInBvaXNzb24iKQ0Kc3VtbWFyeShkcm9uZXMuZ24pDQpkcm9uZXMuZ25iIDwtIGdsbS5uYihkcm9uZXMgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyByb3VuZCArIHFybywgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShkcm9uZXMuZ25iKQ0KQUlDKGRyb25lcy5nbiwgZHJvbmVzLmduYikNCg0KQW5vdmEoZHJvbmVzLmduYikNCg0KZHJvbmUuZW0gPC0gZW1tZWFucyhkcm9uZXMuZ25iLCAidHJlYXRtZW50IikNCnBhaXJzKGRyb25lLmVtKQ0KDQpkcm9uZXMuc3VtIDwtIGJyb29kICU+JQ0KICBncm91cF9ieSh0cmVhdG1lbnQpICU+JQ0Kc3VtbWFyaXNlKG0uZCA9IG1lYW4oZHJvbmVzKSwgDQogICAgICAgICAgc2QuZCA9IHNkKGRyb25lcykpDQoNCmRyb25lcy5zdW0NCg0KYGBgDQoNCg==