load libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(ggpubr)
#library(reshape2)
#library(hrbrthemes)
library("gridExtra")
#library(cowplot)
library(plotly)
library(scales) # to calculate percentages, and put into dataframe
library(multcompView) # to correct for multiple comparisons after Tukey HSD
knitr::opts_chunk$set(echo = TRUE)

load data

data <- read.csv("/Users/nuriteliash/Documents/GitHub/OA_Biogaz/data/OA_Yosef.csv")

df = data %>% 
  #dplyr::mutate(fall_norm = fall/falling_time) %>%
  #dplyr::mutate(date = as.Date(date, format = "%d/%m/%Y")) %>% # Convert to Date object 
    dplyr::mutate(date = factor(date, level = c("March10th", "May7th"))) %>%
dplyr::mutate(varroa = as.numeric(varroa)) %>%
  dplyr::mutate(hive = as.character(hive)) %>%
  dplyr::mutate(honey_weight_net = as.numeric(honey_weight_net)) %>%
  dplyr::mutate(bees_strength = as.numeric(bees_strength)) %>%
  dplyr::mutate(honey_stores = as.numeric(honey_stores))

Honey yield

filtered_df <- df %>%
  filter(date == "May7th") %>%
  na.omit()

p_honey = filtered_df %>% 
  ggplot( aes(x = treatment, y = honey_weight_net, color = treatment, text = paste("Hive:", hive, "<br>Weight:", honey_weight_net))) +
  geom_boxplot() + 
  geom_jitter(width = 0.2) +
  theme(legend.position = "none") +
  ggtitle("Honey yield, May7th") +
  facet_wrap(~ yard) 

# Convert the ggplot object to a plotly object and add hover information
p_honey_plotly <- ggplotly(p_honey, tooltip = "text")

# Print the plotly object
p_honey_plotly
# Run the ANOVA
anova_result <- aov(honey_weight_net ~ treatment * yard, data = filtered_df) %>%
  na.omit()

# Summarize the ANOVA results
summary(anova_result)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment       2    401   200.4   2.566   0.0827 .  
## yard            2   3874  1937.2  24.806 3.12e-09 ***
## treatment:yard  3    128    42.6   0.545   0.6525    
## Residuals      86   6716    78.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Diagnostic plots for ANOVA
#par(mfrow = c(2, 2)) # Set up the plotting area for 2x2 plots
#plot(anova_result)

# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for treatment
tukey_treatment <- TukeyHSD(anova_result, "treatment")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = honey_weight_net ~ treatment * yard, data = filtered_df)
## 
## $treatment
##                         diff        lwr       upr     p adj
## Control-cardboards  2.815421  -2.091936 7.7227780 0.3619976
## Towels-cardboards  -2.689936  -8.486582 3.1067104 0.5125506
## Towels-Control     -5.505357 -11.413220 0.4025059 0.0730716
# Extract the Tukey HSD results for yard
tukey_yard <- TukeyHSD(anova_result, "yard")
print(tukey_yard)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = honey_weight_net ~ treatment * yard, data = filtered_df)
## 
## $yard
##                    diff      lwr       upr     p adj
## Mujea-Biogaz  13.659144  8.57972 18.738567 0.0000000
## Zeytim-Biogaz 11.369619  5.78125 16.957989 0.0000160
## Zeytim-Mujea  -2.289524 -7.77661  3.197562 0.5819204
# Extract the Tukey HSD results for treatment:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "treatment:yard")
#print(tukey_interaction)

Conclusions, honey yield

  • The biggest difference was between yards: Biogaz had the lowest honey yields, regardless of treatment.
  • No difference in honey yield between the treatments.

Colony strength

filtered_df <- df %>%
  filter(date == "May7th") %>%
  na.omit()

p_colony = filtered_df %>% 
  ggplot( aes(x = treatment, y = bees_strength, color = treatment, text = paste("Hive:", hive, "<br>Strength:", bees_strength))) +
  geom_boxplot() + 
  geom_jitter(width = 0.2) +
  theme(legend.position = "none") +
  ggtitle("Colony strength, May7th") +
  facet_wrap(~ yard) 

# Convert the ggplot object to a plotly object and add hover information
p_colony_plotly <- ggplotly(p_colony, tooltip = "text")

# Print the plotly object
p_colony_plotly
# Run the ANOVA
anova_result <- aov(bees_strength ~ treatment * yard, data = filtered_df) %>%
  na.omit()

# Summarize the ANOVA results
summary(anova_result)
##                Df Sum Sq Mean Sq F value Pr(>F)
## treatment       2   0.62  0.3120   0.464  0.631
## yard            2   2.09  1.0429   1.550  0.218
## treatment:yard  3   0.65  0.2180   0.324  0.808
## Residuals      86  57.87  0.6729
# Diagnostic plots for ANOVA
#par(mfrow = c(2, 2)) # Set up the plotting area for 2x2 plots
#plot(anova_result)

# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for treatment
tukey_treatment <- TukeyHSD(anova_result, "treatment")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = bees_strength ~ treatment * yard, data = filtered_df)
## 
## $treatment
##                           diff        lwr       upr     p adj
## Control-cardboards  0.18388278 -0.2716390 0.6394046 0.6022891
## Towels-cardboards   0.09102564 -0.4470438 0.6290951 0.9142759
## Towels-Control     -0.09285714 -0.6412502 0.4555360 0.9141258
# Extract the Tukey HSD results for yard
tukey_yard <- TukeyHSD(anova_result, "yard")
print(tukey_yard)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = bees_strength ~ treatment * yard, data = filtered_df)
## 
## $yard
##                     diff        lwr       upr     p adj
## Mujea-Biogaz  0.29657380 -0.1749200 0.7680676 0.2959241
## Zeytim-Biogaz 0.31381108 -0.2049253 0.8325474 0.3237290
## Zeytim-Mujea  0.01723728 -0.4920975 0.5265720 0.9964149
# Extract the Tukey HSD results for treatment:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "treatment:yard")
#print(tukey_interaction)

Conclusions, colony strength

  • No difference in colony strength between the treatments, and between the yards.

Queen problems?

# Function to replace empty strings with NA
replace_empty <- function(x) {
  if (is.character(x)) {
    ifelse(x == "", NA, x)
  } else {
    x
  }
}

df %>% 
    dplyr::filter(date ==  "May7th") %>%
mutate(across(everything(), ~replace_empty(.))) %>% 
    filter(!is.na(Queen_problems))%>% 
  ggplot(aes(x = treatment, fill = factor(Queen_problems), y = ..count..)) +
  geom_bar(position = "dodge") +
  facet_wrap(~yard) +
  labs(title = "Queen condition ",
       x = "Treatment",
       y = "Hive count",
       fill = "Queen_problems") +
 # theme_minimal() +
  theme(legend.position = "bottom")

Conclusions, Queen problems

  • No difference in the number of hives with queen problems between the treatments, and between the yards.

Varroa infestation level

Ethanol shake

p_varroa_bio = df %>% 
  filter(yard == "Biogaz") %>%
  ggplot( aes(x = date, y = varroa,  text = paste("Hive:", hive, "<br>Varroa:", varroa))) +
  geom_boxplot() + 
  geom_jitter(width = 0.2) +
  theme(legend.position = "none") +
  ggtitle("Varroa level, Biogaz") +
  facet_wrap(~ treatment) 

p_varroa_zyt = df %>% 
  filter(yard == "Zeytim") %>%
  ggplot( aes(x = date, y = varroa,  text = paste("Hive:", hive, "<br>Varroa:", varroa))) +
  geom_boxplot() + 
  geom_jitter(width = 0.2) +
  theme(legend.position = "none") +
  ggtitle("Varroa level, Zeytim") +
  facet_wrap(~ treatment)

p_varroa_muj = df %>% 
  filter(yard == "Mujea") %>%
  ggplot( aes(x = date, y = varroa,  text = paste("Hive:", hive, "<br>Varroa:", varroa))) +
  geom_boxplot() + 
  geom_jitter(width = 0.2) +
  theme(legend.position = "none") +
  ggtitle("Varroa level, Mujea") +
  facet_wrap(~ treatment)

# Convert the ggplot object to a plotly object and add hover information
p_varroa_bio_plotly <- ggplotly(p_varroa_bio, tooltip = "text")
p_varroa_zyt_plotly <- ggplotly(p_varroa_zyt, tooltip = "text")
p_varroa_muj_plotly <- ggplotly(p_varroa_muj, tooltip = "text")

# Print the plotly object
p_varroa_bio_plotly
p_varroa_zyt_plotly
p_varroa_muj_plotly

An increase in mite levels was observed in all yards. and in all treatments.

Among the 3 treatments, the lowest increase was observed in the towels treatment

 df %>%
  group_by(treatment, date) %>%
  summarize(mean = mean(varroa, na.rm = TRUE), sd = sd(varroa, na.rm = TRUE)) 
## # A tibble: 6 × 4
## # Groups:   treatment [3]
##   treatment  date       mean    sd
##   <chr>      <fct>     <dbl> <dbl>
## 1 Control    March10th 0.128 0.339
## 2 Control    May7th    3.72  4.51 
## 3 Towels     March10th 0.308 0.838
## 4 Towels     May7th    1.29  1.65 
## 5 cardboards March10th 1.53  6.20 
## 6 cardboards May7th    5.69  6.62

Among the 3 yards, the lowest increase was observed in Zeytim

 df %>%
  group_by(yard, date) %>%
  summarize(mean = mean(varroa, na.rm = TRUE), sd = sd(varroa, na.rm = TRUE)) 
## # A tibble: 6 × 4
## # Groups:   yard [3]
##   yard   date       mean    sd
##   <chr>  <fct>     <dbl> <dbl>
## 1 Biogaz March10th 2.08  7.10 
## 2 Biogaz May7th    5.19  5.46 
## 3 Mujea  March10th 0.179 0.389
## 4 Mujea  May7th    4.59  6.33 
## 5 Zeytim March10th 0.105 0.311
## 6 Zeytim May7th    1.52  2.10

Is there a difference between treatments in varroa levels? post treatment

filtered_df <- df %>%
  filter(date == "May7th") 

# Run the ANOVA
anova_result <- aov(varroa ~ treatment * yard, data = filtered_df) %>%
  na.omit()

# Summarize the ANOVA results
summary(anova_result)
##                Df Sum Sq Mean Sq F value  Pr(>F)   
## treatment       2  277.4  138.68   5.635 0.00493 **
## yard            2  182.6   91.31   3.710 0.02825 * 
## treatment:yard  3  140.2   46.72   1.898 0.13544   
## Residuals      91 2239.7   24.61                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 19 observations deleted due to missingness
# Diagnostic plots for ANOVA
#par(mfrow = c(2, 2)) # Set up the plotting area for 2x2 plots
#plot(anova_result)

# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for treatment
tukey_treatment <- TukeyHSD(anova_result, "treatment")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ treatment * yard, data = filtered_df)
## 
## $treatment
##                         diff       lwr        upr     p adj
## Control-cardboards -1.968254 -4.653024  0.7165156 0.1936001
## Towels-cardboards  -4.404762 -7.563922 -1.2456021 0.0036523
## Towels-Control     -2.436508 -5.682236  0.8092202 0.1791334
# Extract the Tukey HSD results for yard
tukey_yard <- TukeyHSD(anova_result, "yard")
print(tukey_yard)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ treatment * yard, data = filtered_df)
## 
## $yard
##                    diff       lwr        upr     p adj
## Mujea-Biogaz  -1.544402 -4.292610  1.2038074 0.3774609
## Zeytim-Biogaz -3.492896 -6.553172 -0.4326198 0.0211735
## Zeytim-Mujea  -1.948494 -5.008770  1.1117817 0.2877416
# Extract the Tukey HSD results for treatment:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "treatment:yard")
#print(tukey_interaction)
  • Yes, the most significant difference in varroa level were found between towels and cardboards treatment, in Mujea and in Biogaz:
    Towels:Zeytim-cardboards:Biogaz p adj = 0.055
    Towels:Zeytim-cardboards:Mujea p adj = 0.04.
    —

Is there a difference between dates in varroa levels, in the different treatments?

Control
control_df <- df %>%
    filter(hive != "94") %>%
    filter(treatment == "Control")

# Run the ANOVA
anova_result <- aov(varroa ~ date * yard, data = control_df) 

# Summarize the ANOVA results
summary(anova_result)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## date         1  241.8  241.81  27.424 1.68e-06 ***
## yard         2   57.1   28.54   3.236   0.0453 *  
## date:yard    2   52.1   26.06   2.956   0.0587 .  
## Residuals   69  608.4    8.82                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 5 observations deleted due to missingness
# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for yard
tukey_treatment <- TukeyHSD(anova_result, "yard")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = control_df)
## 
## $yard
##                     diff       lwr         upr     p adj
## Mujea-Biogaz  -2.0605263 -4.025395 -0.09565809 0.0377057
## Zeytim-Biogaz -1.6799397 -4.026262  0.66638244 0.2068947
## Zeytim-Mujea   0.3805866 -1.694763  2.45593649 0.8992700
# Extract the Tukey HSD results for date
tukey_date <- TukeyHSD(anova_result, "date")
print(tukey_date)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = control_df)
## 
## $date
##                      diff      lwr      upr   p adj
## May7th-March10th 3.594017 2.224895 4.963139 1.7e-06
# Extract the Tukey HSD results for date:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "date:yard")
print(tukey_interaction)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = control_df)
## 
## $`date:yard`
##                                           diff         lwr        upr     p adj
## May7th:Biogaz-March10th:Biogaz     6.300000000   2.4073772 10.1926228 0.0001547
## March10th:Mujea-March10th:Biogaz  -0.094736842  -3.4952897  3.3058160 0.9999995
## May7th:Mujea-March10th:Biogaz      2.273684211  -1.1268686  5.6742370 0.3756977
## March10th:Zeytim-March10th:Biogaz -0.100000000  -3.9926228  3.7926228 0.9999996
## May7th:Zeytim-March10th:Biogaz     2.942857143  -1.3466043  7.2323186 0.3468080
## March10th:Mujea-May7th:Biogaz     -6.394736842  -9.7952897 -2.9941840 0.0000083
## May7th:Mujea-May7th:Biogaz        -4.026315789  -7.4268686 -0.6257630 0.0111539
## March10th:Zeytim-May7th:Biogaz    -6.400000000 -10.2926228 -2.5073772 0.0001172
## May7th:Zeytim-May7th:Biogaz       -3.357142857  -7.6466043  0.9323186 0.2105104
## May7th:Mujea-March10th:Mujea       2.368421053  -0.4555844  5.1924265 0.1511725
## March10th:Zeytim-March10th:Mujea  -0.005263158  -3.4058160  3.3952897 1.0000000
## May7th:Zeytim-March10th:Mujea      3.037593985  -0.8108767  6.8860647 0.2027099
## March10th:Zeytim-May7th:Mujea     -2.373684211  -5.7742370  1.0268686 0.3276898
## May7th:Zeytim-May7th:Mujea         0.669172932  -3.1792977  4.5176436 0.9956398
## May7th:Zeytim-March10th:Zeytim     3.042857143  -1.2466043  7.3323186 0.3100994
  • Yes, in the control treatment, an increase in mite level was obsereved in all yards. the increase was significant in Biogaz yard.

May7th:Biogaz-March10th:Biogaz p adj = 0.0001
May7th:Mujea-March10th:Mujea p adj = 0.15
May7th:Zeytim-March10th:Zeytim p adj = 0.310

Cardboards
card_df <- df %>%
    filter(hive != "94") %>%
    filter(treatment == "cardboards")

# Run the ANOVA
anova_result <- aov(varroa ~ date * yard, data = card_df) 

# Summarize the ANOVA results
summary(anova_result)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## date         1  501.5   501.5  21.618 1.24e-05 ***
## yard         2  132.4    66.2   2.853   0.0633 .  
## date:yard    2   93.5    46.7   2.015   0.1398    
## Residuals   83 1925.4    23.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness
# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for yard
tukey_treatment <- TukeyHSD(anova_result, "yard")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = card_df)
## 
## $yard
##                     diff       lwr       upr     p adj
## Mujea-Biogaz  -0.8254766 -3.718554 2.0676006 0.7752302
## Zeytim-Biogaz -3.0989442 -6.323535 0.1256466 0.0622736
## Zeytim-Mujea  -2.2734676 -5.270390 0.7234545 0.1724837
# Extract the Tukey HSD results for date
tukey_date <- TukeyHSD(anova_result, "date")
print(tukey_date)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = card_df)
## 
## $date
##                      diff      lwr      upr    p adj
## May7th-March10th 4.762195 2.725019 6.799372 1.24e-05
# Extract the Tukey HSD results for date:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "date:yard")
print(tukey_interaction)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = card_df)
## 
## $`date:yard`
##                                         diff         lwr        upr     p adj
## May7th:Biogaz-March10th:Biogaz     4.0494505  -1.3625343  9.4614354 0.2568432
## March10th:Mujea-March10th:Biogaz  -2.0576923  -7.0635752  2.9481906 0.8360403
## May7th:Mujea-March10th:Biogaz      4.5256410  -0.5886249  9.6399069 0.1134101
## March10th:Zeytim-March10th:Biogaz -2.2410256  -7.5654431  3.0833918 0.8219704
## May7th:Zeytim-March10th:Biogaz    -0.7521368  -6.8451033  5.3408298 0.9991781
## March10th:Mujea-May7th:Biogaz     -6.1071429 -11.0034745 -1.2108112 0.0061230
## May7th:Mujea-May7th:Biogaz         0.4761905  -4.5308962  5.4832771 0.9997700
## March10th:Zeytim-May7th:Biogaz    -6.2904762 -11.5120298 -1.0689226 0.0090525
## May7th:Zeytim-May7th:Biogaz       -4.8015873 -10.8048733  1.2016987 0.1926688
## May7th:Mujea-March10th:Mujea       6.5833333   2.0182266 11.1484401 0.0008962
## March10th:Zeytim-March10th:Mujea  -0.1833333  -4.9826984  4.6160317 0.9999975
## May7th:Zeytim-March10th:Mujea      1.3055556  -4.3343669  6.9454780 0.9841968
## March10th:Zeytim-May7th:Mujea     -6.7666667 -11.6789734 -1.8543600 0.0017337
## May7th:Zeytim-May7th:Mujea        -5.2777778 -11.0141162  0.4585606 0.0893717
## May7th:Zeytim-March10th:Zeytim     1.4888889  -4.4355759  7.4133537 0.9772355
  • Yes, in the cardboards treatment, an increase in mite level was observed in all yards. the increase was significant in Mujea yard.

May7th:Biogaz-March10th:Biogaz p adj = 0.26
May7th:Mujea-March10th:Mujea p adj = 0.0009 ** significant increase in mite level after treatment with cardboards!
May7th:Zeytim-March10th:Zeytim p adj = 0.98

Towels
towel_df <- df %>%
    filter(treatment == "Towels")

# Run the ANOVA
anova_result <- aov(varroa ~ date * yard, data = towel_df) 

# Summarize the ANOVA results
summary(anova_result)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## date         1  11.11  11.112   8.949 0.00458 **
## yard         1  11.50  11.500   9.261 0.00398 **
## date:yard    1   6.93   6.928   5.580 0.02277 * 
## Residuals   43  53.40   1.242                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for yard
tukey_treatment <- TukeyHSD(anova_result, "yard")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = towel_df)
## 
## $yard
##                     diff       lwr       upr     p adj
## Zeytim-Biogaz -0.9888312 -1.645769 -0.331893 0.0040662
# Extract the Tukey HSD results for date
tukey_date <- TukeyHSD(anova_result, "date")
print(tukey_date)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = towel_df)
## 
## $date
##                      diff       lwr      upr     p adj
## May7th-March10th 0.978022 0.3186818 1.637362 0.0045834
# Extract the Tukey HSD results for date:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "date:yard")
print(tukey_interaction)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ date * yard, data = towel_df)
## 
## $`date:yard`
##                                          diff        lwr        upr     p adj
## May7th:Biogaz-March10th:Biogaz     1.62179487  0.4296472  2.8139426 0.0039602
## March10th:Zeytim-March10th:Biogaz -0.30769231 -1.4757537  0.8603691 0.8949600
## May7th:Zeytim-March10th:Biogaz    -0.23931624 -1.5306567  1.0520242 0.9597192
## March10th:Zeytim-May7th:Biogaz    -1.92948718 -3.1216349 -0.7373395 0.0005000
## May7th:Zeytim-May7th:Biogaz       -1.86111111 -3.1742786 -0.5479436 0.0025465
## May7th:Zeytim-March10th:Zeytim     0.06837607 -1.2229644  1.3597166 0.9989730
  • Increase in mite level in Biogaz yard only
    May7th:Biogaz-March10th:Biogaz p adj = 0.0039602.
    May7th:Zeytim-March10th:Zeytim p adj = 0.9989730.

check that the varroa levels were similar in the begining of the experiments, prior treatment.

March_df <- df %>%
      filter(date == "March10th")

# Run the ANOVA
anova_result <- aov(varroa ~ treatment * yard, data = March_df) 

# Summarize the ANOVA results
summary(anova_result)
##                 Df Sum Sq Mean Sq F value Pr(>F)  
## treatment        2   50.0   24.98   1.616 0.2035  
## yard             2  106.5   53.26   3.445 0.0355 *
## treatment:yard   3  120.8   40.25   2.604 0.0557 .
## Residuals      106 1638.8   15.46                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 4 observations deleted due to missingness
# Perform post hoc tests using Tukey's Honest Significant Difference
tukey_result <- TukeyHSD(anova_result)

# Extract the Tukey HSD results for treatment
tukey_treatment <- TukeyHSD(anova_result, "treatment")
print(tukey_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ treatment * yard, data = March_df)
## 
## $treatment
##                          diff       lwr       upr     p adj
## Control-cardboards -1.4024071 -3.408120 0.6033056 0.2246295
## Towels-cardboards  -1.2229199 -3.490714 1.0448744 0.4085143
## Towels-Control      0.1794872 -2.186953 2.5459269 0.9822423
# Extract the Tukey HSD results for yard
tukey_date <- TukeyHSD(anova_result, "yard")
print(tukey_date)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = varroa ~ treatment * yard, data = March_df)
## 
## $yard
##                      diff       lwr       upr     p adj
## Mujea-Biogaz  -2.02707344 -4.172094 0.1179467 0.0681542
## Zeytim-Biogaz -1.99709959 -4.155817 0.1616174 0.0759453
## Zeytim-Mujea   0.02997385 -2.100514 2.1604614 0.9993836
# Extract the Tukey HSD results for date:yard interaction
tukey_interaction <- TukeyHSD(anova_result, "treatment:yard")
#print(tukey_interaction)
  • no significant differnece in varroa level, between the hives prior treatment