---
title: "Script_4_Spatial_model_phenotypic_resistance_data"
author: "Hylke Kortenbosch"
date: "2024-09-26"
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
  "MASS",        # For stepwise model selection
  "stars",       # For converting coordinate systems
  "tidyverse",   # For data manipulation
  "terra",       # For making spatial raster
  "grid",        # For gridded plotting
  "gridExtra",   # For arranging 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)
```


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


#How many samples are negative for the different types of phenotypic resistance?
```{r}
Schimmelradar360spat_focal$RF.ITR
Schimmelradar360spat_focal$RF.VOR

# Count zeros in each column
num_zeros_RF_VOR <- sum(Schimmelradar360spat_focal$RF.VOR == 0, na.rm = TRUE)
num_zeros_RF_ITR <- sum(Schimmelradar360spat_focal$RF.ITR == 0, na.rm = TRUE)

# Count rows where both columns are zero
num_both_zero <- sum(Schimmelradar360spat_focal$RF.VOR == 0 & Schimmelradar360spat_focal$RF.ITR == 0, na.rm = TRUE)
```

#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 resistance fraction points on a map (Figure 1B)
```{r}
#Set coordinates to amersfoort system for plotting
Schimmelradar360spat_focal_amersfoort <- st_transform(Schimmelradar360spat_focal, 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 custom breaks for the color scale

true_breaks <- c(0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.30)  # Adjust as needed

letter_grob <- textGrob("B", x = 0.02, y = 0.98, 
                        hjust = 1, vjust = 0, 
                        gp = gpar(fontsize = 5, fontface = "bold"))  # Reduced fontsize

# Updated plot for Itraconazole RF with alpha transparency and no legend
Itraconazole_RF_points_NL_alpha <- ggplot() + 
  geom_sf(data = provincesNL_transformed, fill = "white", linewidth = 0.05) +  # Reduced line width
  coord_sf(datum = st_crs(28992)) +
  geom_point(
    data = as.data.frame(Schimmelradar360spat_focal_amersfoort), 
    aes(x = long, y = lat, alpha = RF.ITR),  
    color = "black",  # Black outline for points
    size = 0.2  # Reduced point size
  ) + 
  # Fix alpha scale with adjusted breaks
  scale_alpha_continuous(
    name = "Fraction",  
    breaks = true_breaks,  
    labels = true_breaks,  
    range = c(0.05, 1),  # Ensure useful alpha range
    guide = guide_legend(
      reverse = TRUE, 
      keywidth = unit(0.15, "cm"),  
      keyheight = unit(0.15, "cm"),
    limits = c(0, 0.35) # Ensure one legend will work for all res plots
    )
  ) +
  labs(
    title = "Itraconazole Resistance",
    x = NULL,  # Remove x-axis label
    y = NULL   # Remove y-axis label
  ) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  # Maintain aspect ratio
    plot.title = element_text(size = 5, face = "bold", hjust = 0.5),  # Reduced title size
    axis.title = element_blank(),  # Remove axis titles
    axis.text = element_blank(),   # Remove axis text
    axis.ticks = element_blank(), 
    legend.position = "none",  # Remove the legend   
     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)


# Display the updated plot
Itraconazole_RF_points_NL_alpha


ggsave("Itraconazole_resistance_points_NL.pdf", plot = Itraconazole_RF_points_NL_alpha, width = 7, height = 7, dpi = 300)
```

Model selection for spatial Itraconazole resistance land-use model (Tables S6 and S7)
```{r}
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_ITRRES_naive <- NULL
model_counter <- 0  # Initialize counter for number of models run

# Initialize an empty DataFrame to store all models
saved_models_ITRRES_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(RES.ITR, Tot.ITR - RES.ITR) ~", 
                                 paste(predictors, collapse = " + ")))
    # Fit model with fixed dispersion and no zero-inflation
    model <- glmmTMB(formula, 
                     family = betabinomial(), 
                     dispformula = ~1, 
                     ziformula = ~0, 
                     data = as.data.frame(Schimmelradar360spat_focal))
    
    aic <- AIC(model)
    
    if (aic < best_aic) {
        best_aic <- aic
        best_model_ITRRES_naive <- formula
    }

    
    # Store all models in the DataFrame
    saved_models_ITRRES_naive_df <- rbind(saved_models_ITRRES_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_ITRES_naive_df <- saved_models_ITRRES_naive_df[saved_models_ITRRES_naive_df$AIC <= best_aic + 2, ]

write.csv(qualified_models_ITRES_naive_df,"qualified_models_ITRRES_naive_df.csv", row.names = FALSE)
qualified_models_ITRRES_naive_df_cleaned <- read.csv("qualified_models_ITRRES_naive_df.csv")


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

# Second loop: Generate only the qualified models
for (i in 1:nrow(qualified_models_ITRRES_naive_df_cleaned)) {
    # Recreate the model formula from the saved string
    model_formula <- as.formula(paste("cbind(RES.ITR, Tot.ITR - RES.ITR) ~", qualified_models_ITRRES_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)
    
    # Store the qualified model in the list
    qualified_model_ITRRES_naive_df_list[[i]] <- model
}

#generate avaraged model
AVGITRRESmodels_naive <- model.avg(qualified_model_ITRRES_naive_df_list)
Summary_AVGITRRESmodels_naive <- summary(AVGITRRESmodels_naive)

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


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

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

Output_summary_ITRRES_naive_model  <- signif(Conditional_summary_ITRRES_naive_model, digits = 3)

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


null_betamodelITRRES <- glmmTMB(cbind(RES.ITR, Tot.ITR - RES.ITR) ~ 1, family = betabinomial, data = Schimmelradar360spat_focal.df, dispformula = ~1, ziformula = ~0)

# use best fitting model for plotting prediction
best_betamodelITRRES <- glmmTMB(cbind(RES.ITR, Tot.ITR - RES.ITR) ~ bulb22fra + forest22fra + maize22fra + Onion22fra, family = betabinomial, data = Schimmelradar360spat_focal.df, dispformula = ~1, ziformula = ~0)
check_collinearity(best_betamodelITRRES) #check VIFS

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

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

Semivariogram of phenotyic itraconazole resistance model 
```{r}
vg.ITR <- variogram(resid(best_betamodelITRRES, type = "response")~1, Schimmelradar360spat_focal, boundaries = 1:20 * 5e3)
plot(vg.ITR) 

vg.NUll_ITR <- variogram(resid(null_betamodelITRRES, type = "response")~1, Schimmelradar360spat_focal, boundaries = 1:20 * 5e3)
plot(vg.ITR) # no spatial structure in the data


vg.ITR$dist_km <- vg.ITR$dist / 1000
vg.NUll_ITR$dist_km <- vg.NUll_ITR$dist / 1000

max_range <- max(max(vg.ITR$dist_km), max(vg.NUll_ITR$dist_km))


# Compute the overall range of the semi-variance values
overall_range <- range(c(0, vg.ITR$gamma, vg.NUll_ITR$gamma))
dev.off()
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.ITR$dist_km, vg.ITR$gamma, col = "blue", pch = 16)
text(vg.ITR$dist_km, vg.ITR$gamma, labels = vg.ITR$n, pos = 2, cex = 0.7, col = "blue")

# Plot the points for the second variogram
points(vg.NUll_ITR$dist_km, vg.NUll_ITR$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("ITR_variogram_plot.tiff", width = 7, height = 7, units = "in", res = 300)
dev.print(device = tiff, filename = "ITR_variogram_plot.tiff")
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 resistance through space and plotting the prediction (Figure 1H)
```{r}
Predictor_raster_NL <- rast("Predictor_raster_NL.tif")
plot(Predictor_raster_NL[[11]])
str(Predictor_raster_NL)

Itraconazole_naive_RF_NL <- terra::predict(Predictor_raster_NL, best_betamodelITRRES, filename = "Itraconazole_naive_RF_NL.tif", type = "response", overwrite = T)

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

Itraconazole_naive_RF_NL_df <- as.data.frame(Itraconazole_naive_RF_NL, xy = TRUE)

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

# Define consistent breaks for both plots (true range: 0 to 0.25)
true_breaks <- c(0, 0.02, 0.05, 0.1, 0.15)

# Letter annotation
letter_grob <- textGrob("H", x = 0.02, y = 0.98, 
                        hjust = 1, vjust = 0, 
                        gp = gpar(fontsize = 5, fontface = "bold"))

hist(Itraconazole_naive_RF_NL_df$lyr1)
Itraconazole_naive_RF_NL_df$capped_lyr1 <- pmin(Itraconazole_naive_RF_NL_df$lyr1, 0.15)

# Create the plot
Itraconazole_naive_resistance_fraction_NL <- ggplot(Itraconazole_naive_RF_NL_df, aes(x = x, y = y, fill = capped_lyr1)) +
  geom_raster() +
 # Apply custom white-to-blue gradient
  scale_fill_gradientn(
    name = "Resistance Fraction",
    colors = custom_colors,
    breaks = true_breaks,
    labels = true_breaks,
    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  # High values at the top
    )
  ) +
  labs(
    title = "Itraconazole Resistance",
    x = "Easting (meters)",
    y = "Northing (meters)"
  ) +
  scale_x_continuous(
    breaks = seq(0, max(Itraconazole_naive_RF_NL_df$x), by = 50000)
  ) +
  scale_y_continuous(
    breaks = seq(0, max(Itraconazole_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,  
    plot.title = element_text(size = 5, face = "bold", hjust = 0.5),  
    axis.title = element_blank(),  
    axis.text = element_blank(),   
    axis.ticks = element_blank(),
    legend.title = element_text(size = 3.5),  
    legend.text = element_text(size = 2.5),    
    legend.box.margin = margin(0, 0, 0, 0),  
    legend.spacing.y = unit(0.05, 'cm'),      
    legend.spacing.x = unit(0.05, 'cm'),
    legend.position = "none",  # Remove the legend   
     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)

# Display the plot
Itraconazole_naive_resistance_fraction_NL


# Save the plot without legend
ggsave("Itraconazole_resistance_fraction_no_legend.jpeg",
       plot = Itraconazole_naive_resistance_fraction_NL,
       width = 1.2,  # Main plot width
       height = 1.5, # Height adjusted for aspect ratio
       dpi = 600)

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

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


for(i in 0:4){

 

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

 

  train <- as.data.frame(Schimmelradar360spat_focal)[-ids,]

 

  test  <- as.data.frame(Schimmelradar360spat_focal)[ids,]

 

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

             family = betabinomial(),

             data = train)

 

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

 

res <- c(res, test$RF.ITR - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_NULL_RMSE_ITRRES  <- sqrt(mean(res^2))
five_fold_RMSE_ITRRES  <- 0.03768774
```

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


for(i in 0:4){

 

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

 

  train <- as.data.frame(Schimmelradar360spat_focal)[-ids,]

 

  test  <- as.data.frame(Schimmelradar360spat_focal)[ids,]

 

  mdl <- glmmTMB(cbind(RES.ITR, Tot.ITR - RES.ITR) ~ forest22fra + maize22fra +  
    Onion22fra + bulb22fra,

             family = betabinomial(),

             data = train)

 

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

 

res <- c(res, test$RF.ITR - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_RMSE_ITRRES  <- sqrt(mean(res^2))
five_fold_RMSE_ITRRES  <- 0.03578144
```

#Make plot voriconazole resistance fraction points on a map (Figure 1C)
```{r}
#Set coordinates to amersfoort system for plotting
Schimmelradar360spat_focal_amersfoort <- st_transform(Schimmelradar360spat_focal, 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 Itraconazole)
true_breaks <- c(0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.30)  

# Letter annotation for Voriconazole plot
letter_grob <- textGrob("C", x = 0.02, y = 0.98, 
                        hjust = 1, vjust = 0, 
                        gp = gpar(fontsize = 5, fontface = "bold"))  

# Use continuous alpha scaling instead of binning
Voriconazole_RF_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_amersfoort), 
    aes(x = long, y = lat, alpha = RF.VOR),  # Use actual continuous resistance values
    color = "black",  
    size = 0.2  
  ) + 
  # Apply continuous alpha scale with limits to unify the legend
  scale_alpha_continuous(
    name = "Fraction",  
    breaks = true_breaks,  
    labels = true_breaks,  
    range = c(0.05, 1),  # Ensures meaningful transparency variation
    limits = c(0, 0.35),  # Match Itraconazole’s legend for comparability
    guide = guide_legend(
      reverse = TRUE, 
      keywidth = unit(0.15, "cm"),  
      keyheight = unit(0.15, "cm")
    )
  ) +
  labs(
    title = "Voriconazole Resistance",
    x = NULL,
    y = NULL
  ) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    aspect.ratio = 1.25,  
    plot.title = element_text(size = 5, face = "bold", hjust = 0.5),  
    axis.title = element_blank(),  
    axis.text = element_blank(),   
    axis.ticks = element_blank(),
    legend.title = element_text(size = 3.5),  
    legend.text = element_text(size = 2.5),    
    legend.box.margin = margin(0, 0, 0, 0),  
    legend.spacing.y = unit(0.05, 'cm'),      
    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)

# Display the fixed plot
Voriconazole_RF_points_NL_alpha


#save plot adjust to desired format, dimensions and resolution
ggsave("Voriconazole_resistance_points_NL.pdf", plot = Voriconazole_RF_points_NL_alpha, width = 7, height = 7, dpi = 600)
```

#Model selection for spatial Voriconaozle resistance land-use model (Tables S8 and S9)
```{r}
# fit null betabinomial model
null_betamodelVORRES_naive <- glmmTMB(cbind(RES.VOR, Tot.VOR - RES.VOR) ~ 1, family = betabinomial, data = as.data.frame(Schimmelradar360spat_focal))
summary(null_betamodelVORRES)

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

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_VORRES_naive <- NULL
model_counter <- 0  # Initialize counter for number of models run

# Initialize an empty DataFrame to store all models
saved_models_VORRES_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(RES.VOR, Tot.VOR - RES.VOR) ~", 
                                 paste(predictors, collapse = " + ")))
    # Fit model with fixed dispersion and no zero-inflation
    model <- glmmTMB(formula, 
                     family = betabinomial(), 
                     dispformula = ~1, 
                     ziformula = ~0, 
                     data = Schimmelradar360spat_focal.df)
    
    aic <- AIC(model)
    
    if (aic < best_aic) {
        best_aic <- aic
        best_model_VORRES_naive <- formula
    }

    # Store all models in the DataFrame
    saved_models_VORRES_naive_df <- rbind(saved_models_VORRES_naive_df, 
                              data.frame(Model = paste(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_VORRES_naive_df <- saved_models_VORRES_naive_df[saved_models_VORRES_naive_df$AIC <= best_aic + 2, ]
View
View(qualified_models_VORRES_naive_df)

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

# Make all models with an AIC of <2 higher than the best model and summarise them


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

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

#generate avaraged model
AVGVORRESmodels_naive <- model.avg(qualified_model_VORRES_naive_df_list)
# Generate summary table of conitional avarage for use in supplement
Summary_AVGVORRESmodels_naive <- summary(AVGVORRESmodels_naive)

Sum_of_weights <- as.data.frame(Summary_AVGVORRESmodels_naive$sw)
colnames(Sum_of_weights) <- "Sum of Weights"
Conditional_summary_VORRES_naive_model <- merge(Summary_AVGVORRESmodels_naive$coefmat.subset[,-2], Sum_of_weights, by = "row.names")
Conditional_summary_VORRES_naive_model <- Conditional_summary_VORRES_naive_model[match(rownames(Summary_AVGVORRESmodels_naive$coefmat.subset[,]), Conditional_summary_VORRES_naive_model$Row.names), ]
Conditional_summary_VORRES_naive_model[1,] <- c(as.character("Intercept"), Summary_AVGVORRESmodels_naive$coefmat.subset[1,-2], 1)


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

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

Output_summary_VORRES_naive_model  <- signif(Conditional_summary_VORRES_naive_model, digits = 3)

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

best_betamodelVORRES <- glmmTMB(cbind(RES.VOR, Tot.VOR - RES.VOR) ~ bulb22fra + buildings22fra + grass22fra + greenhouse22fra +  maize22fra + potato22fra  +  water22fra, family = betabinomial, data = Schimmelradar360spat_focal.df, dispformula = ~1, ziformula = ~0)
summary(best_betamodelVORRES)
check_collinearity(best_betamodelVORRES) #check VIFS

#Generate and export summary table best VORRES model (LGN + Onions)
Output_best_VORRES_naive_model <- as.data.frame(signif(summary(best_betamodelVORRES)$ coefficients$cond, digits = 3))

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

Semivariogram of residuals of voriconazole resistance 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_VOR <- variogram(resid(null_betamodelVORRES_naive, type = "response")~1, Schimmelradar360spat_focal, boundaries = 1:20 * 5e3)
plot(vg.NULL_VOR, plot.numbers = F)

vg.VOR <- variogram(resid(best_betamodelVORRES, type = "response")~1, Schimmelradar360spat_focal, boundaries = 1:20 * 5e3)
plot(vg.VOR, plot.numbers = F)

vg.VOR$dist_km <- vg.VOR$dist / 1000
vg.NULL_VOR$dist_km <- vg.NULL_VOR$dist / 1000

max_range <- max(max(vg.VOR$dist_km), max(vg.NULL_VOR$dist_km))


# Compute the overall range of the semi-variance values
overall_range <- range(c(0, vg.VOR$gamma, vg.NULL_VOR$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", cex.main = 1)
# Plot the points for the first variogram
points(vg.VOR$dist_km, vg.VOR$gamma, col = "blue", pch = 16)
text(vg.VOR$dist_km, vg.VOR$gamma, labels = vg.VOR$n, pos = 2, cex = 0.7, col = "blue")
# Plot the points for the second variogram
points(vg.NULL_VOR$dist_km, vg.NULL_VOR$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("VOR_variogram_plot.tiff", width = 7, height = 7, units = "in", res = 300)
dev.off()

# There is no spatial structure in the unexplained variance and therefore we cannot fit a meaningful semivariogram model
```
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 Voriconazole resistance through space and plotting the prediction (Figure 1I)
```{r}
Predictor_raster_NL <- rast("Predictor_raster_NL.tif")
summary(Naive_selected_modelVORRES)
Voriconazole_naive_RF_NL <- terra::predict(Predictor_raster_NL, Naive_selected_modelVORRES, filename = "Voriconazole_naive_RF_NL.tif", type = "response", overwrite = T)

Voriconazole_naive_RF_NL <- rast("Voriconazole_naive_RF_NL.tif")
# convert to dataframe
Voriconazole_naive_RF_NL_df <- as.data.frame(Voriconazole_naive_RF_NL, xy = TRUE)
# Define meaningful breaks for the legend (true range: 0 to 0.25)
true_breaks <- c(0, 0.02, 0.05, 0.1, 0.15)

# Keep the values continuous
Voriconazole_naive_RF_NL_df$lyr1_transformed <- Voriconazole_naive_RF_NL_df$lyr1

str(Voriconazole_naive_RF_NL_df$lyr1)

Voriconazole_naive_RF_NL_df$capped_lyr1 <- pmin(Voriconazole_naive_RF_NL_df$lyr1, 0.15)

# Letter annotation
letter_grob <- textGrob("I", x = 0.02, y = 0.98, 
                        hjust = 1, vjust = 0, 
                        gp = gpar(fontsize = 5, fontface = "bold"))

# Updated plot with continuous scale
Voriconazole_naive_resistance_fraction_NL <- ggplot(Voriconazole_naive_RF_NL_df, aes(x = x, y = y, fill = capped_lyr1)) +
  geom_raster() +
 # Apply custom white-to-blue gradient
  scale_fill_gradientn(
    name = "Fraction",
    colors = custom_colors,
    breaks = true_breaks,
    labels = true_breaks,
    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(
    title = "Voriconazole Resistance",
    x = "Easting (meters)",
    y = "Northing (meters)"
  ) +
  scale_x_continuous(
    breaks = seq(0, max(Voriconazole_naive_RF_NL_df$x), by = 50000)
  ) +
  scale_y_continuous(
    breaks = seq(0, max(Voriconazole_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,  
    plot.title = element_text(size = 5, face = "bold", hjust = 0.5),  
    axis.title = element_blank(),  
    axis.text = element_blank(),   
    axis.ticks = element_blank(),
    legend.title = element_text(size = 3.5),  
    legend.text = element_text(size = 2.5),    
    legend.box.margin = margin(0, 0, 0, 0),  
    legend.spacing.y = unit(0.05, 'cm'),      
    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)

# Display the plot
Voriconazole_naive_resistance_fraction_NL

# Save the plot without legend
ggsave("Voriconazole_resistance_fraction_no_legend.jpeg",
       plot = Voriconazole_naive_resistance_fraction_NL,
       width = 1.2,  # Main plot width
       height = 1.5, # Height adjusted for aspect ratio
       dpi = 600)


ggsave("Voriconazole_naive_resistance_fraction_NL_blackbg.pdf", plot = Voriconazole_naive_resistance_fraction_NL, width = 7, height = 7, dpi = 300)
```

#Cross validation of voriconazole resistance fraction nullmodel
```{r}
# 5(k)-fold cross-validation

set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal))


for(i in 0:4){

 

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

 

  train <- as.data.frame(Schimmelradar360spat_focal)[-ids,]

 

  test  <- as.data.frame(Schimmelradar360spat_focal)[ids,]

 

  mdl <- glmmTMB(cbind(RES.VOR, Tot.VOR - RES.VOR) ~ 1,

             family = betabinomial(),

             data = train)


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

res <- c(res, test$RF.VOR - prd)

 
}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_null_RMSE_VORRES  <- sqrt(mean(res^2))
five_fold_null_RMSE_VORRES <- 0.03890009
```

#Cross validation of best voriconazole resistance fraction model
```{r}
# 5(k)-fold cross-validation

set.seed(123) # make reproducible

res <- c()


idx <- sample(1:nrow(Schimmelradar360spat_focal))


for(i in 0:4){

 

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

 

  train <- as.data.frame(Schimmelradar360spat_focal)[-ids,]

 

  test  <- as.data.frame(Schimmelradar360spat_focal)[ids,]

 

  mdl <- glmmTMB(cbind(RES.VOR, Tot.VOR - RES.VOR) ~ buildings22fra + maize22fra +  
    bulb22fra + greenhouse22fra + water22fra + potato22fra +  grass22fra,

             family = betabinomial(),

             data = train)

 

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

 

res <- c(res, test$RF.VOR - prd)

 

}

min(res) 
max(res)

ME    <- mean(res)

MAE   <- mean(abs(res))

five_fold_RMSE_VORRES  <- sqrt(mean(res^2))
five_fold_RMSE_VORRES  <- 0.03509
```

# Neutral map of all sample locations (Used in Figure 1F)
```{r}
#Set coordinates to amersfoort system for plotting
Schimmelradar360spat_focal_amersfoort <- st_transform(Schimmelradar360spat_focal, 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)


Sample_points_NL_neutral <-  ggplot() + 
  geom_sf(data = provincesNL_transformed, fill = "white", linewidth = 0.2) +  # Reduced line width
  coord_sf(datum = st_crs(28992)) +
  geom_point(
    data = as.data.frame(Schimmelradar360spat_focal_amersfoort), 
    aes(x = long, y = lat),  # Removed color aesthetic
    color = "darkred",  # Black dots
    size = 0.2  # Reduced point size
  ) + 
  labs(
    x = NULL,  # Removed x-axis title
    y = NULL   # Removed y-axis title
  ) +
  theme(
    plot.margin = margin(0, 0, 0, 0),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
     panel.background = element_rect(fill = "white", color = NA),  # Set panel background to white
    plot.background = element_rect(fill = "white", color = NA),   # Set plot background to white
    aspect.ratio = 1.25,  # Maintain aspect ratio
    plot.title = element_blank(),  # Removed plot title
    axis.title = element_blank(),  # Removed axis titles
    axis.text = element_blank(),   # Removed axis text
    axis.ticks = element_blank(),  # Removed axis ticks
    legend.position = "none"       # Removed legend
  )

# Save the neutral map as a PDF
ggsave("Sample_points_NL_neutral.svg", 
       plot = Sample_points_NL_neutral, 
       width = 2, height = 2, dpi = 600)
```

Generate lollipop figure of voriconazole resistance summary model (Figure 2)
```{r}
str(Output_summary_VORRES_naive_model)
Output_summary_VORRES_naive_model_lolly <- Output_summary_VORRES_naive_model[-1,]



Output_summary_VORRES_naive_model_lolly$Variable <- rownames(Output_summary_VORRES_naive_model[-1,])
Output_summary_VORRES_naive_model_lolly$Variable[2] <- "Urbanization"

# Add significance labels
Output_summary_VORRES_naive_model_lolly$Significance <- ifelse(
  Output_summary_VORRES_naive_model_lolly$`Pr(>|z|)` < 0.001, "***",
  ifelse(Output_summary_VORRES_naive_model_lolly$`Pr(>|z|)` < 0.01, "**",
         ifelse(Output_summary_VORRES_naive_model_lolly$`Pr(>|z|)` < 0.05, "*", ""))
)

Output_summary_VORRES_naive_model_lolly_sorted <- Output_summary_VORRES_naive_model_lolly[order(-Output_summary_VORRES_naive_model_lolly$Estimate), ]

# Reorder Variable factor by Estimate in descending order
Output_summary_VORRES_naive_model_lolly_sorted$Variable <- factor(
  Output_summary_VORRES_naive_model_lolly_sorted$Variable,
  levels = rev(Output_summary_VORRES_naive_model_lolly_sorted$Variable[order(-Output_summary_VORRES_naive_model_lolly_sorted$Estimate)]
))


# Create the lollipop plot
Summary_VORRES_naive_model_lolly_plot <- ggplot(Output_summary_VORRES_naive_model_lolly_sorted, aes(x = Variable, y = Estimate)) +
  # Lollipop sticks
  geom_segment(aes(x = Variable, xend = Variable, y = 0, yend = Estimate), color = "black", linewidth = 1) +  
  # Lollipop heads with size and color linked to `Sum of Weights`
  geom_point(aes(size = `Sum of Weights`, color = `Sum of Weights`)) + 
  # Significance text
  geom_text(
    aes(label = Significance), 
    size = 5, 
    color = "red", 
    nudge_x = 0.4  # Adjust position based on estimate
  ) +
  # Horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 2) +
  # Set y-axis breaks in increments of 5
  scale_y_continuous(
    name = "Effect size",
    breaks = seq(-5, max(Output_summary_VORRES_naive_model_lolly_sorted$Estimate, na.rm = TRUE), by = 5)
  ) +
  # Viridis color scale for Sum of Weights (linked to size as well)
  scale_color_viridis_c(name = "Sum of Weights", option = "viridis") +
  scale_size_continuous(name = "Sum of Weights") +
  # Combine legends for size and color
  guides(
    color = guide_legend(order = 1, reverse = T),
    size = guide_legend(order = 1, reverse = T)
  ) +
  # Minimal theme and axis customizations
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major = element_blank(),   # Remove major grid lines
    panel.grid.minor = element_blank(),   # Remove minor grid lines
    axis.text.x = element_text(size = 12), # No angle for labels
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    legend.position = "right"  # Unified legend on the right
  ) +
  # Labels for axes
  labs(
    x = "Land-use classes",
    y = "Estimate"
  ) +
  # Flip the coordinates
  coord_flip()

Summary_VORRES_naive_model_lolly_plot

# Save the plot without legend
ggsave("Summary_VORRES_model_lolly_plot.pdf",
       plot = Summary_VORRES_naive_model_lolly_plot,
       width = 7,  # Main plot width
       height = 5, # Height adjusted for aspect ratio
       dpi = 600)
```




Combining raw data and prediction plots into 2x2 grid
```{r}
Figure_naive_resistance <- grid.arrange(Itraconazole_RF_points_NL, Itraconazole_naive_resistance_fraction_NL, Voriconazole_RF_points_NL, Voriconazole_naive_resistance_fraction_NL,  nrow = 2, ncol = 2)

ggsave("Figure_naive_resistance.pdf", plot = Figure_naive_resistance, width = 6, height = 4, dpi = 600)
ggsave("Figure_naive_resistance.jpeg", plot = Figure_naive_resistance, width = 6, height = 4, dpi = 300)
```

