Greenhouse Data Analysis

6 Point Data Analysis

gh = factor, 6 unique IDs, each greenhouse coded as an individual and seperate entity source = koppert or biobest 6 data points for all analyses averaged across greenhouse

input greenhouse data

greenhouse_analysis <- read_csv("greenhouse_analysis.csv", 
    col_types = cols(year = col_factor(levels = c("2022", 
        "2024")), gh = col_factor(levels = c("1", 
        "2", "3", "4", "5", "6"))))

greenhouse_analysis <- greenhouse_analysis %>%
  mutate(across(c(
    coragen_L, lalstop_k61_L, previcurflex_L, pylon_L, lunatranquility_L,
    rootshielplus_L, milstop_L, quadristop_L, azaguard_L, botanigard22wp_L,
    botanigardes_L, captivaprime_L, fontelis_L, mpede_L, nofly_L, veneratecg_L,
    beleaf50sg_L, entrustsc_L, grotto_L
  ), as.logical))

greenhouse_analysis$source <- as.factor(greenhouse_analysis$source)

look at general data structure

library(ggplot2)
library(tidyr)
library(dplyr)

# List of your logical treatment variables
treat_vars <- c(
  "coragen_L", "lalstop_k61_L", "previcurflex_L", "pylon_L", 
  "lunatranquility_L", "rootshielplus_L", "milstop_L", "quadristop_L", 
  "azaguard_L", "botanigard22wp_L", "botanigardes_L", "captivaprime_L", 
  "fontelis_L", "mpede_L", "nofly_L", "veneratecg_L", "beleaf50sg_L", 
  "entrustsc_L", "grotto_L"
)

# Reshape the data to long format
greenhouse_long <- greenhouse_analysis %>%
  pivot_longer(cols = all_of(treat_vars),
               names_to = "treatment",
               values_to = "applied")

# Create boxplots of crith_prop vs each treatment
ggplot(greenhouse_long, aes(x = as.factor(applied), y = crith_prop)) +
  geom_boxplot(fill = "gray90", color = "black", outlier.shape = 21) +
  facet_wrap(~ treatment, scales = "free_y", ncol = 4) +
  labs(
    x = "Applied (TRUE = Treated, FALSE = Control)",
    y = "Crithidia infection proportion",
    title = "Effect of each treatment on Crithidia infection"
  ) +
  theme_bw(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 0, hjust = 0.5)
  )

ggplot(greenhouse_long, aes(x = as.factor(applied), y = mean_pol)) +
  geom_boxplot(fill = "gray90", color = "black", outlier.shape = 21) +
  facet_wrap(~ treatment, scales = "free_y", ncol = 4) +
  labs(
    x = "Applied (TRUE = Treated, FALSE = Control)",
    y = "Mean anther bruising",
    title = "Effect of each treatment on mean anther bruising"
  ) +
  theme_bw(base_size = 14) +
  theme(
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 0, hjust = 0.5)
  )

ggplot(greenhouse_analysis, aes(x=crith_prop, y = prop_1)) +
  geom_point(size = 3) + 
  geom_smooth(method = "glm", se = TRUE, size = 2)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `geom_smooth()` using formula = 'y ~ x'

ggplot(greenhouse_analysis, aes(x=total_application, y = prop_1)) +
  geom_point(size = 3) + 
  geom_smooth(method = "glm", se = TRUE, size = 2)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(greenhouse_analysis, aes(x=total_application, y = prop_2)) +
  geom_point(size = 3) + 
  geom_smooth(method = "glm", se = TRUE, size = 2)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(greenhouse_analysis, aes(x=total_application, y = prop_3)) +
  geom_point(size = 3) + 
  geom_smooth(method = "glm", se = TRUE, size = 2)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(greenhouse_analysis, aes(x=total_application, y = mean_pol)) +
  geom_point(size = 3) + 
  geom_smooth(method = "glm", se = TRUE, size = 2)

## `geom_smooth()` using formula = 'y ~ x'

ggplot(greenhouse_analysis, aes(x=source, y = mean_pol, fill = source)) +
  geom_boxplot()

ggplot(greenhouse_analysis, aes(x=source, y = crith_prop, fill = source)) +
  geom_boxplot()

Crithidia vs. individual applied pesticide presence/absence

library(ggplot2)
library(dplyr)

# Summarize the data
greenhouse_summary <- greenhouse_long %>%
  group_by(applied, treatment) %>%
  summarise(
    mean_crith = mean(crith_prop, na.rm = TRUE),
    se_crith = sd(crith_prop, na.rm = TRUE)/sqrt(n()),
    .groups = "drop"
  )

# Bar plot
ggplot(greenhouse_summary, aes(x = treatment, y = mean_crith, fill = applied)) +
  geom_col(position = position_dodge(width = 0.8), color = "black") +
  geom_errorbar(aes(ymin = mean_crith - se_crith, ymax = mean_crith + se_crith),
                width = 0.2, position = position_dodge(width = 0.8)) +
  labs(
    x = "Pesticide",
    y = "Mean Crithidia infection proportion",
    title = "Effect of each treatment on Crithidia infection"
  ) +
  scale_fill_manual(values = c("FALSE" = "gray70", "TRUE" = "skyblue"), 
                    name = "Applied") +
  theme_bw(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    legend.position = "top"
  )

Pollination analysis

Pollination score vs. Crithidia

library(ggplot2)
library(dplyr)
library(lme4)
library(cowplot)

# Fit model
crith_mod1 <- glm(mean_pol ~ crith_prop + source + year, data = greenhouse_analysis)
Anova(crith_mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: mean_pol
##            LR Chisq Df Pr(>Chisq)  
## crith_prop   2.9585  1    0.08543 .
## source       0.6591  1    0.41688  
## year         6.2548  1    0.01239 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(crith_mod1)

## 
## Call:
## glm(formula = mean_pol ~ crith_prop + source + year, data = greenhouse_analysis)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  2.27761    0.12593  18.086  0.00304 **
## crith_prop  -0.27087    0.15748  -1.720  0.22757   
## sourcekop   -0.06171    0.07602  -0.812  0.50214   
## year2024     0.15049    0.06017   2.501  0.12953   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.003709939)
## 
##     Null deviance: 0.0851667  on 5  degrees of freedom
## Residual deviance: 0.0074199  on 2  degrees of freedom
## AIC: -13.145
## 
## Number of Fisher Scoring iterations: 2

AIC(crith_mod1)

## [1] -13.14485

sim_resgh_cm <- simulateResiduals(fittedModel = crith_mod1)
plot(sim_resgh_cm)

testDispersion(crith_mod1)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.40046, p-value = 0.296
## alternative hypothesis: two.sided

# Create a grid of crith_prop values
crith_grid <- data.frame(
  crith_prop = seq(min(greenhouse_analysis$crith_prop, na.rm = TRUE),
                   max(greenhouse_analysis$crith_prop, na.rm = TRUE),
                   length.out = 100),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year   = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict with standard errors (fixed effects only)
pred <- predict(crith_mod1, newdata = crith_grid, re.form = NA, se.fit = TRUE)
crith_grid$predicted <- pred$fit
crith_grid$lower <- pred$fit - 1.96 * pred$se.fit
crith_grid$upper <- pred$fit + 1.96 * pred$se.fit

# Plot with CI ribbon
ggplot(greenhouse_analysis, aes(x = crith_prop, y = mean_pol)) +
  geom_point(shape = 21, size = 6, fill = "darkgreen", color = "black") +
  geom_ribbon(data = crith_grid,
              aes(x = crith_prop, ymin = lower, ymax = upper),
              alpha = 0.2, fill = "darkgray",
              inherit.aes = FALSE) +  # <-- prevent inheriting y = mean_pol
  geom_line(data = crith_grid,
            aes(x = crith_prop, y = predicted),
            color = "black", size = 2,
            inherit.aes = FALSE) +
  theme_cowplot() +
  labs(
    y = "Average Anther Bruising",
    x = "Average Proportion Crithidia"
  ) +
  annotate(
    geom = "text",
    x = 0.2,
    y = 2.5,
    label = "P = 0.08",
    size = 12
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  )

Pollination vs. total chemicals applied

library(ggplot2)
library(dplyr)
library(lme4)
library(cowplot)

# Fit the model
total_applied_mod1 <- glm(mean_pol ~ total_applied + source + year, data = greenhouse_analysis)
Anova(total_applied_mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: mean_pol
##               LR Chisq Df Pr(>Chisq)   
## total_applied  10.1418  1   0.001449 **
## source          3.0744  1   0.079532 . 
## year            0.0212  1   0.884245   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(total_applied_mod1)

## 
## Call:
## glm(formula = mean_pol ~ total_applied + source + year, data = greenhouse_analysis)
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    2.439e+00  1.181e-01  20.652  0.00234 **
## total_applied -2.379e-06  7.471e-07  -3.185  0.08606 . 
## sourcekop     -8.941e-02  5.099e-02  -1.753  0.22162   
## year2024       1.019e-02  6.998e-02   0.146  0.89759   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.001515072)
## 
##     Null deviance: 0.0851667  on 5  degrees of freedom
## Residual deviance: 0.0030301  on 2  degrees of freedom
## AIC: -18.518
## 
## Number of Fisher Scoring iterations: 2

sim_resgh_tt <- simulateResiduals(fittedModel = total_applied_mod1)
plot(sim_resgh_tt)

testDispersion(total_applied_mod1)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.40046, p-value = 0.296
## alternative hypothesis: two.sided

# Create a grid of total_applied values
total_grid <- data.frame(
  total_applied = seq(min(greenhouse_analysis$total_applied, na.rm = TRUE),
                      max(greenhouse_analysis$total_applied, na.rm = TRUE),
                      length.out = 100),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year   = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict with standard errors (fixed effects only)
pred_total <- predict(total_applied_mod1, newdata = total_grid, re.form = NA, se.fit = TRUE)
total_grid$predicted <- pred_total$fit
total_grid$lower <- pred_total$fit - 1.96 * pred_total$se.fit
total_grid$upper <- pred_total$fit + 1.96 * pred_total$se.fit

# Plot with CI ribbon
ggplot(greenhouse_analysis, aes(x = total_applied, y = mean_pol)) +
  geom_point(shape = 21, size = 6, fill = "darkgreen", color = "black") +
  geom_ribbon(data = total_grid,
              aes(x = total_applied, ymin = lower, ymax = upper),
              alpha = 0.2, fill = "seagreen",
              inherit.aes = FALSE) +
  geom_line(data = total_grid,
            aes(x = total_applied, y = predicted),
            color = "black", size = 2,
            inherit.aes = FALSE) +
  theme_cowplot() +
  labs(
    y = "Average Anther Bruising",
    x = "Total Pesticides Applied (ml)"
  ) +
  annotate(
    geom = "text",
    x = min(greenhouse_analysis$total_applied) + 0.1 * diff(range(greenhouse_analysis$total_applied)),
    y = 2.7,
    label = "P < 0.01",  # Replace with your actual P-value if you want
    size = 12
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  )

Pollination vs. fungicide

total_fung_mod1 <- glm(mean_pol ~ fungicide_b + source + year, data = greenhouse_analysis)
Anova(total_fung_mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: mean_pol
##             LR Chisq Df Pr(>Chisq)    
## fungicide_b   28.224  1  1.081e-07 ***
## source         0.030  1     0.8627    
## year         133.624  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

total_fung_mod1

## 
## Call:  glm(formula = mean_pol ~ fungicide_b + source + year, data = greenhouse_analysis)
## 
## Coefficients:
## (Intercept)  fungicide_b    sourcekop     year2024  
##   2.134e+00   -2.097e-06    3.812e-03    2.481e-01  
## 
## Degrees of Freedom: 5 Total (i.e. Null);  2 Residual
## Null Deviance:       0.08517 
## Residual Deviance: 0.001217  AIC: -23.99

summary(total_fung_mod1)

## 
## Call:
## glm(formula = mean_pol ~ fungicide_b + source + year, data = greenhouse_analysis)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.134e+00  1.973e-02 108.139 8.55e-05 ***
## fungicide_b -2.097e-06  3.948e-07  -5.313   0.0337 *  
## sourcekop    3.812e-03  2.203e-02   0.173   0.8786    
## year2024     2.481e-01  2.146e-02  11.560   0.0074 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0006086432)
## 
##     Null deviance: 0.0851667  on 5  degrees of freedom
## Residual deviance: 0.0012173  on 2  degrees of freedom
## AIC: -23.99
## 
## Number of Fisher Scoring iterations: 2

AIC(total_fung_mod1)

## [1] -23.99008

sim_resgh_fung <- simulateResiduals(fittedModel = total_fung_mod1)
qqnorm(resid(total_fung_mod1));qqline(resid(total_fung_mod1))

testUniformity(sim_resgh_fung)

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.236, p-value = 0.8238
## alternative hypothesis: two-sided

testDispersion(sim_resgh_fung)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.40046, p-value = 0.296
## alternative hypothesis: two.sided

# Create a grid of total_applied values
total_grid <- data.frame(
  fungicide_b = seq(min(greenhouse_analysis$fungicide_b, na.rm = TRUE),
                      max(greenhouse_analysis$fungicide_b, na.rm = TRUE),
                      length.out = 100),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year   = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict with standard errors (fixed effects only)
pred_total <- predict(total_fung_mod1, newdata = total_grid, re.form = NA, se.fit = TRUE)
total_grid$predicted <- pred_total$fit
total_grid$lower <- pred_total$fit - 1.96 * pred_total$se.fit
total_grid$upper <- pred_total$fit + 1.96 * pred_total$se.fit

# Plot with CI ribbon
ggplot(greenhouse_analysis, aes(x = fungicide_b, y = mean_pol)) +
  geom_point(shape = 21, size = 6, fill = "darkgreen", color = "black") +
  geom_ribbon(data = total_grid,
              aes(x = fungicide_b, ymin = lower, ymax = upper),
              alpha = 0.2, fill = "seagreen",
              inherit.aes = FALSE) +
  geom_line(data = total_grid,
            aes(x = fungicide_b, y = predicted),
            color = "black", size = 2,
            inherit.aes = FALSE) +
  theme_cowplot() +
  labs(
    y = "Average Anther Bruising",
    x = "Total Fungicides Applied (ml)"
  ) +
  annotate(
    geom = "text",
    x = 10000,
    y = 2.6,
    label = "P < 0.01",  # Replace with your actual P-value if you want
    size = 9
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  )

Pollination vs. insecticide

total_insec_mod1 <- glm(mean_pol ~ insecticide_1 + source + year, data = greenhouse_analysis)
Anova(total_insec_mod1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: mean_pol
##               LR Chisq Df Pr(>Chisq)
## insecticide_1  0.09021  1     0.7639
## source         0.24497  1     0.6206
## year           1.15623  1     0.2822

total_insec_mod1

## 
## Call:  glm(formula = mean_pol ~ insecticide_1 + source + year, data = greenhouse_analysis)
## 
## Coefficients:
##   (Intercept)  insecticide_1      sourcekop       year2024  
##     1.985e+00      6.986e-07      5.665e-02      2.849e-01  
## 
## Degrees of Freedom: 5 Total (i.e. Null);  2 Residual
## Null Deviance:       0.08517 
## Residual Deviance: 0.0176    AIC: -7.962

summary(total_insec_mod1)

## 
## Call:
## glm(formula = mean_pol ~ insecticide_1 + source + year, data = greenhouse_analysis)
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   1.985e+00  2.947e-01   6.736   0.0213 *
## insecticide_1 6.986e-07  2.326e-06   0.300   0.7923  
## sourcekop     5.665e-02  1.145e-01   0.495   0.6697  
## year2024      2.849e-01  2.650e-01   1.075   0.3947  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.008800879)
## 
##     Null deviance: 0.085167  on 5  degrees of freedom
## Residual deviance: 0.017602  on 2  degrees of freedom
## AIC: -7.9618
## 
## Number of Fisher Scoring iterations: 2

AIC(total_insec_mod1)

## [1] -7.961834

sim_resgh_insec <- simulateResiduals(fittedModel = total_insec_mod1)
qqnorm(resid(total_insec_mod1));qqline(resid(total_insec_mod1))

testUniformity(sim_resgh_insec)

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.188, p-value = 0.9563
## alternative hypothesis: two-sided

testDispersion(sim_resgh_insec)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.40046, p-value = 0.296
## alternative hypothesis: two.sided

# Create a grid of insecticide_1 values for prediction
total_grid_insec <- data.frame(
  insecticide_1 = seq(
    min(greenhouse_analysis$insecticide_1, na.rm = TRUE),
    max(greenhouse_analysis$insecticide_1, na.rm = TRUE),
    length.out = 100
  ),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year   = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict with standard errors (fixed effects only)
pred_total_insec <- predict(total_insec_mod1, newdata = total_grid_insec, se.fit = TRUE)
total_grid_insec$predicted <- pred_total_insec$fit
total_grid_insec$lower <- pred_total_insec$fit - 1.96 * pred_total_insec$se.fit
total_grid_insec$upper <- pred_total_insec$fit + 1.96 * pred_total_insec$se.fit

# Plot model with CI ribbon
ggplot(greenhouse_analysis, aes(x = insecticide_1, y = mean_pol)) +
  geom_point(shape = 21, size = 6, fill = "darkgreen", color = "black") +
  geom_ribbon(
    data = total_grid_insec,
    aes(x = insecticide_1, ymin = lower, ymax = upper),
    alpha = 0.2, fill = "darkgray",
    inherit.aes = FALSE
  ) +
  geom_line(
    data = total_grid_insec,
    aes(x = insecticide_1, y = predicted),
    color = "black", size = 2,
    inherit.aes = FALSE
  ) +
  scale_x_continuous(labels = scales::comma_format(accuracy = 1)) +  # avoid scientific notation
  theme_cowplot() +
  labs(
    y = "Average Anther Bruising",
    x = "Total Insecticides Applied (ml)"
  ) +
  annotate(
    geom = "text",
    x = 17000,
    y = 2.6,
    label = "P = 0.8",  # adjust if desired
    size = 10
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  )

Crithidia analysis

Crithidia vs. total chemicals applied

# Fit the binomial GLM using plain numeric total_applied
gh_ins_crith <- glm(
  cbind(total_crith, crith_neg) ~ total_applied + year + source,
  data = greenhouse_analysis,
  family = binomial("logit")
)

Anova(gh_ins_crith)

## Analysis of Deviance Table (Type II tests)
## 
## Response: cbind(total_crith, crith_neg)
##               LR Chisq Df Pr(>Chisq)    
## total_applied   49.182  1  2.332e-12 ***
## year            11.915  1  0.0005569 ***
## source           0.263  1  0.6080041    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Create prediction grid
total_grid <- data.frame(
  total_applied = seq(
    min(greenhouse_analysis$total_applied, na.rm = TRUE),
    max(greenhouse_analysis$total_applied, na.rm = TRUE),
    length.out = 100
  ),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict on the link scale with SE
pred <- predict(
  gh_ins_crith,
  newdata = total_grid,
  type = "link",
  se.fit = TRUE
)

# Compute 95% CI on link scale and transform back to probability
inv_logit <- function(x) exp(x) / (1 + exp(x))
total_grid <- total_grid %>%
  mutate(
    fit_link = pred$fit,
    lower_link = fit_link - 1.96 * pred$se.fit,
    upper_link = fit_link + 1.96 * pred$se.fit,
    predicted = inv_logit(fit_link),
    lower = inv_logit(lower_link),
    upper = inv_logit(upper_link)
  )

# Plot
ggplot(greenhouse_analysis, aes(x = total_applied, y = total_crith / (total_crith + crith_neg))) +
  geom_point(shape = 21, size = 5, fill = "darkorchid", color = "black") +
  geom_ribbon(data = total_grid,
              aes(x = total_applied, ymin = lower, ymax = upper),
              alpha = 0.2, fill = "orchid",
              inherit.aes = FALSE) +
  geom_line(data = total_grid,
            aes(x = total_applied, y = predicted),
            color = "black", size = 1.5,
            inherit.aes = FALSE) +
  coord_cartesian(ylim = c(0, 1)) +  # set y-axis from 0 to 1
  theme_cowplot() +
  labs(
    y = "Average Crithidia Prevelance",
    x = "Total Pesticides Applied (ml)"
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  ) +
  annotate(geom = "text", 
           x = 45000,
           y = 1,
           label = "P < 0.001",
           size = 12)

Crithidia vs. fungicide

library(ggplot2)
library(dplyr)
library(cowplot)

greenhouse_analysis <- greenhouse_analysis %>%
  mutate(
    fungicide_1_sc = scale(fungicide_1),   # rescaled version
    # keep source and year as factors
    source = factor(source),
    year = factor(year)
  )

# Refit model (from before)
gh_ins_crith1_rescaled <- glmer(
  cbind(total_crith, crith_neg) ~ fungicide_1_sc + source + year + (1 | gh),
  data = greenhouse_analysis,
  family = binomial("logit"),
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)

simres <- simulateResiduals(gh_ins_crith1_rescaled)
testDispersion(simres)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98085, p-value = 0.952
## alternative hypothesis: two.sided

# Create prediction grid using the scaled variable
total_grid <- data.frame(
  fungicide_1_sc = seq(
    min(greenhouse_analysis$fungicide_1_sc, na.rm = TRUE),
    max(greenhouse_analysis$fungicide_1_sc, na.rm = TRUE),
    length.out = 100
  ),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict on link scale with SE
pred <- predict(
  gh_ins_crith1_rescaled,
  newdata = total_grid,
  type = "link",
  se.fit = TRUE,
  re.form = NA  # marginal prediction (average over random effects)
)

## Warning in predict.merMod(gh_ins_crith1_rescaled, newdata = total_grid, :
## se.fit computation uses an approximation to estimate the sampling distribution
## of the parameters

# Back-transform the scaled x-axis for plotting
fung_mean <- attr(greenhouse_analysis$fungicide_1_sc, "scaled:center")
fung_sd <- attr(greenhouse_analysis$fungicide_1_sc, "scaled:scale")
total_grid$fungicide_1 <- total_grid$fungicide_1_sc * fung_sd + fung_mean

# Compute 95% CI on link scale and transform back to probability
inv_logit <- function(x) exp(x) / (1 + exp(x))
total_grid <- total_grid %>%
  mutate(
    fit_link = pred$fit,
    lower_link = fit_link - 1.96 * pred$se.fit,
    upper_link = fit_link + 1.96 * pred$se.fit,
    predicted = inv_logit(fit_link),
    lower = inv_logit(lower_link),
    upper = inv_logit(upper_link)
  )

# Plot the model fit with CI
ggplot(greenhouse_analysis, aes(x = fungicide_1, y = total_crith / (total_crith + crith_neg))) +
  geom_point(shape = 21, size = 6, fill = "darkorchid", color = "black") +
  geom_ribbon(
    data = total_grid,
    aes(x = fungicide_1, ymin = lower, ymax = upper),
    alpha = 0.2, fill = "orchid",
    inherit.aes = FALSE
  ) +
  geom_line(
    data = total_grid,
    aes(x = fungicide_1, y = predicted),
    color = "black", size = 2,
    inherit.aes = FALSE
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  theme_cowplot() +
  labs(
    y = "Average Crithidia Prevalence",
    x = "Total Fungicides Applied (ml)"
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  ) +
  annotate(
    geom = "text",
    x = 17000,
    y = 0.95,
    label = "P < 0.001",
    size = 10
  )

Crithidia vs. insecticide

# Rescale and refit model
greenhouse_analysis <- greenhouse_analysis %>%
  mutate(
    insecticide_1_sc = scale(insecticide_1),
    source = factor(source),
    year = factor(year)
  )

# Fit GLMM with scaled predictor
gh_ins_crith1_rescaled <- glmer(
  cbind(total_crith, crith_neg) ~ insecticide_1_sc + source + year + (1 | gh),
  data = greenhouse_analysis,
  family = binomial("logit"),
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
)

summary(gh_ins_crith1_rescaled)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(total_crith, crith_neg) ~ insecticide_1_sc + source + year +  
##     (1 | gh)
##    Data: greenhouse_analysis
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##     60.8     59.7    -25.4     50.8        1 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -0.29961 -0.16038  0.03638  0.15290  0.55952 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  gh     (Intercept) 0.4713   0.6865  
## Number of obs: 6, groups:  gh, 6
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        3.4224     1.2917   2.650  0.00806 **
## insecticide_1_sc  -1.9379     1.1738  -1.651  0.09873 . 
## sourcekop         -2.9814     0.9208  -3.238  0.00120 **
## year2024          -4.5702     2.0754  -2.202  0.02766 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ins_1_ sorckp
## insctcd_1_s -0.933              
## sourcekop   -0.794  0.728       
## year2024    -0.960  0.959  0.706

Anova(gh_ins_crith1_rescaled)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: cbind(total_crith, crith_neg)
##                    Chisq Df Pr(>Chisq)   
## insecticide_1_sc  2.7259  1   0.098730 . 
## source           10.4843  1   0.001204 **
## year              4.8494  1   0.027655 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Model diagnostics
sim_res_gh_in1 <- simulateResiduals(fittedModel = gh_ins_crith1_rescaled)
plot(sim_res_gh_in1)

testDispersion(sim_res_gh_in1)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.012, p-value = 0.896
## alternative hypothesis: two.sided

# Create prediction grid using scaled variable
total_grid <- data.frame(
  insecticide_1_sc = seq(
    min(greenhouse_analysis$insecticide_1_sc, na.rm = TRUE),
    max(greenhouse_analysis$insecticide_1_sc, na.rm = TRUE),
    length.out = 100
  ),
  source = factor(levels(greenhouse_analysis$source)[1], levels = levels(greenhouse_analysis$source)),
  year = factor(levels(greenhouse_analysis$year)[1], levels = levels(greenhouse_analysis$year))
)

# Predict on link scale with SE
pred <- predict(
  gh_ins_crith1_rescaled,
  newdata = total_grid,
  type = "link",
  se.fit = TRUE,
  re.form = NA  # marginal predictions
)

## Warning in predict.merMod(gh_ins_crith1_rescaled, newdata = total_grid, :
## se.fit computation uses an approximation to estimate the sampling distribution
## of the parameters

# Back-transform scaled variable to original units
insect_mean <- attr(greenhouse_analysis$insecticide_1_sc, "scaled:center")
insect_sd <- attr(greenhouse_analysis$insecticide_1_sc, "scaled:scale")
total_grid$insecticide_1 <- total_grid$insecticide_1_sc * insect_sd + insect_mean

# Compute 95% CI on link scale and back-transform to probability
inv_logit <- function(x) exp(x) / (1 + exp(x))
total_grid <- total_grid %>%
  mutate(
    fit_link = pred$fit,
    lower_link = fit_link - 1.96 * pred$se.fit,
    upper_link = fit_link + 1.96 * pred$se.fit,
    predicted = inv_logit(fit_link),
    lower = inv_logit(lower_link),
    upper = inv_logit(upper_link)
  )

# Final plot on ORIGINAL insecticide scale
ggplot(greenhouse_analysis, aes(x = insecticide_1, y = total_crith / (total_crith + crith_neg))) +
  geom_point(shape = 21, size = 6, fill = "darkorchid", color = "black") +
  geom_ribbon(
    data = total_grid,
    aes(x = insecticide_1, ymin = lower, ymax = upper),
    alpha = 0.2, fill = "darkgray",
    inherit.aes = FALSE
  ) +
  geom_line(
    data = total_grid,
    aes(x = insecticide_1, y = predicted),
    color = "black", size = 2,
    inherit.aes = FALSE
  ) +
  scale_x_continuous(
    labels = scales::comma_format(accuracy = 1),  # Show regular numbers with commas
    name = "Total Insecticides Applied (ml)"
  ) +
  coord_cartesian(ylim = c(0, 1.2)) +
  theme_cowplot() +
  labs(
    y = "Average Crithidia Prevalence"
  ) +
  theme(
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    text = element_text(size = 18)
  ) +
  annotate(
    geom = "text",
    x = 20000,
    y = 1.2,
    label = "P = 0.1",
    size = 12
  )