---
title: "Script_6_Spatial_model_itraconazole_resistant_cyp51_TR_data_science"
author: "Hylke Kortenbosch"
date: "2024-09-30"
output: html_document
---

Load required packages
```{r}
packages <- c(
  "sf",          # For spatial objects
  "glmmTMB",     # For beta-binomial models
  "MuMIn",       # For model selection
  "ggplot2",     # For plotting
  "tmap",        # For polygon of the Netherlands
  "gstat",       # For geostatistical analyses
  "ggspatial",   # For spatial plotting
  "MASS",        # For stepwise model selection
  "stars",       # For converting coordinate systems
  "tidyverse",   # For data manipulation
  "terra",       # For making spatial raster
  "grid",        # For grid graphics
  "gridExtra",   # For arranging multiple plots
  "performance", # For checking model VIFs
  "svglite"      # For generating SVG files
)

# Function to check, install (if needed), and load packages
install_and_load <- function(pkg) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}

# Apply the function to each package
lapply(packages, install_and_load)
```

# Key session info
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default


attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] svglite_2.1.3      performance_0.12.4 gridExtra_2.3      terra_1.7-71       lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
[12] tibble_3.2.1       tidyverse_2.0.0    stars_0.6-5        abind_1.4-5        MASS_7.3-60.2      gstat_2.1-1        tmap_3.3-4         ggspatial_1.1.9    ggplot2_3.5.1      MuMIn_1.47.5       glmmTMB_1.1.9     
[23] sf_1.0-16  


# Load sf object schimmelradar data
```{r}
load("Schimmelradar360spat_focal_TR_ITRA.RData")
```

#Generate spatial polygon data frame of the Netherlands for use in plotting
```{r}
data(NLD_prov)
provincesNL <- sf::st_sf(NLD_prov)

# transform to appropriate coordinate system
st_crs(provincesNL) <- st_crs(provincesNL)
provincesNL_GPS <- st_transform(provincesNL, crs = 4326)
provincesNL_transformed <- st_transform(provincesNL, crs = 28992)
```



#Make plot Itraconazole resistant TR34 fraction points on a map
```{r}
Schimmelradar360spat_focal_amersfoort <- st_transform(Schimmelradar360spat_focal_TR_ITRA, crs = 28992)

# Extract the new coordinates and update the longitude and latitude columns
Schimmelradar360spat_focal_amersfoort$long <- st_coordinates(Schimmelradar360spat_focal_amersfoort)[, 1]  # X (eastings)
Schimmelradar360spat_focal_amersfoort$lat <- st_coordinates(Schimmelradar360spat_focal_amersfoort)[, 2]   # Y (northings)

# Define meaningful breaks for the alpha scale (consistent with other plots)
true_breaks <- c(0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.30)  

letter_grob <- textGrob("A", x = 0.1, y = 0.9,  
                        hjust = 1, vjust = 0, gp = gpar(fontsize = 8, fontface = "bold"))

# Use continuous alpha scaling instead of binning
ITR34frac_points_NL_alpha <- ggplot() + 
  geom_sf(data = provincesNL_transformed, fill = "white", linewidth = 0.05) +  # Thin province outlines
  coord_sf(datum = st_crs(28992)) +
  geom_point(
    data = as.data.frame(Schimmelradar360spat_focal_TR_VORI_plot),
    aes(x = long, y = lat, alpha = TR34VORIfrac),  # Use actual continuous resistance values
    color = "black",  
    size = 0.2  
  ) +
  # Apply continuous alpha scale with limits to unify the legend
  scale_alpha_continuous(
    name = "TR34 Fraction",  
    breaks = true_breaks,  
    labels = true_breaks,  
    range = c(0.05, 1),  # Ensures meaningful transparency variation
    limits = c(0, 0.35),  # Match other resistance plots for comparability
    guide = guide_legend(
      reverse = TRUE, 
      keywidth = unit(0.15, "cm"),  
      keyheight = unit(0.15, "cm")
    )
  ) +
  labs(
    x = "Easting (meters)",
    y = "Northing (meters)"
  ) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  # Adjust aspect ratio for compactness
    axis.title = element_blank(),  # Smaller axis titles
    axis.text = element_blank(),   # Smaller axis text
    axis.ticks = element_blank(), 
    legend.title = element_text(size = 6),  # Smaller legend title
    legend.text = element_text(size = 5),     # Smaller legend text
    legend.box.margin = margin(0, 0, 0, 0),   # Compact legend margin
    legend.spacing.y = unit(0.05, 'cm'),      # Reduced legend spacing
    legend.spacing.x = unit(0.05, 'cm'),       # Reduced legend spacing 
     panel.background = element_rect(fill = "white", color = NA),  # White panel background
    plot.background = element_rect(fill = "white", color = NA)  # White plot background    
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      keywidth = unit(0.2, "cm"),  # Compact legend key width
      keyheight = unit(0.2, "cm")  # Compact legend key height
    )) +
  annotation_scale(
    location = "bl",  # Bottom-left corner
   # width_hint = 0.4, # Adjust the width of the scalebar
   # bar_cols = c("black", "white"),  # Colors for the scalebar     style = "ticks",  # Tick-style scale bar
   style = "ticks",
   line_width = 0.5, # Thickness of ticks
   tick_height = 0.5, # Height of ticks
    text_cex = 0.5  # Text size for the scale label
  ) +
  annotation_custom(grob = letter_grob)


ggsave("ITRA_TR34_points_NL.pdf", plot = ITRATR34frac_points_NL, width = 7, height = 7, dpi = 600)
```

Model selection for spatial Itraconazole resistant TR34 land-use model (Tables S14 and S15)
```{r}
cor(as.data.frame(Schimmelradar360spat_focal_TR_ITRA)[,c(27:40,55:60,62,63)]) # strong correlation between beetroot and potato, only use potato because most informative of the two

binmodelTR34potato <- glm(cbind(TR34, Tot.ITR - TR34) ~ bulb22fra + fruit22fra + tree_farming22fra + greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = binomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(binmodelTR34potato) # model is strongly overdispersed, go for betabinomial instead

#specify global model
betamodelITRTR34 <- glmmTMB(cbind(TR34, Tot.ITR - TR34) ~ bulb22fra + fruit22fra + tree_farming22fra +  greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(betamodelVORTR34)

#Specify null model
null_betamodelITRTR34 <- glmmTMB(cbind(TR34, Tot.ITR - TR34) ~ 1, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(null_betamodelITRTR34)

Schimmelradar360spat_focal_TR_ITRA.df <- as.data.frame(Schimmelradar360spat_focal_TR_ITRA)

variable_names <- c("bulb22fra", "buildings22fra", "forest22fra", "fruit22fra", "grain22fra", "grass22fra", "greenhouse22fra","maize22fra", "Onion22fra",  "potato22fra", "tree_farming22fra", "water22fra") 

num_vars <- length(variable_names)
best_aic <- Inf  # Initialize with a large value
best_model_TR34ITRRES_naive <- NULL
model_counter <- 0  # Initialize counter for number of models run

# Initialize an empty DataFrame to store all models
saved_models_TR34ITRRES_naive_df <- data.frame(Model = character(), AIC = numeric(), stringsAsFactors = FALSE)


for (i in 1:num_vars) {
  subsets <- combn(num_vars, i)
  for (j in 1:ncol(subsets)) {
    predictors <- variable_names[subsets[, j]]
    formula <- as.formula(paste("cbind(TR34, Tot.ITR - TR34) ~", 
                                 paste(predictors, collapse = " + ")))
    # Fit model with fixed dispersion and no zero-inflation
    model <- glmmTMB(formula, 
                     family = betabinomial(), 
                     dispformula = ~1, 
                     ziformula = ~0, 
                     data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    aic <- AIC(model)
    
    if (aic < best_aic) {
        best_aic <- aic
        best_model_ITRRES_naive <- formula
    }

    
    # Store all models in the DataFrame
    saved_models_TR34ITRRES_naive_df <- rbind(saved_models_TR34ITRRES_naive_df, 
                              data.frame(Model = as.character(formula), AIC = aic))
    
    # Increment model counter
    model_counter <- model_counter + 1
    
    # Print progress (model combination number and its AIC)
    cat("Tested model:", model_counter, "AIC:", aic, "with predictors:", predictors, "\n")
  }
}

# Print the best model formula and AIC value
cat("Best model found: \n", best_model, "\nBest AIC:", best_aic, "\n")

# Filter the DataFrame to keep only models with AIC <= best AIC + 2
qualified_models_TR34ITRRES_naive_df <- saved_models_TR34ITRRES_naive_df[saved_models_TR34ITRRES_naive_df$AIC <= best_aic + 2, ]

# Remove rows where any value starts with "~" or "cbind"
qualified_models_TR34ITRRES_naive_df_cleaned <- qualified_models_TR34ITRRES_naive_df[!apply(qualified_models_TR34ITRRES_naive_df, 1, function(row) {
  any(grepl("^~", row) | grepl("^cbind", row))
}), ]

write.csv(qualified_models_TR34ITRRES_naive_df_cleaned,"qualified_models_TR34ITRRES_naive_df.csv", row.names = FALSE)
qualified_models_TR34ITRRES_naive_df_cleaned <- read.csv("qualified_models_TR34ITRRES_naive_df.csv")


# Initialize a list to store the actual model objects of qualified models
qualified_model_TR34ITRRES_naive_df_list <- list()

# Second loop: Generate only the qualified models
for (i in 1:nrow(qualified_models_TR34ITRRES_naive_df_cleaned)) {
    # Recreate the model formula from the saved string
    model_formula <- as.formula(paste("cbind(TR34, Tot.ITR - TR34) ~", qualified_models_TR34ITRRES_naive_df_cleaned$Model[i]))
    
    # Fit the model again (only for qualified models)
    model <- glmmTMB(model_formula, family = betabinomial(), dispformula = ~1, ziformula = ~0, data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    # Store the qualified model in the list
    qualified_model_TR34ITRRES_naive_df_list[[i]] <- model
}

#generate avaraged model
AVGTR34ITRRESmodels_naive <- model.avg(qualified_model_TR34ITRRES_naive_df_list)
Summary_AVGTR34ITRRESmodels_naive <- summary(AVGTR34ITRRESmodels_naive)

# Generate output table of avaraged model model
Sum_of_weights <- as.data.frame(Summary_AVGTR34ITRRESmodels_naive$sw)
colnames(Sum_of_weights) <- "Sum of Weights"
Conditional_summary_TR34ITRRES_naive_model <- merge(Summary_AVGTR34ITRRESmodels_naive$coefmat.subset[,-2], Sum_of_weights, by = "row.names")
Conditional_summary_TR34ITRRES_naive_model <- Conditional_summary_TR34ITRRES_naive_model[match(rownames(Summary_AVGTR34ITRRESmodels_naive$coefmat.subset[,]), Conditional_summary_TR34ITRRES_naive_model$Row.names), ]
Conditional_summary_TR34ITRRES_naive_model[1,] <- c(as.character("Intercept"), Summary_AVGTR34ITRRESmodels_naive$coefmat.subset[1,-2], 1)

Conditional_summary_TR34ITRRES_naive_model <- Conditional_summary_TR34ITRRES_naive_model[,-1]
rownames(Conditional_summary_TR34ITRRES_naive_model) <- c("Intercept" , "Forest" , "Grains" , "Grassland", "Flowerbulbs", "Water" ,"Potatoes", "Greenhouses" , "Maize",  "Onions" ,  "Fruit" ,  "Tree Farming", "Urbanisation") 

#Turn dataframe numeric again
Conditional_summary_TR34ITRRES_naive_model[] <- lapply(Conditional_summary_TR34ITRRES_naive_model, function(x) {
  if (is.character(x)) as.numeric(x) else x
})

Output_summary_TR34ITRRES_naive_model  <- signif(Conditional_summary_TR34ITRRES_naive_model, digits = 3)

# Check if the file already exists
if(file.exists("Output_summary_TR34ITRRES_naive_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_summary_TR34ITRRES_naive_model.csv")
}
write.csv(Output_summary_TR34ITRRES_naive_model, "Output_summary_TR34ITRRES_naive_model.csv", row.names = T, quote = F) # Table S13

# use best fitting model for plotting prediction
best_betamodelTR34ITRRES <- glmmTMB(cbind(TR34, Tot.ITR - TR34) ~ grain22fra, family = betabinomial, data = Schimmelradar360spat_focal_TR_ITRA.df, dispformula = ~1, ziformula = ~0)
#Only one variable in minimal AIC model, no need to check for VIF

Output_best_TR34ITRRES_model <- as.data.frame(signif(summary(best_betamodelTR34ITRRES)$ coefficients$cond, digits = 3))

# Check if the file already exists
if(file.exists("Output_best_TR34ITRRES_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_best_TR34ITRRES_model.csv")
}
write.csv(Output_best_TR34ITRRES_model, "Output_best_TR34ITRRES_model.csv", row.names = T, quote = F) # Table S14
```

#Experimental variogram of residuals of Itraconazole resistant TR34 fraction 
```{r}
# experimental variogram log transformed data; it is recommended using GLM if 
# data are counts but let us first check this approach, which is quite common.
vg.NULL_TR34ITR <- variogram(resid(null_betamodelITRTR34, type = "response") ~ 1, Schimmelradar360spat_focal_TR_ITRA, boundaries = 1:20 * 5e3)
plot(vg.NULL_TR34ITR, plot.numbers = T)

vg.TR34ITR <- variogram(resid(best_betamodelTR34ITRRES, type = "response") ~ 1, Schimmelradar360spat_focal_TR_ITRA, boundaries = 1:20 * 5e3)
plot(vg.TR34ITR, plot.numbers = T)

#set distance to km
vg.TR34ITR$dist_km <- vg.TR34ITR$dist / 1000
vg.NULL_TR34ITR$dist_km <- vg.NULL_TR34ITR$dist / 1000

max_range <- max(max(vg.TR34ITR$dist_km), max(vg.NULL_TR34ITR$dist_km))


# Compute the overall range of the semi-variance values
overall_range <- range(c(0, vg.TR34ITR$gamma, vg.NULL_TR34ITR$gamma))
dev.off()
# Create an empty plot with appropriate range
plot(1, type = "n", xlim = c(0, max_range), ylim = overall_range,
     xlab = "Distance (km)", ylab = "Semi-variance")
# Plot the points for the first variogram
points(vg.TR34ITR$dist_km, vg.TR34ITR$gamma, col = "blue", pch = 16)
text(vg.TR34ITR$dist_km, vg.TR34ITR$gamma, labels = vg.TR34ITR$n, pos = 2, cex = 0.7, col = "blue")

# Plot the points for the second variogram
points(vg.NULL_TR34ITR$dist_km, vg.NULL_TR34ITR$gamma, col = "red", pch = 16)


legend("bottomright", legend = c("Minimal AIC Model Residuals", "Null Model Residuals"),
       col = c("blue", "red"), pch = 16, cex = 1,  
       x.intersp = 1.5, y.intersp = 1.5, box.lwd = 2)
# Save the plot as TIFF with DPI 300
tiff("TR34ITR_variogram_plot.tiff", width = 7, height = 7, units = "in", res = 300)
dev.off()
```
Note on semivariogram: 
The numbers represent the number of pairwise distances per distance bin. The null model is included to visualize how the land use model reduces the residual semivariance through space. The semivariance is the variance between a subgroup of samples that are a certain distance apart, for example all samples that are 5-10 km apart. Here, there is no evidence of spatial correlation of the residuals for the land use model here and the independence assumption is reasonable.

#Predicting Itraconazole resistant TR34 through space and plotting the prediction
```{r}
Predictor_raster_NL <- rast("Predictor_raster_NL.tif")

Itraconazole_resistant_TR34_naive_RF_NL <- terra::predict(Predictor_raster_NL, best_betamodelTR34ITRRES, filename = "Itraconazole_resistant_TR34_naive_NL.tif", type = "response", overwrite = T)

Itraconazole_resistant_TR34_naive_RF_NL <- rast("Itraconazole_resistant_TR34_naive_NL.tif")
# Change to standard coordinate system for familiarity readers


Itraconazole_resistant_TR34_naive_RF_NL_df <- as.data.frame(Itraconazole_resistant_TR34_naive_RF_NL, xy = TRUE)

#cap at 0.15
Itraconazole_resistant_TR34_naive_RF_NL_df$capped_lyr1 <- pmin(Itraconazole_resistant_TR34_naive_RF_NL_df$lyr1, 0.15)


custom_colors <- c("white", "darkblue") 


# Define consistent breaks for all resistance fraction plots
true_breaks <- c(0, 0.02, 0.05, 0.1, 0.15)

letter_grob <- textGrob("B", 
                        x = 0.1, 
                        y = 0.9, 
                        hjust = 1, 
                        vjust = 0, 
                        gp = gpar(fontsize = 8, fontface = "bold"))

Itraconazole_resistant_TR34_naive_resistance_fraction_NL <- ggplot(Itraconazole_resistant_TR34_naive_RF_NL_df, aes(x = x, y = y, fill = capped_lyr1)) +
  geom_raster() +
 # Apply custom white-to-blue gradient
  scale_fill_gradientn(
    name = "TR34 Fraction",
    colors = custom_colors,
    breaks = true_breaks,
    labels = c("0", "0.02", "0.05", "0.10", expression("\u2265 0.15")),
    limits = c(0, 0.15)  # Explicitly set range
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = unit(0.15, "cm"),  
      barheight = unit(2, "cm"),
      direction = "vertical",
      reverse = FALSE  
    )
  ) +
  labs(
    x = "Easting (meters)",
    y = "Northing (meters)"
  ) +
  scale_x_continuous(breaks = seq(0, max(Itraconazole_resistant_TR34_naive_RF_NL_df$x), by = 50000)) +
  scale_y_continuous(breaks = seq(0, max(Itraconazole_resistant_TR34_naive_RF_NL_df$y), by = 50000)) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  # Adjust aspect ratio for compactness
    axis.title = element_blank(),  # Smaller axis titles
    axis.text = element_blank(),   # Smaller axis text
    axis.ticks = element_blank(), 
    legend.title = element_text(size = 6),  # Smaller legend title
    legend.text = element_text(size = 5),     # Smaller legend text
    legend.box.margin = margin(0, 0, 0, 0),   # Compact legend margin
    legend.spacing.y = unit(0.05, 'cm'),      # Reduced legend spacing
    legend.spacing.x = unit(0.05, 'cm'),      # Reduced legend spacing
     panel.background = element_rect(fill = "white", color = NA),  # White panel background
    plot.background = element_rect(fill = "white", color = NA)  # White plot background         
  )  +
  annotation_custom(grob = letter_grob)

ggsave("Itraconazole_resistant_TR34_naive_resistance_fraction_NL.eps", plot = Itraconazole_resistant_TR34_naive_resistance_fraction_NL, width = 7, height = 7, dpi = 300)
```

#K-fold cross validation of itraconazole resistant TR34 fraction model
```{r}
# 5(k)-fold cross-validation
set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){
  
  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(TR34, Tot.ITR - TR34) ~ 1,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$TR34frac - prd)

 
}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_NULL_RMSE_TR34ITRRES  <- sqrt(mean(res^2))
five_fold_NULL_RMSE_TR34ITRRES  <- 0.02038664
```

#Leave one out cross validation of itraconazole resistance fraction model
```{r}
# 5(k)-fold cross-validation

set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){

 

  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(TR34, Tot.ITR - TR34) ~ grain22fra,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$TR34frac - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_RMSE_TR34ITRRES   <- sqrt(mean(res^2))
five_fold_RMSE_TR34ITRRES   <- 0.0203696
```

#Make plot Itraconazole resistant TR46 fraction points on a map
```{r}
Schimmelradar360spat_focal_amersfoort <- st_transform(Schimmelradar360spat_focal_TR_ITRA, crs = 28992)

# Extract the new coordinates and update the longitude and latitude columns
Schimmelradar360spat_focal_amersfoort$long <- st_coordinates(Schimmelradar360spat_focal_amersfoort)[, 1]  # X (eastings)
Schimmelradar360spat_focal_amersfoort$lat <- st_coordinates(Schimmelradar360spat_focal_amersfoort)[, 2]   # Y (northings)

# Define meaningful breaks for the alpha scale (consistent with other plots)
true_breaks <- c(0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.30)  

# Adjusted letter annotation
letter_grob <- textGrob("C", x = 0.1, y = 0.9, 
                        hjust = 1, vjust = 0, gp = gpar(fontsize = 8, fontface = "bold"))

# Plot with legend retained and compact adjustments
ITRATR46frac_points_NL_alpha <- ggplot() + 
  geom_sf(data = provincesNL_transformed, fill = "white", linewidth = 0.05) +
  coord_sf(datum = st_crs(28992)) +
  geom_point(
    data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA_plot), 
    aes(x = long, y = lat, alpha = TR46ITRfrac), size = 0.2, color =
    "black"
  )  + 
  # Apply continuous alpha scale with limits to unify the legend
  scale_alpha_continuous(
    name = "TR46 Fraction",  
    breaks = true_breaks,  
    labels = true_breaks,  
    range = c(0.05, 1),  # Ensures meaningful transparency variation
    limits = c(0, 0.35),  # Match other resistance plots for comparability
    guide = guide_legend(
      reverse = TRUE, 
      keywidth = unit(0.15, "cm"),  
      keyheight = unit(0.15, "cm")
    )
  ) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  # Adjust aspect ratio for compactness
    axis.title = element_blank(),  # Smaller axis titles
    axis.text = element_blank(),   # Smaller axis text
    axis.ticks = element_blank(), 
    legend.title = element_text(size = 6),  # Smaller legend title
    legend.text = element_text(size = 5),     # Smaller legend text
    legend.box.margin = margin(0, 0, 0, 0),   # Compact legend margin
    legend.spacing.y = unit(0.05, 'cm'),      # Reduced legend spacing
    legend.spacing.x = unit(0.05, 'cm'),       # Reduced legend spacing 
     panel.background = element_rect(fill = "white", color = NA),  # White panel background
    plot.background = element_rect(fill = "white", color = NA)  # White plot background 
  ) +
  guides(
    color = guide_legend(
      reverse = TRUE,
      keywidth = unit(0.2, "cm"),  # Compact legend key width
      keyheight = unit(0.2, "cm")  # Compact legend key height
    )) +
  annotation_custom(grob = letter_grob)

ggsave("ITRA_TR46_points_NL.pdf", plot = ITRATR46frac_points_NL, width = 7, height = 7, dpi = 600)
```

#Model selection for spatial Itraconazole-resistant TR46 land-use model (Tables S16 and S17)
```{r}
binmodelTR46potato <- glm(cbind(TR46, Tot.ITR - TR46) ~ bulb22fra + fruit22fra + tree_farming22fra + greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = binomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(binmodelTR46potato) # model is strongly overdispersed, go for betabinomial instead

#specify global model
betamodelITRTR46 <- glmmTMB(cbind(TR46, Tot.ITR - TR46) ~ bulb22fra + fruit22fra + tree_farming22fra +  greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(betamodelITRTR46)

#Specify null model
null_betamodelITRTR46 <- glmmTMB(cbind(TR46, Tot.ITR - TR46) ~ 1, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(null_betamodelITRTR46)

Schimmelradar360spat_focal_TR_ITRA.df <- as.data.frame(Schimmelradar360spat_focal_TR_ITRA)


variable_names <- c("bulb22fra", "buildings22fra", "forest22fra", "fruit22fra", "grain22fra", "grass22fra", "greenhouse22fra","maize22fra", "Onion22fra",  "potato22fra", "tree_farming22fra", "water22fra") 

num_vars <- length(variable_names)
best_aic <- Inf  # Initialize with a large value
best_model_TR46ITRRES_naive <- NULL
model_counter <- 0  # Initialize counter for number of models run

# Initialize an empty DataFrame to store all models
saved_models_TR46ITRRES_naive_df <- data.frame(Model = character(), AIC = numeric(), stringsAsFactors = FALSE)


for (i in 1:num_vars) {
  subsets <- combn(num_vars, i)
  for (j in 1:ncol(subsets)) {
    predictors <- variable_names[subsets[, j]]
    formula <- as.formula(paste("cbind(TR46, Tot.ITR - TR46) ~", 
                                 paste(predictors, collapse = " + ")))
    # Fit model with fixed dispersion and no zero-inflation
    model <- glmmTMB(formula, 
                     family = betabinomial(), 
                     dispformula = ~1, 
                     ziformula = ~0, 
                     data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    aic <- AIC(model)
    
    if (aic < best_aic) {
        best_aic <- aic
        best_model_ITRRES_naive <- formula
    }

    
    # Store all models in the DataFrame
    saved_models_TR46ITRRES_naive_df <- rbind(saved_models_TR46ITRRES_naive_df, 
                              data.frame(Model = as.character(formula), AIC = aic))
    
    # Increment model counter
    model_counter <- model_counter + 1
    
    # Print progress (model combination number and its AIC)
    cat("Tested model:", model_counter, "AIC:", aic, "with predictors:", predictors, "\n")
  }
}

# Print the best model formula and AIC value
cat("Best model found: \n", best_model, "\nBest AIC:", best_aic, "\n")

# Filter the DataFrame to keep only models with AIC <= best AIC + 2
qualified_models_TR46ITRRES_naive_df <- saved_models_TR46ITRRES_naive_df[saved_models_TR46ITRRES_naive_df$AIC <= best_aic + 2, ]

# Remove rows where any value starts with "~" or "cbind"
qualified_models_TR46ITRRES_naive_df_cleaned <- qualified_models_TR46ITRRES_naive_df[!apply(qualified_models_TR46ITRRES_naive_df, 1, function(row) {
  any(grepl("^~", row) | grepl("^cbind", row))
}), ]

write.csv(qualified_models_TR46ITRRES_naive_df_cleaned,"qualified_models_TR46ITRRES_naive_df.csv", row.names = FALSE)
qualified_models_TR46ITRRES_naive_df_cleaned <- read.csv("qualified_models_TR46ITRRES_naive_df.csv")


# Initialize a list to store the actual model objects of qualified models
qualified_model_TR46ITRRES_naive_df_list <- list()

# Second loop: Generate only the qualified models
for (i in 1:nrow(qualified_models_TR46ITRRES_naive_df_cleaned)) {
    # Recreate the model formula from the saved string
    model_formula <- as.formula(paste("cbind(TR46, Tot.ITR - TR46) ~", qualified_models_TR46ITRRES_naive_df_cleaned$Model[i]))
    
    # Fit the model again (only for qualified models)
    model <- glmmTMB(model_formula, family = betabinomial(), dispformula = ~1, ziformula = ~0, data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    # Store the qualified model in the list
    qualified_model_TR46ITRRES_naive_df_list[[i]] <- model
}

#generate avaraged model
AVGTR46ITRRESmodels_naive <- model.avg(qualified_model_TR46ITRRES_naive_df_list)
Summary_AVGTR46ITRRESmodels_naive <- summary(AVGTR46ITRRESmodels_naive)

# Generate output table of avaraged model model
Sum_of_weights <- as.data.frame(Summary_AVGTR46ITRRESmodels_naive$sw)
colnames(Sum_of_weights) <- "Sum of Weights"
Conditional_summary_TR46ITRRES_naive_model <- merge(Summary_AVGTR46ITRRESmodels_naive$coefmat.subset[,-2], Sum_of_weights, by = "row.names")
Conditional_summary_TR46ITRRES_naive_model <- Conditional_summary_TR46ITRRES_naive_model[match(rownames(Summary_AVGTR46ITRRESmodels_naive$coefmat.subset[,]), Conditional_summary_TR46ITRRES_naive_model$Row.names), ]
Conditional_summary_TR46ITRRES_naive_model[1,] <- c(as.character("Intercept"), Summary_AVGTR46ITRRESmodels_naive$coefmat.subset[1,-2], 1)

Conditional_summary_TR46ITRRES_naive_model <- Conditional_summary_TR46ITRRES_naive_model[,-1]
rownames(Conditional_summary_TR46ITRRES_naive_model) <- c("Intercept" , "Flowerbulbs",     "Urbanisation", "Forest" ,"Maize", "Water" ,  "Grassland","Onions" , "Greenhouses" , "Fruit", "Tree Farming", "Grains" , "Potatoes") 

#Turn dataframe numeric again
Conditional_summary_TR46ITRRES_naive_model[] <- lapply(Conditional_summary_TR46ITRRES_naive_model, function(x) {
  if (is.character(x)) as.numeric(x) else x
})

Output_summary_TR46ITRRES_naive_model  <- signif(Conditional_summary_TR46ITRRES_naive_model, digits = 3)

# Check if the file already exists
if(file.exists("Output_summary_TR46ITRRES_naive_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_summary_TR46ITRRES_naive_model.csv")
}
write.csv(Output_summary_TR46ITRRES_naive_model, "Output_summary_TR46ITRRES_naive_model.csv", row.names = T, quote = F)

# use best fitting model for plotting prediction
best_betamodelTR46ITRRES <- glmmTMB(cbind(TR46, Tot.ITR - TR46) ~ bulb22fra + buildings22fra + forest22fra + maize22fra + water22fra, family = betabinomial, data = Schimmelradar360spat_focal_TR_ITRA.df, dispformula = ~1, ziformula = ~0)
check_collinearity(best_betamodelTR46ITRRES) # Check VIFs

Output_best_TR46ITRRES_model <- as.data.frame(signif(summary(best_betamodelTR46ITRRES)$ coefficients$cond, digits = 3))

# Check if the file already exists
if(file.exists("Output_best_TR46ITRRES_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_best_TR46ITRRES_model.csv")
}
write.csv(Output_best_TR46ITRRES_model, "Output_best_TR46ITRRES_model.csv", row.names = T, quote = F)
```

#Experimental variogram of residuals of Itraconazole resistant TR46 fraction
```{r}
# experimental variogram log transformed data; it is recommended using GLM if 
# data are counts but let us first check this approach, which is quite common.
vg.NULL_TR46ITR <- variogram(resid(null_betamodelITRTR46, type = "response") ~ 1, Schimmelradar360spat_focal_TR_ITRA, boundaries = 1:20 * 5e3)
plot(vg.NULL_TR46ITR, plot.numbers = T)

vg.TR46ITR <- variogram(resid(best_betamodelTR46ITRRES, type = "response") ~ 1, Schimmelradar360spat_focal_TR_ITRA, boundaries = 1:20 * 5e3)
plot(vg.TR46ITR, plot.numbers = T)

vg.NULL_TR46ITR$dist_km <- vg.NULL_TR46ITR$dist / 1000
vg.TR46ITR$dist_km <- vg.TR46ITR$dist / 1000

max_range <- max(max(vg.TR46ITR$dist_km), max(vg.NULL_TR46ITR$dist_km))


# Compute the overall range of the semi-variance values
overall_range <- range(c(0, vg.TR46ITR$gamma, vg.NULL_TR46ITR$gamma))
dev.off()
# Create an empty plot with appropriate range
plot(1, type = "n", xlim = c(0, max_range), ylim = overall_range,
     xlab = "Distance (km)", ylab = "Semi-variance")
# Plot the points for the first variogram
points(vg.TR46ITR$dist_km, vg.TR46ITR$gamma, col = "blue", pch = 16)
text(vg.TR46ITR$dist_km, vg.TR46ITR$gamma, labels = vg.TR46ITR$n, pos = 2, cex = 0.7, col = "blue")

# Plot the points for the second variogram
points(vg.NULL_TR46ITR$dist_km, vg.NULL_TR46ITR$gamma, col = "red", pch = 16)


legend("bottomright", legend = c("Minimal AIC Model Residuals", "Null Model Residuals"),
       col = c("blue", "red"), pch = 16, cex = 1,  
       x.intersp = 1.5, y.intersp = 1.5, box.lwd = 2)
# Save the plot as TIFF with DPI 300
tiff("TR46ITR_variogram_plot.tiff", width = 7, height = 7, units = "in", res = 300)
dev.off()
```
Note on semivariogram:
The numbers represent the number of pairwise distances per distance bin. Due to our sampling strategy, there are relatively few pairwise distances between samples of less than 5 km. The null model is included to visualize how the land use model reduces the residual variance through space.  Here, here the residuals show evidence of limited spatial autocorrelation up to approximately 5 km distance. However half of the residual spatial variance occurs at distances below the shortest distance between sample pairs. Furthermore, due to our sampling strategy, there are relatively few pairwise distances between samples of less than 5 km. The independence assumption in our model is therefore reasonable. 

#Predicting Itraconazole resistant TR46 through space and plotting the prediction
```{r}
Predictor_raster_NL <- rast("Predictor_raster_NL.tif")

Itraconazole_resistant_TR46_naive_RF_NL <- terra::predict(Predictor_raster_NL, best_betamodelTR46ITRRES, filename = "Itraconazole_resistant_TR46_naive_NL.tif", type = "response", overwrite = T)

Itraconazole_resistant_TR46_naive_RF_NL <- rast("Itraconazole_resistant_TR46_naive_NL.tif")
# Change to standard coordinate system for familiarity readers


Itraconazole_resistant_TR46_naive_RF_NL_df <- as.data.frame(Itraconazole_resistant_TR46_naive_RF_NL, xy = TRUE)

#cap at 0.15
Itraconazole_resistant_TR46_naive_RF_NL_df$capped_lyr1 <- pmin(Itraconazole_resistant_TR46_naive_RF_NL_df$lyr1, 0.15)


# Define consistent breaks for all resistance fraction plots
true_breaks <- c(0, 0.02, 0.05, 0.1, 0.15)

letter_grob <- textGrob("D", x = 0.1, y = 0.9, 
                        hjust = 1, vjust = 0, gp = gpar(fontsize = 8, fontface = "bold"))

# Plot using ggplot2
Itraconazole_resistant_TR46_naive_resistance_fraction_NL <- ggplot(Itraconazole_resistant_TR46_naive_RF_NL_df, aes(x = x, y = y, fill = capped_lyr1)) +
  geom_raster() +
 # Apply custom white-to-blue gradient
  scale_fill_gradientn(
    name = "TR46 Fraction",
    colors = custom_colors,
    breaks = true_breaks,
    labels = c("0", "0.02", "0.05", "0.10", expression("\u2265 0.15")),
    limits = c(0, 0.15)  # Explicitly set range
  ) +
  guides(
    fill = guide_colorbar(
      barwidth = unit(0.15, "cm"),  
      barheight = unit(2, "cm"),
      direction = "vertical",
      reverse = FALSE  
    )
  ) +
  labs(
    x = "Easting (meters)",
    y = "Northing (meters)"
  ) +
  scale_x_continuous(breaks = seq(0, max(Itraconazole_resistant_TR46_naive_RF_NL_df$x), by = 50000)) +
scale_y_continuous(breaks = seq(0, max(Itraconazole_resistant_TR46_naive_RF_NL_df$y), by = 50000)) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  # Adjust aspect ratio for compactness
    axis.title = element_blank(),  # Smaller axis titles
    axis.text = element_blank(),   # Smaller axis text
    axis.ticks = element_blank(), 
    legend.title = element_text(size = 6),  # Smaller legend title
    legend.text = element_text(size = 5),     # Smaller legend text
    legend.box.margin = margin(0, 0, 0, 0),   # Compact legend margin
    legend.spacing.y = unit(0.05, 'cm'),      # Reduced legend spacing
    legend.spacing.x = unit(0.05, 'cm'), 
     panel.background = element_rect(fill = "white", color = NA),  # White panel background
    plot.background = element_rect(fill = "white", color = NA)  # White plot background         
  ) +
  annotation_custom(grob = letter_grob)

ggsave("Itraconazole_resistant_TR46_naive_resistance_fraction_NL.eps", plot = Itraconazole_resistant_TR46_naive_resistance_fraction_NL, width = 7, height = 7, dpi = 300)
```

#K-fold cross validation of itraconazole resistant TR46 fraction null model
```{r}
# 5(k)-fold cross-validation
set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){
  
  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(TR46, Tot.ITR - TR46) ~ 1,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$TR46ITRfrac - prd)

 
}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_NULL_RMSE_TR46ITRRES  <- sqrt(mean(res^2))
five_fold_NULL_RMSE_TR46ITRRES  <- 0.02309723
```

#K-fold cross validation of itraconazole resistant TR46 fraction land use model
```{r}
# 5(k)-fold cross-validation

set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){

 

  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(TR46, Tot.ITR - TR46) ~ bulb22fra + buildings22fra + forest22fra + maize22fra + water22fra,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$TR46ITRfrac - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_RMSE_TR46ITRRES   <- sqrt(mean(res^2))
five_fold_RMSE_TR46ITRRES   <- 0.02078558
```



Combining raw data and prediction plot into 2x2 grid
```{r}
Figure_naive_itraconazole_resistance_TR <- grid.arrange(ITR34frac_points_NL_alpha, Itraconazole_resistant_TR34_naive_resistance_fraction_NL, ITRATR46frac_points_NL_alpha, Itraconazole_resistant_TR46_naive_resistance_fraction_NL,  nrow = 2, ncol = 2)

ggsave("Figure_naive_itraconazole_resistance_TR.svg", plot = Figure_naive_itraconazole_resistance_TR, width = 6, height = 4, dpi = 600)
ggsave("Figure_naive_itraconazole_resistance_TR.jpeg", plot = Figure_naive_itraconazole_resistance_TR, width = 6, height = 4, dpi = 600)
```

#Model selection for spatial Itraconazole-resistant WT land-use model (Tables S18 and S19)
```{r}
#check for correlation explanatory variables
cor(as.data.frame(Schimmelradar360spat_focal_TR_ITRA)[,c(27:40,55:60,62,63)]) # strong correlation between beetroot and potato, only use potato because most informative of the two

binmodelwt <- glm(cbind(wt, Tot.ITR - wt) ~ bulb22fra + fruit22fra + tree_farming22fra + greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = binomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(binmodelwt) # model is strongly overdispersed, go for betabinomial instead

#specify global model
betamodelITRwt <- glmmTMB(cbind(wt, Tot.ITR - wt) ~ bulb22fra + fruit22fra + tree_farming22fra +  greenhouse22fra + grain22fra + forest22fra + Onion22fra + buildings22fra + grass22fra + maize22fra + potato22fra + water22fra, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(betamodelITRWT)


#Specify null model
null_betamodelITRWT <- glmmTMB(cbind(wt, Tot.ITR - wt) ~ 1, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal_TR_ITRA))
summary(null_betamodelITRWT)

Schimmelradar360spat_focal_TR_ITRA.df <- as.data.frame(Schimmelradar360spat_focal_TR_ITRA)

variable_names <- c("bulb22fra", "buildings22fra", "forest22fra", "fruit22fra", "grain22fra", "grass22fra", "greenhouse22fra","maize22fra", "Onion22fra",  "potato22fra", "tree_farming22fra", "water22fra") 

num_vars <- length(variable_names)
best_aic <- Inf  # Initialize with a large value
best_model_wtITRRES_naive <- NULL
model_counter <- 0  # Initialize counter for number of models run

# Initialize an empty DataFrame to store all models
saved_models_wtITRRES_naive_df <- data.frame(Model = character(), AIC = numeric(), stringsAsFactors = FALSE)


for (i in 1:num_vars) {
  subsets <- combn(num_vars, i)
  for (j in 1:ncol(subsets)) {
    predictors <- variable_names[subsets[, j]]
    formula <- as.formula(paste("cbind(wt, Tot.ITR - wt) ~", 
                                 paste(predictors, collapse = " + ")))
    # Fit model with fixed dispersion and no zero-inflation
    model <- glmmTMB(formula, 
                     family = betabinomial(), 
                     dispformula = ~1, 
                     ziformula = ~0, 
                     data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    aic <- AIC(model)
    
    if (aic < best_aic) {
        best_aic <- aic
        best_model_wtITRRES_naive <- formula
    }

    
    # Store all models in the DataFrame
    saved_models_wtITRRES_naive_df <- rbind(saved_models_wtITRRES_naive_df, 
                              data.frame(Model = as.character(formula), AIC = aic))
    
    # Increment model counter
    model_counter <- model_counter + 1
    
    # Print progress (model combination number and its AIC)
    cat("Tested model:", model_counter, "AIC:", aic, "with predictors:", predictors, "\n")
  }
}


# Filter the DataFrame to keep only models with AIC <= best AIC + 2
qualified_models_wtITRRES_naive_df <- saved_models_wtITRRES_naive_df[saved_models_wtITRRES_naive_df$AIC <= best_aic + 2, ]

# Remove rows where any value starts with "~" or "cbind"
qualified_models_wtITRRES_naive_df_cleaned <- qualified_models_wtITRRES_naive_df[!apply(qualified_models_wtITRRES_naive_df, 1, function(row) {
  any(grepl("^~", row) | grepl("^cbind", row))
}), ]

write.csv(qualified_models_wtITRRES_naive_df_cleaned,"qualified_models_wtITRRES_naive_df.csv", row.names = FALSE)
qualified_models_wtITRRES_naive_df_cleaned <- read.csv("qualified_models_wtITRRES_naive_df.csv")


# Initialize a list to store the actual model objects of qualified models
qualified_model_wtITRRES_naive_df_list <- list()

# Second loop: Generate only the qualified models
for (i in 1:nrow(qualified_models_wtITRRES_naive_df_cleaned)) {
    # Recreate the model formula from the saved string
    model_formula <- as.formula(paste("cbind(wt, Tot.ITR - wt) ~", qualified_models_wtITRRES_naive_df_cleaned$Model[i]))
    
    # Fit the model again (only for qualified models)
    model <- glmmTMB(model_formula, family = betabinomial(), dispformula = ~1, ziformula = ~0, data = Schimmelradar360spat_focal_TR_ITRA.df)
    
    # Store the qualified model in the list
    qualified_model_wtITRRES_naive_df_list[[i]] <- model
}

#generate avaraged model
AVGwtITRRESmodels_naive <- model.avg(qualified_model_wtITRRES_naive_df_list)
Summary_AVGwtITRRESmodels_naive <- summary(AVGwtITRRESmodels_naive)

# Generate output table of avaraged model model
Sum_of_weights <- as.data.frame(Summary_AVGwtITRRESmodels_naive$sw)
colnames(Sum_of_weights) <- "Sum of Weights"
Conditional_summary_wtITRRES_naive_model <- merge(Summary_AVGwtITRRESmodels_naive$coefmat.subset[,-2], Sum_of_weights, by = "row.names")
Conditional_summary_wtITRRES_naive_model <- Conditional_summary_wtITRRES_naive_model[match(rownames(Summary_AVGwtITRRESmodels_naive$coefmat.subset[,]), Conditional_summary_wtITRRES_naive_model$Row.names), ]
Conditional_summary_wtITRRES_naive_model[1,] <- c(as.character("Intercept"), Summary_AVGwtITRRESmodels_naive$coefmat.subset[1,-2], 1)

Conditional_summary_wtITRRES_naive_model <- Conditional_summary_wtITRRES_naive_model[,-1]
rownames(Conditional_summary_wtITRRES_naive_model) <- c("Intercept" ,"Flowerbulbs" , "Grassland",   "Grains" ,  "Tree Farming",  "Forest" , "Fruit", "Urbanisation", "Onions" , "Greenhouses" , "Potatoes", "Water", "Maize") 

#Turn dataframe numeric again
Conditional_summary_wtITRRES_naive_model[] <- lapply(Conditional_summary_wtITRRES_naive_model, function(x) {
  if (is.character(x)) as.numeric(x) else x
})

Output_summary_wtITRRES_naive_model  <- signif(Conditional_summary_wtITRRES_naive_model, digits = 3)

# Check if the file already exists
if(file.exists("Output_summary_wtITRRES_naive_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_summary_wtITRRES_naive_model.csv")
}
write.csv(Output_summary_wtITRRES_naive_model, "Output_summary_wtITRRES_naive_model.csv", row.names = T, quote = F)

# use best fitting model for plotting prediction
best_betamodelWTITRRES <- glmmTMB(cbind(wt, Tot.ITR - wt) ~ bulb22fra, family = betabinomial, data = Schimmelradar360spat_focal_TR_ITRA.df, dispformula = ~1, ziformula = ~0)

Output_best_WTITRRES_model <- as.data.frame(signif(summary(best_betamodelWTITRRES)$ coefficients$cond, digits = 3))

# Check if the file already exists
if(file.exists("Output_best_WTITRRES_model.csv")) {
  # If the file exists, remove it
  file.remove("Output_best_WTITRRES_model.csv")
}
write.csv(Output_best_WTITRRES_model, "Output_best_WTITRRES_model.csv", row.names = T, quote = F)
```

#K-fold cross validation of itraconazole resistant WT fraction model
```{r}
# 5(k)-fold cross-validation
set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){
  
  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(wt, Tot.ITR - wt) ~ 1,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$WTITRfrac - prd)

 
}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_NULL_RMSE_WTITRRES  <- sqrt(mean(res^2))
five_fold_NULL_RMSE_WTITRRES  <- 0.01323007
```

#Leave one out cross validation of itraconazole resistant wt fraction model
```{r}
# 5(k)-fold cross-validation

set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal_TR_ITRA.df))


for(i in 0:4){

 

  ids <- idx[(i*71 + 1):((i+1)*71)]

 

  train <- Schimmelradar360spat_focal_TR_ITRA.df[-ids,]

 

  test  <- Schimmelradar360spat_focal_TR_ITRA.df[ids,]

 

  mdl <- glmmTMB(cbind(wt, Tot.ITR - wt) ~ bulb22fra,

             family = betabinomial(),

             data = train)

 

  prd <- predict(mdl, test, type="response")

 

res <- c(res, test$WTITRfrac - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_RMSE_WTITRRES   <- sqrt(mean(res^2))
five_fold_RMSE_WTITRRES   <- 0.01320854
```
