Effect of Pristine on Bumble Bee Microcolonies - Brood Production
Emily Runnion
Data Collected 2022
Bees are frequently exposed to fungicides in agricultural landscapes, and while these chemicals are generally not considered to be harmful to insect pollinators, the sublethal effects of fungicides are not well understood. We investigated the non-target effects of exposure to field-realistic concentrations of the fungicide, Pristine®, for the common eastern bumble bee (Bombus impatiens). We fed queenless microcolonies of 5 bumble bees pollen patties contaminated with Pristine® ranging from 0 ppb-150,000 ppb, and evaluated the effects of chronic fungicide exposure on brood production, colony weight change, pollen consumption, and drone body quality. We found that while microcolonies fed pollen containing 1,500 ppb consumed the most pollen and produced the most brood, microcolonies exposed to 15,000 ppb gained weight at significantly lower rates. Fungicide exposure also negatively impacted brood development time and drone body quality. Our findings indicate that Pristine® has sublethal effects on bumble bee microcolony growth and brood development, highlighting the potential negative impacts of chronic fungicide exposure to bee health. In the face of the global decline of pollinator species, it is important for us to broaden our understanding of how bees react towards exposure of often overlooked stressors, such as fungicides, to account for potential consequences for pollinator health and pollination services.
Input relevent data files
### Figure out average pollen consumption by treatment
pollen <- read_csv("pollen1.csv", col_types = cols(round = col_factor(levels = c("1",
"2")), treatment = col_factor(levels = c("1",
"2", "3", "4", "5")), replicate = col_factor(levels = c("1",
"2", "3", "4", "5", "6", "7", "9", "11",
"12")), start_date = col_date(format = "%m/%d/%Y"),
start_time = col_character(), end_date = col_date(format = "%m/%d/%Y"),
end_time = col_character()))
pollen$colony <- as.factor(pollen$colony)
pollen$pid <- as.factor(pollen$pid)
pollen$count <- as.factor(pollen$count)
pollen <- subset(pollen, pollen$round != 1)
pollen <- na.omit(pollen)
range(pollen$difference)## [1] -0.30646 1.56542# get rid of negative numbers
pollen$difference[pollen$difference < 0] <- NA
pollen <- na.omit(pollen)
range(pollen$difference)## [1] 0.002715 1.565420# add queenright original colony column
qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)
pollen <- merge(pollen, qro, by.x = "colony")
pollen$pid <- as.factor(pollen$pid)
pollen$dose_consumed <- as.numeric(pollen$dose_consumed)Average pollen consumption per colony data input
pollen$whole_dif <- as.numeric(pollen$difference)
wd.1 <- lm(difference ~ treatment, data = pollen)
summary(wd.1)##
## Call:
## lm(formula = difference ~ treatment, data = pollen)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4853 -0.2416 -0.1504 0.1576 1.0638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.429541 0.024170 17.772 <2e-16 ***
## treatment2 0.072099 0.034886 2.067 0.0390 *
## treatment3 0.078078 0.034406 2.269 0.0235 *
## treatment4 0.058490 0.034988 1.672 0.0949 .
## treatment5 0.004982 0.035040 0.142 0.8870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3375 on 915 degrees of freedom
## Multiple R-squared: 0.009927, Adjusted R-squared: 0.005599
## F-statistic: 2.294 on 4 and 915 DF, p-value: 0.05773wd.emm <- emmeans(wd.1, "treatment")
summary(wd.emm)## treatment emmean SE df lower.CL upper.CL
## 1 0.430 0.0242 915 0.382 0.477
## 2 0.502 0.0252 915 0.452 0.551
## 3 0.508 0.0245 915 0.460 0.556
## 4 0.488 0.0253 915 0.438 0.538
## 5 0.435 0.0254 915 0.385 0.484
##
## Confidence level used: 0.95wd.summary <- pollen %>%
group_by(colony) %>%
summarize(whole.mean = mean(difference),
mean.dose = mean(dose_consumed))
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 mean.dose
## 1 1.11R2 0.2829509 0.00000
## 2 1.12R2 0.1697964 0.00000
## 3 1.1R2 0.5213458 0.00000
## 4 1.2R2 0.3374200 0.00000
## 5 1.3R2 0.4512891 0.00000
## 6 1.4R2 0.6063016 0.00000
## 7 1.5R2 0.7079516 0.00000
## 8 1.7R2 0.7400381 0.00000
## 9 1.9R2 0.2240081 0.00000
## 10 2.11R2 0.4178270 35.99958
## 11 2.12R2 0.4035568 35.62453
## 12 2.1R2 0.6101589 55.13368
## 13 2.2R2 0.5109234 47.65174
## 14 2.3R2 0.3280036 24.29676
## 15 2.4R2 0.3881394 32.86886
## 16 2.5R2 0.7194915 73.77696
## 17 2.7R2 0.5299685 46.45851
## 18 2.9R2 0.5857152 52.49618
## 19 3.11R2 0.4216410 357.65779
## 20 3.12R2 0.3390993 281.21803
## 21 3.1R2 0.3711948 307.88606
## 22 3.2R2 0.3461010 310.28287
## 23 3.3R2 0.8465806 919.09823
## 24 3.4R2 0.4120433 406.20875
## 25 3.5R2 0.8906211 969.89690
## 26 3.7R2 0.6266680 720.47811
## 27 3.9R2 0.4409511 417.98849
## 28 4.11R2 0.3456011 3697.92253
## 29 4.12R2 0.2496295 2085.61694
## 30 4.1R2 0.7074755 8936.65823
## 31 4.2R2 0.3871742 4380.60484
## 32 4.3R2 0.5800074 6846.91060
## 33 4.4R2 0.8356247 10230.80238
## 34 4.5R2 0.2335967 1757.69297
## 35 4.7R2 0.6237400 7030.06345
## 36 4.9R2 0.5322950 5682.26011
## 37 5.11R2 0.4113960 43070.09972
## 38 5.12R2 0.2077741 15383.81515
## 39 5.1R2 0.3782286 37367.48316
## 40 5.2R2 0.4912468 52985.85612
## 41 5.3R2 0.2128778 16364.60071
## 42 5.4R2 0.4563045 48624.01684
## 43 5.5R2 0.7095479 82144.93157
## 44 5.7R2 0.6113189 67415.74103
## 45 5.9R2 0.4188290 42601.27230Worker Mortality data input
mortality <- read_csv("surviving workers per colony.csv")
mortality$colony <- as.factor(mortality$colony)Begin data analysis
brood <- read_csv("brood.csv", col_types = cols(round = col_factor(levels = c("1",
"2")), treatment = col_factor(levels = c("1",
"2", "3", "4", "5")), replicate = col_factor(levels = c("1",
"2", "3", "4", "5", "7", "9", "11", "12"))))
brood$colony <- as.factor(brood$colony)
brood <- subset(brood, brood$round != 1)
brood <- merge(wd.summary, brood, by.x = "colony")
brood <- merge(brood, mortality, by.x = "colony")
brood$qro <- as.factor(brood$qro)Total Brood Cells
ggplot(brood, aes(x=brood_cells, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 5 ,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Total Brood Cells") +
labs(y = "Count", x = "Brood Cells")
descdist(brood$brood_cells, discrete = TRUE)
## summary statistics
## ------
## min: 0 max: 96
## median: 34
## mean: 36.48889
## estimated sd: 19.61635
## estimated skewness: 0.4372852
## estimated kurtosis: 3.842208Step 1 - Check for colinearity
brood.model <- lm(brood_cells~ treatment + whole.mean + alive + duration + qro, data = brood)
vif(brood.model)## GVIF Df GVIF^(1/(2*Df))
## treatment 1.289595 4 1.032302
## whole.mean 1.746024 1 1.321372
## alive 1.833697 1 1.354141
## duration 1.300597 1 1.140437
## qro 1.705126 3 1.093015Step 2 - Work down from most complicated model
brood.pois <- glm(brood_cells ~ treatment + whole.mean + alive + qro + duration, data = brood, family = "poisson")
summary(brood.pois) #overdispersed ##
## Call:
## glm(formula = brood_cells ~ treatment + whole.mean + alive +
## qro + duration, family = "poisson", data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.4129 -0.7154 0.1717 0.9328 2.5346
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.166903 0.180000 12.038 < 2e-16 ***
## treatment2 0.031131 0.084180 0.370 0.71152
## treatment3 0.097683 0.079475 1.229 0.21903
## treatment4 -0.120553 0.083167 -1.450 0.14719
## treatment5 -0.105504 0.091388 -1.154 0.24831
## whole.mean 2.594752 0.181938 14.262 < 2e-16 ***
## alive 0.089117 0.027771 3.209 0.00133 **
## qroB3 -0.219413 0.095896 -2.288 0.02214 *
## qroB4 -0.112401 0.091838 -1.224 0.22099
## qroB5 0.188688 0.070510 2.676 0.00745 **
## duration -0.007129 0.003543 -2.012 0.04417 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 526.03 on 44 degrees of freedom
## Residual deviance: 112.91 on 34 degrees of freedom
## AIC: 366.75
##
## Number of Fisher Scoring iterations: 5brood.int <- glm.nb(brood_cells ~ treatment*whole.mean + alive + qro + duration, data = brood, control = glm.control(maxit = 50)) #Keep this model
summary(brood.int)##
## Call:
## glm.nb(formula = brood_cells ~ treatment * whole.mean + alive +
## qro + duration, data = brood, control = glm.control(maxit = 50),
## init.theta = 42.40994585, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.56845 -0.81143 0.08355 0.56407 1.69786
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.492306 0.359861 4.147 3.37e-05 ***
## treatment2 0.726403 0.437787 1.659 0.097063 .
## treatment3 0.093813 0.355250 0.264 0.791722
## treatment4 1.057873 0.359654 2.941 0.003268 **
## treatment5 -0.172584 0.403621 -0.428 0.668951
## whole.mean 3.553586 0.472099 7.527 5.18e-14 ***
## alive 0.157933 0.039878 3.960 7.48e-05 ***
## qroB3 -0.272755 0.136724 -1.995 0.046050 *
## qroB4 -0.252943 0.142877 -1.770 0.076670 .
## qroB5 0.376043 0.104587 3.596 0.000324 ***
## duration -0.010259 0.005594 -1.834 0.066666 .
## treatment2:whole.mean -1.312705 0.812055 -1.617 0.105981
## treatment3:whole.mean -0.130133 0.618945 -0.210 0.833473
## treatment4:whole.mean -2.134235 0.640504 -3.332 0.000862 ***
## treatment5:whole.mean 0.082161 0.763115 0.108 0.914262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(42.4099) family taken to be 1)
##
## Null deviance: 311.811 on 44 degrees of freedom
## Residual deviance: 51.348 on 30 degrees of freedom
## AIC: 341.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 42.4
## Std. Err.: 18.9
##
## 2 x log-likelihood: -309.802brood.g1 <- glm.nb(brood_cells ~ treatment + whole.mean + alive + qro + duration, data = brood, control = glm.control(maxit = 50))
summary(brood.int)##
## Call:
## glm.nb(formula = brood_cells ~ treatment * whole.mean + alive +
## qro + duration, data = brood, control = glm.control(maxit = 50),
## init.theta = 42.40994585, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.56845 -0.81143 0.08355 0.56407 1.69786
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.492306 0.359861 4.147 3.37e-05 ***
## treatment2 0.726403 0.437787 1.659 0.097063 .
## treatment3 0.093813 0.355250 0.264 0.791722
## treatment4 1.057873 0.359654 2.941 0.003268 **
## treatment5 -0.172584 0.403621 -0.428 0.668951
## whole.mean 3.553586 0.472099 7.527 5.18e-14 ***
## alive 0.157933 0.039878 3.960 7.48e-05 ***
## qroB3 -0.272755 0.136724 -1.995 0.046050 *
## qroB4 -0.252943 0.142877 -1.770 0.076670 .
## qroB5 0.376043 0.104587 3.596 0.000324 ***
## duration -0.010259 0.005594 -1.834 0.066666 .
## treatment2:whole.mean -1.312705 0.812055 -1.617 0.105981
## treatment3:whole.mean -0.130133 0.618945 -0.210 0.833473
## treatment4:whole.mean -2.134235 0.640504 -3.332 0.000862 ***
## treatment5:whole.mean 0.082161 0.763115 0.108 0.914262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(42.4099) family taken to be 1)
##
## Null deviance: 311.811 on 44 degrees of freedom
## Residual deviance: 51.348 on 30 degrees of freedom
## AIC: 341.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 42.4
## Std. Err.: 18.9
##
## 2 x log-likelihood: -309.802brood.g2 <- glm.nb(brood_cells~ treatment + alive + qro + duration, data = brood)
summary(brood.g2)##
## Call:
## glm.nb(formula = brood_cells ~ treatment + alive + qro + duration,
## data = brood, init.theta = 4.697292782, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4708 -0.7907 -0.3366 0.5161 2.3254
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.836399 0.561580 3.270 0.001075 **
## treatment2 -0.037368 0.246291 -0.152 0.879405
## treatment3 0.037539 0.243169 0.154 0.877317
## treatment4 -0.023171 0.237151 -0.098 0.922165
## treatment5 -0.530436 0.252218 -2.103 0.035458 *
## alive 0.345855 0.059390 5.823 5.77e-09 ***
## qroB3 -0.174541 0.262515 -0.665 0.506128
## qroB4 0.684829 0.250440 2.734 0.006248 **
## qroB5 0.733897 0.196371 3.737 0.000186 ***
## duration 0.003106 0.011657 0.266 0.789919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.6973) family taken to be 1)
##
## Null deviance: 85.306 on 44 degrees of freedom
## Residual deviance: 51.007 on 35 degrees of freedom
## AIC: 396.94
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.70
## Std. Err.: 1.19
##
## 2 x log-likelihood: -374.94brood.g3 <- glm.nb(brood_cells~ whole.mean + alive + qro + duration, data = brood)
summary(brood.g3)##
## Call:
## glm.nb(formula = brood_cells ~ whole.mean + alive + qro + duration,
## data = brood, init.theta = 22.88343634, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3713 -0.2900 0.0137 0.4515 1.6161
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.847750 0.304716 6.064 1.33e-09 ***
## whole.mean 2.959131 0.285828 10.353 < 2e-16 ***
## alive 0.130280 0.039136 3.329 0.000872 ***
## qroB3 -0.223824 0.149938 -1.493 0.135495
## qroB4 -0.156689 0.155823 -1.006 0.314629
## qroB5 0.327225 0.112740 2.902 0.003702 **
## duration -0.009228 0.006037 -1.529 0.126355
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(22.8834) family taken to be 1)
##
## Null deviance: 237.182 on 44 degrees of freedom
## Residual deviance: 54.705 on 38 degrees of freedom
## AIC: 343.01
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 22.88
## Std. Err.: 8.02
##
## 2 x log-likelihood: -327.014anova(brood.int, brood.g1)## Likelihood ratio tests of Negative Binomial Models
##
## Response: brood_cells
## Model theta Resid. df
## 1 treatment + whole.mean + alive + qro + duration 25.88373 34
## 2 treatment * whole.mean + alive + qro + duration 42.40995 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -323.5991
## 2 -309.8017 1 vs 2 4 13.79742 0.007970486AIC(brood.int, brood.g1)## df AIC
## brood.int 16 341.8017
## brood.g1 12 347.5991anova(brood.g2, brood.g3)## Likelihood ratio tests of Negative Binomial Models
##
## Response: brood_cells
## Model theta Resid. df 2 x log-lik.
## 1 whole.mean + alive + qro + duration 22.883436 38 -327.0140
## 2 treatment + alive + qro + duration 4.697293 35 -374.9404
## Test df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 3 -47.9264 1AIC(brood.g2, brood.g3)## df AIC
## brood.g2 11 396.9404
## brood.g3 8 343.0140anova(brood.g1, brood.g3)## Likelihood ratio tests of Negative Binomial Models
##
## Response: brood_cells
## Model theta Resid. df
## 1 whole.mean + alive + qro + duration 22.88344 38
## 2 treatment + whole.mean + alive + qro + duration 25.88373 34
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -327.0140
## 2 -323.5991 1 vs 2 4 3.414869 0.4909402anova(brood.int, brood.g3)## Likelihood ratio tests of Negative Binomial Models
##
## Response: brood_cells
## Model theta Resid. df
## 1 whole.mean + alive + qro + duration 22.88344 38
## 2 treatment * whole.mean + alive + qro + duration 42.40995 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -327.0140
## 2 -309.8017 1 vs 2 8 17.21229 0.02797292Step 2a - Drop 1
drop1(brood.int, test = "Chisq") #try getting rid of duration## Single term deletions
##
## Model:
## brood_cells ~ treatment * whole.mean + alive + qro + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.348 339.80
## alive 1 68.194 354.65 16.8451 4.056e-05 ***
## qro 3 78.139 360.59 26.7901 6.515e-06 ***
## duration 1 54.742 341.20 3.3937 0.065448 .
## treatment:whole.mean 4 66.475 346.93 15.1264 0.004446 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1brood.int.update <- update(brood.int, .~. -duration)
anova(brood.int, brood.int.update)## Likelihood ratio tests of Negative Binomial Models
##
## Response: brood_cells
## Model theta
## 1 treatment + whole.mean + alive + qro + treatment:whole.mean 36.82998
## 2 treatment * whole.mean + alive + qro + duration 42.40995
## Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 31 -313.1059
## 2 30 -309.8017 1 vs 2 1 3.304188 0.06910348AIC(brood.int, brood.int.update) #keep everything ## df AIC
## brood.int 16 341.8017
## brood.int.update 15 343.1059Step 3 - Check model fit
qqnorm(resid(brood.int));qqline(resid(brood.int))
#not terrible summary(brood.int)##
## Call:
## glm.nb(formula = brood_cells ~ treatment * whole.mean + alive +
## qro + duration, data = brood, control = glm.control(maxit = 50),
## init.theta = 42.40994585, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.56845 -0.81143 0.08355 0.56407 1.69786
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.492306 0.359861 4.147 3.37e-05 ***
## treatment2 0.726403 0.437787 1.659 0.097063 .
## treatment3 0.093813 0.355250 0.264 0.791722
## treatment4 1.057873 0.359654 2.941 0.003268 **
## treatment5 -0.172584 0.403621 -0.428 0.668951
## whole.mean 3.553586 0.472099 7.527 5.18e-14 ***
## alive 0.157933 0.039878 3.960 7.48e-05 ***
## qroB3 -0.272755 0.136724 -1.995 0.046050 *
## qroB4 -0.252943 0.142877 -1.770 0.076670 .
## qroB5 0.376043 0.104587 3.596 0.000324 ***
## duration -0.010259 0.005594 -1.834 0.066666 .
## treatment2:whole.mean -1.312705 0.812055 -1.617 0.105981
## treatment3:whole.mean -0.130133 0.618945 -0.210 0.833473
## treatment4:whole.mean -2.134235 0.640504 -3.332 0.000862 ***
## treatment5:whole.mean 0.082161 0.763115 0.108 0.914262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(42.4099) family taken to be 1)
##
## Null deviance: 311.811 on 44 degrees of freedom
## Residual deviance: 51.348 on 30 degrees of freedom
## AIC: 341.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 42.4
## Std. Err.: 18.9
##
## 2 x log-likelihood: -309.802Step 4 - Results
Anova(brood.int)## Analysis of Deviance Table (Type II tests)
##
## Response: brood_cells
## LR Chisq Df Pr(>Chisq)
## treatment 4.462 4 0.347063
## whole.mean 121.197 1 < 2.2e-16 ***
## alive 16.845 1 4.056e-05 ***
## qro 26.790 3 6.515e-06 ***
## duration 3.394 1 0.065448 .
## treatment:whole.mean 15.126 4 0.004446 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BA <- as.data.frame(Anova(brood.int))
BA <- setDT(BA)
BA## LR Chisq Df Pr(>Chisq)
## 1: 4.462117 4 3.470631e-01
## 2: 121.196532 1 3.460908e-28
## 3: 16.845059 1 4.055878e-05
## 4: 26.790101 3 6.515060e-06
## 5: 3.393650 1 6.544790e-02
## 6: 15.126447 4 4.446017e-03brood_emm <- emmeans(brood.int, pairwise ~ treatment | whole.mean)
pairs(brood_emm)## whole.mean = 0.48:
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.09565 0.122 Inf -0.783 0.9356
## treatment1 - treatment3 -0.03128 0.121 Inf -0.259 0.9990
## treatment1 - treatment4 -0.03238 0.119 Inf -0.272 0.9988
## treatment1 - treatment5 0.13311 0.128 Inf 1.040 0.8370
## treatment2 - treatment3 0.06437 0.117 Inf 0.549 0.9821
## treatment2 - treatment4 0.06327 0.117 Inf 0.539 0.9833
## treatment2 - treatment5 0.22876 0.120 Inf 1.913 0.3103
## treatment3 - treatment4 -0.00109 0.119 Inf -0.009 1.0000
## treatment3 - treatment5 0.16439 0.117 Inf 1.402 0.6261
## treatment4 - treatment5 0.16548 0.125 Inf 1.328 0.6738
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatesplot(brood_emm, comparisons = TRUE)
plot(brood$treatment, brood$brood_cells)
brood.sum <- brood %>%
group_by(treatment) %>%
summarise(brood.m = mean(brood_cells),
brood.sd = sd(brood_cells),
brood.n = length(brood_cells)) %>%
mutate(brood.se = brood.sd/sqrt(brood.n))
brood.sum## # A tibble: 5 × 5
## treatment brood.m brood.sd brood.n brood.se
## <fct> <dbl> <dbl> <int> <dbl>
## 1 1 33.8 22.6 9 7.53
## 2 2 36.9 11.2 9 3.74
## 3 3 45.6 26.2 9 8.73
## 4 4 36.7 18.3 9 6.09
## 5 5 29.6 17.5 9 5.82Step 5 - Visualize
broodtuk.means <- emmeans(object = brood.int,
specs = "treatment",
adjust = "Tukey",
type = "response")
broodtuk.means <- emmeans(object = brood.int,
specs = "treatment",
adjust = "sidak",
type = "response")
brood.cld.model <- cld(object = broodtuk.means,
adjust = "Tukey",
Letters = letters,
alpha = 0.05)
brood.cld.model## treatment response SE df asymp.LCL asymp.UCL .group
## 5 25.5 2.39 Inf 20.0 32.4 a
## 1 29.1 2.67 Inf 23.0 36.8 a
## 3 30.0 2.70 Inf 23.8 37.8 a
## 4 30.1 2.72 Inf 23.8 37.9 a
## 2 32.0 2.84 Inf 25.5 40.2 a
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 5 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.broodtuk.means <- as.data.frame(broodtuk.means)
brood_max <- brood %>%
group_by(treatment) %>%
summarize(maxbrood = max((brood_cells)))
brood_for_plotting <- full_join(broodtuk.means, brood_max,
by="treatment")
brood_for_plotting$maxb <- brood_for_plotting$response + brood_for_plotting$SEggplot(data = brood_for_plotting, aes(x = treatment, y = response, fill = treatment)) +
geom_col(col = "black")+
coord_cartesian(ylim=c(20, 40)) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
geom_errorbar(aes(ymin = response-SE,
ymax = response+SE),
position = position_dodge(0.9), width = 0.8) +
labs(y = "Days",) +
ggtitle("Average Colony Duration per Treatment") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot()+
theme_classic(base_size = 30) +
theme(legend.position = "none") +
annotate(geom = "text",
x = 1, y = 40,
label = "P > 0.05",
size = 12) +
annotate(geom = "text",
x = c(1, 2, 3, 4, 5),
y = c(brood_for_plotting$maxb + 1),
label = c("a", "a", "a", "a", "a"),
size = 12) +
theme(legend.position = "none")
Brood Cells Summary Results
b_A <- setDT(as.data.frame(Anova(brood.int)))
b_A## LR Chisq Df Pr(>Chisq)
## 1: 4.462117 4 3.470631e-01
## 2: 121.196532 1 3.460908e-28
## 3: 16.845059 1 4.055878e-05
## 4: 26.790101 3 6.515060e-06
## 5: 3.393650 1 6.544790e-02
## 6: 15.126447 4 4.446017e-03b_emm <- emmeans(brood.int, pairwise ~ treatment | whole.mean)
b_contrasts <- setDT(as.data.frame(b_emm$contrasts))
b_means <- as.data.frame(b_emm$emmeans)
b_means <- setDT(b_means)
b_means## treatment whole.mean emmean SE df asymp.LCL asymp.UCL
## 1: 1 0.480499 3.370747 0.09169486 Inf 3.191028 3.550466
## 2: 2 0.480499 3.466397 0.08859853 Inf 3.292747 3.640046
## 3: 3 0.480499 3.402031 0.08982339 Inf 3.225981 3.578082
## 4: 4 0.480499 3.403123 0.09060210 Inf 3.225546 3.580700
## 5: 5 0.480499 3.237641 0.09386184 Inf 3.053676 3.421607b_contrasts## contrast whole.mean estimate SE df z.ratio
## 1: treatment1 - treatment2 0.480499 -0.095649482 0.1221553 Inf -0.783015242
## 2: treatment1 - treatment3 0.480499 -0.031284274 0.1209296 Inf -0.258698325
## 3: treatment1 - treatment4 0.480499 -0.032375882 0.1190949 Inf -0.271849424
## 4: treatment1 - treatment5 0.480499 0.133105716 0.1280396 Inf 1.039567198
## 5: treatment2 - treatment3 0.480499 0.064365208 0.1172839 Inf 0.548798554
## 6: treatment2 - treatment4 0.480499 0.063273600 0.1174068 Inf 0.538926216
## 7: treatment2 - treatment5 0.480499 0.228755198 0.1195785 Inf 1.913013276
## 8: treatment3 - treatment4 0.480499 -0.001091608 0.1187291 Inf -0.009194109
## 9: treatment3 - treatment5 0.480499 0.164389990 0.1172214 Inf 1.402388712
## 10: treatment4 - treatment5 0.480499 0.165481598 0.1246195 Inf 1.327894655
## p.value
## 1: 0.9356128
## 2: 0.9990164
## 3: 0.9988045
## 4: 0.8369782
## 5: 0.9821157
## 6: 0.9832866
## 7: 0.3102514
## 8: 1.0000000
## 9: 0.6261134
## 10: 0.6737628Eggs
brood_egg <- subset(brood, select = c(eggs, treatment, whole.mean, alive, duration, qro, colony))
ggplot(brood_egg, aes(x=eggs, fill = treatment)) +
geom_histogram(position = "identity", binwidth =1 ,col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Eggs per Treatment") +
labs(y = "Count", x = "Egg Count")
Step 1 - Check for co-linearity
egg.model <- lm(eggs~ treatment + whole.mean + alive + duration + qro, data = brood_egg)
vif(egg.model)## GVIF Df GVIF^(1/(2*Df))
## treatment 1.289595 4 1.032302
## whole.mean 1.746024 1 1.321372
## alive 1.833697 1 1.354141
## duration 1.300597 1 1.140437
## qro 1.705126 3 1.093015Step 2 - Work down from most complicated model
egg.int <- glm(eggs ~ treatment*whole.mean + alive + duration + qro, data = brood_egg, family = "poisson")
summary(egg.int) #overdispersed##
## Call:
## glm(formula = eggs ~ treatment * whole.mean + alive + duration +
## qro, family = "poisson", data = brood_egg)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.701 -2.094 -1.117 1.874 5.284
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.239698 0.657513 0.365 0.7154
## treatment2 2.873774 0.567212 5.066 4.05e-07 ***
## treatment3 0.448156 0.653289 0.686 0.4927
## treatment4 -0.911984 0.658684 -1.385 0.1662
## treatment5 0.316782 0.663177 0.478 0.6329
## whole.mean 3.931099 0.732667 5.365 8.07e-08 ***
## alive -0.074626 0.051996 -1.435 0.1512
## duration -0.004335 0.011863 -0.365 0.7148
## qroB3 0.879233 0.195251 4.503 6.70e-06 ***
## qroB4 1.642314 0.190086 8.640 < 2e-16 ***
## qroB5 0.777832 0.190750 4.078 4.55e-05 ***
## treatment2:whole.mean -6.020061 1.017987 -5.914 3.35e-09 ***
## treatment3:whole.mean -2.504669 1.025819 -2.442 0.0146 *
## treatment4:whole.mean 0.415824 1.081746 0.384 0.7007
## treatment5:whole.mean -2.377079 1.135097 -2.094 0.0362 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 624.61 on 44 degrees of freedom
## Residual deviance: 243.98 on 30 degrees of freedom
## AIC: 394.37
##
## Number of Fisher Scoring iterations: 7egg.m1 <- glm(eggs ~ treatment + whole.mean + alive + duration + qro, data = brood_egg, family = "poisson")
summary(egg.m1) #overdispersed##
## Call:
## glm(formula = eggs ~ treatment + whole.mean + alive + duration +
## qro, family = "poisson", data = brood_egg)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.030 -2.502 -1.118 1.264 6.167
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.146417 0.425846 2.692 0.00710 **
## treatment2 -0.409144 0.153405 -2.667 0.00765 **
## treatment3 -1.154044 0.204435 -5.645 1.65e-08 ***
## treatment4 -0.759351 0.162709 -4.667 3.06e-06 ***
## treatment5 -1.030736 0.205385 -5.019 5.21e-07 ***
## whole.mean 2.721246 0.430958 6.314 2.71e-10 ***
## alive -0.056929 0.048742 -1.168 0.24282
## duration -0.009316 0.009847 -0.946 0.34409
## qroB3 1.026369 0.187233 5.482 4.21e-08 ***
## qroB4 1.391845 0.178009 7.819 5.33e-15 ***
## qroB5 0.961607 0.166765 5.766 8.11e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 624.61 on 44 degrees of freedom
## Residual deviance: 289.92 on 34 degrees of freedom
## AIC: 432.31
##
## Number of Fisher Scoring iterations: 6egg.int2 <- glm.nb(eggs ~ treatment*whole.mean + qro + duration, data = brood_egg) #if I include the other variables here the model won't converge
summary(egg.int2) #not overdispersed anymore##
## Call:
## glm.nb(formula = eggs ~ treatment * whole.mean + qro + duration,
## data = brood_egg, init.theta = 0.9604946898, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1808 -1.2604 -0.1925 0.4184 1.7803
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.607936 1.533573 0.396 0.6918
## treatment2 1.759648 1.899392 0.926 0.3542
## treatment3 -0.302662 1.523702 -0.199 0.8425
## treatment4 -2.144087 1.571473 -1.364 0.1724
## treatment5 -2.826548 1.752008 -1.613 0.1067
## whole.mean 2.727897 1.955357 1.395 0.1630
## qroB3 0.841177 0.608696 1.382 0.1670
## qroB4 1.454739 0.658342 2.210 0.0271 *
## qroB5 0.964937 0.428774 2.250 0.0244 *
## duration -0.005203 0.029154 -0.178 0.8584
## treatment2:whole.mean -4.170295 3.729881 -1.118 0.2635
## treatment3:whole.mean -1.017972 2.864754 -0.355 0.7223
## treatment4:whole.mean 2.839330 2.947353 0.963 0.3354
## treatment5:whole.mean 4.274859 3.476783 1.230 0.2189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9605) family taken to be 1)
##
## Null deviance: 88.408 on 44 degrees of freedom
## Residual deviance: 52.051 on 31 degrees of freedom
## AIC: 275.64
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.960
## Std. Err.: 0.290
##
## 2 x log-likelihood: -245.644egg.m2 <- glm.nb(eggs ~ treatment + whole.mean + qro + duration, data = brood_egg) #stick with this model
anova(egg.int2, egg.m2)## Likelihood ratio tests of Negative Binomial Models
##
## Response: eggs
## Model theta Resid. df 2 x log-lik.
## 1 treatment + whole.mean + qro + duration 0.8672220 35 -248.8184
## 2 treatment * whole.mean + qro + duration 0.9604947 31 -245.6442
## Test df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 4 3.174184 0.5291107AIC(egg.int2, egg.m2)## df AIC
## egg.int2 15 275.6442
## egg.m2 11 270.8184Step 2.a - Drop 1
drop1(egg.m2, test = "Chisq")## Single term deletions
##
## Model:
## eggs ~ treatment + whole.mean + qro + duration
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.984 268.82
## treatment 4 56.333 265.17 4.3484 0.360909
## whole.mean 1 60.181 275.01 8.1969 0.004196 **
## qro 3 59.491 270.32 7.5070 0.057379 .
## duration 1 52.905 267.74 0.9209 0.337249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1egg.m2.update <- update(egg.m2, .~. -duration)
anova(egg.m2, egg.m2.update)## Likelihood ratio tests of Negative Binomial Models
##
## Response: eggs
## Model theta Resid. df 2 x log-lik.
## 1 treatment + whole.mean + qro 0.855366 36 -249.7371
## 2 treatment + whole.mean + qro + duration 0.867222 35 -248.8184
## Test df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 1 0.9187158 0.3378124drop1(egg.m2.update, test = "Chisq")## Single term deletions
##
## Model:
## eggs ~ treatment + whole.mean + qro
## Df Deviance AIC LRT Pr(>Chi)
## <none> 52.460 267.74
## treatment 4 56.735 264.01 4.2749 0.370083
## whole.mean 1 59.659 272.94 7.1995 0.007292 **
## qro 3 59.293 268.57 6.8329 0.077417 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Step 3 - Check model fit
qqnorm(resid(egg.m2.update));qqline(resid(egg.m2.update))
plot(resid(egg.m2.update)) +
abline(h=0, col="red", lwd=2) 
## integer(0)Step 4 - Results
Anova(egg.m2.update)## Analysis of Deviance Table (Type II tests)
##
## Response: eggs
## LR Chisq Df Pr(>Chisq)
## treatment 4.2749 4 0.370083
## whole.mean 7.1995 1 0.007292 **
## qro 6.8329 3 0.077417 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(egg.m2.update)##
## Call:
## glm.nb(formula = eggs ~ treatment + whole.mean + qro, data = brood_egg,
## init.theta = 0.8553660464, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1528 -1.3208 -0.3127 0.2674 2.0186
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2479 0.6239 0.397 0.69114
## treatment2 -0.4207 0.5452 -0.772 0.44035
## treatment3 -0.9663 0.5628 -1.717 0.08600 .
## treatment4 -0.7371 0.5549 -1.328 0.18412
## treatment5 -1.0130 0.5662 -1.789 0.07359 .
## whole.mean 3.3460 1.0673 3.135 0.00172 **
## qroB3 0.6962 0.5768 1.207 0.22745
## qroB4 0.9635 0.6004 1.605 0.10858
## qroB5 1.1071 0.4358 2.540 0.01108 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.8554) family taken to be 1)
##
## Null deviance: 81.501 on 44 degrees of freedom
## Residual deviance: 52.460 on 36 degrees of freedom
## AIC: 269.74
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.855
## Std. Err.: 0.253
##
## 2 x log-likelihood: -249.737egg.em <- emmeans(egg.m2.update, "treatment")
pairs(egg.em)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 0.4207 0.545 Inf 0.772 0.9388
## treatment1 - treatment3 0.9663 0.563 Inf 1.717 0.4234
## treatment1 - treatment4 0.7371 0.555 Inf 1.328 0.6736
## treatment1 - treatment5 1.0130 0.566 Inf 1.789 0.3798
## treatment2 - treatment3 0.5456 0.560 Inf 0.975 0.8666
## treatment2 - treatment4 0.3163 0.553 Inf 0.572 0.9792
## treatment2 - treatment5 0.5923 0.570 Inf 1.040 0.8369
## treatment3 - treatment4 -0.2293 0.566 Inf -0.405 0.9944
## treatment3 - treatment5 0.0467 0.588 Inf 0.079 1.0000
## treatment4 - treatment5 0.2760 0.580 Inf 0.476 0.9895
##
## Results are averaged over the levels of: qro
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatesStep 5 - Visualize
plot(egg.em, comparisons = TRUE)
Results Summary
egg.m2.update##
## Call: glm.nb(formula = eggs ~ treatment + whole.mean + qro, data = brood_egg,
## init.theta = 0.8553660464, link = log)
##
## Coefficients:
## (Intercept) treatment2 treatment3 treatment4 treatment5 whole.mean
## 0.2479 -0.4207 -0.9663 -0.7371 -1.0130 3.3460
## qroB3 qroB4 qroB5
## 0.6962 0.9635 1.1071
##
## Degrees of Freedom: 44 Total (i.e. Null); 36 Residual
## Null Deviance: 81.5
## Residual Deviance: 52.46 AIC: 269.7eggs_A <- as.data.frame(Anova(egg.m2.update))
eggs_A <- setDT(eggs_A)
eggs_A## LR Chisq Df Pr(>Chisq)
## 1: 4.274866 4 0.370083013
## 2: 7.199511 1 0.007292346
## 3: 6.832948 3 0.077417264egg_emm <- emmeans(egg.m2.update, pairwise ~ treatment)
egg_contrasts <- as.data.frame(egg_emm$contrasts)
egg_contrasts <- setDT(egg_contrasts)
egg_contrasts## contrast estimate SE df z.ratio p.value
## 1: treatment1 - treatment2 0.42071755 0.5452495 Inf 0.77160561 0.9388022
## 2: treatment1 - treatment3 0.96632181 0.5628396 Inf 1.71686877 0.4234423
## 3: treatment1 - treatment4 0.73705732 0.5549398 Inf 1.32817536 0.6735858
## 4: treatment1 - treatment5 1.01301709 0.5661951 Inf 1.78916616 0.3798409
## 5: treatment2 - treatment3 0.54560425 0.5595746 Inf 0.97503404 0.8665549
## 6: treatment2 - treatment4 0.31633976 0.5531408 Inf 0.57189738 0.9791567
## 7: treatment2 - treatment5 0.59229954 0.5696577 Inf 1.03974637 0.8368919
## 8: treatment3 - treatment4 -0.22926449 0.5664141 Inf -0.40476479 0.9943637
## 9: treatment3 - treatment5 0.04669529 0.5875714 Inf 0.07947168 0.9999910
## 10: treatment4 - treatment5 0.27595978 0.5795266 Inf 0.47618140 0.9895118eggs_means <- as.data.frame(egg_emm$emmeans)
eggs_means <- setDT(eggs_means)
eggs_means## treatment emmean SE df asymp.LCL asymp.UCL
## 1: 1 2.547329 0.4011395 Inf 1.7611096 3.333548
## 2: 2 2.126611 0.4076368 Inf 1.3276576 2.925564
## 3: 3 1.581007 0.4281296 Inf 0.7418882 2.420125
## 4: 4 1.810271 0.4206514 Inf 0.9858098 2.634733
## 5: 5 1.534312 0.4307688 Inf 0.6900202 2.378603Honey Pots
brood_hp <- subset(brood, select = c(honey_pot, treatment, whole.mean, alive, duration, qro, colony))
ggplot(brood_hp, aes(x=honey_pot, fill = treatment)) +
geom_histogram(position = "identity", binwidth = 0.5, col=I("black")) +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
name = "Pristine Level",
labels = c("Treatment 1 (control)", "Treatment 2",
"Treatment 3", "Treatment 4", "Treatment 5")) +
ggtitle("Honey Pots per Treatment") +
labs(y = "Count", x = "HP Count")
Step 1 - Check for co-linearity
hp.model <- lm(honey_pot ~ treatment + whole.mean + alive + duration + qro, data = brood_hp)
vif(hp.model)## GVIF Df GVIF^(1/(2*Df))
## treatment 1.289595 4 1.032302
## whole.mean 1.746024 1 1.321372
## alive 1.833697 1 1.354141
## duration 1.300597 1 1.140437
## qro 1.705126 3 1.093015Step 2 - Work down from most complicated model
hp.int <- glm(honey_pot ~ treatment*whole.mean + alive + duration + qro, data = brood_hp, family = "poisson")
summary(hp.int) #OD##
## Call:
## glm(formula = honey_pot ~ treatment * whole.mean + alive + duration +
## qro, family = "poisson", data = brood_hp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8385 -1.0952 -0.2029 0.7360 2.4498
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.207199 0.689885 -0.300 0.7639
## treatment2 0.182784 0.761564 0.240 0.8103
## treatment3 0.030914 0.680473 0.045 0.9638
## treatment4 0.962727 0.636482 1.513 0.1304
## treatment5 -0.154971 0.694137 -0.223 0.8233
## whole.mean 1.430035 0.920802 1.553 0.1204
## alive 0.190088 0.074132 2.564 0.0103 *
## duration 0.005145 0.010980 0.469 0.6394
## qroB3 -0.100110 0.231679 -0.432 0.6657
## qroB4 -0.058380 0.253005 -0.231 0.8175
## qroB5 0.233847 0.180926 1.293 0.1962
## treatment2:whole.mean 0.502602 1.397351 0.360 0.7191
## treatment3:whole.mean -0.055874 1.206928 -0.046 0.9631
## treatment4:whole.mean -1.151006 1.147432 -1.003 0.3158
## treatment5:whole.mean 0.816948 1.310028 0.624 0.5329
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 112.51 on 44 degrees of freedom
## Residual deviance: 65.86 on 30 degrees of freedom
## AIC: 248.07
##
## Number of Fisher Scoring iterations: 5hp.m1 <- glm(honey_pot ~ treatment + whole.mean + alive + duration + qro, data = brood_hp, family = "poisson")
summary(hp.m1) #OD##
## Call:
## glm(formula = honey_pot ~ treatment + whole.mean + alive + duration +
## qro, family = "poisson", data = brood_hp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7140 -1.1077 0.1059 0.5110 2.4649
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.002638 0.487579 -0.005 0.99568
## treatment2 0.458112 0.212700 2.154 0.03126 *
## treatment3 0.024516 0.224193 0.109 0.91292
## treatment4 0.351763 0.210800 1.669 0.09518 .
## treatment5 0.241896 0.226352 1.069 0.28522
## whole.mean 1.151562 0.454945 2.531 0.01137 *
## alive 0.177660 0.068287 2.602 0.00928 **
## duration 0.004613 0.010005 0.461 0.64474
## qroB3 -0.126997 0.219477 -0.579 0.56284
## qroB4 0.098654 0.236156 0.418 0.67613
## qroB5 0.181534 0.167930 1.081 0.27969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 112.507 on 44 degrees of freedom
## Residual deviance: 68.751 on 34 degrees of freedom
## AIC: 242.96
##
## Number of Fisher Scoring iterations: 5hp.int2 <- glm.nb(honey_pot ~ treatment*whole.mean + alive + duration + qro, data = brood_hp)
summary(hp.int2) #not OD ##
## Call:
## glm.nb(formula = honey_pot ~ treatment * whole.mean + alive +
## duration + qro, data = brood_hp, init.theta = 23.2488872,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6310 -0.9602 -0.2100 0.6821 2.1208
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.165242 0.752625 -0.220 0.8262
## treatment2 0.176578 0.849680 0.208 0.8354
## treatment3 -0.007041 0.741812 -0.009 0.9924
## treatment4 0.957068 0.698971 1.369 0.1709
## treatment5 -0.221696 0.763112 -0.291 0.7714
## whole.mean 1.427621 1.007733 1.417 0.1566
## alive 0.193531 0.080537 2.403 0.0163 *
## duration 0.003925 0.012205 0.322 0.7478
## qroB3 -0.128467 0.262985 -0.488 0.6252
## qroB4 -0.069338 0.289485 -0.240 0.8107
## qroB5 0.245354 0.203477 1.206 0.2279
## treatment2:whole.mean 0.512522 1.579259 0.325 0.7455
## treatment3:whole.mean 0.011416 1.328816 0.009 0.9931
## treatment4:whole.mean -1.137667 1.276957 -0.891 0.3730
## treatment5:whole.mean 0.946490 1.461589 0.648 0.5173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(23.2489) family taken to be 1)
##
## Null deviance: 92.563 on 44 degrees of freedom
## Residual deviance: 54.556 on 30 degrees of freedom
## AIC: 249.09
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 23.2
## Std. Err.: 26.8
##
## 2 x log-likelihood: -217.092hp.m2 <- glm.nb(honey_pot ~ treatment + whole.mean + alive + duration + qro, data = brood_hp)
summary(hp.m2) #keep this model ##
## Call:
## glm.nb(formula = honey_pot ~ treatment + whole.mean + alive +
## duration + qro, data = brood_hp, init.theta = 18.15285757,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4708 -1.0404 0.1104 0.4711 2.0842
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.001796 0.559486 0.003 0.9974
## treatment2 0.457373 0.245041 1.867 0.0620 .
## treatment3 0.021826 0.253828 0.086 0.9315
## treatment4 0.370384 0.240900 1.537 0.1242
## treatment5 0.231539 0.258995 0.894 0.3713
## whole.mean 1.208730 0.526496 2.296 0.0217 *
## alive 0.182876 0.076220 2.399 0.0164 *
## duration 0.003274 0.011574 0.283 0.7773
## qroB3 -0.161323 0.256370 -0.629 0.5292
## qroB4 0.089812 0.278004 0.323 0.7466
## qroB5 0.204025 0.196376 1.039 0.2988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(18.1529) family taken to be 1)
##
## Null deviance: 88.334 on 44 degrees of freedom
## Residual deviance: 54.400 on 34 degrees of freedom
## AIC: 243.39
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 18.2
## Std. Err.: 17.1
##
## 2 x log-likelihood: -219.385anova(hp.int2, hp.m2)## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pot
## Model theta Resid. df
## 1 treatment + whole.mean + alive + duration + qro 18.15286 34
## 2 treatment * whole.mean + alive + duration + qro 23.24889 30
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -219.3853
## 2 -217.0919 1 vs 2 4 2.293383 0.6819741Step 2.a - Drop 1
drop1(hp.m2, test = "Chisq")## Single term deletions
##
## Model:
## honey_pot ~ treatment + whole.mean + alive + duration + qro
## Df Deviance AIC LRT Pr(>Chi)
## <none> 54.400 241.38
## treatment 4 60.772 239.76 6.3720 0.17304
## whole.mean 1 59.494 244.48 5.0936 0.02401 *
## alive 1 60.506 245.49 6.1060 0.01347 *
## duration 1 54.479 239.46 0.0791 0.77854
## qro 3 56.343 237.33 1.9423 0.58447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hp.m3 <- glm.nb(honey_pot ~ treatment + whole.mean + alive + qro, data = brood_hp)
hp.m4 <- glm.nb(honey_pot ~ treatment + whole.mean + alive + duration, data = brood_hp)
hp.m5 <- glm.nb(honey_pot ~ treatment + whole.mean + alive, data = brood_hp)
anova(hp.m2, hp.m3)## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pot
## Model theta Resid. df
## 1 treatment + whole.mean + alive + qro 17.58049 35
## 2 treatment + whole.mean + alive + duration + qro 18.15286 34
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -219.4632
## 2 -219.3853 1 vs 2 1 0.07794046 0.7801081anova(hp.m3, hp.m4)## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pot
## Model theta Resid. df 2 x log-lik.
## 1 treatment + whole.mean + alive + duration 17.52117 37 -221.3263
## 2 treatment + whole.mean + alive + qro 17.58049 35 -219.4632
## Test df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 2 1.863058 0.393951drop1(hp.m3, test = "Chisq")## Single term deletions
##
## Model:
## honey_pot ~ treatment + whole.mean + alive + qro
## Df Deviance AIC LRT Pr(>Chi)
## <none> 54.128 239.46
## treatment 4 60.941 238.28 6.8140 0.14605
## whole.mean 1 60.034 243.37 5.9063 0.01509 *
## alive 1 60.166 243.50 6.0383 0.01400 *
## qro 3 56.012 235.35 1.8847 0.59669
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(hp.m4, hp.m5)## Likelihood ratio tests of Negative Binomial Models
##
## Response: honey_pot
## Model theta Resid. df 2 x log-lik.
## 1 treatment + whole.mean + alive 17.17565 38 -221.3473
## 2 treatment + whole.mean + alive + duration 17.52117 37 -221.3263
## Test df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 1 0.02101159 0.8847474AIC(hp.m2, hp.m3, hp.m4, hp.m5)## df AIC
## hp.m2 12 243.3853
## hp.m3 11 241.4632
## hp.m4 9 239.3263
## hp.m5 8 237.3473Step 3 - Check model fit
qqnorm(resid(hp.m5));qqline(resid(hp.m5))
plot(resid(hp.m5)) +
abline(h=0, col="red", lwd=2) 
## integer(0)Step 4 - Results
Anova(hp.m5)## Analysis of Deviance Table (Type II tests)
##
## Response: honey_pot
## LR Chisq Df Pr(>Chisq)
## treatment 6.9398 4 0.139103
## whole.mean 9.7702 1 0.001774 **
## alive 5.5759 1 0.018209 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(hp.m5)##
## Call:
## glm.nb(formula = honey_pot ~ treatment + whole.mean + alive,
## data = brood_hp, init.theta = 17.17564907, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.27773 -1.05832 0.08925 0.62513 2.39932
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.25431 0.33133 0.768 0.44276
## treatment2 0.48891 0.23517 2.079 0.03762 *
## treatment3 0.04221 0.25015 0.169 0.86601
## treatment4 0.37416 0.23896 1.566 0.11740
## treatment5 0.28792 0.24734 1.164 0.24441
## whole.mean 1.34428 0.42736 3.146 0.00166 **
## alive 0.14517 0.06298 2.305 0.02116 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(17.1756) family taken to be 1)
##
## Null deviance: 87.307 on 44 degrees of freedom
## Residual deviance: 55.751 on 38 degrees of freedom
## AIC: 237.35
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 17.2
## Std. Err.: 16.1
##
## 2 x log-likelihood: -221.347hp.em <- emmeans(hp.m5, "treatment")
pairs(hp.em)## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 -0.4889 0.235 Inf -2.079 0.2292
## treatment1 - treatment3 -0.0422 0.250 Inf -0.169 0.9998
## treatment1 - treatment4 -0.3742 0.239 Inf -1.566 0.5195
## treatment1 - treatment5 -0.2879 0.247 Inf -1.164 0.7719
## treatment2 - treatment3 0.4467 0.219 Inf 2.035 0.2491
## treatment2 - treatment4 0.1147 0.210 Inf 0.546 0.9825
## treatment2 - treatment5 0.2010 0.216 Inf 0.931 0.8849
## treatment3 - treatment4 -0.3319 0.224 Inf -1.480 0.5758
## treatment3 - treatment5 -0.2457 0.231 Inf -1.061 0.8263
## treatment4 - treatment5 0.0862 0.224 Inf 0.384 0.9954
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatesStep 5 - Visualize
plot(hp.em, comparisons = TRUE)
#### Results Summary
hp.m5##
## Call: glm.nb(formula = honey_pot ~ treatment + whole.mean + alive,
## data = brood_hp, init.theta = 17.17564907, link = log)
##
## Coefficients:
## (Intercept) treatment2 treatment3 treatment4 treatment5 whole.mean
## 0.25431 0.48891 0.04221 0.37416 0.28792 1.34428
## alive
## 0.14517
##
## Degrees of Freedom: 44 Total (i.e. Null); 38 Residual
## Null Deviance: 87.31
## Residual Deviance: 55.75 AIC: 237.3hp_A <- as.data.frame(Anova(hp.m5))
hp_A <- setDT(hp_A)
hp_A## LR Chisq Df Pr(>Chisq)
## 1: 6.939816 4 0.13910307
## 2: 9.770212 1 0.00177362
## 3: 5.575943 1 0.01820885hp_emm <- emmeans(hp.m5, pairwise ~ treatment)
hp_contrasts <- as.data.frame(hp_emm$contrasts)
hp_contrasts <- setDT(hp_contrasts)
hp_contrasts## contrast estimate SE df z.ratio p.value
## 1: treatment1 - treatment2 -0.48890518 0.2351688 Inf -2.0789539 0.2291565
## 2: treatment1 - treatment3 -0.04220913 0.2501519 Inf -0.1687340 0.9998188
## 3: treatment1 - treatment4 -0.37415824 0.2389592 Inf -1.5657826 0.5194963
## 4: treatment1 - treatment5 -0.28791680 0.2473415 Inf -1.1640455 0.7719242
## 5: treatment2 - treatment3 0.44669605 0.2194870 Inf 2.0351819 0.2490958
## 6: treatment2 - treatment4 0.11474695 0.2103100 Inf 0.5456086 0.9825001
## 7: treatment2 - treatment5 0.20098838 0.2158809 Inf 0.9310152 0.8849405
## 8: treatment3 - treatment4 -0.33194911 0.2243466 Inf -1.4796265 0.5757829
## 9: treatment3 - treatment5 -0.24570767 0.2314810 Inf -1.0614592 0.8262678
## 10: treatment4 - treatment5 0.08624144 0.2243847 Inf 0.3843464 0.9953836hp_means <- as.data.frame(hp_emm$emmeans)
hp_means <- setDT(hp_means)
hp_means## treatment emmean SE df asymp.LCL asymp.UCL
## 1: 1 1.497060 0.1836461 Inf 1.137120 1.856999
## 2: 2 1.985965 0.1460636 Inf 1.699685 2.272244
## 3: 3 1.539269 0.1694482 Inf 1.207157 1.871381
## 4: 4 1.871218 0.1545253 Inf 1.568354 2.174082
## 5: 5 1.784977 0.1622127 Inf 1.467045 2.102908Proportion of Living Larvae
brood_l <- subset(brood, select = c(dead_larvae, live_larvae, treatment, whole.mean, alive, duration, qro, colony))Step 1 - Check for co-linearity
dl.model <- lm(dead_larvae ~ treatment + whole.mean + alive + duration + qro, data = brood_l)
vif(dl.model) #why am I just now realizing this is always going to be same answer for all response variables ## GVIF Df GVIF^(1/(2*Df))
## treatment 1.289595 4 1.032302
## whole.mean 1.746024 1 1.321372
## alive 1.833697 1 1.354141
## duration 1.300597 1 1.140437
## qro 1.705126 3 1.093015Step 2 - Model
larvae_bind <- glm(cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + alive + duration + qro, data = brood_l, family = binomial("logit"))
summary(larvae_bind)##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ treatment + whole.mean +
## alive + duration + qro, family = binomial("logit"), data = brood_l)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3783 -1.6202 0.2547 1.6725 4.2477
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.61656 0.70409 6.557 5.50e-11 ***
## treatment2 -1.29401 0.34956 -3.702 0.000214 ***
## treatment3 -0.82538 0.32215 -2.562 0.010404 *
## treatment4 -1.42772 0.34922 -4.088 4.34e-05 ***
## treatment5 0.04359 0.41910 0.104 0.917156
## whole.mean -0.36806 0.68021 -0.541 0.588442
## alive -0.42076 0.12715 -3.309 0.000935 ***
## duration 0.01530 0.01160 1.318 0.187420
## qroB3 -0.70932 0.35668 -1.989 0.046735 *
## qroB4 -1.05099 0.34160 -3.077 0.002093 **
## qroB5 -1.17233 0.24750 -4.737 2.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 329.89 on 42 degrees of freedom
## Residual deviance: 248.19 on 32 degrees of freedom
## AIC: 341.18
##
## Number of Fisher Scoring iterations: 5larvae_bind.int <- glm(cbind(live_larvae, dead_larvae) ~ treatment*whole.mean + alive + duration + qro, data = brood_l, family = binomial("logit"))
summary(larvae_bind.int)##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ treatment * whole.mean +
## alive + duration + qro, family = binomial("logit"), data = brood_l)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.144 -1.546 0.000 1.605 3.739
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.10019 1.31013 0.840 0.40104
## treatment2 -0.55785 1.54707 -0.361 0.71841
## treatment3 1.62555 1.07879 1.507 0.13186
## treatment4 3.33573 1.56370 2.133 0.03291 *
## treatment5 -0.12798 1.54365 -0.083 0.93393
## whole.mean 5.35298 2.03633 2.629 0.00857 **
## alive -0.33173 0.14786 -2.244 0.02486 *
## duration 0.01836 0.01343 1.367 0.17154
## qroB3 -0.93805 0.37836 -2.479 0.01317 *
## qroB4 -1.69871 0.43619 -3.894 9.84e-05 ***
## qroB5 -0.75880 0.33014 -2.298 0.02154 *
## treatment2:whole.mean -1.53975 2.68372 -0.574 0.56614
## treatment3:whole.mean -4.71595 1.93895 -2.432 0.01501 *
## treatment4:whole.mean -8.55810 2.70125 -3.168 0.00153 **
## treatment5:whole.mean 0.40336 2.99327 0.135 0.89281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 329.89 on 42 degrees of freedom
## Residual deviance: 230.73 on 28 degrees of freedom
## AIC: 331.72
##
## Number of Fisher Scoring iterations: 5AIC(larvae_bind, larvae_bind.int)## df AIC
## larvae_bind 11 341.1843
## larvae_bind.int 15 331.7192Anova(larvae_bind)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## treatment 34.467 4 5.977e-07 ***
## whole.mean 0.294 1 0.5874330
## alive 13.381 1 0.0002542 ***
## duration 1.684 1 0.1943949
## qro 23.229 3 3.618e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(larvae_bind, test = "Chisq")## Single term deletions
##
## Model:
## cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + alive +
## duration + qro
## Df Deviance AIC LRT Pr(>Chi)
## <none> 248.19 341.18
## treatment 4 282.66 367.65 34.467 5.977e-07 ***
## whole.mean 1 248.49 339.48 0.294 0.5874330
## alive 1 261.57 352.57 13.381 0.0002542 ***
## duration 1 249.87 340.87 1.684 0.1943949
## qro 3 271.42 358.41 23.229 3.618e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1larvae_bind1 <- glm(cbind(live_larvae, dead_larvae) ~ treatment + alive + qro, data = brood_l, family = binomial("logit"))
anova(larvae_bind, larvae_bind1)## Analysis of Deviance Table
##
## Model 1: cbind(live_larvae, dead_larvae) ~ treatment + whole.mean + alive +
## duration + qro
## Model 2: cbind(live_larvae, dead_larvae) ~ treatment + alive + qro
## Resid. Df Resid. Dev Df Deviance
## 1 32 248.19
## 2 34 249.90 -2 -1.712AIC(larvae_bind, larvae_bind1)## df AIC
## larvae_bind 11 341.1843
## larvae_bind1 9 338.8962summary(larvae_bind1)##
## Call:
## glm(formula = cbind(live_larvae, dead_larvae) ~ treatment + alive +
## qro, family = binomial("logit"), data = brood_l)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3333 -1.6428 0.1578 1.7751 4.4306
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.9851 0.6097 8.177 2.92e-16 ***
## treatment2 -1.1840 0.3379 -3.504 0.000459 ***
## treatment3 -0.7199 0.3105 -2.319 0.020412 *
## treatment4 -1.4046 0.3292 -4.267 1.98e-05 ***
## treatment5 0.1646 0.4063 0.405 0.685502
## alive -0.4274 0.1232 -3.470 0.000520 ***
## qroB3 -0.5639 0.3134 -1.799 0.071961 .
## qroB4 -1.0480 0.2632 -3.982 6.84e-05 ***
## qroB5 -1.1037 0.2326 -4.746 2.08e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 329.89 on 42 degrees of freedom
## Residual deviance: 249.90 on 34 degrees of freedom
## AIC: 338.9
##
## Number of Fisher Scoring iterations: 5qqnorm(resid(larvae_bind1));qqline(resid(larvae_bind1)) #Not bad! 
plot(resid(larvae_bind1)) +
abline(h=0, col="red", lwd=2) 
## integer(0)lm.emm <- emmeans(larvae_bind1, pairwise ~ treatment, type = "response")
pairs(lm.emm)## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.267 1.1040 Inf 1 3.504 0.0042
## treatment1 / treatment3 2.054 0.6378 Inf 1 2.319 0.1389
## treatment1 / treatment4 4.074 1.3410 Inf 1 4.267 0.0002
## treatment1 / treatment5 0.848 0.3447 Inf 1 -0.405 0.9944
## treatment2 / treatment3 0.629 0.1686 Inf 1 -1.730 0.4151
## treatment2 / treatment4 1.247 0.3577 Inf 1 0.769 0.9394
## treatment2 / treatment5 0.260 0.0953 Inf 1 -3.672 0.0022
## treatment3 / treatment4 1.983 0.4714 Inf 1 2.881 0.0324
## treatment3 / treatment5 0.413 0.1377 Inf 1 -2.653 0.0611
## treatment4 / treatment5 0.208 0.0688 Inf 1 -4.749 <.0001
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleStep 3 - Visualize and Results
lm.emm <- emmeans(larvae_bind1, pairwise ~ treatment, type = "response")
pairs(lm.emm)## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.267 1.1040 Inf 1 3.504 0.0042
## treatment1 / treatment3 2.054 0.6378 Inf 1 2.319 0.1389
## treatment1 / treatment4 4.074 1.3410 Inf 1 4.267 0.0002
## treatment1 / treatment5 0.848 0.3447 Inf 1 -0.405 0.9944
## treatment2 / treatment3 0.629 0.1686 Inf 1 -1.730 0.4151
## treatment2 / treatment4 1.247 0.3577 Inf 1 0.769 0.9394
## treatment2 / treatment5 0.260 0.0953 Inf 1 -3.672 0.0022
## treatment3 / treatment4 1.983 0.4714 Inf 1 2.881 0.0324
## treatment3 / treatment5 0.413 0.1377 Inf 1 -2.653 0.0611
## treatment4 / treatment5 0.208 0.0688 Inf 1 -4.749 <.0001
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scalesummary(lm.emm)## $emmeans
## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.928 0.0178 Inf 0.884 0.955
## 2 0.797 0.0362 Inf 0.717 0.859
## 3 0.862 0.0193 Inf 0.819 0.895
## 4 0.759 0.0369 Inf 0.679 0.823
## 5 0.938 0.0183 Inf 0.891 0.965
##
## Results are averaged over the levels of: qro
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## treatment1 / treatment2 3.267 1.1040 Inf 1 3.504 0.0042
## treatment1 / treatment3 2.054 0.6378 Inf 1 2.319 0.1389
## treatment1 / treatment4 4.074 1.3410 Inf 1 4.267 0.0002
## treatment1 / treatment5 0.848 0.3447 Inf 1 -0.405 0.9944
## treatment2 / treatment3 0.629 0.1686 Inf 1 -1.730 0.4151
## treatment2 / treatment4 1.247 0.3577 Inf 1 0.769 0.9394
## treatment2 / treatment5 0.260 0.0953 Inf 1 -3.672 0.0022
## treatment3 / treatment4 1.983 0.4714 Inf 1 2.881 0.0324
## treatment3 / treatment5 0.413 0.1377 Inf 1 -2.653 0.0611
## treatment4 / treatment5 0.208 0.0688 Inf 1 -4.749 <.0001
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio scaleAnova(larvae_bind1)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_larvae, dead_larvae)
## LR Chisq Df Pr(>Chisq)
## treatment 41.740 4 1.888e-08 ***
## alive 15.094 1 0.0001023 ***
## qro 27.616 3 4.372e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cldl <- cld(object = lm.emm,
adjust = "Tukey",
alpha = 0.05,
Letters = letters)
cldl <- as.data.frame(cldl)
cldl## treatment prob SE df asymp.LCL asymp.UCL .group
## 4 4 0.7585118 0.03689167 Inf 0.6518489 0.8404928 a
## 2 2 0.7966044 0.03619343 Inf 0.6881304 0.8742441 ab
## 3 3 0.8616744 0.01931800 Inf 0.8042284 0.9042703 bc
## 1 1 0.9275176 0.01775261 Inf 0.8665569 0.9618555 c
## 5 5 0.9378315 0.01833481 Inf 0.8705600 0.9712941 cemm3 <- emmeans(larvae_bind1 , ~ treatment, type = "response")
emm3## treatment prob SE df asymp.LCL asymp.UCL
## 1 0.928 0.0178 Inf 0.884 0.955
## 2 0.797 0.0362 Inf 0.717 0.859
## 3 0.862 0.0193 Inf 0.819 0.895
## 4 0.759 0.0369 Inf 0.679 0.823
## 5 0.938 0.0183 Inf 0.891 0.965
##
## Results are averaged over the levels of: qro
## 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 3.267 1.1040 Inf 1 3.504 0.0042
## treatment1 / treatment3 2.054 0.6378 Inf 1 2.319 0.1389
## treatment1 / treatment4 4.074 1.3410 Inf 1 4.267 0.0002
## treatment1 / treatment5 0.848 0.3447 Inf 1 -0.405 0.9944
## treatment2 / treatment3 0.629 0.1686 Inf 1 -1.730 0.4151
## treatment2 / treatment4 1.247 0.3577 Inf 1 0.769 0.9394
## treatment2 / treatment5 0.260 0.0953 Inf 1 -3.672 0.0022
## treatment3 / treatment4 1.983 0.4714 Inf 1 2.881 0.0324
## treatment3 / treatment5 0.413 0.1377 Inf 1 -2.653 0.0611
## treatment4 / treatment5 0.208 0.0688 Inf 1 -4.749 <.0001
##
## Results are averaged over the levels of: qro
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log odds ratio 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.5,1))
p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
labs(y = "Probability of Survival",) +
ggtitle("Larvae probability of being alive upon colony dissection") +
scale_x_discrete(name = "Treatment",
labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
theme_cowplot() +
theme_classic(base_size = 12)+
annotate(geom = "text",
x = 2, y = 1,
label = "P < 0.05",
size = 5) +
annotate(geom = "text",
x = c(4, 2, 3, 1, 5),
y = c(0.87, 0.89, 0.92, 0.98, 0.99),
label = c("a", "ab", "bc", "c", "c"),
size = 5) +
theme(legend.position = "none")
Results Summary
larvae_bind1##
## Call: glm(formula = cbind(live_larvae, dead_larvae) ~ treatment + alive +
## qro, family = binomial("logit"), data = brood_l)
##
## Coefficients:
## (Intercept) treatment2 treatment3 treatment4 treatment5 alive
## 4.9851 -1.1840 -0.7199 -1.4046 0.1646 -0.4274
## qroB3 qroB4 qroB5
## -0.5639 -1.0480 -1.1037
##
## Degrees of Freedom: 42 Total (i.e. Null); 34 Residual
## Null Deviance: 329.9
## Residual Deviance: 249.9 AIC: 338.9lar_A <- as.data.frame(Anova(larvae_bind1))
lar_A <- setDT(lar_A)
lar_A## LR Chisq Df Pr(>Chisq)
## 1: 41.74011 4 1.888427e-08
## 2: 15.09402 1 1.022869e-04
## 3: 27.61603 3 4.372447e-06lar_emm <- emmeans(larvae_bind1, pairwise ~ treatment)
lar_contrasts <- as.data.frame(lar_emm$contrasts)
lar_contrasts <- setDT(lar_contrasts)
lar_means <- as.data.frame(lar_emm$emmeans)
lar_means <- setDT(lar_means)
lar_means## treatment emmean SE df asymp.LCL asymp.UCL
## 1: 1 2.549168 0.2640629 Inf 2.0316140 3.066722
## 2: 2 1.365205 0.2233806 Inf 0.9273871 1.803023
## 3: 3 1.829267 0.1620752 Inf 1.5116058 2.146929
## 4: 4 1.144538 0.2014049 Inf 0.7497913 1.539284
## 5: 5 2.713722 0.3144716 Inf 2.0973694 3.330076lar_contrasts## contrast estimate SE df z.ratio p.value
## 1: treatment1 - treatment2 1.1839628 0.3378995 Inf 3.5038898 4.186710e-03
## 2: treatment1 - treatment3 0.7199005 0.3104781 Inf 2.3186839 1.388642e-01
## 3: treatment1 - treatment4 1.4046301 0.3291624 Inf 4.2672863 1.921229e-04
## 4: treatment1 - treatment5 -0.1645546 0.4063417 Inf -0.4049661 9.943529e-01
## 5: treatment2 - treatment3 -0.4640623 0.2681777 Inf -1.7304286 4.151211e-01
## 6: treatment2 - treatment4 0.2206674 0.2868448 Inf 0.7692917 9.394366e-01
## 7: treatment2 - treatment5 -1.3485174 0.3672080 Inf -3.6723525 2.236525e-03
## 8: treatment3 - treatment4 0.6847297 0.2376688 Inf 2.8810244 3.236834e-02
## 9: treatment3 - treatment5 -0.8844551 0.3333431 Inf -2.6532876 6.114583e-02
## 10: treatment4 - treatment5 -1.5691848 0.3304342 Inf -4.7488567 2.018816e-05Proportion of Living Pupae
brood_p <- subset(brood, select = c(dead_pupae, live_pupae, treatment, whole.mean, alive, duration, qro, colony))Step 1 - Model
pupae_bind <- glm(cbind(live_pupae, dead_pupae) ~ treatment + whole.mean + alive + duration + qro, data = brood_p, family = binomial("logit"))
summary(pupae_bind)##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ treatment + whole.mean +
## alive + duration + qro, family = binomial("logit"), data = brood_p)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1267 -1.4858 0.2463 1.2575 2.5505
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.00525 1.02143 1.963 0.049626 *
## treatment2 -1.44833 0.37385 -3.874 0.000107 ***
## treatment3 -1.71440 0.37957 -4.517 6.28e-06 ***
## treatment4 -0.78717 0.42115 -1.869 0.061609 .
## treatment5 -1.30054 0.45335 -2.869 0.004121 **
## whole.mean -0.54630 0.89542 -0.610 0.541792
## alive -0.36419 0.16580 -2.197 0.028051 *
## duration 0.02100 0.01717 1.223 0.221303
## qroB3 -0.78239 0.45515 -1.719 0.085624 .
## qroB4 0.38729 0.41808 0.926 0.354255
## qroB5 0.05733 0.45289 0.127 0.899259
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.65 on 40 degrees of freedom
## Residual deviance: 120.44 on 30 degrees of freedom
## AIC: 212.21
##
## Number of Fisher Scoring iterations: 5pupae_bind.int <- glm(cbind(live_pupae, dead_pupae) ~ treatment*whole.mean + alive + duration + qro, data = brood_p, family = binomial("logit"))
summary(pupae_bind.int)##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ treatment * whole.mean +
## alive + duration + qro, family = binomial("logit"), data = brood_p)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.576 -1.198 0.000 1.157 2.788
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.30229 1.56676 0.193 0.8470
## treatment2 -1.80387 1.84143 -0.980 0.3273
## treatment3 0.17761 1.43896 0.123 0.9018
## treatment4 1.76762 1.58768 1.113 0.2656
## treatment5 2.87056 1.73514 1.654 0.0981 .
## whole.mean 2.24436 1.82750 1.228 0.2194
## alive -0.41386 0.18412 -2.248 0.0246 *
## duration 0.02928 0.02223 1.317 0.1879
## qroB3 -0.61244 0.48549 -1.262 0.2071
## qroB4 0.40805 0.52970 0.770 0.4411
## qroB5 -0.03983 0.49472 -0.081 0.9358
## treatment2:whole.mean 0.61909 3.21872 0.192 0.8475
## treatment3:whole.mean -3.37593 2.49154 -1.355 0.1754
## treatment4:whole.mean -4.36209 2.64390 -1.650 0.0990 .
## treatment5:whole.mean -7.62423 3.14317 -2.426 0.0153 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.65 on 40 degrees of freedom
## Residual deviance: 109.20 on 26 degrees of freedom
## AIC: 208.97
##
## Number of Fisher Scoring iterations: 5AIC(pupae_bind, pupae_bind.int)## df AIC
## pupae_bind 11 212.2052
## pupae_bind.int 15 208.9660Anova(pupae_bind)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## treatment 25.4719 4 4.043e-05 ***
## whole.mean 0.3732 1 0.54125
## alive 5.3536 1 0.02068 *
## duration 1.5622 1 0.21134
## qro 5.4534 3 0.14145
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(pupae_bind, pupae_bind.int)## Analysis of Deviance Table
##
## Model 1: cbind(live_pupae, dead_pupae) ~ treatment + whole.mean + alive +
## duration + qro
## Model 2: cbind(live_pupae, dead_pupae) ~ treatment * whole.mean + alive +
## duration + qro
## Resid. Df Resid. Dev Df Deviance
## 1 30 120.44
## 2 26 109.20 4 11.239drop1(pupae_bind.int, test = "Chisq")## Single term deletions
##
## Model:
## cbind(live_pupae, dead_pupae) ~ treatment * whole.mean + alive +
## duration + qro
## Df Deviance AIC LRT Pr(>Chi)
## <none> 109.20 208.97
## alive 1 114.61 212.38 5.4087 0.02004 *
## duration 1 111.13 208.89 1.9256 0.16524
## qro 3 113.12 206.88 3.9174 0.27052
## treatment:whole.mean 4 120.44 212.21 11.2393 0.02400 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pupae_int1 <- glm(cbind(live_pupae, dead_pupae) ~ treatment*whole.mean + alive, data = brood_p, family = binomial("logit"))
anova(pupae_bind.int, pupae_int1)## Analysis of Deviance Table
##
## Model 1: cbind(live_pupae, dead_pupae) ~ treatment * whole.mean + alive +
## duration + qro
## Model 2: cbind(live_pupae, dead_pupae) ~ treatment * whole.mean + alive
## Resid. Df Resid. Dev Df Deviance
## 1 26 109.20
## 2 30 114.42 -4 -5.2214AIC(pupae_bind.int, pupae_int1)## df AIC
## pupae_bind.int 15 208.9660
## pupae_int1 11 206.1873summary(pupae_int1)##
## Call:
## glm(formula = cbind(live_pupae, dead_pupae) ~ treatment * whole.mean +
## alive, family = binomial("logit"), data = brood_p)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.471 -1.207 0.000 1.307 2.617
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5150 1.1622 1.304 0.1924
## treatment2 -2.1224 1.5646 -1.357 0.1749
## treatment3 -0.2798 1.2029 -0.233 0.8161
## treatment4 1.8973 1.5512 1.223 0.2213
## treatment5 2.7287 1.5571 1.752 0.0797 .
## whole.mean 2.4931 1.7797 1.401 0.1613
## alive -0.4725 0.1152 -4.101 4.11e-05 ***
## treatment2:whole.mean 1.6668 2.7280 0.611 0.5412
## treatment3:whole.mean -2.2551 1.9972 -1.129 0.2588
## treatment4:whole.mean -4.4230 2.5848 -1.711 0.0871 .
## treatment5:whole.mean -6.8751 2.8288 -2.430 0.0151 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 172.65 on 40 degrees of freedom
## Residual deviance: 114.42 on 30 degrees of freedom
## AIC: 206.19
##
## Number of Fisher Scoring iterations: 5qqnorm(resid(pupae_int1));qqline(resid(pupae_int1)) #Not bad! 
plot(resid(pupae_int1)) +
abline(h=0, col="red", lwd=2) 
## integer(0)pm.emm <- emmeans(pupae_int1, pairwise ~ treatment | whole.mean)
pairs(pm.emm)## whole.mean = 0.48:
## contrast estimate SE df z.ratio p.value
## treatment1 - treatment2 1.3215 0.415 Inf 3.186 0.0126
## treatment1 - treatment3 1.3633 0.414 Inf 3.297 0.0087
## treatment1 - treatment4 0.2279 0.489 Inf 0.466 0.9903
## treatment1 - treatment5 0.5748 0.430 Inf 1.338 0.6675
## treatment2 - treatment3 0.0418 0.369 Inf 0.113 1.0000
## treatment2 - treatment4 -1.0936 0.451 Inf -2.426 0.1085
## treatment2 - treatment5 -0.7468 0.388 Inf -1.924 0.3043
## treatment3 - treatment4 -1.1355 0.450 Inf -2.523 0.0856
## treatment3 - treatment5 -0.7886 0.385 Inf -2.048 0.2433
## treatment4 - treatment5 0.3469 0.465 Inf 0.746 0.9456
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimatesResults Summary
pupae_int1##
## Call: glm(formula = cbind(live_pupae, dead_pupae) ~ treatment * whole.mean +
## alive, family = binomial("logit"), data = brood_p)
##
## Coefficients:
## (Intercept) treatment2 treatment3
## 1.5150 -2.1224 -0.2798
## treatment4 treatment5 whole.mean
## 1.8973 2.7287 2.4931
## alive treatment2:whole.mean treatment3:whole.mean
## -0.4725 1.6668 -2.2551
## treatment4:whole.mean treatment5:whole.mean
## -4.4230 -6.8751
##
## Degrees of Freedom: 40 Total (i.e. Null); 30 Residual
## Null Deviance: 172.6
## Residual Deviance: 114.4 AIC: 206.2Anova(pupae_int1)## Analysis of Deviance Table (Type II tests)
##
## Response: cbind(live_pupae, dead_pupae)
## LR Chisq Df Pr(>Chisq)
## treatment 22.1077 4 0.0001908 ***
## whole.mean 0.1242 1 0.7245039
## alive 19.7012 1 9.055e-06 ***
## treatment:whole.mean 11.8444 4 0.0185468 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pup_A <- as.data.frame(Anova(pupae_int1))
pup_A <- setDT(pup_A)
pup_A## LR Chisq Df Pr(>Chisq)
## 1: 22.1076956 4 1.907658e-04
## 2: 0.1242181 1 7.245039e-01
## 3: 19.7011670 1 9.054606e-06
## 4: 11.8443827 4 1.854678e-02pup_emm <- emmeans(pupae_int1, pairwise ~ treatment)
pup_contrasts <- as.data.frame(pup_emm$contrasts)
pup_contrasts <- setDT(pup_contrasts)
pup_means <- as.data.frame(pup_emm$emmeans)
pup_means <- setDT(pup_means)
pup_means## treatment emmean SE df asymp.LCL asymp.UCL
## 1: 1 0.7703572 0.3270155 Inf 0.1294187 1.41129573
## 2: 2 -0.5511688 0.2532684 Inf -1.0475657 -0.05477182
## 3: 3 -0.5929809 0.2667758 Inf -1.1158519 -0.07010987
## 4: 4 0.5424762 0.3717617 Inf -0.1861633 1.27111560
## 5: 5 0.1956066 0.2921055 Inf -0.3769096 0.76812278pup_contrasts## contrast estimate SE df z.ratio p.value
## 1: treatment1 - treatment2 1.32152598 0.4148233 Inf 3.1857562 0.012568922
## 2: treatment1 - treatment3 1.36333813 0.4135072 Inf 3.2970115 0.008667262
## 3: treatment1 - treatment4 0.22788107 0.4886199 Inf 0.4663770 0.990308234
## 4: treatment1 - treatment5 0.57475061 0.4296308 Inf 1.3377778 0.667518151
## 5: treatment2 - treatment3 0.04181215 0.3692871 Inf 0.1132240 0.999962988
## 6: treatment2 - treatment4 -1.09364491 0.4508910 Inf -2.4255195 0.108494556
## 7: treatment2 - treatment5 -0.74677537 0.3880885 Inf -1.9242397 0.304305037
## 8: treatment3 - treatment4 -1.13545706 0.4500733 Inf -2.5228267 0.085585660
## 9: treatment3 - treatment5 -0.78858752 0.3851184 Inf -2.0476496 0.243308077
## 10: treatment4 - treatment5 0.34686954 0.4649569 Inf 0.7460252 0.945586806Proportion of Surviving Larvae and Pupae
pl_bind <- glm(cbind(alive_lp, dead_lp) ~ treatment + whole.mean + qro + duration + alive, data = brood, family = binomial("logit"))
summary(pl_bind )##
## Call:
## glm(formula = cbind(alive_lp, dead_lp) ~ treatment + whole.mean +
## qro + duration + alive, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.4463 -1.4406 0.0016 2.1953 4.4558
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.703903 0.522750 7.085 1.39e-12 ***
## treatment2 -1.574526 0.232860 -6.762 1.36e-11 ***
## treatment3 -1.000069 0.217843 -4.591 4.42e-06 ***
## treatment4 -0.988606 0.243484 -4.060 4.90e-05 ***
## treatment5 -0.589486 0.265445 -2.221 0.0264 *
## whole.mean -0.615161 0.489011 -1.258 0.2084
## qroB3 -0.439782 0.250719 -1.754 0.0794 .
## qroB4 -0.135490 0.236338 -0.573 0.5664
## qroB5 -0.381178 0.188001 -2.028 0.0426 *
## duration 0.012090 0.007963 1.518 0.1290
## alive -0.391459 0.091394 -4.283 1.84e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 411.27 on 42 degrees of freedom
## Residual deviance: 307.21 on 32 degrees of freedom
## AIC: 450.3
##
## Number of Fisher Scoring iterations: 4pl_bind.int <- glm(cbind(alive_lp, dead_lp) ~ treatment*whole.mean + qro + duration + alive, data = brood, family = binomial("logit"))
summary(pl_bind.int)##
## Call:
## glm(formula = cbind(alive_lp, dead_lp) ~ treatment * whole.mean +
## qro + duration + alive, family = binomial("logit"), data = brood)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.158 -1.587 0.000 2.063 3.952
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.918986 0.899447 1.022 0.30691
## treatment2 -0.711330 1.074012 -0.662 0.50777
## treatment3 1.059043 0.801206 1.322 0.18623
## treatment4 3.065617 1.048575 2.924 0.00346 **
## treatment5 2.227395 1.026291 2.170 0.02998 *
## whole.mean 3.670974 1.325106 2.770 0.00560 **
## qroB3 -0.435193 0.264866 -1.643 0.10037
## qroB4 -0.242295 0.285810 -0.848 0.39658
## qroB5 -0.012576 0.229488 -0.055 0.95630
## duration 0.010924 0.008764 1.247 0.21258
## alive -0.297905 0.100233 -2.972 0.00296 **
## treatment2:whole.mean -1.563899 1.904546 -0.821 0.41157
## treatment3:whole.mean -3.802220 1.422744 -2.672 0.00753 **
## treatment4:whole.mean -6.975948 1.772700 -3.935 8.31e-05 ***
## treatment5:whole.mean -5.303407 1.910637 -2.776 0.00551 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 411.27 on 42 degrees of freedom
## Residual deviance: 285.19 on 28 degrees of freedom
## AIC: 436.29
##
## Number of Fisher Scoring iterations: 5anova(pl_bind.int, pl_bind)## Analysis of Deviance Table
##
## Model 1: cbind(alive_lp, dead_lp) ~ treatment * whole.mean + qro + duration +
## alive
## Model 2: cbind(alive_lp, dead_lp) ~ treatment + whole.mean + qro + duration +
## alive
## Resid. Df Resid. Dev Df Deviance
## 1 28 285.19
## 2 32 307.21 -4 -22.016AIC(pl_bind, pl_bind.int)## df AIC
## pl_bind 11 450.3017
## pl_bind.int 15 436.2854drop1(pl_bind.int, test = "Chisq")## Single term deletions
##
## Model:
## cbind(alive_lp, dead_lp) ~ treatment * whole.mean + qro + duration +
## alive
## Df Deviance AIC LRT Pr(>Chi)
## <none> 285.19 436.29
## qro 3 287.90 433.00 2.7126 0.4380932
## duration 1 286.73 435.82 1.5359 0.2152266
## alive 1 294.72 443.81 9.5282 0.0020234 **
## treatment:whole.mean 4 307.21 450.30 22.0163 0.0001989 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pl_update <- glm(cbind(alive_lp, dead_lp) ~ treatment*whole.mean + duration + alive, data = brood, family = binomial("logit"))
pl_update2 <- glm(cbind(alive_lp, dead_lp) ~ treatment*whole.mean + alive, data = brood, family = binomial("logit"))
anova(pl_bind.int, pl_update, pl_update2)## Analysis of Deviance Table
##
## Model 1: cbind(alive_lp, dead_lp) ~ treatment * whole.mean + qro + duration +
## alive
## Model 2: cbind(alive_lp, dead_lp) ~ treatment * whole.mean + duration +
## alive
## Model 3: cbind(alive_lp, dead_lp) ~ treatment * whole.mean + alive
## Resid. Df Resid. Dev Df Deviance
## 1 28 285.19
## 2 31 287.90 -3 -2.71258
## 3 32 288.58 -1 -0.67913AIC(pl_bind.int, pl_update, pl_update2)## df AIC
## pl_bind.int 15 436.2854
## pl_update 12 432.9980
## pl_update2 11 431.6771qqnorm(resid(pl_update));qqline(resid(pl_update))
LS0tDQp0aXRsZTogIkVmZmVjdCBvZiBQcmlzdGluZSBvbiBCdW1ibGUgQmVlIE1pY3JvY29sb25pZXMgLSBCcm9vZCBQcm9kdWN0aW9uIg0KYXV0aG9yOiAiRW1pbHkgUnVubmlvbiINCmRhdGU6ICJEYXRhIENvbGxlY3RlZCAyMDIyIg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIHRvYzogdHJ1ZQ0KICAgIHRvY19kZXB0aDogNA0KICAgIG51bWJlcl9zZWN0aW9uczogZmFsc2UNCiAgICB0b2NfZmxvYXQ6IHRydWUNCiAgICB0aGVtZTogam91cm5hbA0KICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCi0tLQ0KDQpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0NCmtuaXRyOjpvcHRzX2NodW5rJHNldChtZXNzYWdlID0gRkFMU0UpDQpgYGANCg0KYGBge3IgbG9hZCBsaWJyYXJpZXMsIGluY2x1ZGU9RkFMU0V9DQpsaWJyYXJ5KHJlYWRyKQ0KbGlicmFyeShrYWJsZUV4dHJhKQ0KbGlicmFyeShzdGF0cykNCmxpYnJhcnkoZ2dwbG90MikNCmxpYnJhcnkoY2FyKQ0KbGlicmFyeShlbW1lYW5zKQ0KbGlicmFyeShNQVNTKQ0KbGlicmFyeShsbWU0KQ0KbGlicmFyeShibG1lY28pDQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkoZHBseXIpDQpsaWJyYXJ5KGNvd3Bsb3QpDQpsaWJyYXJ5KGJlc3ROb3JtYWxpemUpDQpsaWJyYXJ5KHBsb3RseSkNCmxpYnJhcnkoYWdyaWNvbGFlKSANCmxpYnJhcnkoZ2dwdWJyKQ0KbGlicmFyeShnbHVlKQ0KbGlicmFyeShtdWx0Y29tcCkNCmxpYnJhcnkobXVsdGNvbXBWaWV3KQ0KbGlicmFyeShnbG1tVE1CKQ0KbGlicmFyeShyc3RhdGl4KQ0KbGlicmFyeShmaXRkaXN0cnBsdXMpDQpsaWJyYXJ5KGxvZ3NwbGluZSkNCmxpYnJhcnkob2xzcnIpDQpsaWJyYXJ5KEdHYWxseSkNCmxpYnJhcnkoZGF0YS50YWJsZSkNCmBgYA0KDQpCZWVzIGFyZSBmcmVxdWVudGx5IGV4cG9zZWQgdG8gZnVuZ2ljaWRlcyBpbiBhZ3JpY3VsdHVyYWwgbGFuZHNjYXBlcywgYW5kIHdoaWxlIHRoZXNlIGNoZW1pY2FscyBhcmUgZ2VuZXJhbGx5IG5vdCBjb25zaWRlcmVkIHRvIGJlIGhhcm1mdWwgdG8gaW5zZWN0IHBvbGxpbmF0b3JzLCB0aGUgc3VibGV0aGFsIGVmZmVjdHMgb2YgZnVuZ2ljaWRlcyBhcmUgbm90IHdlbGwgdW5kZXJzdG9vZC4gV2UgaW52ZXN0aWdhdGVkIHRoZSBub24tdGFyZ2V0IGVmZmVjdHMgb2YgZXhwb3N1cmUgdG8gZmllbGQtcmVhbGlzdGljIGNvbmNlbnRyYXRpb25zIG9mIHRoZSBmdW5naWNpZGUsIFByaXN0aW5lwq4sIGZvciB0aGUgY29tbW9uIGVhc3Rlcm4gYnVtYmxlIGJlZSAoKkJvbWJ1cyBpbXBhdGllbnMqKS4gV2UgZmVkIHF1ZWVubGVzcyBtaWNyb2NvbG9uaWVzIG9mIDUgYnVtYmxlIGJlZXMgcG9sbGVuIHBhdHRpZXMgY29udGFtaW5hdGVkIHdpdGggUHJpc3RpbmXCriByYW5naW5nIGZyb20gMCBwcGItMTUwLDAwMCBwcGIsIGFuZCBldmFsdWF0ZWQgdGhlIGVmZmVjdHMgb2YgY2hyb25pYyBmdW5naWNpZGUgZXhwb3N1cmUgb24gYnJvb2QgcHJvZHVjdGlvbiwgY29sb255IHdlaWdodCBjaGFuZ2UsIHBvbGxlbiBjb25zdW1wdGlvbiwgYW5kIGRyb25lIGJvZHkgcXVhbGl0eS4gV2UgZm91bmQgdGhhdCB3aGlsZSBtaWNyb2NvbG9uaWVzIGZlZCBwb2xsZW4gY29udGFpbmluZyAxLDUwMCBwcGIgY29uc3VtZWQgdGhlIG1vc3QgcG9sbGVuIGFuZCBwcm9kdWNlZCB0aGUgbW9zdCBicm9vZCwgbWljcm9jb2xvbmllcyBleHBvc2VkIHRvIDE1LDAwMCBwcGIgZ2FpbmVkIHdlaWdodCBhdCBzaWduaWZpY2FudGx5IGxvd2VyIHJhdGVzLiBGdW5naWNpZGUgZXhwb3N1cmUgYWxzbyBuZWdhdGl2ZWx5IGltcGFjdGVkIGJyb29kIGRldmVsb3BtZW50IHRpbWUgYW5kIGRyb25lIGJvZHkgcXVhbGl0eS4gT3VyIGZpbmRpbmdzIGluZGljYXRlIHRoYXQgUHJpc3RpbmXCriBoYXMgc3VibGV0aGFsIGVmZmVjdHMgb24gYnVtYmxlIGJlZSBtaWNyb2NvbG9ueSBncm93dGggYW5kIGJyb29kIGRldmVsb3BtZW50LCBoaWdobGlnaHRpbmcgdGhlIHBvdGVudGlhbCBuZWdhdGl2ZSBpbXBhY3RzIG9mIGNocm9uaWMgZnVuZ2ljaWRlIGV4cG9zdXJlIHRvIGJlZSBoZWFsdGguIEluIHRoZSBmYWNlIG9mIHRoZSBnbG9iYWwgZGVjbGluZSBvZiBwb2xsaW5hdG9yIHNwZWNpZXMsIGl0IGlzIGltcG9ydGFudCBmb3IgdXMgdG8gYnJvYWRlbiBvdXIgdW5kZXJzdGFuZGluZyBvZiBob3cgYmVlcyByZWFjdCB0b3dhcmRzIGV4cG9zdXJlIG9mIG9mdGVuIG92ZXJsb29rZWQgc3RyZXNzb3JzLCBzdWNoIGFzIGZ1bmdpY2lkZXMsIHRvIGFjY291bnQgZm9yIHBvdGVudGlhbCBjb25zZXF1ZW5jZXMgZm9yIHBvbGxpbmF0b3IgaGVhbHRoIGFuZCBwb2xsaW5hdGlvbiBzZXJ2aWNlcy4NCg0KIyMjIElucHV0IHJlbGV2ZW50IGRhdGEgZmlsZXMgDQoNCmBgYHtyfQ0KIyMjIEZpZ3VyZSBvdXQgYXZlcmFnZSBwb2xsZW4gY29uc3VtcHRpb24gYnkgdHJlYXRtZW50IA0KDQoNCnBvbGxlbiA8LSByZWFkX2NzdigicG9sbGVuMS5jc3YiLCBjb2xfdHlwZXMgPSBjb2xzKHJvdW5kID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCIxIiwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiMiIpKSwgdHJlYXRtZW50ID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCIxIiwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICIyIiwgIjMiLCAiNCIsICI1IikpLCByZXBsaWNhdGUgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjEiLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiMiIsICIzIiwgIjQiLCAiNSIsICI2IiwgIjciLCAiOSIsICIxMSIsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICIxMiIpKSwgc3RhcnRfZGF0ZSA9IGNvbF9kYXRlKGZvcm1hdCA9ICIlbS8lZC8lWSIpLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0YXJ0X3RpbWUgPSBjb2xfY2hhcmFjdGVyKCksIGVuZF9kYXRlID0gY29sX2RhdGUoZm9ybWF0ID0gIiVtLyVkLyVZIiksIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZW5kX3RpbWUgPSBjb2xfY2hhcmFjdGVyKCkpKQ0KDQoNCnBvbGxlbiRjb2xvbnkgPC0gYXMuZmFjdG9yKHBvbGxlbiRjb2xvbnkpDQpwb2xsZW4kcGlkIDwtIGFzLmZhY3Rvcihwb2xsZW4kcGlkKQ0KcG9sbGVuJGNvdW50IDwtIGFzLmZhY3Rvcihwb2xsZW4kY291bnQpDQoNCnBvbGxlbiA8LSBzdWJzZXQocG9sbGVuLCBwb2xsZW4kcm91bmQgIT0gMSkNCg0KcG9sbGVuIDwtIG5hLm9taXQocG9sbGVuKQ0KDQpyYW5nZShwb2xsZW4kZGlmZmVyZW5jZSkNCg0KIyBnZXQgcmlkIG9mIG5lZ2F0aXZlIG51bWJlcnMNCnBvbGxlbiRkaWZmZXJlbmNlW3BvbGxlbiRkaWZmZXJlbmNlIDwgMF0gPC0gTkENCnBvbGxlbiA8LSBuYS5vbWl0KHBvbGxlbikNCnJhbmdlKHBvbGxlbiRkaWZmZXJlbmNlKQ0KDQoNCiMgYWRkIHF1ZWVucmlnaHQgb3JpZ2luYWwgY29sb255IGNvbHVtbiANCnFybyA8LSByZWFkX2NzdigicXJvLmNzdiIpDQpxcm8kY29sb255IDwtIGFzLmZhY3Rvcihxcm8kY29sb255KQ0KcXJvJHFybyA8LSBhcy5mYWN0b3IocXJvJHFybykNCg0KcG9sbGVuIDwtIG1lcmdlKHBvbGxlbiwgcXJvLCBieS54ID0gImNvbG9ueSIpDQoNCnBvbGxlbiRwaWQgPC0gYXMuZmFjdG9yKHBvbGxlbiRwaWQpDQpwb2xsZW4kZG9zZV9jb25zdW1lZCA8LSBhcy5udW1lcmljKHBvbGxlbiRkb3NlX2NvbnN1bWVkKQ0KDQpgYGANCg0KIyMjIEF2ZXJhZ2UgcG9sbGVuIGNvbnN1bXB0aW9uIHBlciBjb2xvbnkgZGF0YSBpbnB1dA0KDQpgYGB7cn0NCnBvbGxlbiR3aG9sZV9kaWYgPC0gYXMubnVtZXJpYyhwb2xsZW4kZGlmZmVyZW5jZSkNCndkLjEgPC0gbG0oZGlmZmVyZW5jZSB+IHRyZWF0bWVudCwgZGF0YSA9IHBvbGxlbikNCnN1bW1hcnkod2QuMSkNCndkLmVtbSA8LSBlbW1lYW5zKHdkLjEsICJ0cmVhdG1lbnQiKQ0Kc3VtbWFyeSh3ZC5lbW0pDQoNCndkLnN1bW1hcnkgPC0gcG9sbGVuICU+JSANCiAgZ3JvdXBfYnkoY29sb255KSAlPiUNCiAgc3VtbWFyaXplKHdob2xlLm1lYW4gPSBtZWFuKGRpZmZlcmVuY2UpLA0KICAgICAgICAgICAgbWVhbi5kb3NlID0gbWVhbihkb3NlX2NvbnN1bWVkKSkNCg0KYXMuZGF0YS5mcmFtZSh3ZC5zdW1tYXJ5KSAgIyB0aGlzIGlzIHRoZSBkYXRhIGZyYW1lIEkgd2lsbCBtZXJnZSB3aXRoIHN1YnNlcXVlbnQgZGF0YSBmcmFtZXMgdG8gY29udGFpbiBhdmVyYWdlIHBvbGxlbiBjb25zdW1wdGlvbiBwZXIgY29sb255IGFzIGEgbmV3IGNvbHVtbiAgDQoNCmBgYA0KDQojIyMgV29ya2VyIE1vcnRhbGl0eSBkYXRhIGlucHV0IA0KDQpgYGB7cn0NCm1vcnRhbGl0eSAgPC0gcmVhZF9jc3YoInN1cnZpdmluZyB3b3JrZXJzIHBlciBjb2xvbnkuY3N2IikNCg0KbW9ydGFsaXR5JGNvbG9ueSA8LSBhcy5mYWN0b3IobW9ydGFsaXR5JGNvbG9ueSkNCmBgYA0KDQoNCiMgQmVnaW4gZGF0YSBhbmFseXNpcyANCg0KYGBge3IgaW5wdXQgYnJvb2QgZGF0YX0NCmJyb29kIDwtIHJlYWRfY3N2KCJicm9vZC5jc3YiLCBjb2xfdHlwZXMgPSBjb2xzKHJvdW5kID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCIxIiwgDQogICAgIjIiKSksIHRyZWF0bWVudCA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMSIsIA0KICAgICIyIiwgIjMiLCAiNCIsICI1IikpLCByZXBsaWNhdGUgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjEiLCANCiAgICAiMiIsICIzIiwgIjQiLCAiNSIsICI3IiwgIjkiLCAiMTEiLCAiMTIiKSkpKQ0KDQpicm9vZCRjb2xvbnkgPC0gYXMuZmFjdG9yKGJyb29kJGNvbG9ueSkNCg0KYnJvb2QgPC0gc3Vic2V0KGJyb29kLCBicm9vZCRyb3VuZCAhPSAxKQ0KDQpicm9vZCA8LSBtZXJnZSh3ZC5zdW1tYXJ5LCBicm9vZCwgYnkueCA9ICJjb2xvbnkiKQ0KDQpicm9vZCA8LSBtZXJnZShicm9vZCwgbW9ydGFsaXR5LCBieS54ID0gImNvbG9ueSIpDQoNCmJyb29kJHFybyA8LSBhcy5mYWN0b3IoYnJvb2QkcXJvKQ0KDQpgYGANCg0KDQojIyMgVG90YWwgQnJvb2QgQ2VsbHMNCg0KYGBge3J9DQoNCmdncGxvdChicm9vZCwgYWVzKHg9YnJvb2RfY2VsbHMsIGZpbGwgPSB0cmVhdG1lbnQpKSArDQogIGdlb21faGlzdG9ncmFtKHBvc2l0aW9uID0gImlkZW50aXR5IiwgYmlud2lkdGggPSA1ICxjb2w9SSgiYmxhY2siKSkgKw0KICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjKCJwZWFjaHB1ZmYzIiwgImRhcmtzZWFncmVlbiIsICJsaWdodHNlYWdyZWVuIiwgImRhcmtvbGl2ZWdyZWVuMyIsICJkYXJrb2xpdmVncmVlbjQiKSwNCiAgICAgICAgICAgICAgICAgICAgbmFtZSA9ICJQcmlzdGluZSBMZXZlbCIsDQogICAgICAgICAgICAgICAgICAgIGxhYmVscyA9IGMoIlRyZWF0bWVudCAxIChjb250cm9sKSIsICJUcmVhdG1lbnQgMiIsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJUcmVhdG1lbnQgMyIsICJUcmVhdG1lbnQgNCIsICJUcmVhdG1lbnQgNSIpKSArDQogIGdndGl0bGUoIlRvdGFsIEJyb29kIENlbGxzIikgKw0KICBsYWJzKHkgPSAiQ291bnQiLCB4ID0gIkJyb29kIENlbGxzIikNCmBgYA0KDQpgYGB7cn0NCmRlc2NkaXN0KGJyb29kJGJyb29kX2NlbGxzLCBkaXNjcmV0ZSA9IFRSVUUpDQpgYGANCg0KDQojIyMjIFN0ZXAgMSAtIENoZWNrIGZvciBjb2xpbmVhcml0eSANCg0KYGBge3J9DQpicm9vZC5tb2RlbCA8LSBsbShicm9vZF9jZWxsc34gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gKyBxcm8sIGRhdGEgPSBicm9vZCkNCnZpZihicm9vZC5tb2RlbCkNCmBgYA0KDQoNCiMjIyMgU3RlcCAyIC0gV29yayBkb3duIGZyb20gbW9zdCBjb21wbGljYXRlZCBtb2RlbA0KDQpgYGB7cn0NCmJyb29kLnBvaXMgPC0gZ2xtKGJyb29kX2NlbGxzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgcXJvICsgZHVyYXRpb24sIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gInBvaXNzb24iKQ0Kc3VtbWFyeShicm9vZC5wb2lzKSAgICNvdmVyZGlzcGVyc2VkIA0KDQpicm9vZC5pbnQgPC0gZ2xtLm5iKGJyb29kX2NlbGxzIH4gdHJlYXRtZW50Kndob2xlLm1lYW4gKyBhbGl2ZSArIHFybyArIGR1cmF0aW9uLCBkYXRhID0gYnJvb2QsIGNvbnRyb2wgPSBnbG0uY29udHJvbChtYXhpdCA9IDUwKSkgI0tlZXAgdGhpcyBtb2RlbCANCnN1bW1hcnkoYnJvb2QuaW50KQ0KDQpicm9vZC5nMSA8LSBnbG0ubmIoYnJvb2RfY2VsbHMgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyBxcm8gKyBkdXJhdGlvbiwgZGF0YSA9IGJyb29kLCBjb250cm9sID0gZ2xtLmNvbnRyb2wobWF4aXQgPSA1MCkpDQpzdW1tYXJ5KGJyb29kLmludCkNCg0KYnJvb2QuZzIgPC0gZ2xtLm5iKGJyb29kX2NlbGxzfiB0cmVhdG1lbnQgKyBhbGl2ZSAgKyBxcm8gKyBkdXJhdGlvbiwgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShicm9vZC5nMikNCg0KYnJvb2QuZzMgPC0gZ2xtLm5iKGJyb29kX2NlbGxzfiB3aG9sZS5tZWFuICsgYWxpdmUgKyBxcm8gKyBkdXJhdGlvbiwgZGF0YSA9IGJyb29kKQ0Kc3VtbWFyeShicm9vZC5nMykNCmBgYA0KDQpgYGB7cn0NCmFub3ZhKGJyb29kLmludCwgYnJvb2QuZzEpDQpBSUMoYnJvb2QuaW50LCBicm9vZC5nMSkNCg0KYW5vdmEoYnJvb2QuZzIsIGJyb29kLmczKQ0KQUlDKGJyb29kLmcyLCBicm9vZC5nMykNCg0KYW5vdmEoYnJvb2QuZzEsIGJyb29kLmczKQ0KDQphbm92YShicm9vZC5pbnQsIGJyb29kLmczKQ0KYGBgDQojIyMjIFN0ZXAgMmEgLSBEcm9wIDENCg0KYGBge3J9DQpkcm9wMShicm9vZC5pbnQsIHRlc3QgPSAiQ2hpc3EiKSAgICN0cnkgZ2V0dGluZyByaWQgb2YgZHVyYXRpb24NCg0KYnJvb2QuaW50LnVwZGF0ZSA8LSB1cGRhdGUoYnJvb2QuaW50LCAufi4gLWR1cmF0aW9uKQ0KYW5vdmEoYnJvb2QuaW50LCBicm9vZC5pbnQudXBkYXRlKQ0KQUlDKGJyb29kLmludCwgYnJvb2QuaW50LnVwZGF0ZSkgICAgI2tlZXAgZXZlcnl0aGluZyANCg0KYGBgDQoNCg0KDQojIyMjIFN0ZXAgMyAtIENoZWNrIG1vZGVsIGZpdA0KDQpgYGB7cn0NCnFxbm9ybShyZXNpZChicm9vZC5pbnQpKTtxcWxpbmUocmVzaWQoYnJvb2QuaW50KSkNCiNub3QgdGVycmlibGUgDQpgYGANCg0KYGBge3J9DQpzdW1tYXJ5KGJyb29kLmludCkNCmBgYA0KDQoNCiMjIyMgU3RlcCA0IC0gUmVzdWx0cyANCg0KYGBge3J9DQpBbm92YShicm9vZC5pbnQpDQpCQSA8LSBhcy5kYXRhLmZyYW1lKEFub3ZhKGJyb29kLmludCkpDQpCQSA8LSBzZXREVChCQSkNCkJBDQpgYGANCg0KDQpgYGB7cn0NCmJyb29kX2VtbSA8LSBlbW1lYW5zKGJyb29kLmludCwgcGFpcndpc2UgfiB0cmVhdG1lbnQgfCB3aG9sZS5tZWFuKQ0KcGFpcnMoYnJvb2RfZW1tKQ0KDQpwbG90KGJyb29kX2VtbSwgY29tcGFyaXNvbnMgPSBUUlVFKQ0KDQpwbG90KGJyb29kJHRyZWF0bWVudCwgYnJvb2QkYnJvb2RfY2VsbHMpDQpgYGANCg0KYGBge3J9DQpicm9vZC5zdW0gPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcmlzZShicm9vZC5tID0gbWVhbihicm9vZF9jZWxscyksIA0KICAgICAgICAgICAgYnJvb2Quc2QgPSBzZChicm9vZF9jZWxscyksDQogICAgICAgICAgICBicm9vZC5uID0gbGVuZ3RoKGJyb29kX2NlbGxzKSkgJT4lDQogIG11dGF0ZShicm9vZC5zZSA9IGJyb29kLnNkL3NxcnQoYnJvb2QubikpDQpicm9vZC5zdW0NCmBgYA0KDQoNCiMjIyMgU3RlcCA1IC0gVmlzdWFsaXplDQoNCmBgYHtyfQ0KICAgICAgICANCmJyb29kdHVrLm1lYW5zIDwtIGVtbWVhbnMob2JqZWN0ID0gYnJvb2QuaW50LA0KICAgICAgICAgICAgICAgICAgICAgICAgc3BlY3MgPSAidHJlYXRtZW50IiwNCiAgICAgICAgICAgICAgICAgICAgICAgIGFkanVzdCA9ICJUdWtleSIsDQogICAgICAgICAgICAgICAgICAgICAgICB0eXBlID0gInJlc3BvbnNlIikNCg0KYnJvb2R0dWsubWVhbnMgPC0gZW1tZWFucyhvYmplY3QgPSBicm9vZC5pbnQsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICBzcGVjcyA9ICJ0cmVhdG1lbnQiLA0KICAgICAgICAgICAgICAgICAgICAgICAgICBhZGp1c3QgPSAic2lkYWsiLA0KICAgICAgICAgICAgICAgICAgICAgICAgICB0eXBlID0gInJlc3BvbnNlIikNCg0KYnJvb2QuY2xkLm1vZGVsIDwtIGNsZChvYmplY3QgPSBicm9vZHR1ay5tZWFucywNCiAgICAgICAgICAgICAgICAgICAgIGFkanVzdCA9ICJUdWtleSIsDQogICAgICAgICAgICAgICAgICAgICBMZXR0ZXJzID0gbGV0dGVycywNCiAgICAgICAgICAgICAgICAgICAgIGFscGhhID0gMC4wNSkNCmJyb29kLmNsZC5tb2RlbA0KDQpicm9vZHR1ay5tZWFucyA8LSBhcy5kYXRhLmZyYW1lKGJyb29kdHVrLm1lYW5zKQ0KDQpicm9vZF9tYXggPC0gYnJvb2QgJT4lDQogIGdyb3VwX2J5KHRyZWF0bWVudCkgJT4lDQogIHN1bW1hcml6ZShtYXhicm9vZCA9IG1heCgoYnJvb2RfY2VsbHMpKSkNCg0KDQpicm9vZF9mb3JfcGxvdHRpbmcgPC0gZnVsbF9qb2luKGJyb29kdHVrLm1lYW5zLCBicm9vZF9tYXgsDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBieT0idHJlYXRtZW50IikNCg0KYnJvb2RfZm9yX3Bsb3R0aW5nJG1heGIgPC0gYnJvb2RfZm9yX3Bsb3R0aW5nJHJlc3BvbnNlICsgYnJvb2RfZm9yX3Bsb3R0aW5nJFNFDQpgYGANCg0KDQoNCmBgYHtyLCBmaWcud2lkdGg9IDEwLCBmaWcuaGVpZ2h0PSAxMn0NCmdncGxvdChkYXRhID0gYnJvb2RfZm9yX3Bsb3R0aW5nLCBhZXMoeCA9IHRyZWF0bWVudCwgeSA9IHJlc3BvbnNlLCBmaWxsID0gdHJlYXRtZW50KSkgKw0KICBnZW9tX2NvbChjb2wgPSAiYmxhY2siKSsNCiAgY29vcmRfY2FydGVzaWFuKHlsaW09YygyMCwgNDApKSArDQogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGMoInBlYWNocHVmZjMiLCAiZGFya3NlYWdyZWVuIiwgImxpZ2h0c2VhZ3JlZW4iLCAiZGFya29saXZlZ3JlZW4zIiwgImRhcmtvbGl2ZWdyZWVuNCIpLA0KICAgICAgICAgICAgICAgICAgICBuYW1lID0gIlByaXN0aW5lIExldmVsIiwNCiAgICAgICAgICAgICAgICAgICAgbGFiZWxzID0gYygiVHJlYXRtZW50IDEgKGNvbnRyb2wpIiwgIlRyZWF0bWVudCAyIiwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIlRyZWF0bWVudCAzIiwgIlRyZWF0bWVudCA0IiwgIlRyZWF0bWVudCA1IikpICsgDQogIGdlb21fZXJyb3JiYXIoYWVzKHltaW4gPSByZXNwb25zZS1TRSwgDQogICAgICAgICAgICAgICAgICAgIHltYXggPSByZXNwb25zZStTRSksDQogICAgICAgICAgICAgICAgcG9zaXRpb24gPSBwb3NpdGlvbl9kb2RnZSgwLjkpLCB3aWR0aCA9IDAuOCkgKw0KICBsYWJzKHkgPSAiRGF5cyIsKSArDQogIGdndGl0bGUoIkF2ZXJhZ2UgQ29sb255IER1cmF0aW9uIHBlciBUcmVhdG1lbnQiKSArDQogIHNjYWxlX3hfZGlzY3JldGUobmFtZSA9ICJUcmVhdG1lbnQiLCANCiAgICAgICAgICAgICAgICAgICBsYWJlbHMgPSBjKCIwIFBQQiIsICIxNTAgUFBCIiwgIjEsNTAwIFBQQiIsICIxNSwwMDAgUFBCIiwgIjE1MCwwMDAgUFBCIikpICsNCiAgdGhlbWVfY293cGxvdCgpKw0KICB0aGVtZV9jbGFzc2ljKGJhc2Vfc2l6ZSA9IDMwKSArDQogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKw0KICBhbm5vdGF0ZShnZW9tID0gInRleHQiLCANCiAgICAgICAgICB4ID0gMSwgeSA9IDQwLA0KICAgICAgICAgIGxhYmVsID0gIlAgPiAwLjA1IiwNCiAgICAgICAgICBzaXplID0gMTIpICsNCiAgYW5ub3RhdGUoZ2VvbSA9ICJ0ZXh0IiwNCiAgICAgICAgICAgeCA9IGMoMSwgMiwgMywgNCwgNSksDQogICAgICAgICAgIHkgPSBjKGJyb29kX2Zvcl9wbG90dGluZyRtYXhiICsgMSksDQogICAgICAgICAgIGxhYmVsID0gYygiYSIsICJhIiwgImEiLCAiYSIsICJhIiksDQogICAgICAgICAgIHNpemUgPSAxMikgKw0KICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAgIm5vbmUiKQ0KYGBgDQoNCg0KIyMjIyBCcm9vZCBDZWxscyBTdW1tYXJ5IFJlc3VsdHMgDQoNCmBgYHtyfQ0KDQpiX0EgPC0gc2V0RFQoYXMuZGF0YS5mcmFtZShBbm92YShicm9vZC5pbnQpKSkNCmJfQQ0KYl9lbW0gPC0gZW1tZWFucyhicm9vZC5pbnQsIHBhaXJ3aXNlIH4gdHJlYXRtZW50IHwgd2hvbGUubWVhbikNCmJfY29udHJhc3RzIDwtIHNldERUKGFzLmRhdGEuZnJhbWUoYl9lbW0kY29udHJhc3RzKSkNCmJfbWVhbnMgPC0gYXMuZGF0YS5mcmFtZShiX2VtbSRlbW1lYW5zKQ0KYl9tZWFucyA8LSBzZXREVChiX21lYW5zKQ0KYl9tZWFucw0KYl9jb250cmFzdHMNCg0KYGBgDQoNCg0KDQojIyMgRWdncw0KDQpgYGB7cn0NCg0KYnJvb2RfZWdnIDwtIHN1YnNldChicm9vZCwgc2VsZWN0ID0gYyhlZ2dzLCB0cmVhdG1lbnQsIHdob2xlLm1lYW4sIGFsaXZlLCBkdXJhdGlvbiwgcXJvLCBjb2xvbnkpKQ0KDQpnZ3Bsb3QoYnJvb2RfZWdnLCBhZXMoeD1lZ2dzLCBmaWxsID0gdHJlYXRtZW50KSkgKw0KICBnZW9tX2hpc3RvZ3JhbShwb3NpdGlvbiA9ICJpZGVudGl0eSIsIGJpbndpZHRoID0xICxjb2w9SSgiYmxhY2siKSkgKw0KICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjKCJwZWFjaHB1ZmYzIiwgImRhcmtzZWFncmVlbiIsICJsaWdodHNlYWdyZWVuIiwgImRhcmtvbGl2ZWdyZWVuMyIsICJkYXJrb2xpdmVncmVlbjQiKSwNCiAgICAgICAgICAgICAgICAgICAgbmFtZSA9ICJQcmlzdGluZSBMZXZlbCIsDQogICAgICAgICAgICAgICAgICAgIGxhYmVscyA9IGMoIlRyZWF0bWVudCAxIChjb250cm9sKSIsICJUcmVhdG1lbnQgMiIsIA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJUcmVhdG1lbnQgMyIsICJUcmVhdG1lbnQgNCIsICJUcmVhdG1lbnQgNSIpKSArDQogIGdndGl0bGUoIkVnZ3MgcGVyIFRyZWF0bWVudCIpICsNCiAgbGFicyh5ID0gIkNvdW50IiwgeCA9ICJFZ2cgQ291bnQiKQ0KDQpgYGANCg0KDQojIyMjIFN0ZXAgMSAtIENoZWNrIGZvciBjby1saW5lYXJpdHkNCg0KYGBge3J9DQplZ2cubW9kZWwgPC0gbG0oZWdnc34gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gICsgcXJvLCBkYXRhID0gYnJvb2RfZWdnKQ0KdmlmKGVnZy5tb2RlbCkNCmBgYA0KDQoNCiMjIyMgU3RlcCAyIC0gV29yayBkb3duIGZyb20gbW9zdCBjb21wbGljYXRlZCBtb2RlbCANCg0KDQpgYGB7cn0NCmVnZy5pbnQgPC0gZ2xtKGVnZ3MgfiB0cmVhdG1lbnQqd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gKyBxcm8sIGRhdGEgPSBicm9vZF9lZ2csIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoZWdnLmludCkgI292ZXJkaXNwZXJzZWQNCg0KZWdnLm0xIDwtICBnbG0oZWdncyB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIGR1cmF0aW9uICsgcXJvLCBkYXRhID0gYnJvb2RfZWdnLCBmYW1pbHkgPSAicG9pc3NvbiIpDQpzdW1tYXJ5KGVnZy5tMSkgI292ZXJkaXNwZXJzZWQNCg0KZWdnLmludDIgPC0gZ2xtLm5iKGVnZ3MgfiB0cmVhdG1lbnQqd2hvbGUubWVhbiArIHFybyArIGR1cmF0aW9uLCBkYXRhID0gYnJvb2RfZWdnKSAgICNpZiBJIGluY2x1ZGUgdGhlIG90aGVyIHZhcmlhYmxlcyBoZXJlIHRoZSBtb2RlbCB3b24ndCBjb252ZXJnZSANCnN1bW1hcnkoZWdnLmludDIpICAgI25vdCBvdmVyZGlzcGVyc2VkIGFueW1vcmUNCg0KZWdnLm0yIDwtIGdsbS5uYihlZ2dzIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIHFybyArIGR1cmF0aW9uLCBkYXRhID0gYnJvb2RfZWdnKSAgICNzdGljayB3aXRoIHRoaXMgbW9kZWwgDQoNCg0KYW5vdmEoZWdnLmludDIsIGVnZy5tMikNCkFJQyhlZ2cuaW50MiwgZWdnLm0yKQ0KDQpgYGANCg0KDQojIyMjIFN0ZXAgMi5hIC0gRHJvcCAxDQoNCmBgYHtyfQ0KZHJvcDEoZWdnLm0yLCB0ZXN0ID0gIkNoaXNxIikNCg0KZWdnLm0yLnVwZGF0ZSA8LSB1cGRhdGUoZWdnLm0yLCAufi4gLWR1cmF0aW9uKQ0KYW5vdmEoZWdnLm0yLCBlZ2cubTIudXBkYXRlKQ0KDQpkcm9wMShlZ2cubTIudXBkYXRlLCB0ZXN0ID0gIkNoaXNxIikNCg0KYGBgDQoNCg0KIyMjIyBTdGVwIDMgLSBDaGVjayBtb2RlbCBmaXQgDQoNCmBgYHtyfQ0KcXFub3JtKHJlc2lkKGVnZy5tMi51cGRhdGUpKTtxcWxpbmUocmVzaWQoZWdnLm0yLnVwZGF0ZSkpDQpgYGANCg0KYGBge3J9DQoNCnBsb3QocmVzaWQoZWdnLm0yLnVwZGF0ZSkpICsgDQogIGFibGluZShoPTAsIGNvbD0icmVkIiwgbHdkPTIpIA0KDQpgYGANCg0KDQojIyMjIFN0ZXAgNCAtIFJlc3VsdHMgDQoNCmBgYHtyfQ0KQW5vdmEoZWdnLm0yLnVwZGF0ZSkNCnN1bW1hcnkoZWdnLm0yLnVwZGF0ZSkNCmVnZy5lbSA8LSBlbW1lYW5zKGVnZy5tMi51cGRhdGUsICJ0cmVhdG1lbnQiKQ0KcGFpcnMoZWdnLmVtKQ0KYGBgDQoNCiMjIyMgU3RlcCA1IC0gVmlzdWFsaXplIA0KDQpgYGB7cn0NCnBsb3QoZWdnLmVtLCBjb21wYXJpc29ucyA9IFRSVUUpDQpgYGANCg0KDQpSZXN1bHRzIFN1bW1hcnkgDQoNCmBgYHtyfQ0KZWdnLm0yLnVwZGF0ZQ0KZWdnc19BIDwtIGFzLmRhdGEuZnJhbWUoQW5vdmEoZWdnLm0yLnVwZGF0ZSkpDQplZ2dzX0EgPC0gc2V0RFQoZWdnc19BKQ0KZWdnc19BDQoNCg0KZWdnX2VtbSA8LSBlbW1lYW5zKGVnZy5tMi51cGRhdGUsIHBhaXJ3aXNlIH4gdHJlYXRtZW50KQ0KZWdnX2NvbnRyYXN0cyA8LSBhcy5kYXRhLmZyYW1lKGVnZ19lbW0kY29udHJhc3RzKQ0KZWdnX2NvbnRyYXN0cyA8LSBzZXREVChlZ2dfY29udHJhc3RzKQ0KZWdnX2NvbnRyYXN0cw0KDQoNCmVnZ3NfbWVhbnMgPC0gYXMuZGF0YS5mcmFtZShlZ2dfZW1tJGVtbWVhbnMpDQplZ2dzX21lYW5zIDwtIHNldERUKGVnZ3NfbWVhbnMpDQplZ2dzX21lYW5zDQoNCmBgYA0KDQoNCg0KIyMjIEhvbmV5IFBvdHMNCg0KDQpgYGB7cn0NCg0KYnJvb2RfaHAgPC0gc3Vic2V0KGJyb29kLCBzZWxlY3QgPSBjKGhvbmV5X3BvdCwgdHJlYXRtZW50LCB3aG9sZS5tZWFuLCBhbGl2ZSwgZHVyYXRpb24sIHFybywgY29sb255KSkNCg0KZ2dwbG90KGJyb29kX2hwLCBhZXMoeD1ob25leV9wb3QsIGZpbGwgPSB0cmVhdG1lbnQpKSArDQogIGdlb21faGlzdG9ncmFtKHBvc2l0aW9uID0gImlkZW50aXR5IiwgYmlud2lkdGggPSAwLjUsIGNvbD1JKCJibGFjayIpKSArDQogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGMoInBlYWNocHVmZjMiLCAiZGFya3NlYWdyZWVuIiwgImxpZ2h0c2VhZ3JlZW4iLCAiZGFya29saXZlZ3JlZW4zIiwgImRhcmtvbGl2ZWdyZWVuNCIpLA0KICAgICAgICAgICAgICAgICAgICBuYW1lID0gIlByaXN0aW5lIExldmVsIiwNCiAgICAgICAgICAgICAgICAgICAgbGFiZWxzID0gYygiVHJlYXRtZW50IDEgKGNvbnRyb2wpIiwgIlRyZWF0bWVudCAyIiwgDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIlRyZWF0bWVudCAzIiwgIlRyZWF0bWVudCA0IiwgIlRyZWF0bWVudCA1IikpICsNCiAgZ2d0aXRsZSgiSG9uZXkgUG90cyBwZXIgVHJlYXRtZW50IikgKw0KICBsYWJzKHkgPSAiQ291bnQiLCB4ID0gIkhQIENvdW50IikNCg0KYGBgDQoNCg0KIyMjIyBTdGVwIDEgLSBDaGVjayBmb3IgY28tbGluZWFyaXR5DQoNCmBgYHtyfQ0KaHAubW9kZWwgPC0gbG0oaG9uZXlfcG90IH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gICsgcXJvLCBkYXRhID0gYnJvb2RfaHApDQp2aWYoaHAubW9kZWwpDQpgYGANCg0KIyMjIyBTdGVwIDIgLSBXb3JrIGRvd24gZnJvbSBtb3N0IGNvbXBsaWNhdGVkIG1vZGVsIA0KDQpgYGB7cn0NCmhwLmludCA8LSBnbG0oaG9uZXlfcG90IH4gdHJlYXRtZW50Kndob2xlLm1lYW4gKyBhbGl2ZSArIGR1cmF0aW9uICsgcXJvLCBkYXRhID0gYnJvb2RfaHAsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoaHAuaW50KSAgI09EDQoNCg0KaHAubTEgPC0gZ2xtKGhvbmV5X3BvdCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIGR1cmF0aW9uICsgcXJvLCBkYXRhID0gYnJvb2RfaHAsIGZhbWlseSA9ICJwb2lzc29uIikNCnN1bW1hcnkoaHAubTEpICAgI09EDQoNCmhwLmludDIgPC0gZ2xtLm5iKGhvbmV5X3BvdCB+IHRyZWF0bWVudCp3aG9sZS5tZWFuICsgYWxpdmUgKyBkdXJhdGlvbiArIHFybywgZGF0YSA9IGJyb29kX2hwKQ0Kc3VtbWFyeShocC5pbnQyKSAgICNub3QgT0QgDQoNCmhwLm0yIDwtIGdsbS5uYihob25leV9wb3QgfiB0cmVhdG1lbnQgKyB3aG9sZS5tZWFuICsgYWxpdmUgKyBkdXJhdGlvbiArIHFybywgZGF0YSA9IGJyb29kX2hwKQ0Kc3VtbWFyeShocC5tMikgICAja2VlcCB0aGlzIG1vZGVsIA0KDQphbm92YShocC5pbnQyLCBocC5tMikNCg0KYGBgDQoNCg0KIyMjIyBTdGVwIDIuYSAtIERyb3AgMQ0KDQpgYGB7cn0NCmRyb3AxKGhwLm0yLCB0ZXN0ID0gIkNoaXNxIikNCg0KaHAubTMgPC0gZ2xtLm5iKGhvbmV5X3BvdCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIHFybywgZGF0YSA9IGJyb29kX2hwKQ0KDQpocC5tNCA8LSBnbG0ubmIoaG9uZXlfcG90IH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24sIGRhdGEgPSBicm9vZF9ocCkNCg0KaHAubTUgPC0gZ2xtLm5iKGhvbmV5X3BvdCB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSwgZGF0YSA9IGJyb29kX2hwKQ0KDQphbm92YShocC5tMiwgaHAubTMpDQphbm92YShocC5tMywgaHAubTQpDQoNCmRyb3AxKGhwLm0zLCB0ZXN0ID0gIkNoaXNxIikNCg0KYW5vdmEoaHAubTQsIGhwLm01KQ0KDQpBSUMoaHAubTIsIGhwLm0zLCBocC5tNCwgaHAubTUpDQoNCmBgYA0KDQoNCiMjIyMgU3RlcCAzIC0gQ2hlY2sgbW9kZWwgZml0IA0KDQpgYGB7cn0NCnFxbm9ybShyZXNpZChocC5tNSkpO3FxbGluZShyZXNpZChocC5tNSkpDQoNCnBsb3QocmVzaWQoaHAubTUpKSArIA0KICBhYmxpbmUoaD0wLCBjb2w9InJlZCIsIGx3ZD0yKSANCmBgYA0KDQoNCg0KIyMjIyBTdGVwIDQgLSBSZXN1bHRzIA0KDQpgYGB7cn0NCkFub3ZhKGhwLm01KQ0Kc3VtbWFyeShocC5tNSkNCg0KaHAuZW0gPC0gZW1tZWFucyhocC5tNSwgInRyZWF0bWVudCIpDQpwYWlycyhocC5lbSkNCmBgYA0KDQoNCg0KIyMjIyBTdGVwIDUgLSBWaXN1YWxpemUgDQoNCmBgYHtyfQ0KcGxvdChocC5lbSwgY29tcGFyaXNvbnMgPSBUUlVFKQ0KDQpgYGANCg0KDQoNCiAjIyMjIFJlc3VsdHMgU3VtbWFyeSANCg0KDQpgYGB7cn0NCmhwLm01DQoNCmhwX0EgPC0gYXMuZGF0YS5mcmFtZShBbm92YShocC5tNSkpDQpocF9BIDwtIHNldERUKGhwX0EpDQpocF9BDQoNCg0KaHBfZW1tIDwtIGVtbWVhbnMoaHAubTUsIHBhaXJ3aXNlIH4gdHJlYXRtZW50KQ0KaHBfY29udHJhc3RzIDwtIGFzLmRhdGEuZnJhbWUoaHBfZW1tJGNvbnRyYXN0cykNCmhwX2NvbnRyYXN0cyA8LSBzZXREVChocF9jb250cmFzdHMpDQpocF9jb250cmFzdHMNCg0KaHBfbWVhbnMgPC0gYXMuZGF0YS5mcmFtZShocF9lbW0kZW1tZWFucykNCmhwX21lYW5zIDwtIHNldERUKGhwX21lYW5zKQ0KaHBfbWVhbnMNCg0KYGBgDQoNCg0KDQojIyMjIFByb3BvcnRpb24gb2YgTGl2aW5nIExhcnZhZQ0KDQpgYGB7cn0NCmJyb29kX2wgPC0gc3Vic2V0KGJyb29kLCBzZWxlY3QgPSBjKGRlYWRfbGFydmFlLCBsaXZlX2xhcnZhZSwgdHJlYXRtZW50LCB3aG9sZS5tZWFuLCBhbGl2ZSwgZHVyYXRpb24sIHFybywgY29sb255KSkNCmBgYA0KDQoNCiMjIyMgU3RlcCAxIC0gQ2hlY2sgZm9yIGNvLWxpbmVhcml0eQ0KDQpgYGB7cn0NCmRsLm1vZGVsIDwtIGxtKGRlYWRfbGFydmFlIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gICsgcXJvLCBkYXRhID0gYnJvb2RfbCkNCnZpZihkbC5tb2RlbCkgICAjd2h5IGFtIEkganVzdCBub3cgcmVhbGl6aW5nIHRoaXMgaXMgYWx3YXlzIGdvaW5nIHRvIGJlIHNhbWUgYW5zd2VyIGZvciBhbGwgcmVzcG9uc2UgdmFyaWFibGVzIA0KYGBgDQoNCg0KIyMjIyBTdGVwIDIgLSBNb2RlbCANCg0KYGBge3J9DQpsYXJ2YWVfYmluZCA8LSBnbG0oY2JpbmQobGl2ZV9sYXJ2YWUsIGRlYWRfbGFydmFlKSB+IHRyZWF0bWVudCArIHdob2xlLm1lYW4gKyBhbGl2ZSArIGR1cmF0aW9uICsgcXJvLCBkYXRhID0gYnJvb2RfbCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQpzdW1tYXJ5KGxhcnZhZV9iaW5kKQ0KDQpsYXJ2YWVfYmluZC5pbnQgPC0gZ2xtKGNiaW5kKGxpdmVfbGFydmFlLCBkZWFkX2xhcnZhZSkgfiB0cmVhdG1lbnQqd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gKyBxcm8sIGRhdGEgPSBicm9vZF9sLCBmYW1pbHkgPSBiaW5vbWlhbCgibG9naXQiKSkNCnN1bW1hcnkobGFydmFlX2JpbmQuaW50KQ0KDQpBSUMobGFydmFlX2JpbmQsIGxhcnZhZV9iaW5kLmludCkNCg0KQW5vdmEobGFydmFlX2JpbmQpDQoNCmRyb3AxKGxhcnZhZV9iaW5kLCB0ZXN0ID0gIkNoaXNxIikNCg0KbGFydmFlX2JpbmQxIDwtIGdsbShjYmluZChsaXZlX2xhcnZhZSwgZGVhZF9sYXJ2YWUpIH4gdHJlYXRtZW50ICsgYWxpdmUgKyBxcm8sIGRhdGEgPSBicm9vZF9sLCBmYW1pbHkgPSBiaW5vbWlhbCgibG9naXQiKSkNCg0KYW5vdmEobGFydmFlX2JpbmQsIGxhcnZhZV9iaW5kMSkNCg0KQUlDKGxhcnZhZV9iaW5kLCBsYXJ2YWVfYmluZDEpDQoNCnN1bW1hcnkobGFydmFlX2JpbmQxKQ0KDQpxcW5vcm0ocmVzaWQobGFydmFlX2JpbmQxKSk7cXFsaW5lKHJlc2lkKGxhcnZhZV9iaW5kMSkpICAjTm90IGJhZCEgDQoNCnBsb3QocmVzaWQobGFydmFlX2JpbmQxKSkgKyANCiAgYWJsaW5lKGg9MCwgY29sPSJyZWQiLCBsd2Q9MikgDQoNCmxtLmVtbSA8LSBlbW1lYW5zKGxhcnZhZV9iaW5kMSwgcGFpcndpc2UgfiAgdHJlYXRtZW50LCB0eXBlID0gInJlc3BvbnNlIikNCnBhaXJzKGxtLmVtbSkNCmBgYA0KDQoNCiMjIyMgU3RlcCAzIC0gVmlzdWFsaXplIGFuZCBSZXN1bHRzIA0KDQpgYGB7cn0NCg0KbG0uZW1tIDwtIGVtbWVhbnMobGFydmFlX2JpbmQxLCBwYWlyd2lzZSB+ICB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KcGFpcnMobG0uZW1tKQ0Kc3VtbWFyeShsbS5lbW0pDQpBbm92YShsYXJ2YWVfYmluZDEpDQoNCmNsZGwgPC0gY2xkKG9iamVjdCA9IGxtLmVtbSwNCiAgICAgICAgICAgICBhZGp1c3QgPSAiVHVrZXkiLCANCiAgICAgICAgICAgICBhbHBoYSA9IDAuMDUsIA0KICAgICAgICAgICAgIExldHRlcnMgPSBsZXR0ZXJzKQ0KDQpjbGRsIDwtIGFzLmRhdGEuZnJhbWUoY2xkbCkNCmNsZGwNCg0KDQplbW0zIDwtIGVtbWVhbnMobGFydmFlX2JpbmQxICwgfiB0cmVhdG1lbnQsIHR5cGUgPSAicmVzcG9uc2UiKQ0KZW1tMw0KcGFpcnMoZW1tMykNCg0KZW1tZGYzIDwtIGFzLmRhdGEuZnJhbWUoZW1tMykNCg0KcCA8LSBnZ3Bsb3QoZGF0YSA9IGVtbWRmMywgYWVzKHg9dHJlYXRtZW50LCB5PXByb2IsIGZpbGw9dHJlYXRtZW50KSkgKyANCiAgZ2VvbV9jb2wocG9zaXRpb24gPSAiZG9kZ2UiLCBjb2xvciA9ICJibGFjayIpICsNCiAgY29vcmRfY2FydGVzaWFuKHlsaW0gPSBjKDAuNSwxKSkNCg0KcCArIGdlb21fZXJyb3JiYXIocG9zaXRpb249cG9zaXRpb25fZG9kZ2UoLjkpLHdpZHRoPS4yNSwgYWVzKHltYXg9YXN5bXAuVUNMLCB5bWluPWFzeW1wLkxDTCksYWxwaGE9LjYpKw0KICBsYWJzKHg9IlRyZWF0bWVudCIsIHk9IlByb2JhYmlsaXR5IG9mIFN1cnZpdmFsIiwgZmlsbCA9ICJ0cmVhdG1lbnQiKSArIHRoZW1lX2J3KCkgKw0KICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjKCJwZWFjaHB1ZmYzIiwgImRhcmtzZWFncmVlbiIsICJsaWdodHNlYWdyZWVuIiwgImRhcmtvbGl2ZWdyZWVuMyIsICJkYXJrb2xpdmVncmVlbjQiKSkgKw0KICAgbGFicyh5ID0gIlByb2JhYmlsaXR5IG9mIFN1cnZpdmFsIiwpICsNCiAgZ2d0aXRsZSgiTGFydmFlIHByb2JhYmlsaXR5IG9mIGJlaW5nIGFsaXZlIHVwb24gY29sb255IGRpc3NlY3Rpb24iKSArDQogIHNjYWxlX3hfZGlzY3JldGUobmFtZSA9ICJUcmVhdG1lbnQiLCANCiAgICAgICAgICAgICAgICAgICBsYWJlbHMgPSBjKCIwIFBQQiIsICIxNTAgUFBCIiwgIjEsNTAwIFBQQiIsICIxNSwwMDAgUFBCIiwgIjE1MCwwMDAgUFBCIikpICsNCiAgdGhlbWVfY293cGxvdCgpICsNCiAgdGhlbWVfY2xhc3NpYyhiYXNlX3NpemUgPSAxMikrDQogIGFubm90YXRlKGdlb20gPSAidGV4dCIsDQogICAgICAgICAgIHggPSAyLCB5ID0gMSwgDQogICAgICAgICAgIGxhYmVsID0gIlAgPCAwLjA1IiwgDQogICAgICAgICAgIHNpemUgPSA1KSArDQogIGFubm90YXRlKGdlb20gPSAidGV4dCIsIA0KICAgICAgICAgICB4ID0gYyg0LCAyLCAzLCAxLCA1KSwNCiAgICAgICAgICAgeSA9IGMoMC44NywgMC44OSwgMC45MiwgMC45OCwgMC45OSksIA0KICAgICAgICAgICBsYWJlbCA9IGMoImEiLCAiYWIiLCAiYmMiLCAiYyIsICJjIiksIA0KICAgICAgICAgICBzaXplID0gNSkgKw0KICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpDQoNCmBgYA0KDQoNCg0KIyMjIyBSZXN1bHRzIFN1bW1hcnkgDQoNCmBgYHtyfQ0KDQpsYXJ2YWVfYmluZDENCg0KbGFyX0EgPC0gYXMuZGF0YS5mcmFtZShBbm92YShsYXJ2YWVfYmluZDEpKQ0KbGFyX0EgPC0gc2V0RFQobGFyX0EpDQpsYXJfQQ0KDQoNCmxhcl9lbW0gPC0gZW1tZWFucyhsYXJ2YWVfYmluZDEsIHBhaXJ3aXNlIH4gdHJlYXRtZW50KQ0KbGFyX2NvbnRyYXN0cyA8LSBhcy5kYXRhLmZyYW1lKGxhcl9lbW0kY29udHJhc3RzKQ0KbGFyX2NvbnRyYXN0cyA8LSBzZXREVChsYXJfY29udHJhc3RzKQ0KbGFyX21lYW5zIDwtIGFzLmRhdGEuZnJhbWUobGFyX2VtbSRlbW1lYW5zKQ0KbGFyX21lYW5zIDwtIHNldERUKGxhcl9tZWFucykNCmxhcl9tZWFucw0KbGFyX2NvbnRyYXN0cw0KYGBgDQoNCg0KDQojIyMgUHJvcG9ydGlvbiBvZiBMaXZpbmcgUHVwYWUgDQoNCmBgYHtyfQ0KYnJvb2RfcCA8LSBzdWJzZXQoYnJvb2QsIHNlbGVjdCA9IGMoZGVhZF9wdXBhZSwgbGl2ZV9wdXBhZSwgdHJlYXRtZW50LCB3aG9sZS5tZWFuLCBhbGl2ZSwgZHVyYXRpb24sIHFybywgY29sb255KSkNCmBgYA0KDQoNCiMjIyMgU3RlcCAxIC0gTW9kZWwNCg0KYGBge3J9DQoNCnB1cGFlX2JpbmQgPC0gZ2xtKGNiaW5kKGxpdmVfcHVwYWUsIGRlYWRfcHVwYWUpIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIGFsaXZlICsgZHVyYXRpb24gKyBxcm8sIGRhdGEgPSBicm9vZF9wLCBmYW1pbHkgPSBiaW5vbWlhbCgibG9naXQiKSkNCnN1bW1hcnkocHVwYWVfYmluZCkNCg0KcHVwYWVfYmluZC5pbnQgPC0gZ2xtKGNiaW5kKGxpdmVfcHVwYWUsIGRlYWRfcHVwYWUpIH4gdHJlYXRtZW50Kndob2xlLm1lYW4gKyBhbGl2ZSArIGR1cmF0aW9uICsgcXJvLCBkYXRhID0gYnJvb2RfcCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQpzdW1tYXJ5KHB1cGFlX2JpbmQuaW50KQ0KDQpBSUMocHVwYWVfYmluZCwgcHVwYWVfYmluZC5pbnQpDQoNCkFub3ZhKHB1cGFlX2JpbmQpDQoNCmFub3ZhKHB1cGFlX2JpbmQsIHB1cGFlX2JpbmQuaW50KQ0KDQpkcm9wMShwdXBhZV9iaW5kLmludCwgdGVzdCA9ICJDaGlzcSIpDQoNCnB1cGFlX2ludDEgPC0gZ2xtKGNiaW5kKGxpdmVfcHVwYWUsIGRlYWRfcHVwYWUpIH4gdHJlYXRtZW50Kndob2xlLm1lYW4gKyBhbGl2ZSwgZGF0YSA9IGJyb29kX3AsIGZhbWlseSA9IGJpbm9taWFsKCJsb2dpdCIpKQ0KDQphbm92YShwdXBhZV9iaW5kLmludCwgcHVwYWVfaW50MSkNCg0KQUlDKHB1cGFlX2JpbmQuaW50LCBwdXBhZV9pbnQxKQ0KDQpzdW1tYXJ5KHB1cGFlX2ludDEpDQoNCnFxbm9ybShyZXNpZChwdXBhZV9pbnQxKSk7cXFsaW5lKHJlc2lkKHB1cGFlX2ludDEpKSAgI05vdCBiYWQhIA0KDQpwbG90KHJlc2lkKHB1cGFlX2ludDEpKSArIA0KICBhYmxpbmUoaD0wLCBjb2w9InJlZCIsIGx3ZD0yKSANCg0KcG0uZW1tIDwtIGVtbWVhbnMocHVwYWVfaW50MSwgcGFpcndpc2UgfiAgdHJlYXRtZW50IHwgd2hvbGUubWVhbikNCnBhaXJzKHBtLmVtbSkNCg0KYGBgDQoNCiMjIyMgUmVzdWx0cyBTdW1tYXJ5IA0KDQpgYGB7cn0NCg0KcHVwYWVfaW50MQ0KQW5vdmEocHVwYWVfaW50MSkNCg0KcHVwX0EgPC0gYXMuZGF0YS5mcmFtZShBbm92YShwdXBhZV9pbnQxKSkNCnB1cF9BIDwtIHNldERUKHB1cF9BKQ0KcHVwX0ENCg0KDQpwdXBfZW1tIDwtIGVtbWVhbnMocHVwYWVfaW50MSwgcGFpcndpc2UgfiB0cmVhdG1lbnQpDQpwdXBfY29udHJhc3RzIDwtIGFzLmRhdGEuZnJhbWUocHVwX2VtbSRjb250cmFzdHMpDQpwdXBfY29udHJhc3RzIDwtIHNldERUKHB1cF9jb250cmFzdHMpDQpwdXBfbWVhbnMgPC0gYXMuZGF0YS5mcmFtZShwdXBfZW1tJGVtbWVhbnMpDQpwdXBfbWVhbnMgPC0gc2V0RFQocHVwX21lYW5zKQ0KcHVwX21lYW5zDQpwdXBfY29udHJhc3RzDQoNCmBgYA0KDQoNCiMjIyBQcm9wb3J0aW9uIG9mIFN1cnZpdmluZyBMYXJ2YWUgYW5kIFB1cGFlIA0KDQpgYGB7cn0NCg0KcGxfYmluZCA8LSBnbG0oY2JpbmQoYWxpdmVfbHAsIGRlYWRfbHApIH4gdHJlYXRtZW50ICsgd2hvbGUubWVhbiArIHFybyArIGR1cmF0aW9uICsgYWxpdmUsIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQpzdW1tYXJ5KHBsX2JpbmQgKQ0KDQpwbF9iaW5kLmludCA8LSBnbG0oY2JpbmQoYWxpdmVfbHAsIGRlYWRfbHApIH4gdHJlYXRtZW50Kndob2xlLm1lYW4gKyBxcm8gKyBkdXJhdGlvbiArIGFsaXZlLCBkYXRhID0gYnJvb2QsIGZhbWlseSA9IGJpbm9taWFsKCJsb2dpdCIpKQ0Kc3VtbWFyeShwbF9iaW5kLmludCkNCg0KYW5vdmEocGxfYmluZC5pbnQsIHBsX2JpbmQpDQpBSUMocGxfYmluZCwgcGxfYmluZC5pbnQpDQoNCmRyb3AxKHBsX2JpbmQuaW50LCB0ZXN0ID0gIkNoaXNxIikNCg0KcGxfdXBkYXRlIDwtIGdsbShjYmluZChhbGl2ZV9scCwgZGVhZF9scCkgfiB0cmVhdG1lbnQqd2hvbGUubWVhbiArIGR1cmF0aW9uICsgYWxpdmUsIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQoNCnBsX3VwZGF0ZTIgPC0gZ2xtKGNiaW5kKGFsaXZlX2xwLCBkZWFkX2xwKSB+IHRyZWF0bWVudCp3aG9sZS5tZWFuICsgYWxpdmUsIGRhdGEgPSBicm9vZCwgZmFtaWx5ID0gYmlub21pYWwoImxvZ2l0IikpDQoNCmFub3ZhKHBsX2JpbmQuaW50LCBwbF91cGRhdGUsIHBsX3VwZGF0ZTIpDQoNCkFJQyhwbF9iaW5kLmludCwgcGxfdXBkYXRlLCBwbF91cGRhdGUyKQ0KDQpgYGANCg0KDQoNCmBgYHtyfQ0KDQpxcW5vcm0ocmVzaWQocGxfdXBkYXRlKSk7cXFsaW5lKHJlc2lkKHBsX3VwZGF0ZSkpDQoNCmBgYA0KDQoNCg==