Workers

1 Input Pollen Data

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

### Figure out average pollen consumption by treatment 


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


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

pollen <- na.omit(pollen)

range(pollen$difference)

## [1] -0.98780  1.56542

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

## [1] 0.002715 1.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.5268 -0.2476 -0.1363  0.1874  1.0576 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.460679   0.021143  21.789  < 2e-16 ***
## treatment2   0.047174   0.030208   1.562 0.118630    
## treatment3   0.100480   0.029931   3.357 0.000812 ***
## treatment4   0.053390   0.029931   1.784 0.074703 .  
## treatment5  -0.001879   0.030272  -0.062 0.950524    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3356 on 1231 degrees of freedom
## Multiple R-squared:  0.01281,    Adjusted R-squared:  0.009604 
## F-statistic: 3.994 on 4 and 1231 DF,  p-value: 0.003177

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

##  treatment emmean     SE   df lower.CL upper.CL
##  1          0.461 0.0211 1231    0.419    0.502
##  2          0.508 0.0216 1231    0.466    0.550
##  3          0.561 0.0212 1231    0.520    0.603
##  4          0.514 0.0212 1231    0.473    0.556
##  5          0.459 0.0217 1231    0.416    0.501
## 
## Confidence level used: 0.95

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

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

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

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$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 <- surviving_workers_per_colony <- 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.954
## 2 2         0.983
## 3 3         0.984
## 4 4         0.940
## 5 5         0.981

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

survive.sum

## # A tibble: 5 × 2
##   treatment    sm
##   <fct>     <dbl>
## 1 1          3.83
## 2 2          4.42
## 3 3          4.75
## 4 4          4.25
## 5 5          4.34

2.2 Workers alive at end of experiment

Data is underdispersed

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

## 
## Call:
## glm(formula = alive ~ treatment + whole.mean + round + qro, family = "poisson", 
##     data = wrpm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6126  -0.1836   0.1281   0.3376   0.9117  
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.17250    0.18262   6.421 1.36e-10 ***
## treatment2   0.12086    0.09033   1.338 0.180897    
## treatment3   0.13648    0.09185   1.486 0.137331    
## treatment4   0.06647    0.09322   0.713 0.475797    
## treatment5   0.13529    0.09228   1.466 0.142627    
## whole.mean   0.69594    0.17191   4.048 5.16e-05 ***
## round2      -0.06427    0.15881  -0.405 0.685708    
## qroB3       -0.16415    0.10891  -1.507 0.131767    
## qroB4       -0.38831    0.11987  -3.240 0.001197 ** 
## qroB5       -0.31357    0.08761  -3.579 0.000345 ***
## qroK1       -0.10975    0.16857  -0.651 0.515011    
## qroK2       -0.09008    0.25302  -0.356 0.721813    
## qroK3             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 203.72  on 296  degrees of freedom
## Residual deviance: 159.08  on 285  degrees of freedom
## AIC: 1129.1
## 
## Number of Fisher Scoring iterations: 5

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

## 
## Call:
## glm.nb(formula = alive ~ treatment + whole.mean + round + qro, 
##     data = wrpm, init.theta = 189319.1899, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6126  -0.1836   0.1281   0.3376   0.9117  
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.17250    0.18262   6.421 1.36e-10 ***
## treatment2   0.12086    0.09033   1.338 0.180902    
## treatment3   0.13648    0.09185   1.486 0.137334    
## treatment4   0.06647    0.09322   0.713 0.475803    
## treatment5   0.13529    0.09228   1.466 0.142631    
## whole.mean   0.69594    0.17191   4.048 5.16e-05 ***
## round2      -0.06427    0.15881  -0.405 0.685711    
## qroB3       -0.16415    0.10892  -1.507 0.131771    
## qroB4       -0.38831    0.11987  -3.240 0.001197 ** 
## qroB5       -0.31357    0.08761  -3.579 0.000345 ***
## qroK1       -0.10975    0.16857  -0.651 0.515017    
## qroK2       -0.09009    0.25302  -0.356 0.721816    
## qroK3             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(189319.2) family taken to be 1)
## 
##     Null deviance: 203.72  on 296  degrees of freedom
## Residual deviance: 159.08  on 285  degrees of freedom
## AIC: 1131.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  189319 
##           Std. Err.:  1487176 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -1105.118

AIC(workers.gn, workers.gnb)

##             df      AIC
## workers.gn  12 1129.113
## workers.gnb 13 1131.118

Anova(workers.gnb)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment    3.3219  4   0.505467    
## whole.mean  16.5063  1  4.849e-05 ***
## round                0               
## qro         20.7861  5   0.000889 ***
## ---
## 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     3.322  4  0.5054567    
## whole.mean   16.506  1  4.848e-05 ***
## round                0               
## qro          20.787  5  0.0008888 ***
## ---
## 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.6126  -0.1836   0.1281   0.3376   0.9117  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.10823    0.10297  10.763  < 2e-16 ***
## treatment2   0.12086    0.09033   1.338 0.180897    
## treatment3   0.13648    0.09185   1.486 0.137331    
## treatment4   0.06647    0.09322   0.713 0.475797    
## treatment5   0.13529    0.09228   1.466 0.142627    
## whole.mean   0.69594    0.17191   4.048 5.16e-05 ***
## qroB3       -0.16415    0.10891  -1.507 0.131767    
## qroB4       -0.38831    0.11987  -3.240 0.001197 ** 
## qroB5       -0.31357    0.08761  -3.579 0.000345 ***
## qroK1       -0.04548    0.07980  -0.570 0.568748    
## qroK2       -0.02582    0.20374  -0.127 0.899160    
## qroK3        0.06427    0.15881   0.405 0.685708    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 203.72  on 296  degrees of freedom
## Residual deviance: 159.08  on 285  degrees of freedom
## AIC: 1129.1
## 
## 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 = 189318.4598, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6126  -0.1836   0.1281   0.3376   0.9117  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.10823    0.10297  10.763  < 2e-16 ***
## treatment2   0.12086    0.09033   1.338 0.180902    
## treatment3   0.13648    0.09185   1.486 0.137334    
## treatment4   0.06647    0.09322   0.713 0.475803    
## treatment5   0.13529    0.09228   1.466 0.142631    
## whole.mean   0.69594    0.17191   4.048 5.16e-05 ***
## qroB3       -0.16415    0.10892  -1.507 0.131771    
## qroB4       -0.38831    0.11987  -3.240 0.001197 ** 
## qroB5       -0.31357    0.08761  -3.579 0.000345 ***
## qroK1       -0.04548    0.07981  -0.570 0.568755    
## qroK2       -0.02582    0.20374  -0.127 0.899162    
## qroK3        0.06427    0.15881   0.405 0.685711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(189318.5) family taken to be 1)
## 
##     Null deviance: 203.72  on 296  degrees of freedom
## Residual deviance: 159.08  on 285  degrees of freedom
## AIC: 1131.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  189318 
##           Std. Err.:  1487183 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -1105.118

AIC(workers.gn1, workers.gnb1)

##              df      AIC
## workers.gn1  12 1129.113
## workers.gnb1 13 1131.118

Anova(workers.gn1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: alive
##            LR Chisq Df Pr(>Chisq)    
## treatment     3.322  4  0.5054567    
## whole.mean   16.506  1  4.848e-05 ***
## qro          23.894  6  0.0005463 ***
## ---
## 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    3.3219  4  0.5054666    
## whole.mean  16.5063  1  4.849e-05 ***
## qro         23.8933  6  0.0005464 ***
## ---
## 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.83    60 1.74     60 0.224 
## 2 2          4.42    60 1.39     60 0.180 
## 3 3          4.75    60 0.600    60 0.0775
## 4 4          4.25    59 1.25     59 0.163 
## 5 5          4.34    58 1.41     58 0.185

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 proportion analysis

prop.mort <- glm(cbind(dead, alive) ~ treatment + round + start_wet_weight + qro, data = wrpm, family = binomial("logit"))

summary(prop.mort)

## 
## Call:
## glm(formula = cbind(dead, alive) ~ treatment + round + start_wet_weight + 
##     qro, family = binomial("logit"), data = wrpm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6042  -0.9743  -0.5349  -0.1506   4.7701  
## 
## Coefficients: (1 not defined because of singularities)
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -15.7603   549.5565  -0.029  0.97712    
## treatment2        -1.1161     0.2487  -4.488 7.20e-06 ***
## treatment3        -2.0306     0.3164  -6.417 1.39e-10 ***
## treatment4        -0.7150     0.2363  -3.025  0.00248 ** 
## treatment5        -1.0652     0.2465  -4.321 1.56e-05 ***
## round2            15.3285   549.5565   0.028  0.97775    
## start_wet_weight  -7.4146     1.5949  -4.649 3.34e-06 ***
## qroB3              1.1737     0.2814   4.170 3.04e-05 ***
## qroB4              1.9402     0.2787   6.961 3.37e-12 ***
## qroB5              2.0258     0.2228   9.094  < 2e-16 ***
## qroK1             13.0243   549.5567   0.024  0.98109    
## qroK2             15.8211   549.5567   0.029  0.97703    
## qroK3                  NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 786.31  on 296  degrees of freedom
## Residual deviance: 547.24  on 285  degrees of freedom
## AIC: 706.57
## 
## Number of Fisher Scoring iterations: 15

propdead <- cbind(wrpm$dead, wrpm$alive)
propdead

##        [,1] [,2]
##   [1,]    0    5
##   [2,]    0    5
##   [3,]    0    5
##   [4,]    0    5
##   [5,]    0    5
##   [6,]    5    0
##   [7,]    5    0
##   [8,]    5    0
##   [9,]    5    0
##  [10,]    5    0
##  [11,]    0    5
##  [12,]    0    5
##  [13,]    0    5
##  [14,]    0    5
##  [15,]    0    5
##  [16,]    0    5
##  [17,]    0    5
##  [18,]    0    5
##  [19,]    0    5
##  [20,]    0    5
##  [21,]    0    5
##  [22,]    0    5
##  [23,]    0    5
##  [24,]    0    5
##  [25,]    0    5
##  [26,]    1    4
##  [27,]    1    4
##  [28,]    1    4
##  [29,]    1    4
##  [30,]    1    4
##  [31,]    0    5
##  [32,]    0    5
##  [33,]    0    5
##  [34,]    0    5
##  [35,]    0    5
##  [36,]    0    5
##  [37,]    0    5
##  [38,]    0    5
##  [39,]    0    5
##  [40,]    0    5
##  [41,]    4    1
##  [42,]    4    1
##  [43,]    4    1
##  [44,]    4    1
##  [45,]    4    1
##  [46,]    1    4
##  [47,]    1    4
##  [48,]    1    4
##  [49,]    1    4
##  [50,]    1    4
##  [51,]    0    5
##  [52,]    0    5
##  [53,]    0    5
##  [54,]    0    5
##  [55,]    0    5
##  [56,]    3    2
##  [57,]    3    2
##  [58,]    3    2
##  [59,]    3    2
##  [60,]    3    2
##  [61,]    0    5
##  [62,]    0    5
##  [63,]    0    5
##  [64,]    0    5
##  [65,]    0    5
##  [66,]    0    5
##  [67,]    0    5
##  [68,]    0    5
##  [69,]    0    5
##  [70,]    0    5
##  [71,]    0    5
##  [72,]    0    5
##  [73,]    0    5
##  [74,]    0    5
##  [75,]    0    5
##  [76,]    0    5
##  [77,]    0    5
##  [78,]    0    5
##  [79,]    0    5
##  [80,]    0    5
##  [81,]    0    5
##  [82,]    0    5
##  [83,]    0    5
##  [84,]    0    5
##  [85,]    0    5
##  [86,]    1    4
##  [87,]    1    4
##  [88,]    1    4
##  [89,]    1    4
##  [90,]    1    4
##  [91,]    0    5
##  [92,]    0    5
##  [93,]    0    5
##  [94,]    0    5
##  [95,]    0    5
##  [96,]    0    5
##  [97,]    0    5
##  [98,]    0    5
##  [99,]    0    5
## [100,]    0    5
## [101,]    5    0
## [102,]    5    0
## [103,]    5    0
## [104,]    5    0
## [105,]    5    0
## [106,]    1    4
## [107,]    1    4
## [108,]    1    4
## [109,]    1    4
## [110,]    1    4
## [111,]    0    5
## [112,]    0    5
## [113,]    0    5
## [114,]    0    5
## [115,]    0    5
## [116,]    0    5
## [117,]    0    5
## [118,]    0    5
## [119,]    0    5
## [120,]    0    5
## [121,]    0    5
## [122,]    0    5
## [123,]    0    5
## [124,]    0    5
## [125,]    0    5
## [126,]    0    5
## [127,]    0    5
## [128,]    0    5
## [129,]    0    5
## [130,]    0    5
## [131,]    0    5
## [132,]    0    5
## [133,]    0    5
## [134,]    0    5
## [135,]    0    5
## [136,]    0    5
## [137,]    0    5
## [138,]    0    5
## [139,]    0    5
## [140,]    0    5
## [141,]    0    5
## [142,]    0    5
## [143,]    0    5
## [144,]    0    5
## [145,]    0    5
## [146,]    0    5
## [147,]    0    5
## [148,]    0    5
## [149,]    0    5
## [150,]    0    5
## [151,]    0    5
## [152,]    0    5
## [153,]    0    5
## [154,]    0    5
## [155,]    0    5
## [156,]    0    5
## [157,]    0    5
## [158,]    0    5
## [159,]    0    5
## [160,]    0    5
## [161,]    2    3
## [162,]    2    3
## [163,]    2    3
## [164,]    2    3
## [165,]    2    3
## [166,]    1    4
## [167,]    1    4
## [168,]    1    4
## [169,]    1    4
## [170,]    1    4
## [171,]    0    5
## [172,]    0    5
## [173,]    0    5
## [174,]    0    5
## [175,]    0    5
## [176,]    0    5
## [177,]    0    5
## [178,]    0    5
## [179,]    0    5
## [180,]    0    5
## [181,]    0    5
## [182,]    0    5
## [183,]    0    5
## [184,]    0    5
## [185,]    0    5
## [186,]    2    3
## [187,]    2    3
## [188,]    2    3
## [189,]    2    3
## [190,]    2    3
## [191,]    0    5
## [192,]    0    5
## [193,]    0    5
## [194,]    0    5
## [195,]    0    5
## [196,]    0    5
## [197,]    0    5
## [198,]    0    5
## [199,]    0    5
## [200,]    0    5
## [201,]    0    5
## [202,]    0    5
## [203,]    0    5
## [204,]    0    5
## [205,]    0    5
## [206,]    2    3
## [207,]    2    3
## [208,]    2    3
## [209,]    2    3
## [210,]    2    3
## [211,]    0    5
## [212,]    0    5
## [213,]    0    5
## [214,]    0    5
## [215,]    0    5
## [216,]    0    5
## [217,]    0    5
## [218,]    0    5
## [219,]    0    5
## [220,]    0    5
## [221,]    0    5
## [222,]    0    5
## [223,]    0    5
## [224,]    0    5
## [225,]    0    5
## [226,]    4    1
## [227,]    4    1
## [228,]    4    1
## [229,]    4    1
## [230,]    4    1
## [231,]    0    5
## [232,]    0    5
## [233,]    0    5
## [234,]    0    5
## [235,]    0    5
## [236,]    1    4
## [237,]    1    4
## [238,]    1    4
## [239,]    1    4
## [240,]    0    5
## [241,]    0    5
## [242,]    0    5
## [243,]    0    5
## [244,]    0    5
## [245,]    0    5
## [246,]    0    5
## [247,]    0    5
## [248,]    0    5
## [249,]    0    5
## [250,]    1    4
## [251,]    1    4
## [252,]    1    4
## [253,]    1    4
## [254,]    0    5
## [255,]    0    5
## [256,]    0    5
## [257,]    0    5
## [258,]    0    5
## [259,]    0    5
## [260,]    0    5
## [261,]    0    5
## [262,]    0    5
## [263,]    0    5
## [264,]    1    4
## [265,]    1    4
## [266,]    1    4
## [267,]    1    4
## [268,]    1    4
## [269,]    1    4
## [270,]    1    4
## [271,]    1    4
## [272,]    1    4
## [273,]    5    0
## [274,]    5    0
## [275,]    5    0
## [276,]    5    0
## [277,]    5    0
## [278,]    0    5
## [279,]    0    5
## [280,]    0    5
## [281,]    0    5
## [282,]    0    5
## [283,]    0    5
## [284,]    0    5
## [285,]    0    5
## [286,]    0    5
## [287,]    0    5
## [288,]    0    5
## [289,]    0    5
## [290,]    0    5
## [291,]    0    5
## [292,]    0    5
## [293,]    0    5
## [294,]    0    5
## [295,]    0    5
## [296,]    0    5
## [297,]    0    5

Anova(prop.mort)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cbind(dead, alive)
##                  LR Chisq Df Pr(>Chisq)    
## treatment          56.362  4  1.684e-11 ***
## round                      0               
## start_wet_weight           0               
## qro               117.240  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

## integer(0)

qqnorm(resid(prop.mort));qqline(resid(prop.mort))  # yikes 

mort.emm <- emmeans(prop.mort, "treatment")
pairs(mort.emm)

##  contrast                estimate    SE  df z.ratio p.value
##  treatment1 - treatment2    1.116 0.249 Inf   4.488  0.0001
##  treatment1 - treatment3    2.031 0.316 Inf   6.417  <.0001
##  treatment1 - treatment4    0.715 0.236 Inf   3.025  0.0210
##  treatment1 - treatment5    1.065 0.247 Inf   4.321  0.0002
##  treatment2 - treatment3    0.914 0.335 Inf   2.730  0.0498
##  treatment2 - treatment4   -0.401 0.263 Inf  -1.523  0.5472
##  treatment2 - treatment5   -0.051 0.268 Inf  -0.190  0.9997
##  treatment3 - treatment4   -1.316 0.328 Inf  -4.011  0.0006
##  treatment3 - treatment5   -0.965 0.332 Inf  -2.909  0.0298
##  treatment4 - treatment5    0.350 0.261 Inf   1.341  0.6657
## 
## Results are averaged over the levels of: qro, round 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 5 estimates

summary(mort.emm)

##  treatment emmean   SE  df asymp.LCL asymp.UCL
##  1          -3.85 91.6 Inf      -183       176
##  2          -4.96 91.6 Inf      -184       175
##  3          -5.88 91.6 Inf      -185       174
##  4          -4.56 91.6 Inf      -184       175
##  5          -4.91 91.6 Inf      -184       175
## 
## Results are averaged over the levels of: qro, round 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95

prop.sum <- wrpm %>%
  group_by(treatment) %>%
  summarise(mp = mean(cbind(dead, alive)), 
            sd.p = sd(cbind(dead, alive)), 
            np = length(cbind(dead,alive))) %>%
  mutate(sep = sd.p / sqrt(np))

prop.sum

## # A tibble: 5 × 5
##   treatment    mp  sd.p    np   sep
##   <fct>     <dbl> <dbl> <int> <dbl>
## 1 1           2.5  2.19   120 0.200
## 2 2           2.5  2.37   120 0.217
## 3 3           2.5  2.34   120 0.213
## 4 4           2.5  2.16   118 0.199
## 5 5           2.5  2.32   116 0.216

2.2.2 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.9308, p-value = 1.565e-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.99442, p-value = 0.3498

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

## Linear mixed model fit by REML ['lmerMod']
## Formula: log.dry ~ treatment + +whole.mean + round + qro + (1 | colony)
##    Data: wrpm
## 
## REML criterion at convergence: 167.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9098 -0.6741 -0.0030  0.6681  3.1884 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  colony   (Intercept) 0.0000   0.0000  
##  Residual             0.0916   0.3027  
## Number of obs: 297, groups:  colony, 60
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) -3.680029   0.119934 -30.684
## treatment2  -0.058064   0.055438  -1.047
## treatment3  -0.017046   0.057139  -0.298
## treatment4  -0.008591   0.056498  -0.152
## treatment5   0.059820   0.056546   1.058
## whole.mean   0.801188   0.105995   7.559
## round2       0.077879   0.106895   0.729
## qroB3        0.246031   0.066474   3.701
## qroB4        0.239370   0.069856   3.427
## qroB5        0.246565   0.050810   4.853
## qroK1        0.201241   0.112882   1.783
## qroK2        0.425939   0.164120   2.595
## 
## Correlation of Fixed Effects:
##            (Intr) trtmn2 trtmn3 trtmn4 trtmn5 whl.mn round2 qroB3  qroB4 
## treatment2 -0.205                                                        
## treatment3 -0.256  0.499                                                 
## treatment4 -0.333  0.494  0.514                                          
## treatment5 -0.221  0.495  0.486  0.483                                   
## whole.mean -0.356 -0.070 -0.213 -0.089 -0.020                            
## round2     -0.858  0.000  0.117  0.152 -0.004 -0.001                     
## qroB3       0.022  0.004  0.013  0.002  0.001 -0.060 -0.105              
## qroB4       0.112  0.022  0.067  0.024  0.006 -0.313 -0.099  0.178       
## qroB5       0.024  0.005  0.014  0.001  0.001 -0.065 -0.137  0.223  0.228
## qroK1      -0.761  0.016  0.169  0.179  0.026 -0.173  0.889  0.010  0.054
## qroK2      -0.572 -0.033  0.038  0.095 -0.120  0.051  0.621 -0.004 -0.017
##            qroB5  qroK1 
## treatment2              
## treatment3              
## treatment4              
## treatment5              
## whole.mean              
## round2                  
## qroB3                   
## qroB4                   
## qroB5                   
## qroK1       0.010       
## qroK2      -0.004  0.579
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## 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   4.5054  4     0.3419    
## whole.mean 57.1344  1  4.070e-14 ***
## round       0.5308  1     0.4663    
## qro        41.9998  5  5.891e-08 ***
## ---
## 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.0482 0.0191    60
## 2 2         0.0455 0.0149    60
## 3 3         0.0504 0.0191    60
## 4 4         0.0483 0.0208    59
## 5 5         0.0505 0.0174    58

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 + round + start_wet_weight + qro, data = wrpm, family = binomial("logit"))
summary(mod1)

## 
## Call:
## glm(formula = cbind(alive, dead) ~ treatment + round + start_wet_weight + 
##     qro, family = binomial("logit"), data = wrpm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.7701   0.1506   0.5349   0.9743   2.6042  
## 
## Coefficients: (1 not defined because of singularities)
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       15.7603   549.5565   0.029  0.97712    
## treatment2         1.1161     0.2487   4.488 7.20e-06 ***
## treatment3         2.0306     0.3164   6.417 1.39e-10 ***
## treatment4         0.7150     0.2363   3.025  0.00248 ** 
## treatment5         1.0652     0.2465   4.321 1.56e-05 ***
## round2           -15.3285   549.5565  -0.028  0.97775    
## start_wet_weight   7.4146     1.5949   4.649 3.34e-06 ***
## qroB3             -1.1737     0.2814  -4.170 3.04e-05 ***
## qroB4             -1.9402     0.2787  -6.961 3.37e-12 ***
## qroB5             -2.0258     0.2228  -9.094  < 2e-16 ***
## qroK1            -13.0243   549.5567  -0.024  0.98109    
## qroK2            -15.8211   549.5567  -0.029  0.97703    
## qroK3                  NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 786.31  on 296  degrees of freedom
## Residual deviance: 547.24  on 285  degrees of freedom
## AIC: 706.57
## 
## Number of Fisher Scoring iterations: 15

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                               296     786.31
## treatment         4   46.142       292     740.17
## round             1   73.578       291     666.59
## start_wet_weight  1    2.113       290     664.48
## qro               5  117.240       285     547.24

Anova(mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cbind(alive, dead)
##                  LR Chisq Df Pr(>Chisq)    
## treatment          56.362  4  1.684e-11 ***
## round                      0               
## start_wet_weight           0               
## qro               117.240  5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(mod1)

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

## $emmeans
##  treatment  prob    SE  df asymp.LCL asymp.UCL
##  1         0.979 1.871 Inf  2.22e-16         1
##  2         0.993 0.630 Inf  2.22e-16         1
##  3         0.997 0.255 Inf  2.22e-16         1
##  4         0.990 0.935 Inf  2.22e-16         1
##  5         0.993 0.663 Inf  2.22e-16         1
## 
## Results are averaged over the levels of: qro, round 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.328 0.0815 Inf    1  -4.488  0.0001
##  treatment1 / treatment3      0.131 0.0415 Inf    1  -6.417  <.0001
##  treatment1 / treatment4      0.489 0.1156 Inf    1  -3.025  0.0210
##  treatment1 / treatment5      0.345 0.0850 Inf    1  -4.321  0.0002
##  treatment2 / treatment3      0.401 0.1343 Inf    1  -2.730  0.0498
##  treatment2 / treatment4      1.494 0.3933 Inf    1   1.523  0.5472
##  treatment2 / treatment5      1.052 0.2825 Inf    1   0.190  0.9997
##  treatment3 / treatment4      3.727 1.2223 Inf    1   4.011  0.0006
##  treatment3 / treatment5      2.626 0.8714 Inf    1   2.909  0.0298
##  treatment4 / treatment5      0.705 0.1840 Inf    1  -1.341  0.6657
## 
## Results are averaged over the levels of: qro, round 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

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

##  treatment  prob    SE  df asymp.LCL asymp.UCL
##  1         0.979 1.871 Inf  2.22e-16         1
##  2         0.993 0.630 Inf  2.22e-16         1
##  3         0.997 0.255 Inf  2.22e-16         1
##  4         0.990 0.935 Inf  2.22e-16         1
##  5         0.993 0.663 Inf  2.22e-16         1
## 
## Results are averaged over the levels of: qro, round 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale

pairs(emm1)

##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.328 0.0815 Inf    1  -4.488  0.0001
##  treatment1 / treatment3      0.131 0.0415 Inf    1  -6.417  <.0001
##  treatment1 / treatment4      0.489 0.1156 Inf    1  -3.025  0.0210
##  treatment1 / treatment5      0.345 0.0850 Inf    1  -4.321  0.0002
##  treatment2 / treatment3      0.401 0.1343 Inf    1  -2.730  0.0498
##  treatment2 / treatment4      1.494 0.3933 Inf    1   1.523  0.5472
##  treatment2 / treatment5      1.052 0.2825 Inf    1   0.190  0.9997
##  treatment3 / treatment4      3.727 1.2223 Inf    1   4.011  0.0006
##  treatment3 / treatment5      2.626 0.8714 Inf    1   2.909  0.0298
##  treatment4 / treatment5      0.705 0.1840 Inf    1  -1.341  0.6657
## 
## Results are averaged over the levels of: qro, round 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

emmdf <- as.data.frame(emm1)

p <- ggplot(data = emmdf, aes(x=treatment, y=prob, fill=treatment)) + geom_bar(stat = "identity", position = "dodge", color = "black")
p + 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"))

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.979 1.871 Inf  2.22e-16         1
##  2         0.993 0.630 Inf  2.22e-16         1
##  3         0.997 0.255 Inf  2.22e-16         1
##  4         0.990 0.935 Inf  2.22e-16         1
##  5         0.993 0.663 Inf  2.22e-16         1
## 
## Results are averaged over the levels of: qro, round 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                odds.ratio     SE  df null z.ratio p.value
##  treatment1 / treatment2      0.328 0.0815 Inf    1  -4.488  0.0001
##  treatment1 / treatment3      0.131 0.0415 Inf    1  -6.417  <.0001
##  treatment1 / treatment4      0.489 0.1156 Inf    1  -3.025  0.0210
##  treatment1 / treatment5      0.345 0.0850 Inf    1  -4.321  0.0002
##  treatment2 / treatment3      0.401 0.1343 Inf    1  -2.730  0.0498
##  treatment2 / treatment4      1.494 0.3933 Inf    1   1.523  0.5472
##  treatment2 / treatment5      1.052 0.2825 Inf    1   0.190  0.9997
##  treatment3 / treatment4      3.727 1.2223 Inf    1   4.011  0.0006
##  treatment3 / treatment5      2.626 0.8714 Inf    1   2.909  0.0298
##  treatment4 / treatment5      0.705 0.1840 Inf    1  -1.341  0.6657
## 
## Results are averaged over the levels of: qro, round 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio 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")

3.1 Logit regression (https://stats.oarc.ucla.edu/r/dae/logit-regression/)

I can’t include starting wet weight in this model because I only have it on a per colony basis not on a per bee basis ( I can’t trace starting wet weight to the specific bee who might have died or not died)

4 prop survive

prop <- read_csv("propsurvive.csv")
qro <- read_csv("qro.csv")

prop.q <- merge(prop, qro, by.x = "colony")

prop.q$colony <- as.factor(prop.q$colony)
prop.q$treatment <- as.factor(prop.q$treatment)
prop.q$replicate <- as.factor(prop.q$replicate)
prop.q$round <- as.factor(prop.q$round)
prop.q$qro <- as.factor(prop.q$qro)

mylogit <- glm(alive ~ treatment + qro + round, data = prop.q, family = "binomial")


summary(mylogit)

## 
## Call:
## glm(formula = alive ~ treatment + qro + round, family = "binomial", 
##     data = prop.q)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6822  -0.6810  -0.6471  -0.6033   1.9456  
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.343e+00  3.673e-01  -3.656 0.000256 ***
## treatment2  -1.143e-01  4.710e-01  -0.243 0.808288    
## treatment3   3.623e-03  4.677e-01   0.008 0.993819    
## treatment4   3.623e-03  4.677e-01   0.008 0.993819    
## treatment5  -1.143e-01  4.710e-01  -0.243 0.808288    
## qroB3       -2.721e-01  5.897e-01  -0.461 0.644555    
## qroB4        1.949e-15  5.479e-01   0.000 1.000000    
## qroB5       -1.301e-01  4.308e-01  -0.302 0.762626    
## qroK1       -3.983e-03  4.055e-01  -0.010 0.992163    
## qroK2/K1     7.100e-02  8.440e-01   0.084 0.932959    
## qroK3       -4.691e-02  1.179e+00  -0.040 0.968264    
## qroK3/K2    -4.691e-02  1.179e+00  -0.040 0.968264    
## round2              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 294.61  on 299  degrees of freedom
## Residual deviance: 294.16  on 288  degrees of freedom
## AIC: 318.16
## 
## Number of Fisher Scoring iterations: 4

plot(mylogit)

confint(mylogit)

##                  2.5 %     97.5 %
## (Intercept) -2.1051812 -0.6548088
## treatment2  -1.0511111  0.8121875
## treatment3  -0.9251110  0.9250006
## treatment4  -0.9251110  0.9250006
## treatment5  -1.0511111  0.8121875
## qroB3       -1.5686835  0.7998330
## qroB4       -1.1739816  1.0147447
## qroB5       -1.0191917  0.6865304
## qroK1       -0.8301377  0.7722693
## qroK2/K1    -1.8897724  1.5839543
## qroK3       -3.0841371  2.0145402
## qroK3/K2    -3.0841371  2.0145402
## round2              NA         NA

odds ratio

exp(coef(mylogit))

## (Intercept)  treatment2  treatment3  treatment4  treatment5       qroB3 
##   0.2610585   0.8920051   1.0036294   1.0036294   0.8920051   0.7618067 
##       qroB4       qroB5       qroK1    qroK2/K1       qroK3    qroK3/K2 
##   1.0000000   0.8779922   0.9960247   1.0735809   0.9541765   0.9541765 
##      round2 
##          NA

exp(cbind(OR = coef(mylogit), confint(mylogit)))

##                    OR      2.5 %    97.5 %
## (Intercept) 0.2610585 0.12182360 0.5195414
## treatment2  0.8920051 0.34954914 2.2528307
## treatment3  1.0036294 0.39648741 2.5218697
## treatment4  1.0036294 0.39648741 2.5218697
## treatment5  0.8920051 0.34954914 2.2528307
## qroB3       0.7618067 0.20831925 2.2251693
## qroB4       1.0000000 0.30913363 2.7586589
## qroB5       0.8779922 0.36088651 1.9868100
## qroK1       0.9960247 0.43598923 2.1646730
## qroK2/K1    1.0735809 0.15110619 4.8741919
## qroK3       0.9541765 0.04576951 7.4972797
## qroK3/K2    0.9541765 0.04576951 7.4972797
## round2             NA         NA        NA