---
title: "R_code: manuscript DOI: ..."
author: "Egle Liiskmann, Douglas B. Sponsler, Marika Mänd, Andrea Holzschuh, Jochen Krauss, Ingolf Steffan-Dewenter, Oliver Schweiger, Reet Karise"
date: "2025-09-23"
output: html_document
---
# Packages
```{r}
library(glmmTMB) # model building
library(performance) # model evaluation
library(DHARMa) # model evaluation
library(marginaleffects) # marginal/conditional effects
library(emmeans)
library(ggeffects)
library(sf)
library(dplyr)
library(see) # ggplot themes
library(dagitty) # DAG consistency check
library(ggdag)
library(lavaan)
library(tidyverse) # data handling and plotting
library(readr)
theme_set(theme_lucid())
```

# Load data
```{r}
# Read data
data_WP <- read_rds("../GECCO_files/data_WP.rds")
data_WP_HB <- read_rds("../GECCO_files/data_WP_HB.rds")
```

# DAG
## Causal claims

Our DAG includes the following causal paths:   

1. **forest % (within 2 km radius) -> pollinator richness/abundance (Y):** Forest cover "directly" (within the scope of the DAG) affects pollinator abundance/diversity through its influence on nesting habitat and microclimate    

2. **forest % (within 2 km radius) -> floral resources (landscape-scale):** Forest cover affects landscape-scale floral resources by influencing floristics.

3. **forest % (within 2 km radius) -> hive density:** Forest cover affects hive density because beekeepers usually prefer to place colonies in open areas.

4. **hive density -> floral resources (transect-level):** Hive density influences floral resources availability and seed production through honey bee foraging.

5. **hive density -> floral resource depletion:** Hive density influences floral resource depletion rate through honey bee foraging.

6. **floral resources (landscape-scale) -> floral resources (transect-scale):** Floral resources in landscape influence local floral resources by seed dispersal.

7. **floral resources (transect-scale) -> floral resource depletion:** Landscape-scale floral resources influence floral resource depletion rate by determining floral resource production rate.

8. **floral resource depletion -> pollinator richness/abundance (Y):** Floral resource depletion influences pollinator richness/abundance by modulating competition.   

9. **hive density -> pollinator richness/abundance (Y):** Hive density influences pollinator richness/abundance through disease transmission.

10. **restoration status -> pollinator richness/abundance (Y):** Restoration status (i.e. the restoration of grassland habitat by the removal of forest) influences pollinator richness/abundance through microclimate and nesting substrate.

11.**restoration status -> forest %:** Restoration status influences forest % because it entails the replacement of forest with grassland.

12.**restoration status -> floral resources (landscape):** Restoration status influences landscape-scale floral resources by influencing floristics.

13.**restoration status -> floral resources (transect):** Restoration status influences transect-scale floral resources by influencing floristics.

## Adjustment set

Our DAG includes one causal fork that could bias our estimate of the causal effect of hive density on pollinator richness/abundance: `Y <-- forest % --> hive density`. To obtain an unbiased estimate of the causal effect of vehicle density on pollinator richness/abundance, we need to condition on `forest %`. This adjustment set yields a partial (i.e. unmediated by hive density) causal effect of forest %.

In a second model, we can also obtain an unbiased estimate for the causal effect of restoration status on pollinator richness/abundance by including only restoration status as an explanatory variable.


# Validate DAG

Our DAG implies that hive density should be independent of restoration status, conditional on forest %.

## Visual check

Interestingly, it looks like there are fewer hives in restored landscapes, suggesting that restoration status should be considered a cause of hive density. Moreover, there is no obvious graphical evidence for a causal effect of restoration status on forest %, which is odd given that the restoration process was literally the removal of forest. But maybe the scale of restoration was small relative to the 2 km radius over which we're considering forest cover. There is also no obvious correlation between forest cover and hive density. 

Our conclusion is that the safest decision is to simply add a causal arrow from restoration status to hive density, and therefore condition our model on forest % *and* restoration status to obtain an unbiased effect of hive density on pollinator richness/abundance.

```{r}
ggplot(data_WP, aes(restored, hives_2km)) +
  geom_violin()

ggplot(data_WP, aes(restored, forest_perc)) +
  geom_violin()

ggplot(data_WP, aes(hives_2km, forest_perc)) +
  geom_point()


```

## Call DAG data

For technical reasons in 'DAgitty', we have to convert `restored` to binary.

```{r message=FALSE}
dag_data <- read_rds("../GECCO_files/data_WP.rds")%>%
  mutate(restored = as.integer(if_else(restored == "restored", 1, 0))) %>%
  dplyr::select(pdiv = prop_species_richness, pdens = pollinator_density, 
         flora = plant_species_richness, flower = flower_cover, rest = restored, 
         hd = hives_2km, forest = forest_perc)
```

## Specify DAG, plot
```{r}
dag <- dagitty(
   'dag {
bb="0,0,1,1"
depletion [latent,pos="0.460,0.515"]
flora [latent,pos="0.460,0.090"]
flower [exposure,pos="0.460,0.281"]
forest [adjusted,pos="0.090,0.419"]
hd [exposure,pos="0.090,0.740"]
pdiv [outcome,pos="0.460,0.740"]
rest [adjusted,pos="0.090,0.090"]
depletion -> pdiv
flora -> flower
flower -> depletion
forest -> flora
forest -> hd
forest -> pdiv
hd -> depletion
hd -> flower
hd -> pdiv
rest -> flora
rest -> flower
rest -> forest
rest -> pdiv
}'
)
plot(dag)

dag_renamed <- dagitty(
  'dag {
bb="0,0,1,1"
"Floral resource depletion" [latent,pos="0.460,0.515"]
"Floral resources (landscape)" [latent,pos="0.460,0.074"]
"Floral resources (transect)" [exposure,pos="0.460,0.281"]
"Forest percentage" [adjusted,pos="0.090,0.419"]
"Honey bee colony density" [exposure,pos="0.090,0.740"]
"Restoration status" [adjusted,pos="0.090,0.074"]
Y [outcome,pos="0.460,0.740"]
"Floral resource depletion" -> Y
"Floral resources (landscape)" -> "Floral resources (transect)"
"Floral resources (transect)" -> "Floral resource depletion"
"Forest percentage" -> "Floral resources (landscape)"
"Forest percentage" -> "Honey bee colony density"
"Forest percentage" -> Y
"Honey bee colony density" -> "Floral resource depletion"
"Honey bee colony density" -> "Floral resources (transect)"
"Honey bee colony density" -> Y
"Restoration status" -> "Floral resources (landscape)"
"Restoration status" -> "Floral resources (transect)"
"Restoration status" -> "Forest percentage"
"Restoration status" -> Y
}'
)
plot(dag_renamed)
```
## Validate DAG
```{r}
# to validate our DAG with all conditional independencies implied in our DAG and get estimates for all of them
# dsep in dagitty tests all d-separation implied conditional independencies in our DAG against our data
# by default, dsep returns a summary with only the most significant or some combined output, but not all individual tests.
# to get all estimates (test results) - list of tests with p-values, and variables involved

dag_validate <- dsep(dag,dag_data)
plot_dsep(dag_validate)
print(dag_validate)
```
# Call models

We will use a negative binomial model for pollinator density. For pollinator richness, we will model richness as a proportion of the gamma diversity across all sites, which calls for a beta model. In both models, each continuous explanatory variable will interact with pollinator group.

## mod_10_HB: abundance ~ hive * pollinator_group
```{r}
# Model for pollinator abundance in relation to hive density, forest percentage, restoration status, and floral cover on transects
mod_10_HB <- glmmTMB(pollinator_density ~ 
                    hives_2km + 
                    pollinator_group + 
                    forest_perc + 
                    restored +
                    flower_cover +
                    flower_cover:pollinator_group +
                    restored:pollinator_group +
                    hives_2km:pollinator_group + 
                    forest_perc:pollinator_group,
                  family = "nbinom2",
                  data = data_WP_HB)
```

## mod_11: richness ~ hive * pollinator_group
```{r}
# Model for wild pollinator species richness in relation to hive density, forest percentage, restoration status, and floral cover on transects
mod_11 <- glmmTMB(prop_species_richness ~ 
                    hives_2km + 
                    pollinator_group + 
                    forest_perc + 
                    restored +
                    flower_cover +
                    flower_cover:pollinator_group +
                    restored:pollinator_group +
                    hives_2km:pollinator_group + 
                    forest_perc:pollinator_group,
                  family = beta_family(link = "logit"),
                  data = data_WP)
```

# Validate models

## mod_10_HB, pollinator abundance
```{r}
mod_10_HB.residuals <- simulateResiduals(mod_10_HB)
plot(mod_10_HB.residuals)
check_model(mod_10_HB)
```

## mod_11, pollinator species richness
```{r}
mod_11.residuals <- simulateResiduals(mod_11)
plot(mod_11.residuals)
check_model(mod_11)
```

# Visualize models

## mod_10_HB, pollinator abundance

### Parameter estimates
```{r}
mod_10_HB_hives_trends <- emtrends(mod_10_HB, specs = "pollinator_group", var = "hives_2km")
plot(mod_10_HB_hives_trends)

mod_10_HB_forest_trends <- emtrends(mod_10_HB, specs = "pollinator_group", var = "forest_perc")
plot(mod_10_HB_forest_trends)
```

### Predictions
```{r}
mod_10_HB_restored_pred <- ggpredict(mod_10_HB, terms = c("pollinator_group", "restored"))
plot(mod_10_HB_restored_pred) +
      labs(x = "Pollinator group", y = "Pollinator abundance", 
         color = "Restoration status", fill = "Restoration status",
         title = NULL)

# +
#     ggtitle("pollinator density ~ restoration status")

mod_10_HB_hives_restored_pred <- ggpredict(mod_10_HB, terms = c("hives_2km", "pollinator_group", "restored"))
plot(mod_10_HB_hives_restored_pred) +
      labs(x = "Honey bee colony density", y = "Pollinator abundance", 
         color = "Pollinator group", fill = "Pollinator group",
         title = NULL) 

# +
#     ggtitle("pollinator density ~ colony density")

mod_10_HB_hives_pred <- ggpredict(mod_10_HB, terms = c("hives_2km", "pollinator_group"))
plot(mod_10_HB_hives_pred)+
    labs(x = "Honey bee colony density", y = "Pollinator abundance", 
         color = "Pollinator group", fill = "Pollinator group", title = NULL)
#+
#    ggtitle("pollinator density ~ colony density")

mod_10_HB_forest_pred <- ggpredict(mod_10_HB, terms = c("forest_perc", "pollinator_group"))
plot(mod_10_HB_forest_pred)

mod_10_HB_flower_pred <- ggpredict(mod_10_HB, terms = c("flower_cover", "pollinator_group"))
plot(mod_10_HB_flower_pred)
```

## mod_11, pollinator species richness

### Parameter estimates
```{r}
mod_11_hives_trends <- emtrends(mod_11, specs = "pollinator_group", var = "hives_2km")
plot(mod_11_hives_trends)

mod_11_forest_trends <- emtrends(mod_11, specs = "pollinator_group", var = "forest_perc")
plot(mod_11_forest_trends)
```

### Predictions
```{r}
mod_11_restored_pred <- ggpredict(mod_11, terms = c("pollinator_group", "restored"))
plot(mod_11_restored_pred) +
      labs(x = "Pollinator group", y = "Wild pollinator species richness", 
         color = "Restoration status", fill = "Restoration status",
         title = NULL)

# +
#     ggtitle("pollinator density ~ restoration status")

mod_11_hives_restored_pred <- ggpredict(mod_11, terms = c("hives_2km", "pollinator_group", "restored"))
plot(mod_11_hives_restored_pred) +
      labs(x = "Honey bee colony density", y = "Wild pollinator species richness", 
         color = "Pollinator group", fill = "Pollinator group",
         title = NULL) 

# +
#     ggtitle("pollinator density ~ colony density")

mod_11_hives_pred <- ggpredict(mod_11, terms = c("hives_2km", "pollinator_group"))
plot(mod_11_hives_pred)+
    labs(x = "Honey bee colony density", y = "Wild pollinator species richness", 
         color = "Pollinator group", fill = "Pollinator group",
         title = NULL) 
  
  
  #  ggtitle("pollinator species richness ~ colony density")

mod_11_forest_pred <- ggpredict(mod_11, terms = c("forest_perc", "pollinator_group"))
plot(mod_11_forest_pred)

mod_11_flower_pred <- ggpredict(mod_11, terms = c("flower_cover", "pollinator_group"))
plot(mod_11_flower_pred)
```

# Summaries of models

##  mod_10_HB, pollinator abundance
```{r}
summary(mod_10_HB)
```

## mod_11, pollinator species richness
```{r}
summary(mod_11)
```