Colony Weight Change

1 Input Pollen Data

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


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

pollen <- na.omit(pollen)

range(pollen$difference)

## [1] -0.98780  1.56542

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

## [1] 0.002715 1.565420

# 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")

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 Input Weight Data

w <- read_csv("weights1.csv", col_types = cols(round = col_factor(levels = c("1", 
    "2")), date = col_date(format = "%m/%d/%Y"), 
    treatment = col_factor(levels = c("1", 
        "2", "3", "4", "5"))))

w$number <- as.factor(w$number)
w$replicate <- as.factor(w$replicate)
w$id <- as.factor(w$id)
w$colony <- as.factor(w$colony)

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


qro <- read_csv("qro.csv")
qro$colony <- as.factor(qro$colony)
qro$qro <- as.factor(qro$qro)

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

3 Average Weekly Colony Weights

3.1 Shape of colony weights per week

hist(w$weight)

shapiro.test(w$weight)

## 
##  Shapiro-Wilk normality test
## 
## data:  w$weight
## W = 0.96315, p-value = 1.544e-09

w$logw <- log(w$weight)

shapiro.test(w$logw)

## 
##  Shapiro-Wilk normality test
## 
## data:  w$logw
## W = 0.97882, p-value = 2.124e-06

hist(w$logw)

w$boxw <- bcPower(w$weight, -2, gamma=1)

shapiro.test(w$boxw)

## 
##  Shapiro-Wilk normality test
## 
## data:  w$boxw
## W = 0.99084, p-value = 0.004818

hist(w$boxw)

wmod <- glm(boxw ~ treatment + whole.mean + bees_alive + round + qro + days, data = w)
summary(wmod)

## 
## Call:
## glm(formula = boxw ~ treatment + whole.mean + bees_alive + round + 
##     qro + days, data = w)
## 
## Deviance Residuals: 
##        Min          1Q      Median          3Q         Max  
## -5.164e-05  -8.372e-06   5.930e-07   9.267e-06   4.841e-05  
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  4.997e-01  8.501e-06 58785.055  < 2e-16 ***
## treatment2   1.619e-06  2.191e-06     0.739  0.46034    
## treatment3   2.462e-07  2.214e-06     0.111  0.91153    
## treatment4   1.816e-06  2.191e-06     0.829  0.40769    
## treatment5   2.913e-06  2.213e-06     1.316  0.18874    
## whole.mean   3.865e-05  4.402e-06     8.781  < 2e-16 ***
## bees_alive   3.101e-06  1.297e-06     2.390  0.01726 *  
## round2       3.135e-05  5.516e-06     5.684 2.34e-08 ***
## qroB3       -5.373e-07  2.764e-06    -0.194  0.84594    
## qroB4        2.243e-07  2.909e-06     0.077  0.93856    
## qroB5        1.349e-06  1.972e-06     0.684  0.49415    
## qroK1        1.594e-05  5.673e-06     2.809  0.00518 ** 
## qroK2/K1     1.996e-05  6.731e-06     2.965  0.00318 ** 
## qroK3        5.824e-06  7.748e-06     0.752  0.45261    
## qroK3/K2            NA         NA        NA       NA    
## days         1.497e-06  4.418e-08    33.890  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.183067e-10)
## 
##     Null deviance: 3.8097e-07  on 474  degrees of freedom
## Residual deviance: 1.0042e-07  on 460  degrees of freedom
## AIC: -9201.7
## 
## Number of Fisher Scoring iterations: 2

Anova(wmod)

## Analysis of Deviance Table (Type II tests)
## 
## Response: boxw
##            LR Chisq Df Pr(>Chisq)    
## treatment      2.33  4    0.67516    
## whole.mean    77.10  1    < 2e-16 ***
## bees_alive     5.71  1    0.01686 *  
## round                0               
## qro           12.99  6    0.04327 *  
## days        1148.52  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(wmod)  # these are okay 

wem <- emmeans(wmod, "treatment")
pairs(wem)

##  contrast                 estimate       SE  df t.ratio p.value
##  treatment1 - treatment2 -1.62e-06 2.19e-06 460  -0.739  0.9473
##  treatment1 - treatment3 -2.46e-07 2.21e-06 460  -0.111  1.0000
##  treatment1 - treatment4 -1.82e-06 2.19e-06 460  -0.829  0.9217
##  treatment1 - treatment5 -2.91e-06 2.21e-06 460  -1.316  0.6812
##  treatment2 - treatment3  1.37e-06 2.21e-06 460   0.620  0.9719
##  treatment2 - treatment4 -1.97e-07 2.24e-06 460  -0.088  1.0000
##  treatment2 - treatment5 -1.29e-06 2.20e-06 460  -0.587  0.9770
##  treatment3 - treatment4 -1.57e-06 2.23e-06 460  -0.703  0.9558
##  treatment3 - treatment5 -2.67e-06 2.27e-06 460  -1.173  0.7667
##  treatment4 - treatment5 -1.10e-06 2.27e-06 460  -0.482  0.9890
## 
## Results are averaged over the levels of: qro, round 
## P value adjustment: tukey method for comparing a family of 5 estimates

difmod <- glm(diff ~ treatment + whole.mean + bees_alive + round + qro + days, data = w)
summary(difmod)

## 
## Call:
## glm(formula = diff ~ treatment + whole.mean + bees_alive + round + 
##     qro + days, data = w)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.0143  -1.5313  -0.0843   1.4557  10.3641  
## 
## Coefficients: (1 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.460642   1.470459  -5.754  1.6e-08 ***
## treatment2   0.171997   0.378983   0.454 0.650160    
## treatment3   0.375307   0.382994   0.980 0.327636    
## treatment4   0.476712   0.379069   1.258 0.209179    
## treatment5  -0.056555   0.382806  -0.148 0.882615    
## whole.mean   6.841606   0.761397   8.986  < 2e-16 ***
## bees_alive   0.778262   0.224424   3.468 0.000574 ***
## round2       1.179766   0.954093   1.237 0.216892    
## qroB3       -0.035674   0.478071  -0.075 0.940548    
## qroB4       -0.255983   0.503157  -0.509 0.611168    
## qroB5        0.104310   0.341066   0.306 0.759867    
## qroK1        2.772881   0.981212   2.826 0.004919 ** 
## qroK2/K1     3.006329   1.164346   2.582 0.010131 *  
## qroK3        1.768041   1.340263   1.319 0.187767    
## qroK3/K2           NA         NA      NA       NA    
## days         0.243996   0.007642  31.929  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.531789)
## 
##     Null deviance: 11218.8  on 474  degrees of freedom
## Residual deviance:  3004.6  on 460  degrees of freedom
## AIC: 2256.2
## 
## Number of Fisher Scoring iterations: 2

Anova(difmod)

## Analysis of Deviance Table (Type II tests)
## 
## Response: diff
##            LR Chisq Df Pr(>Chisq)    
## treatment      2.80  4  0.5920500    
## whole.mean    80.74  1  < 2.2e-16 ***
## bees_alive    12.03  1  0.0005247 ***
## round                0               
## qro            9.74  6  0.1359059    
## days        1019.46  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(difmod)

ggplot(w) +
  aes( x = days, y = weight, color = treatment) +
  geom_point(aes(color = treatment)) +
  geom_smooth() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
                    name = "Pristine Level",
                    labels = c("Treatment 1 (control)", "Treatment 2", 
                               "Treatment 3", "Treatment 4", "Treatment 5")) + 
  ylab("Time (days since colony initiation") +
  xlab("Weight (g)")+
  ggtitle("Colony Weight Change Over Time (g)") +
  facet_wrap(vars(treatment))

ggplot(w) +
  aes( x = days, y = diff, color = treatment) +
  geom_point(aes(color = treatment)) +
  geom_smooth() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
                    name = "Pristine Level",
                    labels = c("Treatment 1 (control)", "Treatment 2", 
                               "Treatment 3", "Treatment 4", "Treatment 5")) + 
  ylab("Time (days since colony initiation") +
  xlab("Weight Change(g)")+
  ggtitle("Colony Weight (g) Change Over Time Since First Weigh-in") +
  facet_wrap(vars(treatment))

ggplot(w) +
  aes( x = days, y = change, color = treatment) +
  geom_point(aes(color = treatment)) +
  geom_smooth() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
                    name = "Pristine Level",
                    labels = c("Treatment 1 (control)", "Treatment 2", 
                               "Treatment 3", "Treatment 4", "Treatment 5")) + 
  ylab("Time (days since colony initiation") +
  xlab("Weight Change(g) Since the Previous Week")+
  ggtitle("Colony Weight (g) Change Each Week") +
  facet_wrap(vars(treatment))

colsum <- w %>%
  group_by(treatment) %>%
  summarise( mw = mean(weight),
             sdw = sd(weight),
             nw = length(weight)) %>%
  mutate(sew = sdw/ sqrt(nw))

colsum

## # A tibble: 5 × 5
##   treatment    mw   sdw    nw   sew
##   <fct>     <dbl> <dbl> <int> <dbl>
## 1 1          53.2  4.61    97 0.468
## 2 2          53.5  4.38    93 0.455
## 3 3          53.9  5.25    98 0.530
## 4 4          53.7  5.21    98 0.526
## 5 5          53.4  4.20    89 0.445

3.2 Cut off data at day 47 - last day we have data for all treatments

w47 <- subset(w, days<47)

w47mod <- glm(boxw ~ treatment + whole.mean + bees_alive + round + qro + days, data = w47)
summary(w47mod)

## 
## Call:
## glm(formula = boxw ~ treatment + whole.mean + bees_alive + round + 
##     qro + days, data = w47)
## 
## Deviance Residuals: 
##        Min          1Q      Median          3Q         Max  
## -4.378e-05  -7.791e-06   2.490e-07   7.998e-06   4.412e-05  
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  4.997e-01  8.300e-06 60209.043  < 2e-16 ***
## treatment2   5.046e-07  2.064e-06     0.244  0.80700    
## treatment3  -1.197e-07  2.113e-06    -0.057  0.95483    
## treatment4   1.338e-06  2.103e-06     0.637  0.52475    
## treatment5   3.067e-06  2.102e-06     1.459  0.14521    
## whole.mean   3.198e-05  4.233e-06     7.554 2.50e-13 ***
## bees_alive   2.710e-06  1.260e-06     2.152  0.03198 *  
## round2       3.041e-05  5.449e-06     5.581 4.22e-08 ***
## qroB3       -1.052e-06  2.568e-06    -0.410  0.68220    
## qroB4        2.217e-06  2.750e-06     0.806  0.42053    
## qroB5        1.215e-06  1.878e-06     0.647  0.51808    
## qroK1        1.494e-05  5.596e-06     2.671  0.00786 ** 
## qroK2/K1     1.874e-05  6.645e-06     2.820  0.00502 ** 
## qroK3        2.915e-06  7.440e-06     0.392  0.69539    
## qroK3/K2            NA         NA        NA       NA    
## days         1.639e-06  4.549e-08    36.027  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.877838e-10)
## 
##     Null deviance: 3.5886e-07  on 450  degrees of freedom
## Residual deviance: 8.1874e-08  on 436  degrees of freedom
## AIC: -8803.8
## 
## Number of Fisher Scoring iterations: 2

Anova(w47mod)

## Analysis of Deviance Table (Type II tests)
## 
## Response: boxw
##            LR Chisq Df Pr(>Chisq)    
## treatment      3.05  4    0.54979    
## whole.mean    57.07  1   4.21e-14 ***
## bees_alive     4.63  1    0.03143 *  
## round                0               
## qro           14.46  6    0.02487 *  
## days        1297.97  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(w47mod) 

hist(w47$change)

shapiro.test(w47$change)

## 
##  Shapiro-Wilk normality test
## 
## data:  w47$change
## W = 0.92501, p-value = 3.162e-14

w47$logc <- log((w47$change)+2)
shapiro.test(w47$logc)

## 
##  Shapiro-Wilk normality test
## 
## data:  w47$logc
## W = 0.97236, p-value = 1.632e-07

hist(w47$logc)

w47dif <- glm(change ~ treatment + whole.mean + bees_alive + round + qro + days, data = w47)
summary(w47dif)

## 
## Call:
## glm(formula = change ~ treatment + whole.mean + bees_alive + 
##     round + qro + days, data = w47)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.1248  -1.0261  -0.3313   0.9702   7.2373  
## 
## Coefficients: (1 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.132107   0.957323  -2.227 0.026447 *  
## treatment2  -0.105993   0.238084  -0.445 0.656404    
## treatment3   0.018367   0.243657   0.075 0.939948    
## treatment4   0.132612   0.242507   0.547 0.584771    
## treatment5   0.047297   0.242440   0.195 0.845416    
## whole.mean   1.898744   0.488202   3.889 0.000116 ***
## bees_alive   0.358055   0.145283   2.465 0.014104 *  
## round2       0.348660   0.628432   0.555 0.579309    
## qroB3       -0.138172   0.296152  -0.467 0.641050    
## qroB4        0.156832   0.317131   0.495 0.621178    
## qroB5        0.030872   0.216604   0.143 0.886728    
## qroK1        0.594514   0.645420   0.921 0.357493    
## qroK2/K1     0.666092   0.766454   0.869 0.385294    
## qroK3        0.249402   0.858165   0.291 0.771479    
## qroK3/K2           NA         NA      NA       NA    
## days         0.024203   0.005247   4.613 5.22e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.498135)
## 
##     Null deviance: 1241.0  on 450  degrees of freedom
## Residual deviance: 1089.2  on 436  degrees of freedom
## AIC: 1709.5
## 
## Number of Fisher Scoring iterations: 2

Anova(w47dif)

## Analysis of Deviance Table (Type II tests)
## 
## Response: change
##            LR Chisq Df Pr(>Chisq)    
## treatment    1.0091  4  0.9084204    
## whole.mean  15.1263  1  0.0001006 ***
## bees_alive   6.0739  1  0.0137192 *  
## round                0               
## qro          1.7167  6  0.9438261    
## days        21.2801  1  3.968e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(w47dif) 

ggplot(w47) +
  aes( x = days, y = weight, color = treatment) +
  geom_point(aes(color = treatment)) +
  geom_smooth() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("peachpuff3", "darkseagreen", "lightseagreen", "darkolivegreen3", "darkolivegreen4"),
                    name = "Pristine Level",
                    labels = c("Treatment 1 (control)", "Treatment 2", 
                               "Treatment 3", "Treatment 4", "Treatment 5")) + 
  ylab("Time (days since colony initiation") +
  xlab("Weight (g) ")+
  ggtitle("Colony Weight (g) Each Week") +
  facet_wrap(vars(treatment))