library(openxlsx)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(vegan)
library(corrplot)
library(vcd)
library(car)
library(MASS)
library(performance)
library(DHARMa)
library(patchwork)
library(writexl)Effect of Pavement on the Number of Wildlife Roadkills - General
1 Data Collection Context
Brazilian legislation, through the environmental licensing process, requires a series of studies to be conducted for different types of infrastructure projects, such as highways. In this context, studies are often focused on understanding the spatial patterns of wildlife-vehicle collision hotspots, with the aim of guiding the implementation of mitigation measures.
Within this framework, we recorded wildlife impacted by vehicle collisions along the BR-135 highway, located in western Bahia and northern Minas Gerais, in a transitional area between the Cerrado and Caatinga biomes.
Over a period of five years, we conducted monthly surveys using a vehicle traveling at speeds between 40–60 km/h, covering nearly 370 km of roadway. This stretch includes both paved and unpaved segments. As a result, the dataset includes records of amphibians, birds, reptiles, and mammals affected by vehicle collisions.
To account for how landscape composition might influence roadkill events, we extracted land use and land cover data from the MapBiomas platform (https://brasil.mapbiomas.org/) at three spatial scales from the road axis: 1 km, 3 km, and 5 km. Additionally, traffic volume data were incorporated—field-collected data were available for one segment, while for the remaining segments, traffic volume estimates were obtained through models provided by the Brazilian National Department of Transport Infrastructure (DNIT).
Using the spatialized information, we clipped the polygons derived at the 1 km scale so that each polygon had an area of 2 km². For the 3 km scale, each polygon covered an area of 6 km², and for the 5 km scale, each polygon covered 10 km².
2 Goals
The objective of this document is to prepare the dataset and subsequently investigate the effect of road pavement on the number of wildlife roadkills.
1) Is there an effect of pavement on the number of wildlife roadkills?
Prediction: Considering that pavement allows vehicles to travel at higher speeds and is also associated with increased traffic volume, we expect a higher number of wildlife-vehicle colision (WVC) in paved areas.
2) Does the influence of pavement on the number of roadkill records differ among taxonomic groups?
Prediction: Considering that pavement affects vehicle speed, we expect a more pronounced increase in roadkills of more mobile groups, such as birds and mammals, compared to amphibians and reptiles.
Note
Note: the analyses for this objective are in another script.
3) Does the variation in the spatial scale used to assess landscape components (1 km, 3 km, and 5 km) reflect differences in their relationship with the number of roadkills?
Prediction: At larger spatial scales, considering that the total amount of habitat tends to increase, we expect landscape components to have a stronger influence on the number of roadkills. This influence is expected to occur in both the overall roadkill analyses and those conducted separately for each taxonomic group.
3 Description of the Data Structure
This dataset contains information on roadkill occurrences across different road segments, categorized by taxonomic group, seasonality, and various environmental and infrastructural variables. The dataset includes the following columns:
| Variable | Description |
|---|---|
| polygon | Identifier for the road and landscape segment. |
| taxa | Taxonomic group of the affected species (e.g., mammals, reptiles, amphibians, birds). |
| season | Season in which the observation was recorded (dry or rainy). |
| roadkill count | Number of roadkill occurrences recorded in the segment/taxa/season. |
| road length | Length of the road segment (in kilometers). |
| relative rate | A measure of roadkill occurrences relative to the road length. |
| cult land | Proportion of cultivated land within the segment’s surrounding area. |
| water | Proportion of water bodies within the segment’s surrounding area. |
| shrubby | Proportion of shrubby vegetation within the segment’s surrounding area. |
| pasture | Proportion of pastureland within the segment’s surrounding area. |
| urban area | Proportion of urbanized area within the segment’s surrounding area. |
| forest | Proportion of forested land within the segment’s surrounding area. |
| savanna | Proportion of savanna within the segment’s surrounding area. |
| pavement | Proportion of paved road surface. |
| traffic | A measure of traffic intensity in the road segment. |
4 Packages
5 Datasets
d1km_raw <- read.xlsx("data_1km.xlsx")
d3km_raw <- read.xlsx("data_3km.xlsx")
d5km_raw <- read.xlsx("data_5km.xlsx")5.1 Data organization
5.1.1 General data (without distinguishing between taxa)
Combine the records of roadkill regardless of taxa and season; therefore, we will group by polygon. The roadkill count should be summarized, and for the other columns, we will keep only the first row of each polygon, as as the others are identical copies of the first.
datasets_raw <- list("1km" = d1km_raw, "3km" = d3km_raw, "5km" = d5km_raw)
summarize_dataset <- function(data) {
data %>%
group_by(polygon) %>%
summarise(
total_roadkill_count = sum(roadkill_count),
road_length = first(road_length),
relative_rate = sum (relative_rate),
cult_land = first(cult_land),
water = first(water),
shrubby = first(shrubby),
pasture = first(pasture),
urban_area = first(urban_area),
forest = first(forest),
savanna = first(savanna),
pavement = first(pavement),
traffic = first(traffic),
.groups = "drop"
)
}
sum_datasets <- lapply(datasets_raw, summarize_dataset)
d1km_raw_sum <- sum_datasets[[1]]
d3km_raw_sum <- sum_datasets[[2]]
d5km_raw_sum <- sum_datasets[[3]]5.1.2 Pavement as a binary variable
Since the pavement variable is in proportion but is mostly a binary variable, we will exclude polygons with intermediate values (where the road segment is partially paved and partially unpaved). I want to keep only the rows where the values are exactly 0 or 1, and exclude all rows with intermediate values, i.e., values greater than 0 and less than 1.
d1km_raw_sum <- d1km_raw_sum %>% filter(pavement %in% c(0, 1))
d3km_raw_sum <- d3km_raw_sum %>% filter(pavement %in% c(0, 1))
d5km_raw_sum <- d5km_raw_sum %>% filter(pavement %in% c(0, 1))checking the values
all(d1km_raw_sum$pavement %in% c(0, 1)) [1] TRUEall(d3km_raw_sum$pavement %in% c(0, 1)) [1] TRUEall(d5km_raw_sum$pavement %in% c(0, 1))[1] TRUEnrow(d1km_raw_sum)[1] 370nrow(d3km_raw_sum)[1] 367nrow(d5km_raw_sum)[1] 3645.2 Data Exploration
Now that the data has been reorganized, let’s perform some exploratory analyses and transformations.
Let’s create histograms of the variables to get an idea of their distribution.
combined_raw <- bind_rows(
d1km_raw_sum %>% mutate(scale = "1km"),
d3km_raw_sum %>% mutate(scale = "3km"),
d5km_raw_sum %>% mutate(scale = "5km")
)numeric_vars <- c("total_roadkill_count", "road_length","cult_land", "water", "shrubby", "pasture","urban_area", "forest",
"savanna")
long_numeric <- combined_raw |>
tidyr::pivot_longer(cols = all_of(numeric_vars),
names_to = "variable", values_to = "value") |>
dplyr::select(scale, variable, value)
plot_histograms <- function(scale_level) {
color_map <- c("1km" = "#FFD700", "3km" = "#1f77b4", "5km" = "#9B59B6")
long_numeric %>%
filter(scale == scale_level) %>%
ggplot(aes(x = value)) +
geom_bar(data = ~ filter(., variable == "traffic"),
fill = color_map[[scale_level]], color = "white") +
geom_histogram(data = ~ filter(., variable != "traffic"),
bins = 30, fill = color_map[[scale_level]], color = "white") +
facet_wrap(~ variable, scales = "free", ncol = 5) +
labs(title = paste("Histograms -", scale_level),
x = NULL, y = "frequency") +
theme_minimal()
}
plot_1km <- plot_histograms("1km")
plot_3km <- plot_histograms("3km")
plot_5km <- plot_histograms("5km")
plot_1km
plot_3km
plot_5km
Since the traffic and pavement variables are the same across all subsets, the histograms should look identical for each scale. Let’s generate the histograms to confirm that this information is indeed consistent across the three datasets.
ggplot(combined_raw, aes(x = traffic, fill = interaction(traffic, scale))) +
geom_bar() +
facet_wrap(~ scale) +
theme_minimal() +
labs(
title = "Traffic category frequency by scale",
x = "Traffic category",
y = "Frequency"
) +
scale_fill_manual(values = c(
"low.1km" = "#FFF8B0",
"medium.1km" = "#FFD700",
"high.1km" = "#E6BE00",
"low.3km" = "#AED6F1",
"medium.3km" = "#1f77b4",
"high.3km" = "#154360",
"low.5km" = "#E8DAEF",
"medium.5km" = "#9B59B6",
"high.5km" = "#6C3483"
)) +
theme(legend.position="none")
ggplot(combined_raw, aes(x = pavement, fill = interaction(pavement, scale))) +
geom_bar() +
facet_wrap(~ scale) +
theme_minimal() +
labs(
title = "Pavement category frequency by scale",
x = "Pavement category",
y = "Frequency"
) +
scale_fill_manual(values = c(
"0.1km" = "#FFD700",
"1.1km" = "#E6BE00",
"0.3km" = "#1f77b4",
"1.3km" = "#154360",
"0.5km" = "#9B59B6",
"1.5km" = "#6C3483"
)) +
theme(legend.position = "none")
5.3 Transformations
Due to the heterogeneity in the scales at which the exploratory variables were measured, and aiming to reduce the effect of outliers, the values will be transformed by applying the logarithmic transformation for variables in m² and square root for linear variables. To make variables comparable and ensure they contribute equally to the analysis, all transformed variables were then standardized using z-score normalization (mean = 0, standard deviation = 1). The datasets will be named 1km_sum, 3km_sum, and 5km_sum.
5.3.1 Log
vars_to_log <- c("cult_land", "water", "shrubby", "pasture", "urban_area", "forest",
"savanna")
datasets_sum_t <- lapply(sum_datasets, function(data) {
data %>%
mutate(across(all_of(vars_to_log), ~ log(.x + 1)))
})
d1km_sum <- datasets_sum_t[[1]]
d3km_sum <- datasets_sum_t[[2]]
d5km_sum <- datasets_sum_t[[3]]5.3.2 Square root
vars_to_sqrt <- c("road_length")
datasets_sum_t <- map(datasets_sum_t, ~ .x %>%
mutate(across(all_of(vars_to_sqrt), ~ sqrt(.x))))
d1km_sum <- datasets_sum_t[[1]]
d3km_sum <- datasets_sum_t[[2]]
d5km_sum <- datasets_sum_t[[3]]5.3.3 z-score standardization
vars_to_standardize <- c("cult_land", "water", "shrubby", "pasture",
"urban_area", "forest", "savanna", "road_length")
datasets_sum_t <- lapply(datasets_sum_t, function(data) {
data %>%
mutate(across(all_of(vars_to_standardize), ~ decostand(.x, method = "standardize")))
})
d1km_sum <- datasets_sum_t[[1]]
d3km_sum <- datasets_sum_t[[2]]
d5km_sum <- datasets_sum_t[[3]]5.3.4 Is there collinearity among the covariates?
Collinearity was assessed to avoid redundancy among covariates, which can affect model stability and interpretation.
check_collinearity <- function(data, scale_label) {
landscape_vars <- data %>%
dplyr::select(cult_land, water, shrubby, pasture, urban_area, forest, savanna, pavement)
cor_matrix <- cor(landscape_vars, use = "complete.obs")
# Gradientes simétricos com branco no centro
col_1km <- colorRampPalette(c("#E6BE00", "white", "#E6BE00"))(200)
col_3km <- colorRampPalette(c("#154360", "white", "#154360"))(200)
col_5km <- colorRampPalette(c("#6C3483", "white", "#6C3483"))(200)
col_to_use <- switch(scale_label,
"1 km" = col_1km,
"3 km" = col_3km,
"5 km" = col_5km)
corrplot(cor_matrix,
method = "color",
type = "upper",
tl.col = "black",
tl.srt = 45,
addCoef.col = "black",
col = col_to_use,
cl.lim = c(-1, 1),
cl.pos = "n",
title = paste("Collinearity -", scale_label),
mar = c(0, 0, 2, 0))
}
check_collinearity(d1km_raw_sum, "1 km")
check_collinearity(d3km_raw_sum, "3 km")
check_collinearity(d5km_raw_sum, "5 km")
assocstats(table(d1km_sum$traffic, d1km_sum$pavement))$cramer[1] 0.2383378d1km_sum$traffic <- as.factor(d1km_sum$traffic)
model_1 <- lm(total_roadkill_count ~ cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d1km_sum)
vif(model_1) GVIF Df GVIF^(1/(2*Df))
cult_land 8.070973 1 2.840946
water 1.086858 1 1.042525
shrubby 1.322136 1 1.149842
pasture 10.317503 1 3.212087
urban_area 2.938581 1 1.714229
forest 1.420730 1 1.191944
savanna 13.801790 1 3.715076
pavement 1.125928 1 1.061097
traffic 1.684680 2 1.139277d3km_sum$traffic <- as.factor(d3km_sum$traffic)
model_2 <- lm(total_roadkill_count ~ cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d3km_sum)
vif(model_2) GVIF Df GVIF^(1/(2*Df))
cult_land 5.543270 1 2.354415
water 1.061154 1 1.030123
shrubby 1.218105 1 1.103678
pasture 4.359598 1 2.087965
urban_area 1.598799 1 1.264436
forest 1.763446 1 1.327948
savanna 6.726518 1 2.593553
pavement 1.133821 1 1.064810
traffic 1.964433 2 1.183884d5km_sum$traffic <- as.factor(d5km_sum$traffic)
model_3 <- lm(total_roadkill_count ~ cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d5km_sum)
vif(model_3) GVIF Df GVIF^(1/(2*Df))
cult_land 4.189125 1 2.046735
water 1.117289 1 1.057019
shrubby 1.247112 1 1.116742
pasture 4.107781 1 2.026766
urban_area 1.375354 1 1.172755
forest 2.016615 1 1.420076
savanna 5.734398 1 2.394660
pavement 1.125946 1 1.061106
traffic 1.889809 2 1.172477Pasture and savanna are highly correlated, so we don’t need to include both variables in the model. Therefore, we will redo the VIF analysis without pasture
model_4 <- lm(total_roadkill_count ~ cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d1km_sum)
vif(model_4) GVIF Df GVIF^(1/(2*Df))
cult_land 1.671562 1 1.292889
water 1.086082 1 1.042153
shrubby 1.070989 1 1.034886
urban_area 1.219316 1 1.104226
forest 1.142891 1 1.069061
savanna 1.654020 1 1.286087
pavement 1.122510 1 1.059486
traffic 1.434274 2 1.094355model_5 <- lm(total_roadkill_count ~ cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d3km_sum)
vif(model_5) GVIF Df GVIF^(1/(2*Df))
cult_land 2.162257 1 1.470461
water 1.060612 1 1.029860
shrubby 1.106885 1 1.052086
urban_area 1.322452 1 1.149979
forest 1.564803 1 1.250921
savanna 2.137294 1 1.461949
pavement 1.126024 1 1.061143
traffic 1.560377 2 1.117654model_6 <- lm(total_roadkill_count ~ cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d5km_sum)
vif(model_6) GVIF Df GVIF^(1/(2*Df))
cult_land 2.006813 1 1.416620
water 1.117287 1 1.057018
shrubby 1.141983 1 1.068636
urban_area 1.266050 1 1.125189
forest 1.849897 1 1.360109
savanna 2.115779 1 1.454572
pavement 1.109058 1 1.053118
traffic 1.590907 2 1.123081Apparently, just removing pasture was enough to control the inflation.
6 Statistical Analyses - GLM
It is worth noting that the selected analyses presented here underwent a prior selection process to determine the distribution family that best fit our data. The models used for analyses are referred to as total_Xkm_mx. For example, the global model using the total number of roadkill records at the 1 km scale is named total_1km_m0.
6.1 Total number of roadkill
We want to assess how pavement affects the number of recorded roadkill incidents, regardless of the season or taxonomic group.
VARIABLES
Response variable: Roadkill number
Predictor variable: Pavement as binary
Exploratory variables: Road length; Cultivated land; Water; Shrubby; Pasture*; Urban area; Forest; Savanna and Traffic.
*It were excluded of model due to correlation with Savanna
Since we performed many transformations after initially specifying the nature of the data, it’s better to double-check that everything is still correctly defined
format_data_GLM <- function(df) {
df$polygon <- as.factor(df$polygon)
df$pavement <- factor(df$pavement, levels = c("0", "1"))
df$traffic <- factor(df$traffic, levels = c("low", "medium", "high"), ordered = FALSE)
df$water <- as.numeric(df$water)
df$road_length <- as.numeric(df$road_length)
df$shrubby <- as.numeric(df$shrubby)
df$urban_area <- as.numeric(df$urban_area)
df$cult_land <- as.numeric(df$cult_land)
df$savanna <- as.numeric(df$savanna)
df$forest <- as.numeric(df$forest)
return(df)
}
d1km_sum <- format_data_GLM(d1km_sum)
d3km_sum <- format_data_GLM(d3km_sum)
d5km_sum <- format_data_GLM(d5km_sum)6.1.1 GLM - Negative binomial
We started the analysis by considering traffic as a random effect in a GLM. At first, we tried a Poisson model, but it showed signs of overdispersion. Still keeping traffic as a random effect, we then tried a negative binomial model, which handles overdispersion better. However, it was not possible to fit the model properly. Therefore, we present here a negative binomial GLM, treating traffic as a fixed effect.
6.1.1.1 1km
6.1.1.1.1 Global model
total_1km_m0 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + cult_land + forest+ urban_area + water + shrubby + savanna,
data = d1km_sum)6.1.1.1.2 Overdispersion
(chat <- deviance(total_1km_m0) / df.residual(total_1km_m0))[1] 1.1009576.1.1.1.3 Model selection using drop1
We used the drop1 function in R to compare nested models and identify the most influential predictors. This function tests the effect of removing each term from the full model, one at a time, using likelihood ratio tests. Variables with higher p-values contribute less to the model, so we progressively removed the least significant terms, keeping those with the lowest p-values to improve model parsimony. The same procedure will be applied to the other spacial scales.
Global model
drop1(total_1km_m0, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + cult_land +
forest + urban_area + water + shrubby + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.24 2052.9
traffic 2 568.11 2221.8 172.863 < 2.2e-16 ***
pavement 1 500.46 2156.2 105.219 < 2.2e-16 ***
road_length 1 396.30 2052.0 1.060 0.30313
cult_land 1 398.54 2054.2 3.298 0.06938 .
forest 1 395.26 2050.9 0.018 0.89247
urban_area 1 432.89 2088.6 37.647 8.48e-10 ***
water 1 396.27 2052.0 1.025 0.31142
shrubby 1 397.23 2052.9 1.986 0.15878
savanna 1 411.83 2067.5 16.590 4.64e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without forest
total_1km_m1 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area + water +
cult_land + shrubby + savanna + road_length,data = d1km_sum)
drop1(total_1km_m1, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + water +
cult_land + shrubby + savanna + road_length
Df Deviance AIC LRT Pr(>Chi)
<none> 395.26 2050.9
traffic 2 568.20 2219.9 172.946 < 2.2e-16 ***
pavement 1 500.46 2154.2 105.206 < 2.2e-16 ***
urban_area 1 433.18 2086.9 37.925 7.351e-10 ***
water 1 396.29 2050.0 1.030 0.31010
cult_land 1 398.55 2052.2 3.299 0.06932 .
shrubby 1 397.22 2050.9 1.968 0.16071
savanna 1 412.21 2065.9 16.957 3.824e-05 ***
road_length 1 396.31 2050.0 1.055 0.30426
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
total_1km_m2 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area +
cult_land + shrubby + savanna + road_length,data = d1km_sum)
drop1(total_1km_m2, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + cult_land +
shrubby + savanna + road_length
Df Deviance AIC LRT Pr(>Chi)
<none> 395.04 2050.0
traffic 2 591.37 2242.3 196.327 < 2.2e-16 ***
pavement 1 502.17 2155.1 107.122 < 2.2e-16 ***
urban_area 1 433.51 2086.4 38.465 5.576e-10 ***
cult_land 1 398.72 2051.7 3.672 0.05532 .
shrubby 1 396.85 2049.8 1.810 0.17855
savanna 1 412.77 2065.7 17.730 2.546e-05 ***
road_length 1 396.05 2049.0 1.009 0.31525
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
total_1km_m3 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area +
cult_land + shrubby + savanna, data = d1km_sum)
drop1(total_1km_m3, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + cult_land +
shrubby + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.71 2049.0
traffic 2 594.68 2243.9 198.967 < 2.2e-16 ***
pavement 1 501.84 2153.1 106.128 < 2.2e-16 ***
urban_area 1 433.74 2085.0 38.030 6.966e-10 ***
cult_land 1 399.55 2050.8 3.836 0.05015 .
shrubby 1 397.20 2048.5 1.490 0.22229
savanna 1 413.97 2065.2 18.261 1.926e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
total_1km_m4 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area +
cult_land + savanna, data = d1km_sum)
drop1(total_1km_m4, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + cult_land +
savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 394.16 2048.5
traffic 2 599.21 2249.5 205.049 < 2.2e-16 ***
pavement 1 507.56 2159.9 113.394 < 2.2e-16 ***
urban_area 1 432.05 2084.3 37.890 7.483e-10 ***
cult_land 1 398.45 2050.8 4.289 0.03835 *
savanna 1 414.32 2066.6 20.154 7.147e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without cultivated land
total_1km_m5 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area
+ savanna, data = d1km_sum)
drop1(total_1km_m5, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 393.57 2050.7
traffic 2 605.64 2258.8 212.077 < 2.2e-16 ***
pavement 1 502.10 2157.2 108.538 < 2.2e-16 ***
urban_area 1 427.41 2082.6 33.848 5.960e-09 ***
savanna 1 409.45 2064.6 15.885 6.732e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
total_1km_m6 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area,
data = d1km_sum)
drop1(total_1km_m6, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area
Df Deviance AIC LRT Pr(>Chi)
<none> 393.36 2064.2
traffic 2 600.82 2267.7 207.462 < 2.2e-16 ***
pavement 1 500.09 2168.9 106.735 < 2.2e-16 ***
urban_area 1 418.01 2086.9 24.651 6.869e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
total_1km_m7 <- glm.nb (total_roadkill_count ~ traffic + pavement,
data = d1km_sum)
drop1(total_1km_m7, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 390.44 2085.8
traffic 2 566.42 2257.7 175.980 < 2.2e-16 ***
pavement 1 485.27 2178.6 94.827 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
total_1km_m8 <- glm.nb (total_roadkill_count ~ pavement, data = d1km_sum)
drop1(total_1km_m8, test="Chi")Single term deletions
Model:
total_roadkill_count ~ pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 390.65 2221.2
pavement 1 488.03 2316.6 97.378 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.1.1.2 3km
total_3km_m0 <- glm.nb (total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length + cult_land + shrubby,
data = d3km_sum)6.1.1.2.1 Overdispersion
(chat <- deviance(total_3km_m0) / df.residual(total_3km_m0))[1] 1.1066216.1.1.2.2 Model selection using drop1
drop1(total_3km_m0, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length + cult_land + shrubby
Df Deviance AIC LRT Pr(>Chi)
<none> 393.96 2048.0
traffic 2 552.96 2203.0 159.000 < 2.2e-16 ***
pavement 1 500.47 2152.5 106.518 < 2.2e-16 ***
savanna 1 397.54 2049.6 3.579 0.05851 .
urban_area 1 424.30 2076.4 30.346 3.614e-08 ***
water 1 395.57 2047.6 1.612 0.20417
forest 1 396.94 2049.0 2.984 0.08411 .
road_length 1 397.76 2049.8 3.801 0.05122 .
cult_land 1 394.26 2046.3 0.304 0.58163
shrubby 1 394.15 2046.2 0.197 0.65720
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
total_3km_m1 <- glm.nb (total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length + cult_land,
data = d3km_sum)
drop1(total_3km_m1, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length + cult_land
Df Deviance AIC LRT Pr(>Chi)
<none> 393.46 2046.2
traffic 2 559.61 2208.4 166.152 < 2.2e-16 ***
pavement 1 501.23 2152.0 107.768 < 2.2e-16 ***
savanna 1 397.20 2048.0 3.741 0.05311 .
urban_area 1 424.08 2074.8 30.620 3.138e-08 ***
water 1 394.96 2045.7 1.497 0.22114
forest 1 396.80 2047.5 3.344 0.06745 .
road_length 1 397.17 2047.9 3.707 0.05419 .
cult_land 1 393.79 2044.5 0.329 0.56601
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without cultivated land
total_3km_m2 <- glm.nb (total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length,
data = d3km_sum)
drop1(total_3km_m2, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + savanna + urban_area +
water + forest + road_length
Df Deviance AIC LRT Pr(>Chi)
<none> 393.26 2044.5
traffic 2 610.59 2257.9 217.327 < 2.2e-16 ***
pavement 1 502.09 2151.4 108.833 < 2.2e-16 ***
savanna 1 397.90 2047.2 4.636 0.03130 *
urban_area 1 424.19 2073.5 30.931 2.674e-08 ***
water 1 394.87 2044.1 1.603 0.20548
forest 1 398.60 2047.9 5.338 0.02086 *
road_length 1 396.97 2046.2 3.711 0.05405 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
total_3km_m3 <- glm.nb (total_roadkill_count ~ traffic + pavement + savanna + urban_area +
forest + road_length,
data = d3km_sum)
drop1(total_3km_m3, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + savanna + urban_area +
forest + road_length
Df Deviance AIC LRT Pr(>Chi)
<none> 393.38 2044.1
traffic 2 612.85 2259.6 219.470 < 2.2e-16 ***
pavement 1 503.80 2152.6 110.422 < 2.2e-16 ***
savanna 1 397.85 2046.6 4.465 0.03459 *
urban_area 1 423.07 2071.8 29.687 5.077e-08 ***
forest 1 398.36 2047.1 4.982 0.02561 *
road_length 1 397.10 2045.8 3.715 0.05392 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
total_3km_m4 <- glm.nb (total_roadkill_count ~ traffic + pavement + savanna + urban_area +
forest,
data = d3km_sum)
drop1(total_3km_m4, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + savanna + urban_area +
forest
Df Deviance AIC LRT Pr(>Chi)
<none> 393.77 2045.8
traffic 2 615.40 2263.5 221.634 < 2.2e-16 ***
pavement 1 500.68 2150.7 106.906 < 2.2e-16 ***
savanna 1 398.39 2048.5 4.622 0.03157 *
urban_area 1 422.37 2072.4 28.603 8.884e-08 ***
forest 1 399.47 2049.5 5.705 0.01691 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
total_3km_m5 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area +
forest,
data = d3km_sum)
drop1(total_3km_m5, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area + forest
Df Deviance AIC LRT Pr(>Chi)
<none> 393.24 2048.4
traffic 2 613.45 2264.6 220.211 < 2.2e-16 ***
pavement 1 497.34 2150.5 104.093 < 2.2e-16 ***
urban_area 1 418.17 2071.3 24.929 5.947e-07 ***
forest 1 402.49 2055.7 9.246 0.002361 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without forest
total_3km_m6 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area,
data = d3km_sum)
drop1(total_3km_m6, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area
Df Deviance AIC LRT Pr(>Chi)
<none> 393.50 2055.5
traffic 2 602.17 2260.2 208.673 < 2.2e-16 ***
pavement 1 495.58 2155.6 102.081 < 2.2e-16 ***
urban_area 1 419.34 2079.4 25.834 3.721e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
total_3km_m7 <- glm.nb (total_roadkill_count ~ traffic + pavement,
data = d3km_sum)
drop1(total_3km_m7, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 390.67 2078.2
traffic 2 560.08 2243.6 169.405 < 2.2e-16 ***
pavement 1 481.23 2166.7 90.553 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
total_3km_m8 <- glm.nb (total_roadkill_count ~ pavement,
data = d3km_sum)
drop1(total_3km_m8, test="Chi")Single term deletions
Model:
total_roadkill_count ~ pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 390.55 2209.2
pavement 1 484.43 2301.1 93.877 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.1.1.3 5km
Global model
total_5km_m0 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + forest
+ cult_land + urban_area + water + shrubby + savanna,
data = d5km_sum)6.1.1.3.1 Overdispersion
(chat <- deviance(total_5km_m0) / df.residual(total_5km_m0))[1] 1.1204956.1.1.3.2 Model selection using drop1
drop1(total_5km_m0, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + forest +
cult_land + urban_area + water + shrubby + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.53 2043.2
traffic 2 540.01 2183.7 144.474 < 2.2e-16 ***
pavement 1 500.56 2146.2 105.028 < 2.2e-16 ***
road_length 1 404.76 2050.4 9.229 0.002383 **
forest 1 400.02 2045.7 4.490 0.034094 *
cult_land 1 395.55 2041.2 0.016 0.899454
urban_area 1 420.34 2066.0 24.809 6.331e-07 ***
water 1 397.41 2043.1 1.873 0.171158
shrubby 1 395.96 2041.6 0.424 0.514870
savanna 1 396.90 2042.6 1.370 0.241849
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without cultivated land
total_5km_m1 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + forest
+ urban_area + water + shrubby + savanna,
data = d5km_sum)
drop1(total_5km_m1, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + forest +
urban_area + water + shrubby + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.60 2041.2
traffic 2 598.38 2240.0 202.786 < 2.2e-16 ***
pavement 1 503.34 2147.0 107.745 < 2.2e-16 ***
road_length 1 404.85 2048.5 9.250 0.002355 **
forest 1 401.13 2044.8 5.537 0.018620 *
urban_area 1 421.54 2065.2 25.941 3.521e-07 ***
water 1 397.45 2041.1 1.857 0.172960
shrubby 1 396.01 2039.7 0.417 0.518437
savanna 1 398.44 2042.1 2.845 0.091660 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
total_5km_m2 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + forest
+ urban_area + water + savanna,
data=d5km_sum)
drop1(total_5km_m2, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + forest +
urban_area + water + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.07 2039.7
traffic 2 610.50 2251.1 215.431 < 2.2e-16 ***
pavement 1 504.19 2146.8 109.120 < 2.2e-16 ***
road_length 1 404.33 2046.9 9.261 0.002341 **
forest 1 401.61 2044.2 6.538 0.010558 *
urban_area 1 420.91 2063.5 25.843 3.704e-07 ***
water 1 396.63 2039.2 1.556 0.212208
savanna 1 397.93 2040.5 2.859 0.090843 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
total_5km_m3 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + forest
+ urban_area + savanna,
data=d5km_sum)
drop1(total_5km_m3, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + forest +
urban_area + savanna
Df Deviance AIC LRT Pr(>Chi)
<none> 395.17 2039.2
traffic 2 609.00 2249.0 213.822 < 2.2e-16 ***
pavement 1 506.80 2148.8 111.627 < 2.2e-16 ***
road_length 1 404.35 2046.4 9.172 0.002457 **
forest 1 401.01 2043.0 5.837 0.015696 *
urban_area 1 419.74 2061.8 24.568 7.173e-07 ***
savanna 1 397.99 2040.0 2.815 0.093376 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
total_5km_m4 <- glm.nb (total_roadkill_count ~ traffic + pavement + road_length + forest
+ urban_area,
data=d5km_sum)
drop1(total_5km_m4, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + road_length + forest +
urban_area
Df Deviance AIC LRT Pr(>Chi)
<none> 394.28 2040.0
traffic 2 602.92 2244.6 208.635 < 2.2e-16 ***
pavement 1 503.53 2147.2 109.250 < 2.2e-16 ***
road_length 1 402.84 2046.5 8.558 0.003440 **
forest 1 404.91 2048.6 10.624 0.001116 **
urban_area 1 416.43 2060.1 22.146 2.526e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
total_5km_m5 <- glm.nb (total_roadkill_count ~ traffic + pavement + forest + urban_area,
data=d5km_sum)
drop1(total_5km_m5, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + forest + urban_area
Df Deviance AIC LRT Pr(>Chi)
<none> 395.38 2046.5
traffic 2 604.58 2251.7 209.200 < 2.2e-16 ***
pavement 1 497.77 2146.9 102.387 < 2.2e-16 ***
forest 1 408.87 2058.0 13.489 0.0002399 ***
urban_area 1 415.85 2064.9 20.464 6.075e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without forest
total_5km_m6 <- glm.nb (total_roadkill_count ~ traffic + pavement + urban_area,
data=d5km_sum)
drop1(total_5km_m6, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement + urban_area
Df Deviance AIC LRT Pr(>Chi)
<none> 395.33 2057.7
traffic 2 590.38 2248.7 195.051 < 2.2e-16 ***
pavement 1 494.47 2154.8 99.138 < 2.2e-16 ***
urban_area 1 415.22 2075.6 19.890 8.204e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
total_5km_m7 <- glm.nb (total_roadkill_count ~ traffic + pavement,
data=d5km_sum)
drop1(total_5km_m7, test="Chi")Single term deletions
Model:
total_roadkill_count ~ traffic + pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 392.50 2074.8
traffic 2 556.26 2234.6 163.762 < 2.2e-16 ***
pavement 1 481.42 2161.7 88.921 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
total_5km_m8 <- glm.nb (total_roadkill_count ~ pavement,
data=d5km_sum)
drop1(total_5km_m8, test="Chi")Single term deletions
Model:
total_roadkill_count ~ pavement
Df Deviance AIC LRT Pr(>Chi)
<none> 390.31 2201.5
pavement 1 479.82 2289.0 89.514 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.1.2 AIC Tables
Let’s build the AICc table to understand which model performs best. We tried to use the AICc() function from the MuMIn package, but it didn’t work with our models. The models were fitted using glm.nb (negative binomial), and AICc() gave an error in this case. Because of that, we calculated AICc manually. This allowed us to compare the models properly.
Function to calculate AICc manually
calc_aicc <- function(model) {
k <- length(coef(model)) # Number of parameters
n <- nobs(model) # Number of observations
aic <- AIC(model) # Standard AIC
aicc <- aic + (2 * k * (k + 1)) / (n - k - 1) # AICc formula
return(aicc)
}Lists of all models
# 1km scale models
models_1km <- list(
total_1km_m0, total_1km_m1, total_1km_m2,
total_1km_m3, total_1km_m4, total_1km_m5,
total_1km_m6, total_1km_m7, total_1km_m8
)
# 3km scale models
models_3km <- list(
total_3km_m0, total_3km_m1, total_3km_m2,
total_3km_m3, total_3km_m4, total_3km_m5,
total_3km_m6, total_3km_m7, total_3km_m8
)
# 5km scale models
models_5km <- list(
total_5km_m0, total_5km_m1, total_5km_m2,
total_5km_m3, total_5km_m4, total_5km_m5,
total_5km_m6, total_5km_m7, total_5km_m8
)
# Model names (for identification in the table)
model_names_1km <- paste0("1km_m", 0:8)
model_names_3km <- paste0("3km_m", 0:8)
model_names_5km <- paste0("5km_m", 0:8)Calculate AICc for all models
rank_models <- function(models, model_names) {
results <- data.frame(
Model = model_names,
AICc = sapply(models, calc_aicc),
stringsAsFactors = FALSE
)
results <- results[order(results$AICc), ] # Sort by AICc (best first)
results$delta_AICc <- results$AICc - min(results$AICc)
results$Akaike_weight <- exp(-0.5 * results$delta_AICc) / sum(exp(-0.5 * results$delta_AICc))
return(results)
}Sort by AICc (best model first)
results_1km <- rank_models(models_1km, model_names_1km)
results_3km <- rank_models(models_3km, model_names_3km)
results_5km <- rank_models(models_5km, model_names_5km)
print(results_1km) # Best 1km models (top = best) Model AICc delta_AICc Akaike_weight
5 1km_m4 2050.772 0.0000000 3.529769e-01
4 1km_m3 2051.386 0.6139912 2.596688e-01
3 1km_m2 2052.479 1.7067342 1.503605e-01
6 1km_m5 2052.948 2.1753702 1.189517e-01
2 1km_m1 2053.564 2.7916350 8.740786e-02
1 1km_m0 2055.670 4.8979795 3.049043e-02
7 1km_m6 2066.384 15.6118734 1.437709e-04
8 1km_m7 2087.860 37.0873492 3.121266e-09
9 1km_m8 2223.255 172.4825737 1.240599e-38print(results_3km) # Best 3km models Model AICc delta_AICc Akaike_weight
4 3km_m3 2046.540 0.0000000 3.574868e-01
3 3km_m2 2047.042 0.5023324 2.780865e-01
5 3km_m4 2048.148 1.6081276 1.599777e-01
2 3km_m1 2048.827 2.2871330 1.139240e-01
6 3km_m5 2050.651 4.1110567 4.576731e-02
1 3km_m0 2050.756 4.2166213 4.341425e-02
7 3km_m6 2057.707 11.1678184 1.343380e-03
8 3km_m7 2080.268 33.7284188 1.695223e-08
9 3km_m8 2211.279 164.7390686 6.034284e-37print(results_5km) # Best 5km models Model AICc delta_AICc Akaike_weight
4 5km_m3 2041.608 0.0000000 3.413021e-01
3 5km_m2 2042.158 0.5497347 2.592779e-01
5 5km_m4 2042.311 0.7032999 2.401150e-01
2 5km_m1 2043.857 2.2488427 1.108687e-01
1 5km_m0 2045.968 4.3596518 3.858803e-02
6 5km_m5 2048.707 7.0985222 9.811026e-03
7 5km_m6 2059.857 18.2486354 3.719616e-05
8 5km_m7 2076.925 35.3165524 7.315551e-09
9 5km_m8 2203.496 161.8882747 2.396331e-366.1.3 Nagelkerke’s R²
nagelkerke_r2_nb <- function(model) {
null_model <- update(model, ~1)
LL_model <- as.numeric(logLik(model))
LL_null <- as.numeric(logLik(null_model))
# Calculate Cox-Snell R²
n <- nobs(model)
cox_snell <- 1 - exp((2/n) * (LL_null - LL_model))
# Calculate Nagelkerke's adjustment
max_r2 <- 1 - exp((2/n) * LL_null)
nagelkerke <- cox_snell / max_r2
return(nagelkerke)
}
# ----------------------------------------------
# Apply to all Negative Binomial models
# ----------------------------------------------
# Assuming you have these model lists from previous steps:
# models_1km, models_3km, models_5km
# Calculate Nagelkerke R² for each scale
add_r2_to_results <- function(results_df, models) {
results_df$Nagelkerke_R2 <- sapply(models, nagelkerke_r2_nb)
return(results_df)
}
# Update previous results (from AICc analysis)
results_1km <- add_r2_to_results(results_1km, models_1km)
results_3km <- add_r2_to_results(results_3km, models_3km)
results_5km <- add_r2_to_results(results_5km, models_5km)
# ----------------------------------------------
# Final sorted tables (per scale)
# ----------------------------------------------
# Sort by AICc (best model first) with R² column
results_1km <- results_1km[order(results_1km$AICc), ]
results_3km <- results_3km[order(results_3km$AICc), ]
results_5km <- results_5km[order(results_5km$AICc), ]
# Print results with all metrics
print(results_1km) Model AICc delta_AICc Akaike_weight Nagelkerke_R2
5 1km_m4 2050.772 0.0000000 3.529769e-01 0.5665646
4 1km_m3 2051.386 0.6139912 2.596688e-01 0.5665431
3 1km_m2 2052.479 1.7067342 1.503605e-01 0.5653322
6 1km_m5 2052.948 2.1753702 1.189517e-01 0.5641411
2 1km_m1 2053.564 2.7916350 8.740786e-02 0.5623923
1 1km_m0 2055.670 4.8979795 3.049043e-02 0.5573120
7 1km_m6 2066.384 15.6118734 1.437709e-04 0.5382922
8 1km_m7 2087.860 37.0873492 3.121266e-09 0.5078573
9 1km_m8 2223.255 172.4825737 1.240599e-38 0.2817144print(results_3km) Model AICc delta_AICc Akaike_weight Nagelkerke_R2
4 3km_m3 2046.540 0.0000000 3.574868e-01 0.5604918
3 3km_m2 2047.042 0.5023324 2.780865e-01 0.5602558
5 3km_m4 2048.148 1.6081276 1.599777e-01 0.5598599
2 3km_m1 2048.827 2.2871330 1.139240e-01 0.5579296
6 3km_m5 2050.651 4.1110567 4.576731e-02 0.5534342
1 3km_m0 2050.756 4.2166213 4.341425e-02 0.5478022
7 3km_m6 2057.707 11.1678184 1.343380e-03 0.5363749
8 3km_m7 2080.268 33.7284188 1.695223e-08 0.5040883
9 3km_m8 2211.279 164.7390686 6.034284e-37 0.2826369print(results_5km) Model AICc delta_AICc Akaike_weight Nagelkerke_R2
4 5km_m3 2041.608 0.0000000 3.413021e-01 0.5559657
3 5km_m2 2042.158 0.5497347 2.592779e-01 0.5559461
5 5km_m4 2042.311 0.7032999 2.401150e-01 0.5554368
2 5km_m1 2043.857 2.2488427 1.108687e-01 0.5535285
1 5km_m0 2045.968 4.3596518 3.858803e-02 0.5500746
6 5km_m5 2048.707 7.0985222 9.811026e-03 0.5394359
7 5km_m6 2059.857 18.2486354 3.719616e-05 0.5223404
8 5km_m7 2076.925 35.3165524 7.315551e-09 0.4964806
9 5km_m8 2203.496 161.8882747 2.396331e-36 0.27831776.1.4 % Deviance
deviance_explained <- function(model) {
null_dev <- model$null.deviance # Null deviance (intercept-only model)
resid_dev <- model$deviance # Residual deviance (fitted model)
percent_dev <- (1 - (resid_dev / null_dev)) * 100
return(percent_dev)
}
results_1km$Pct_Deviance_Explained <- sapply(models_1km, deviance_explained)
results_3km$Pct_Deviance_Explained <- sapply(models_3km, deviance_explained)
results_5km$Pct_Deviance_Explained <- sapply(models_5km, deviance_explained)
print(results_1km[order(results_1km$AICc), ]) # Sorted by AICc (best first) Model AICc delta_AICc Akaike_weight Nagelkerke_R2
5 1km_m4 2050.772 0.0000000 3.529769e-01 0.5665646
4 1km_m3 2051.386 0.6139912 2.596688e-01 0.5665431
3 1km_m2 2052.479 1.7067342 1.503605e-01 0.5653322
6 1km_m5 2052.948 2.1753702 1.189517e-01 0.5641411
2 1km_m1 2053.564 2.7916350 8.740786e-02 0.5623923
1 1km_m0 2055.670 4.8979795 3.049043e-02 0.5573120
7 1km_m6 2066.384 15.6118734 1.437709e-04 0.5382922
8 1km_m7 2087.860 37.0873492 3.121266e-09 0.5078573
9 1km_m8 2223.255 172.4825737 1.240599e-38 0.2817144
Pct_Deviance_Explained
5 51.46378
4 51.46144
3 51.31750
6 51.18832
2 50.96341
1 50.36558
7 48.16506
8 44.64964
9 19.953246.1.5 General table
Generate general summary tables for each spatial scale (1km, 3km, 5km)
Extract Model Formulas (Clean Format)
get_formula <- function(model) {
paste(deparse(formula(model)), collapse = " ") |>
stringr::str_squish()
}
compile_results <- function(model_list, scale_name) {
aiccs <- sapply(model_list, calc_aicc)
data.frame(
Scale = rep(scale_name, length(model_list)),
Model = names(model_list),
Formula = sapply(model_list, get_formula),
AICc = aiccs,
Delta_AICc = aiccs - min(aiccs),
Nagelkerke_R2 = sapply(model_list, nagelkerke_r2_nb),
Pct_Deviance = sapply(model_list, deviance_explained),
stringsAsFactors = FALSE
) |>
dplyr::arrange(AICc)
}
names(models_1km) <- paste0("1km_m", 0:8)
names(models_3km) <- paste0("3km_m", 0:8)
names(models_5km) <- paste0("5km_m", 0:8)
# Generate results
final_table <- rbind(
compile_results(models_1km, "1km"),
compile_results(models_3km, "3km"),
compile_results(models_5km, "5km")
)
final_table[,4:7] <- round(final_table[,4:7], 3)
# Add significance flags
final_table$Best <- ifelse(final_table$Delta_AICc == 0, "★", "")
# Save as CSV
write_xlsx(final_table, "final_model_comparison.xlsx")Based on the AIC values, we identified the best-performing models for each spatial scale: model 4 for the 1 km scale, model 3 for the 3 km scale, and again model 3 for the 5 km scale. We will now proceed with model validation.
6.1.6 Basic diagnostics of the model residuals
models_GLM <- list(total_1km_m4, total_3km_m4, total_5km_m4)
model_names <- c("1km", "3km", "5km")
# Plot diagnostics
par(mfrow = c(2, 2))
for (i in 1:length(models_GLM)) {
plot(models_GLM[[i]], main = paste("Model", model_names[i]))
}
par(mfrow = c(1, 3))
chat_values <- sapply(models_GLM, function(m) deviance(m) / df.residual(m))
names(chat_values) <- model_names
chat_values 1km 3km 5km
1.085849 1.093804 1.104432 6.1.7 Advanced diagnostics of the model residuals
par(mfrow = c(1, 3))
simulationOutput <- simulateResiduals(fittedModel = total_1km_m4, plot = TRUE)
simulationOutput <- simulateResiduals(fittedModel = total_3km_m3, plot = TRUE)
simulationOutput <- simulateResiduals(fittedModel = total_5km_m3, plot = TRUE)
6.1.8 Summary
6.1.8.1 1km - Model M4
summary(total_1km_m4)
Call:
glm.nb(formula = total_roadkill_count ~ traffic + pavement +
urban_area + cult_land + savanna, data = d1km_sum, init.theta = 3.355735941,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.29765 0.14558 2.045 0.0409 *
trafficmedium 0.09121 0.08443 1.080 0.2800
traffichigh 1.59690 0.12909 12.371 < 2e-16 ***
pavement1 1.53093 0.15134 10.116 < 2e-16 ***
urban_area -0.27958 0.04248 -6.581 4.67e-11 ***
cult_land -0.09480 0.04476 -2.118 0.0342 *
savanna -0.20591 0.04526 -4.550 5.36e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.3557) family taken to be 1)
Null deviance: 803.81 on 369 degrees of freedom
Residual deviance: 394.16 on 363 degrees of freedom
(6 observations deleted due to missingness)
AIC: 2050.5
Number of Fisher Scoring iterations: 1
Theta: 3.356
Std. Err.: 0.379
2 x log-likelihood: -2034.463 6.1.8.2 3km - Model M3
summary(total_3km_m3)
Call:
glm.nb(formula = total_roadkill_count ~ traffic + pavement +
savanna + urban_area + forest + road_length, data = d3km_sum,
init.theta = 3.265092846, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.31187 0.14860 2.099 0.0358 *
trafficmedium 0.09486 0.08518 1.114 0.2654
traffichigh 1.58522 0.13276 11.940 < 2e-16 ***
pavement1 1.52605 0.15196 10.042 < 2e-16 ***
savanna -0.08441 0.03831 -2.203 0.0276 *
urban_area -0.24907 0.04204 -5.925 3.13e-09 ***
forest 0.09015 0.03943 2.286 0.0222 *
road_length 0.06898 0.03711 1.858 0.0631 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.2651) family taken to be 1)
Null deviance: 784.85 on 366 degrees of freedom
Residual deviance: 393.38 on 359 degrees of freedom
(6 observations deleted due to missingness)
AIC: 2046.1
Number of Fisher Scoring iterations: 1
Theta: 3.265
Std. Err.: 0.366
2 x log-likelihood: -2028.137 6.1.8.3 5km - Model M3
summary(total_5km_m3)
Call:
glm.nb(formula = total_roadkill_count ~ traffic + pavement +
road_length + forest + urban_area + savanna, data = d5km_sum,
init.theta = 3.104198301, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.31794 0.14869 2.138 0.03250 *
trafficmedium 0.05950 0.08701 0.684 0.49407
traffichigh 1.61195 0.13894 11.602 < 2e-16 ***
pavement1 1.55440 0.15158 10.255 < 2e-16 ***
road_length 0.11421 0.03831 2.981 0.00287 **
forest 0.10312 0.04218 2.445 0.01449 *
urban_area -0.22738 0.04166 -5.458 4.81e-08 ***
savanna -0.06984 0.04052 -1.724 0.08478 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(3.1042) family taken to be 1)
Null deviance: 777.05 on 363 degrees of freedom
Residual deviance: 395.17 on 356 degrees of freedom
(6 observations deleted due to missingness)
AIC: 2041.2
Number of Fisher Scoring iterations: 1
Theta: 3.104
Std. Err.: 0.347
2 x log-likelihood: -2023.202 summary_1km_sum <- summary(total_1km_m4)
summary_3km_sum <- summary(total_3km_m3)
summary_5km_sum <- summary(total_5km_m3)
extract_est_p <- function(model_summary) {
coef_table <- model_summary$coefficients
data.frame(
Term = rownames(coef_table),
Estimate = coef_table[, "Estimate"],
P_value = coef_table[, "Pr(>|z|)"],
row.names = NULL
)
}
df_1km_sum <- extract_est_p(summary_1km_sum)
df_3km_sum <- extract_est_p(summary_3km_sum)
df_5km_sum <- extract_est_p(summary_5km_sum)
GLM_scale_comparison <- df_1km_sum %>%
rename(Estimate_1km = Estimate, P_1km = P_value) %>%
full_join(df_3km_sum %>% rename(Estimate_3km = Estimate, P_3km = P_value), by = "Term") %>%
full_join(df_5km_sum %>% rename(Estimate_5km = Estimate, P_5km = P_value), by = "Term")
GLM_scale_comparison Term Estimate_1km P_1km Estimate_3km P_3km
1 (Intercept) 0.29765013 4.090078e-02 0.31186995 3.584338e-02
2 trafficmedium 0.09120898 2.800157e-01 0.09486184 2.654363e-01
3 traffichigh 1.59690467 3.759111e-35 1.58522439 7.283471e-33
4 pavement1 1.53092983 4.703616e-24 1.52604883 9.930132e-24
5 urban_area -0.27958225 4.673442e-11 -0.24907400 3.126391e-09
6 cult_land -0.09480095 3.419241e-02 NA NA
7 savanna -0.20591366 5.362968e-06 -0.08441493 2.756818e-02
8 forest NA NA 0.09014630 2.223214e-02
9 road_length NA NA 0.06897530 6.310264e-02
Estimate_5km P_5km
1 0.31793672 3.249653e-02
2 0.05950344 4.940722e-01
3 1.61194880 4.044971e-31
4 1.55439711 1.127837e-24
5 -0.22737748 4.814821e-08
6 NA NA
7 -0.06983886 8.477630e-02
8 0.10312106 1.448584e-02
9 0.11420893 2.872694e-03write_xlsx(GLM_scale_comparison, "GLM_summary_comparison.xlsx")6.1.9 Summary plot
extract_summary <- function(model, label) {
sum <- summary(model)
coefs <- as.data.frame(sum$coefficients)
coefs$term <- rownames(coefs)
coefs$model <- label
colnames(coefs) <- c("estimate", "std_error", "z_value", "p_value", "term", "model")
return(coefs)
}
coef_1km <- extract_summary(total_1km_m4, "1 km")
coef_3km <- extract_summary(total_3km_m3, "3 km")
coef_5km <- extract_summary(total_5km_m3, "5 km")
coef_data_all <- rbind(coef_1km, coef_3km, coef_5km)
coef_data_all <- coef_data_all %>%
mutate(term = case_when(
str_detect(term, "^trafficmedium") ~ "Traffic (Medium)",
str_detect(term, "^traffichigh") ~ "Traffic (High)",
str_detect(term, "^pavement") ~ "Pavement",
str_detect(term, "^urban_area") ~ "Urban Area",
str_detect(term, "^cult_land") ~ "Cultivated Land",
str_detect(term, "^savanna") ~ "Savanna",
str_detect(term, "^forest") ~ "Forest",
str_detect(term, "^road_length") ~ "Road Length",
TRUE ~ term
)) %>%
filter(term != "(Intercept)") coef_data_all$p_value_num <- as.numeric(gsub("< 0.001", "0.00099", coef_data_all$p_value))
# Calculando intervalo de confiança
coef_data_all$lower <- coef_data_all$estimate - 1.96 * coef_data_all$std_error
coef_data_all$upper <- coef_data_all$estimate + 1.96 * coef_data_all$std_error
# Significância
coef_data_all$significant <- ifelse(coef_data_all$p_value_num < 0.05, "Significant", "Non-significant")coef_data_all$term <- factor(coef_data_all$term,
levels = unique(coef_data_all$term[order(coef_data_all$estimate)]))
cores_modelos <- c("1 km" = "#FFD700", "3 km" = "#1f77b4", "5 km" = "#9b59b6")
coef_1km <- coef_data_all %>% filter(model == "1 km", term != "(Intercept)")
coef_3km <- coef_data_all %>% filter(model == "3 km", term != "(Intercept)")
coef_5km <- coef_data_all %>% filter(model == "5 km", term != "(Intercept)")plot_model <- function(data, cor_model, titulo) {
data <- data %>%
mutate(cor_final = ifelse(p_value_num < 0.05, cor_model, "lightgray"))
ggplot(data, aes(x = estimate, y = term)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
geom_point(aes(color = cor_final), size = 3) +
geom_errorbarh(aes(xmin = lower, xmax = upper, color = cor_final), height = 0.2) +
scale_color_identity() +
labs(title = titulo, x = "Estimated Coefficients", y = "") +
theme_minimal() +
theme(
axis.line = element_line(color = "black"),
legend.position = "none",
plot.title = element_text(hjust = 0.5)
)
}
g1 <- plot_model(coef_1km, "#FFD700", "1 km")
g2 <- plot_model(coef_3km, "#1f77b4", "3 km")
g3 <- plot_model(coef_5km, "#9b59b6", "5 km") # 1. Remover intercepto
coef_data_all <- coef_data_all %>%
filter(term != "(Intercept)")
# 3. Ordenar os termos no eixo y
coef_data_all$term <- factor(coef_data_all$term,
levels = unique(coef_data_all$term[order(coef_data_all$estimate)]))
# 4. Gráfico final
Summary_coeff_GLM <- ggplot(coef_data_all, aes(x = estimate, y = term, color = model, shape = significant)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
geom_point(position = position_dodge(width = 0.6), size = 2.5) +
geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2,
position = position_dodge(width = 0.6)) +
scale_y_discrete(expand = expansion(mult = c(0.1, 0.1)))+
scale_color_manual(
values = c("1 km" = "#FFD700", "3 km" = "#1f77b4", "5 km" = "#9b59b6"),
name = NULL
) +
scale_shape_manual(
values = c("Significant" = 16, "Non-significant" = 1),
labels = c("Non-significant", "Significant"),
name = NULL
) +
labs(
x = "Estimated Coefficients",
y = "",
title = NULL
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black"),
plot.title = element_text(hjust = 0.5),
legend.position = "bottom",
legend.text = element_text(size=13),
axis.title = element_text(size=14),
axis.text = element_text(size = 12),
axis.text.y = element_text(size = 12, margin = margin(r = 10))
)
Summary_coeff_GLM
6.1.10 Relative Rate
rate_pav <- d1km_raw_sum |>
group_by(pavement) |>
summarise(
mean = mean(relative_rate, na.rm = TRUE),
se = sd(relative_rate, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
) |>
mutate(
ci_lower = mean - 1.96 * se,
ci_upper = mean + 1.96 * se
)
plot_rate_pav_point <- ggplot(rate_pav, aes(x = factor(pavement), y = mean, color = factor(pavement))) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.15) +
scale_color_manual(
values = c("0" = "darkgoldenrod3", "1" = "gray40"),
labels = c("Unpaved", "Paved")
) +
scale_x_discrete(labels = c("0" = "Unpaved", "1" = "Paved")) +
labs(
x = "Pavement",
y = expression("WVC rate (individuals km⁻¹)")
) +
theme(
panel.background = element_blank(),
axis.ticks = element_line(),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
axis.text = element_text(size = 10),
axis.line = element_line(),
legend.position = "none"
)
rate_traffic <- d1km_raw_sum |>
mutate(
traffic = factor(
traffic,
levels = c("low", "medium", "high")
)
) |>
group_by(traffic) |>
summarise(
mean = mean(relative_rate, na.rm = TRUE),
se = sd(relative_rate, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
) |>
mutate(
ci_lower = mean - 1.96 * se,
ci_upper = mean + 1.96 * se
)
plot_rate_traffic_point <- ggplot(rate_traffic, aes(x = traffic, y = mean, color = traffic)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.15) +
scale_color_manual(
values = c("low" = "#2ECC71", "medium" = "#F39C12", "high" = "#E63946"),
labels = c("Low", "Medium", "High")
) +
scale_x_discrete(labels = c("low" = "Low", "medium" = "Medium", "high" = "High")) +
labs(
x = "Traffic",
y = expression("Relative_rate")
) +
theme(
panel.background = element_blank(),
axis.ticks = element_line(),
axis.title.x = element_text(size = 12),
axis.title.y = element_blank(),
axis.text = element_text(size = 10),
axis.line = element_line(),
legend.position = "none"
)
combined_rate_point_plot <- wrap_plots(list(plot_rate_pav_point, plot_rate_traffic_point), ncol = 2) +
plot_annotation(
title = "1km",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
)
)
combined_rate_point_plot
6.1.10.1 Landscape
- 1km : Urban area, cultivated land, and savanna
landscape_vars_1km <- c("cult_land", "savanna", "urban_area")
labels <- c("Cultivated land (km²)", "Savanna (km²)", "Urban area (km²)")
landscape_plots <- list()
for (i in seq_along(landscape_vars_1km)) {
var <- landscape_vars_1km[i]
var_seq <- seq(
min(d1km_raw_sum[[var]], na.rm = TRUE),
max(d1km_raw_sum[[var]], na.rm = TRUE),
length.out = 100
)
new_data <- d1km_raw_sum[1, , drop = FALSE]
new_data <- new_data[rep(1, 100), ]
new_data[[var]] <- var_seq
for (v in setdiff(landscape_vars_1km, var)) {
new_data[[v]] <- mean(d1km_raw_sum[[v]], na.rm = TRUE)
}
new_data$pavement <- factor("1", levels = c("0", "1"))
new_data$traffic <- factor("medium", levels = c("low", "medium", "high"))
pred <- predict(total_1km_m4, newdata = new_data, type = "response", se.fit = TRUE)
new_data$fit <- pred$fit
new_data$lower <- pmax(pred$fit - 1.96 * pred$se.fit, 0)
new_data$upper <- pred$fit + 1.96 * pred$se.fit
p <- ggplot() +
geom_point(data = d1km_raw_sum, aes(x = !!sym(var), y = total_roadkill_count), alpha = 0.6) +
geom_line(data = new_data, aes(x = !!sym(var), y = fit), color = "#FFD700", size = 1.2) +
geom_ribbon(data = new_data, aes(x = !!sym(var), ymin = lower, ymax = upper), fill = "blue", alpha = 0.2) +
labs(
x = labels[i],
y = "Roadkill Number") +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(),
aspect.ratio = (1.0))
landscape_plots[[var]] <- p
}
wrap_plots(landscape_plots)+
plot_annotation(
title = "1km",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
))
- 3km: Forest, cultivated land, and savanna
landscape_vars_3km <- c("savanna", "urban_area", "forest")
labels <- c("Savanna (km²)", "Urban area (km²)", "Forest (km²)")
landscape_plots <- list()
for (i in seq_along(landscape_vars_3km)) {
var <- landscape_vars_3km[i]
var_seq <- seq(
min(d3km_raw_sum[[var]], na.rm = TRUE),
max(d3km_raw_sum[[var]], na.rm = TRUE),
length.out = 100
)
new_data <- d3km_raw_sum[1, , drop = FALSE]
new_data <- new_data[rep(1, 100), ]
new_data[[var]] <- var_seq
for (v in setdiff(landscape_vars_3km, var)) {
new_data[[v]] <- mean(d3km_raw_sum[[v]], na.rm = TRUE)
}
new_data$pavement <- factor("1", levels = c("0", "1"))
new_data$traffic <- factor("medium", levels = c("low", "medium", "high"))
pred <- predict(total_3km_m3, newdata = new_data, type = "response", se.fit = TRUE)
new_data$fit <- pred$fit
new_data$lower <- pmax(pred$fit - 1.96 * pred$se.fit, 0)
new_data$upper <- pred$fit + 1.96 * pred$se.fit
p <- ggplot() +
geom_point(data = d3km_raw_sum, aes(x = !!sym(var), y = total_roadkill_count), alpha = 0.3) +
geom_line(data = new_data, aes(x = !!sym(var), y = fit), color = "blue", size = 1.2) +
geom_ribbon(data = new_data, aes(x = !!sym(var), ymin = lower, ymax = upper), fill = "blue", alpha = 0.2) +
labs(
x = labels[i],
y = "Roadkill Number",
) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(),
aspect.ratio = (1.0))
landscape_plots[[var]] <- p
}
wrap_plots(landscape_plots)+
plot_annotation(
title = "3km",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
))
landscape_vars_3km <- c("savanna", "urban_area", "forest")
labels <- c("Savanna (km²)", "Urban area (km²)", "Forest (km²")
landscape_plots <- list()
for (i in seq_along(landscape_vars_3km)) {
var <- landscape_vars_3km[i]
var_seq <- seq(
min(d3km_raw_sum[[var]], na.rm = TRUE),
max(d3km_raw_sum[[var]], na.rm = TRUE),
length.out = 100
)
new_data <- d3km_raw_sum[1, , drop = FALSE]
new_data <- new_data[rep(1, 100), ]
new_data[[var]] <- var_seq
for (v in setdiff(landscape_vars_3km, var)) {
new_data[[v]] <- mean(d3km_raw_sum[[v]], na.rm = TRUE)
}
new_data$pavement <- factor("1", levels = c("0", "1"))
new_data$traffic <- factor("medium", levels = c("low", "medium", "high"))
pred <- predict(total_3km_m3, newdata = new_data, type = "response", se.fit = TRUE)
new_data$fit <- pred$fit
new_data$lower <- pmax(pred$fit - 1.96 * pred$se.fit, 0)
new_data$upper <- pred$fit + 1.96 * pred$se.fit
p <- ggplot() +
geom_point(data = d3km_raw_sum, aes(x = !!sym(var), y = total_roadkill_count), alpha = 0.3) +
geom_line(data = new_data, aes(x = !!sym(var), y = fit), color = "#9B59B6", size = 1.2) +
geom_ribbon(data = new_data, aes(x = !!sym(var), ymin = lower, ymax = upper), fill = "blue", alpha = 0.2) +
labs(
x = labels[i],
y = "Roadkill Number",
) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(),
aspect.ratio = (1.0))
landscape_plots[[var]] <- p
}
wrap_plots(landscape_plots)+
plot_annotation(
title = "5km",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
))