Workers No Round 1

1 Input Pollen Data

Start by imputing pollen data and creating a new data frame with average pollen consumption on a per-colony basis

### Figure out average pollen consumption by treatment 


pollen <- read_csv("pollen1.csv", col_types = cols(round = col_factor(levels = c("1", 
                                                                                 "2")), treatment = col_factor(levels = c("1", 
                                                                                                                          "2", "3", "4", "5")), replicate = col_factor(levels = c("1", 
                                                                                                                                                                                  "2", "3", "4", "5", "6", "7", "9", "11", 
                                                                                                                                                                                  "12")), start_date = col_date(format = "%m/%d/%Y"), 
                                                   start_time = col_character(), end_date = col_date(format = "%m/%d/%Y"), 
                                                   end_time = col_character()))


pollen$colony <- as.factor(pollen$colony)
pollen$pid <- as.factor(pollen$pid)
pollen$count <- as.factor(pollen$count)

pollen <- na.omit(pollen)

range(pollen$difference)

## [1] -0.98780  1.56542

pollen <- subset(pollen, pollen$round != 1)
pollen

## # A tibble: 928 × 15
##    round colony treatment replicate count start_date start_…¹ start…² end_date  
##    <fct> <fct>  <fct>     <fct>     <fct> <date>     <chr>      <dbl> <date>    
##  1 2     1.11R2 1         11        2     2022-08-22 10:30      1.18  2022-08-24
##  2 2     1.11R2 1         11        3     2022-08-24 12:00      1.15  2022-08-26
##  3 2     1.11R2 1         11        4     2022-08-26 10:30      1.09  2022-08-28
##  4 2     1.11R2 1         11        5     2022-08-28 1:00       1.26  2022-08-30
##  5 2     1.11R2 1         11        5     2022-08-30 11:30      1.14  2022-09-01
##  6 2     1.11R2 1         11        21    2022-09-01 10:30      1.07  2022-09-03
##  7 2     1.11R2 1         11        8     2022-09-03 12:20      0.844 2022-09-05
##  8 2     1.11R2 1         11        9     2022-09-05 12:40      1.30  2022-09-07
##  9 2     1.11R2 1         11        10    2022-09-07 10:30      1.22  2022-09-09
## 10 2     1.11R2 1         11        11    2022-09-09 11:30      1.26  2022-09-11
## # … with 918 more rows, 6 more variables: end_time <chr>, end_weight <dbl>,
## #   whole_dif <chr>, difference <dbl>, pid <fct>, bees_alive <dbl>, and
## #   abbreviated variable names ¹​start_time, ²​start_weight

# get rid of negative numbers
pollen$difference[pollen$difference < 0] <- NA
pollen <- na.omit(pollen)
range(pollen$difference)

## [1] 0.002715 1.565420

1.1 Average pollen consumption per colony

pollen$whole_dif <- as.numeric(pollen$difference)
wd.1 <- lm(difference ~ treatment, data = pollen)
summary(wd.1)

## 
## Call:
## lm(formula = difference ~ treatment, data = pollen)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4853 -0.2416 -0.1504  0.1576  1.0638 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.429541   0.024170  17.772   <2e-16 ***
## treatment2  0.072099   0.034886   2.067   0.0390 *  
## treatment3  0.078078   0.034406   2.269   0.0235 *  
## treatment4  0.058490   0.034988   1.672   0.0949 .  
## treatment5  0.004982   0.035040   0.142   0.8870    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3375 on 915 degrees of freedom
## Multiple R-squared:  0.009927,   Adjusted R-squared:  0.005599 
## F-statistic: 2.294 on 4 and 915 DF,  p-value: 0.05773

wd.emm <- emmeans(wd.1, "treatment")
summary(wd.emm)

##  treatment emmean     SE  df lower.CL upper.CL
##  1          0.430 0.0242 915    0.382    0.477
##  2          0.502 0.0252 915    0.452    0.551
##  3          0.508 0.0245 915    0.460    0.556
##  4          0.488 0.0253 915    0.438    0.538
##  5          0.435 0.0254 915    0.385    0.484
## 
## Confidence level used: 0.95

wd.summary <- pollen %>% 
  group_by(colony) %>%
  summarize(whole.mean = mean(difference))

as.data.frame(wd.summary)  # this is the data frame I will merge with subsequent data frames to contain average pollen consumption per colony as a new column  

##    colony whole.mean
## 1  1.11R2  0.2829509
## 2  1.12R2  0.1697964
## 3   1.1R2  0.5213458
## 4   1.2R2  0.3374200
## 5   1.3R2  0.4512891
## 6   1.4R2  0.6063016
## 7   1.5R2  0.7079516
## 8   1.7R2  0.7400381
## 9   1.9R2  0.2240081
## 10 2.11R2  0.4178270
## 11 2.12R2  0.4035568
## 12  2.1R2  0.6101589
## 13  2.2R2  0.5109234
## 14  2.3R2  0.3280036
## 15  2.4R2  0.3881394
## 16  2.5R2  0.7194915
## 17  2.7R2  0.5299685
## 18  2.9R2  0.5857152
## 19 3.11R2  0.4216410
## 20 3.12R2  0.3390993
## 21  3.1R2  0.3711948
## 22  3.2R2  0.3461010
## 23  3.3R2  0.8465806
## 24  3.4R2  0.4120433
## 25  3.5R2  0.8906211
## 26  3.7R2  0.6266680
## 27  3.9R2  0.4409511
## 28 4.11R2  0.3456011
## 29 4.12R2  0.2496295
## 30  4.1R2  0.7074755
## 31  4.2R2  0.3871742
## 32  4.3R2  0.5800074
## 33  4.4R2  0.8356247
## 34  4.5R2  0.2335967
## 35  4.7R2  0.6237400
## 36  4.9R2  0.5322950
## 37 5.11R2  0.4113960
## 38 5.12R2  0.2077741
## 39  5.1R2  0.3782286
## 40  5.2R2  0.4912468
## 41  5.3R2  0.2128778
## 42  5.4R2  0.4563045
## 43  5.5R2  0.7095479
## 44  5.7R2  0.6113189
## 45  5.9R2  0.4188290

2 Worker Data

workers <- read_csv("workers2.csv", col_types = cols(round = col_factor(levels = c("1", 
    "2")), treatment = col_factor(levels = c("1", 
    "2", "3", "4", "5")), replicate = col_factor(levels = c("1", 
    "2", "3", "4", "5", "7", "9", "11", "12")), 
    biobest_origin = col_factor(levels = c("B1", 
        "B2", "B3", "B4", "B5", "K1", "K2", 
        "K3")), start_date = col_skip(), mortality_date = col_skip(), 
    replaced = col_skip(), mortality = col_skip(), `replacement_wet weight` = col_skip(), 
    ending_wet_weight = col_skip()))

workers <- subset(workers, workers$round != 1)


workers$qro <- as.factor(workers$biobest_origin)
workers$dry <- as.double(workers$`dry weight`)
workers$colony <- as.factor(workers$colony)
workers$proportion <- as.double(workers$proportion)

wrfp <- merge(wd.summary, workers, by.x = "colony") 

2.1 mortality

mortality <- read_csv("surviving workers per colony.csv")

mortality$colony <- as.factor(mortality$colony)

mortality 

## # A tibble: 60 × 3
##    colony alive  dead
##    <fct>  <dbl> <dbl>
##  1 1.11R2     5     0
##  2 1.12R2     0     5
##  3 1.1R1      5     0
##  4 1.1R2      5     0
##  5 1.2R1      5     0
##  6 1.2R2      4     1
##  7 1.3R1      5     0
##  8 1.3R2      5     0
##  9 1.4R2      1     4
## 10 1.5R2      4     1
## # … with 50 more rows

wrpm <- merge(wrfp, mortality, by.x = "colony")

wrpm <- na.omit(wrpm)

Proportion of time each worker survived during the experiment per treatment, and average numbers of workers surviving whole experiment per treatment.

Workers in treatment 4 survived the shortest amount of time, but workers in treatment 1 died more on average.

prop.sum <- wrpm %>%
  group_by(treatment) %>%
  summarise(mp = mean(proportion))

prop.sum

## # A tibble: 5 × 2
##   treatment    mp
##   <fct>     <dbl>
## 1 1         0.938
## 2 2         0.978
## 3 3         0.979
## 4 4         0.920
## 5 5         0.988

survive.sum <- wrpm %>%
  group_by(treatment) %>%
  summarise(sm = mean(alive))

survive.sum

## # A tibble: 5 × 2
##   treatment    sm
##   <fct>     <dbl>
## 1 1          3.44
## 2 2          4.22
## 3 3          4.67
## 4 4          4   
## 5 5          4.33

2.2 Workers alive at end of experiment

Data is underdispersed

workers.gn <- glm(alive ~ treatment + whole.mean  + qro, data = wrpm, family = "poisson")
summary(workers.gn)

## 
## Call:
## glm(formula = alive ~ treatment + whole.mean + qro, family = "poisson", 
##     data = wrpm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.55436  -0.24920   0.08058   0.38115   0.93076  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.89623    0.12319   7.275 3.46e-13 ***
## treatment2   0.15526    0.10850   1.431 0.152435    
## treatment3   0.23019    0.10688   2.154 0.031267 *  
## treatment4   0.08540    0.11124   0.768 0.442674    
## treatment5   0.25260    0.10781   2.343 0.019135 *  
## whole.mean   1.03151    0.20305   5.080 3.77e-07 ***
## qroB3       -0.18595    0.10951  -1.698 0.089499 .  
## qroB4       -0.46936    0.12338  -3.804 0.000142 ***
## qroB5       -0.32122    0.08767  -3.664 0.000248 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 195.08  on 223  degrees of freedom
## Residual deviance: 145.73  on 215  degrees of freedom
## AIC: 857.4
## 
## Number of Fisher Scoring iterations: 5

workers.gnb <- glm.nb(alive ~ treatment + whole.mean  + qro, data = wrpm)
summary(workers.gnb)

## 
## Call:
## glm.nb(formula = alive ~ treatment + whole.mean + qro, data = wrpm, 
##     init.theta = 145404.8379, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.55434  -0.24920   0.08058   0.38114   0.93075  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.89623    0.12319   7.275 3.46e-13 ***
## treatment2   0.15526    0.10850   1.431 0.152442    
## treatment3   0.23019    0.10688   2.154 0.031269 *  
## treatment4   0.08540    0.11124   0.768 0.442683    
## treatment5   0.25260    0.10782   2.343 0.019137 *  
## whole.mean   1.03151    0.20305   5.080 3.77e-07 ***
## qroB3       -0.18595    0.10951  -1.698 0.089504 .  
## qroB4       -0.46936    0.12339  -3.804 0.000142 ***
## qroB5       -0.32122    0.08768  -3.664 0.000249 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(145404.8) family taken to be 1)
## 
##     Null deviance: 195.08  on 223  degrees of freedom
## Residual deviance: 145.73  on 215  degrees of freedom
## AIC: 859.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  145405 
##           Std. Err.:  1233053 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -839.402

AIC(workers.gn, workers.gnb)

##             df      AIC
## workers.gn   9 857.3976
## workers.gnb 10 859.4017

Anova(workers.gnb)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment    7.7101  4     0.1028    
## whole.mean  25.9957  1  3.422e-07 ***
## qro         24.7155  3  1.771e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(workers.gn)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment    7.7103  4     0.1028    
## whole.mean  25.9963  1  3.421e-07 ***
## qro         24.7161  3  1.770e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

workers.gn1 <- glm(alive ~ treatment + whole.mean  + qro, data = wrpm, family = "poisson")
summary(workers.gn1)

## 
## Call:
## glm(formula = alive ~ treatment + whole.mean + qro, family = "poisson", 
##     data = wrpm)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.55436  -0.24920   0.08058   0.38115   0.93076  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.89623    0.12319   7.275 3.46e-13 ***
## treatment2   0.15526    0.10850   1.431 0.152435    
## treatment3   0.23019    0.10688   2.154 0.031267 *  
## treatment4   0.08540    0.11124   0.768 0.442674    
## treatment5   0.25260    0.10781   2.343 0.019135 *  
## whole.mean   1.03151    0.20305   5.080 3.77e-07 ***
## qroB3       -0.18595    0.10951  -1.698 0.089499 .  
## qroB4       -0.46936    0.12338  -3.804 0.000142 ***
## qroB5       -0.32122    0.08767  -3.664 0.000248 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 195.08  on 223  degrees of freedom
## Residual deviance: 145.73  on 215  degrees of freedom
## AIC: 857.4
## 
## Number of Fisher Scoring iterations: 5

workers.gnb1 <- glm.nb(alive ~ treatment + whole.mean  + qro, data = wrpm)
summary(workers.gnb1)

## 
## Call:
## glm.nb(formula = alive ~ treatment + whole.mean + qro, data = wrpm, 
##     init.theta = 145404.8379, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.55434  -0.24920   0.08058   0.38114   0.93075  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.89623    0.12319   7.275 3.46e-13 ***
## treatment2   0.15526    0.10850   1.431 0.152442    
## treatment3   0.23019    0.10688   2.154 0.031269 *  
## treatment4   0.08540    0.11124   0.768 0.442683    
## treatment5   0.25260    0.10782   2.343 0.019137 *  
## whole.mean   1.03151    0.20305   5.080 3.77e-07 ***
## qroB3       -0.18595    0.10951  -1.698 0.089504 .  
## qroB4       -0.46936    0.12339  -3.804 0.000142 ***
## qroB5       -0.32122    0.08768  -3.664 0.000249 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(145404.8) family taken to be 1)
## 
##     Null deviance: 195.08  on 223  degrees of freedom
## Residual deviance: 145.73  on 215  degrees of freedom
## AIC: 859.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  145405 
##           Std. Err.:  1233053 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -839.402

AIC(workers.gn1, workers.gnb1)

##              df      AIC
## workers.gn1   9 857.3976
## workers.gnb1 10 859.4017

Anova(workers.gn1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment    7.7103  4     0.1028    
## whole.mean  25.9963  1  3.421e-07 ***
## qro         24.7161  3  1.770e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(workers.gnb1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment    7.7101  4     0.1028    
## whole.mean  25.9957  1  3.422e-07 ***
## qro         24.7155  3  1.771e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(workers.gn)

plot(workers.gn1)

plot(workers.gnb)

plot(workers.gnb1)

plot(wrpm$treatment, wrpm$alive)

alive.sum <- wrpm %>%
  group_by(treatment) %>%
  summarise(a.m= mean(alive), 
            a.n = length(alive), 
            sd.a = sd(alive),
            n.a = length(alive)) %>%
  mutate(sea = sd.a / sqrt(n.a))

alive.sum

## # A tibble: 5 × 6
##   treatment   a.m   a.n  sd.a   n.a   sea
##   <fct>     <dbl> <int> <dbl> <int> <dbl>
## 1 1          3.44    45 1.85     45 0.276
## 2 2          4.22    45 1.57     45 0.233
## 3 3          4.67    45 0.674    45 0.101
## 4 4          4       44 1.36     44 0.206
## 5 5          4.33    45 1.58     45 0.236

ggplot(data = alive.sum, aes(x = treatment, y = a.m, fill = treatment)) +
  geom_col(col = "black")+
  coord_cartesian(ylim=c(3,5)) +
  scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
                    name = "Pristine Level",
                    labels = c("Treatment 1 (control)", "Treatment 2", 
                               "Treatment 3", "Treatment 4", "Treatment 5")) + 
  geom_errorbar(aes(ymin = a.m - sea, 
                    ymax = a.m + sea),
                position = position_dodge(0.9), width = 0.4) +
  labs(y = "Workers Alive at End of Experiment",) +
  ggtitle("Average Workers Alive per Treatment at End of Experiment") +
  scale_x_discrete(name = "Treatment", 
                   labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
  theme_cowplot()+
  theme_classic(base_size = 20)

2.2.1 dry weights

hist(wrpm$dry)

range(wrpm$dry)

## [1] 0.01761 0.12047

shapiro.test(wrpm$dry)

## 
##  Shapiro-Wilk normality test
## 
## data:  wrpm$dry
## W = 0.906, p-value = 1.153e-10

wrpm$log.dry <- log(wrpm$dry)

hist(wrpm$log.dry)

shapiro.test(wrpm$log.dry)    #log actually worked for once heck yeah 

## 
##  Shapiro-Wilk normality test
## 
## data:  wrpm$log.dry
## W = 0.98708, p-value = 0.04033

dry.lmer <- lmer(log.dry ~ treatment + whole.mean  + qro + (1|colony), data = wrpm)
summary(dry.lmer)

## Linear mixed model fit by REML ['lmerMod']
## Formula: log.dry ~ treatment + whole.mean + qro + (1 | colony)
##    Data: wrpm
## 
## REML criterion at convergence: 142.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.79001 -0.67996  0.05816  0.59743  3.15400 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  colony   (Intercept) 0.00000  0.0000  
##  Residual             0.09849  0.3138  
## Number of obs: 224, groups:  colony, 45
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) -3.657569   0.072769 -50.262
## treatment2  -0.054855   0.066468  -0.825
## treatment3  -0.007965   0.066797  -0.119
## treatment4  -0.023149   0.066826  -0.346
## treatment5   0.047289   0.066194   0.714
## whole.mean   0.931825   0.126445   7.369
## qroB3        0.241213   0.068968   3.497
## qroB4        0.212541   0.073568   2.889
## qroB5        0.242613   0.052721   4.602
## 
## Correlation of Fixed Effects:
##            (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn qroB3  qroB4 
## treatment2 -0.384                                                 
## treatment3 -0.351  0.506                                          
## treatment4 -0.382  0.502  0.503                                   
## treatment5 -0.476  0.495  0.491  0.492                            
## whole.mean -0.721 -0.096 -0.138 -0.093  0.030                     
## qroB3      -0.108  0.007  0.009  0.002 -0.002 -0.069              
## qroB4       0.108  0.034  0.049  0.029 -0.011 -0.354  0.181       
## qroB5      -0.153  0.007  0.010  0.002 -0.002 -0.074  0.224  0.231
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Anova(dry.lmer)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: log.dry
##              Chisq Df Pr(>Chisq)    
## treatment   2.4796  4     0.6483    
## whole.mean 54.3085  1  1.714e-13 ***
## qro        29.9820  3  1.392e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(resid(dry.lmer)) + 
  abline(h=0, col="red", lwd=2)

## integer(0)

qqnorm(resid(dry.lmer));qqline(resid(dry.lmer))

wdsum <- wrpm %>%
  group_by(treatment) %>%
  summarise(md = mean(dry),
            sdd = sd(dry), 
            ndry = length(dry))
wdsum

## # A tibble: 5 × 4
##   treatment     md    sdd  ndry
##   <fct>      <dbl>  <dbl> <int>
## 1 1         0.0474 0.0205    45
## 2 2         0.0456 0.0159    45
## 3 3         0.0499 0.0204    45
## 4 4         0.0487 0.0227    44
## 5 5         0.0478 0.0176    45

I’m going to go ahead and say treatment doesn’t influence worker dry weights

plot(wrpm$treatment, wrpm$dry)

3 cbind based on herbivory example

mod1 <- glm(cbind(alive, dead) ~ treatment  + start_wet_weight + qro, data = wrpm, family = binomial("logit"))
summary(mod1)

## 
## Call:
## glm(formula = cbind(alive, dead) ~ treatment + start_wet_weight + 
##     qro, family = binomial("logit"), data = wrpm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.0212   0.1706   0.7094   1.1142   2.7057  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.2649     0.3066   0.864  0.38766    
## treatment2         1.1322     0.2523   4.487 7.21e-06 ***
## treatment3         2.0535     0.3189   6.440 1.20e-10 ***
## treatment4         0.7596     0.2398   3.168  0.00153 ** 
## treatment5         1.3266     0.2614   5.075 3.87e-07 ***
## start_wet_weight   8.2150     1.6432   4.999 5.76e-07 ***
## qroB3             -1.1967     0.2840  -4.213 2.52e-05 ***
## qroB4             -2.0071     0.2824  -7.107 1.19e-12 ***
## qroB5             -2.0798     0.2259  -9.205  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 677.74  on 223  degrees of freedom
## Residual deviance: 515.01  on 215  degrees of freedom
## AIC: 654.06
## 
## Number of Fisher Scoring iterations: 5

anova(mod1)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cbind(alive, dead)
## 
## Terms added sequentially (first to last)
## 
## 
##                  Df Deviance Resid. Df Resid. Dev
## NULL                               223     677.74
## treatment         4   51.925       219     625.82
## start_wet_weight  1    2.897       218     622.92
## qro               3  107.907       215     515.01

Anova(mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cbind(alive, dead)
##                  LR Chisq Df Pr(>Chisq)    
## treatment          59.959  4  2.959e-12 ***
## start_wet_weight   27.744  1  1.385e-07 ***
## qro               107.907  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mod1) # this doesn't look great 

emmeans(mod1, pairwise ~ treatment, type = "response")

## $emmeans
##  treatment  prob     SE  df asymp.LCL asymp.UCL
##  1         0.599 0.0388 Inf     0.521     0.672
##  2         0.822 0.0289 Inf     0.759     0.872
##  3         0.921 0.0205 Inf     0.870     0.953
##  4         0.761 0.0336 Inf     0.690     0.821
##  5         0.849 0.0270 Inf     0.788     0.895
## 
## Results are averaged over the levels of: qro 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.322 0.0813 Inf    1  -4.487  0.0001
##  treatment1 / treatment3      0.128 0.0409 Inf    1  -6.440  <.0001
##  treatment1 / treatment4      0.468 0.1122 Inf    1  -3.168  0.0133
##  treatment1 / treatment5      0.265 0.0694 Inf    1  -5.075  <.0001
##  treatment2 / treatment3      0.398 0.1344 Inf    1  -2.729  0.0499
##  treatment2 / treatment4      1.451 0.3860 Inf    1   1.401  0.6270
##  treatment2 / treatment5      0.823 0.2331 Inf    1  -0.687  0.9594
##  treatment3 / treatment4      3.647 1.2032 Inf    1   3.921  0.0008
##  treatment3 / treatment5      2.069 0.7113 Inf    1   2.114  0.2141
##  treatment4 / treatment5      0.567 0.1555 Inf    1  -2.068  0.2339
## 
## Results are averaged over the levels of: qro 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

emm1 <- emmeans(mod1, ~ treatment, type = "response")
emm1

##  treatment  prob     SE  df asymp.LCL asymp.UCL
##  1         0.599 0.0388 Inf     0.521     0.672
##  2         0.822 0.0289 Inf     0.759     0.872
##  3         0.921 0.0205 Inf     0.870     0.953
##  4         0.761 0.0336 Inf     0.690     0.821
##  5         0.849 0.0270 Inf     0.788     0.895
## 
## Results are averaged over the levels of: qro 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale

pairs(emm1)

##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.322 0.0813 Inf    1  -4.487  0.0001
##  treatment1 / treatment3      0.128 0.0409 Inf    1  -6.440  <.0001
##  treatment1 / treatment4      0.468 0.1122 Inf    1  -3.168  0.0133
##  treatment1 / treatment5      0.265 0.0694 Inf    1  -5.075  <.0001
##  treatment2 / treatment3      0.398 0.1344 Inf    1  -2.729  0.0499
##  treatment2 / treatment4      1.451 0.3860 Inf    1   1.401  0.6270
##  treatment2 / treatment5      0.823 0.2331 Inf    1  -0.687  0.9594
##  treatment3 / treatment4      3.647 1.2032 Inf    1   3.921  0.0008
##  treatment3 / treatment5      2.069 0.7113 Inf    1   2.114  0.2141
##  treatment4 / treatment5      0.567 0.1555 Inf    1  -2.068  0.2339
## 
## Results are averaged over the levels of: qro 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

emmdf <- as.data.frame(emm1)

p <- ggplot(data = emmdf, aes(x=treatment, y=prob, fill=treatment)) + 
  geom_col(position = "dodge", color = "black") +
  coord_cartesian(ylim = c(0.5,1))

p + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=.6)+
  labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
  scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4")) +
   labs(y = "Probability of Survival",) +
  ggtitle("Probability of workers surviving the whole experiment") +
  scale_x_discrete(name = "Treatment", 
                   labels = c("0 PPB", "150 PPB", "1,500 PPB", "15,000 PPB", "150,000 PPB")) +
  theme_cowplot() +
  theme_classic(base_size = 12)

pr <- subset(wrpm, wrpm$round != 1)
pr

mod2 <- glm(cbind(alive, dead) ~ treatment + start_wet_weight + qro, data = pr, family = binomial("logit"))
summary(mod2)

## 
## Call:
## glm(formula = cbind(alive, dead) ~ treatment + start_wet_weight + 
##     qro, family = binomial("logit"), data = pr)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.0212   0.1706   0.7094   1.1142   2.7057  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.2649     0.3066   0.864  0.38766    
## treatment2         1.1322     0.2523   4.487 7.21e-06 ***
## treatment3         2.0535     0.3189   6.440 1.20e-10 ***
## treatment4         0.7596     0.2398   3.168  0.00153 ** 
## treatment5         1.3266     0.2614   5.075 3.87e-07 ***
## start_wet_weight   8.2150     1.6432   4.999 5.76e-07 ***
## qroB3             -1.1967     0.2840  -4.213 2.52e-05 ***
## qroB4             -2.0071     0.2824  -7.107 1.19e-12 ***
## qroB5             -2.0798     0.2259  -9.205  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 677.74  on 223  degrees of freedom
## Residual deviance: 515.01  on 215  degrees of freedom
## AIC: 654.06
## 
## Number of Fisher Scoring iterations: 5

anova(mod2)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: cbind(alive, dead)
## 
## Terms added sequentially (first to last)
## 
## 
##                  Df Deviance Resid. Df Resid. Dev
## NULL                               223     677.74
## treatment         4   51.925       219     625.82
## start_wet_weight  1    2.897       218     622.92
## qro               3  107.907       215     515.01

Anova(mod2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cbind(alive, dead)
##                  LR Chisq Df Pr(>Chisq)    
## treatment          59.959  4  2.959e-12 ***
## start_wet_weight   27.744  1  1.385e-07 ***
## qro               107.907  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mod2)

emmeans(mod1, pairwise ~ treatment, type = "response")

## $emmeans
##  treatment  prob     SE  df asymp.LCL asymp.UCL
##  1         0.599 0.0388 Inf     0.521     0.672
##  2         0.822 0.0289 Inf     0.759     0.872
##  3         0.921 0.0205 Inf     0.870     0.953
##  4         0.761 0.0336 Inf     0.690     0.821
##  5         0.849 0.0270 Inf     0.788     0.895
## 
## Results are averaged over the levels of: qro 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.322 0.0813 Inf    1  -4.487  0.0001
##  treatment1 / treatment3      0.128 0.0409 Inf    1  -6.440  <.0001
##  treatment1 / treatment4      0.468 0.1122 Inf    1  -3.168  0.0133
##  treatment1 / treatment5      0.265 0.0694 Inf    1  -5.075  <.0001
##  treatment2 / treatment3      0.398 0.1344 Inf    1  -2.729  0.0499
##  treatment2 / treatment4      1.451 0.3860 Inf    1   1.401  0.6270
##  treatment2 / treatment5      0.823 0.2331 Inf    1  -0.687  0.9594
##  treatment3 / treatment4      3.647 1.2032 Inf    1   3.921  0.0008
##  treatment3 / treatment5      2.069 0.7113 Inf    1   2.114  0.2141
##  treatment4 / treatment5      0.567 0.1555 Inf    1  -2.068  0.2339
## 
## Results are averaged over the levels of: qro 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

emm2 <- emmeans(mod2, ~ treatment, type = "response")
emm2

##  treatment  prob     SE  df asymp.LCL asymp.UCL
##  1         0.599 0.0388 Inf     0.521     0.672
##  2         0.822 0.0289 Inf     0.759     0.872
##  3         0.921 0.0205 Inf     0.870     0.953
##  4         0.761 0.0336 Inf     0.690     0.821
##  5         0.849 0.0270 Inf     0.788     0.895
## 
## Results are averaged over the levels of: qro 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale

pairs(emm2)

##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.322 0.0813 Inf    1  -4.487  0.0001
##  treatment1 / treatment3      0.128 0.0409 Inf    1  -6.440  <.0001
##  treatment1 / treatment4      0.468 0.1122 Inf    1  -3.168  0.0133
##  treatment1 / treatment5      0.265 0.0694 Inf    1  -5.075  <.0001
##  treatment2 / treatment3      0.398 0.1344 Inf    1  -2.729  0.0499
##  treatment2 / treatment4      1.451 0.3860 Inf    1   1.401  0.6270
##  treatment2 / treatment5      0.823 0.2331 Inf    1  -0.687  0.9594
##  treatment3 / treatment4      3.647 1.2032 Inf    1   3.921  0.0008
##  treatment3 / treatment5      2.069 0.7113 Inf    1   2.114  0.2141
##  treatment4 / treatment5      0.567 0.1555 Inf    1  -2.068  0.2339
## 
## Results are averaged over the levels of: qro 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

emmdf2 <- as.data.frame(emm2)

p2 <- ggplot(data = emmdf2, aes(x=treatment, y=prob, fill=treatment)) + geom_bar(stat = "identity", position = "dodge", color = "black")
p2 + geom_errorbar(position=position_dodge(.9),width=.25, aes(ymax=asymp.UCL, ymin=asymp.LCL),alpha=0.6)+
  labs(x="Treatment", y="Probability of Survival", fill = "treatment") + theme_bw() +
  scale_fill_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"))+
  ggtitle("Probability of a worker surviving the entire experiment")