library(openxlsx)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(vegan)
library(corrplot)
library(vcd)
library(car)
library(glmmTMB)
library(performance)
library(DHARMa)
library(ggeffects)
library(writexl)Effect of Pavement on the Number of Wildlife Roadkills - Taxa
Data Collection Context
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².
Purple points represent records in unpaved areas, and yellow points represent records in paved areas.
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 roadkills in paved areas.
Note
Analyses for this objective are in another script: “Effect of Pavement on the Number of Wildlife Roadkills - General“
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.
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 taxa 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.
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). |
| roadkill count | Number of roadkill occurrences recorded in the polygon/taxa. |
| road length | Length of the road polygon (in kilometers). |
| relative rate | A measure of roadkill occurrences relative to the road length. |
| cult land | Proportion of cultivated land within the polygon’s surrounding area. |
| water | Proportion of water bodies within the polygon’s surrounding area. |
| shrubby | Proportion of shrubby vegetation within the polygon’s surrounding area. |
| pasture | Proportion of pastureland within the polygon’s surrounding area. |
| urban area | Proportion of urbanized area within the polygon’s surrounding area. |
| forest | Proportion of forested land within the polygon’s surrounding area. |
| savanna | Proportion of savanna within the polygon’s surrounding area. |
| pavement | Proportion of paved road surface./Binary (0 and 1) |
| traffic | A measure of traffic intensity in the road segment. |
Packages
Datasets
d1km_raw <- read.xlsx("data_1km.xlsx")
d3km_raw <- read.xlsx("data_3km.xlsx")
d5km_raw <- read.xlsx("data_5km.xlsx")sum(d1km_raw$roadkill_count)[1] 2959Pavement 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 <- d1km_raw %>% filter(pavement %in% c(0, 1))
d3km_raw <- d3km_raw %>% filter(pavement %in% c(0, 1))
d5km_raw <- d5km_raw %>% filter(pavement %in% c(0, 1))checking the values
all(d1km_raw$pavement %in% c(0, 1)) [1] TRUEall(d3km_raw$pavement %in% c(0, 1)) [1] TRUEall(d5km_raw$pavement %in% c(0, 1))[1] TRUEI’ll create a list containing the three datasets to simplify the required data manipulation
datasets_raw <- list("1km" = d1km_raw, "3km" = d3km_raw, "5km" = d5km_raw)Now that all the data has been prepared, we can combine the 1 km, 3 km, and 5 km datasets into a single dataframe, while keeping a column to identify the scale
combined_raw <- bind_rows(
d1km_raw %>% mutate(scale = "1km"),
d3km_raw %>% mutate(scale = "3km"),
d5km_raw %>% mutate(scale = "5km"))Data Exploration
General exploration - dataset with taxa and season
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.
numeric_vars <- c("roadkill_count", "road_length", "relative_rate",
"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 = 4) +
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")
Trasformations
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, disaggregated by taxon and season, that have been transformed and standardized will be referred to as d1km, d3km, and d5km. The datasets containing the taxa number of roadkills (without separating by taxon or season) will be named 1km_sum, 3km_sum, and 5km_sum.
Log
vars_to_log <- c("cult_land", "water", "shrubby", "pasture",
"urban_area", "forest", "savanna")
datasets_t <- lapply(datasets_raw, function(data) {
data %>%
mutate(across(all_of(vars_to_log), ~ log(.x + 1)))
})
d1km <- datasets_t[[1]]
d3km <- datasets_t[[2]]
d5km <- datasets_t[[3]]Square root
vars_to_sqrt <- c("road_length")
datasets_t <- map(datasets_t, ~ .x %>%
mutate(across(all_of(vars_to_sqrt), ~ sqrt(.x))))
d1km <- datasets_t[[1]]
d3km <- datasets_t[[2]]
d5km <- datasets_t[[3]]z-score standardization
vars_to_standardize <- c("cult_land", "water", "shrubby", "pasture",
"urban_area", "forest", "savanna", "road_length")
datasets_t <- lapply(datasets_t, function(data) {
data %>%
mutate(across(all_of(vars_to_standardize), ~ decostand(.x, method = "standardize")))
})
d1km <- datasets_t[[1]]
d3km <- datasets_t[[2]]
d5km <- datasets_t[[3]]Number of records in each category
To understand how many records exist for each combination of taxon, pavement type, regardless of the specific group
nr1 <- d1km %>%
group_by(pavement, taxa) %>%
summarise(total = sum(roadkill_count), .groups = "drop") %>%
arrange(desc(total))
nr3 <-d3km %>%
group_by(pavement, taxa)%>%
summarise(total=sum(roadkill_count), .groups = "drop") %>%
arrange(desc(total))
nr5 <- d5km %>%
group_by(pavement, taxa)%>%
summarise(total=sum(roadkill_count), .groups= "drop") %>%
arrange (desc(total))
nr <- nr1 %>%
rename(total_1km = total) %>%
full_join(nr3 %>% rename(total_3km = total), by = c("pavement", "taxa")) %>%
full_join(nr5 %>% rename(total_5km = total), by = c("pavement", "taxa"))
nr# A tibble: 8 × 5
pavement taxa total_1km total_3km total_5km
<dbl> <chr> <dbl> <dbl> <dbl>
1 1 bird 789 789 791
2 1 amp 783 783 777
3 1 mam 720 715 714
4 1 rep 560 557 557
5 0 amp 21 22 24
6 0 bird 20 20 19
7 0 rep 18 17 18
8 0 mam 13 13 13Is 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")
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),
title = paste("Collinearity -", scale_label),
mar = c(0, 0, 2, 0),
cl.pos = "n")
}
check_collinearity(d1km, "1 km")
check_collinearity(d3km, "3 km")
check_collinearity(d5km, "5 km")
assocstats(table(d1km$taxa, d1km$traffic))$cramer[1] 0assocstats(table(d1km$taxa, d1km$pavement))$cramer[1] 0assocstats(table(d1km$traffic, d1km$pavement))$cramer[1] 0.2072176d1km$traffic <- as.factor(d1km$traffic)
model_1 <- lm(roadkill_count ~ taxa+ cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d1km)
vif(model_1) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 8.140115 1 2.853089
water 1.081052 1 1.039737
shrubby 1.322578 1 1.150034
pasture 10.369777 1 3.220214
urban_area 2.955876 1 1.719266
forest 1.434016 1 1.197504
savanna 13.802385 1 3.715156
pavement 1.123683 1 1.060039
traffic 1.686760 2 1.139629d3km$traffic <- as.factor(d3km$traffic)
model_2 <- lm(roadkill_count ~ taxa+cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d3km)
vif(model_2) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 5.644237 1 2.375760
water 1.061962 1 1.030515
shrubby 1.210668 1 1.100304
pasture 4.382719 1 2.093495
urban_area 1.605080 1 1.266917
forest 1.801305 1 1.342127
savanna 6.847589 1 2.616790
pavement 1.136955 1 1.066281
traffic 1.992440 2 1.188082d5km$traffic <- as.factor(d5km$traffic)
model_3 <- lm(roadkill_count ~ taxa+cult_land + water + shrubby + pasture + urban_area + forest + savanna + pavement + traffic, data = d5km)
vif(model_3) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 4.215958 1 2.053280
water 1.117014 1 1.056889
shrubby 1.251948 1 1.118905
pasture 4.130414 1 2.032342
urban_area 1.376548 1 1.173264
forest 2.048739 1 1.431342
savanna 5.771882 1 2.402474
pavement 1.129891 1 1.062963
traffic 1.923365 2 1.177648Pasture 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(roadkill_count ~ taxa + cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d1km)
vif(model_4) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 1.674769 1 1.294129
water 1.080373 1 1.039410
shrubby 1.071651 1 1.035206
urban_area 1.215980 1 1.102715
forest 1.145829 1 1.070434
savanna 1.650706 1 1.284798
pavement 1.119555 1 1.058090
traffic 1.438187 2 1.095100model_5 <- lm(roadkill_count ~ taxa + cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d3km)
vif(model_5) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 2.196601 1 1.482094
water 1.061471 1 1.030277
shrubby 1.100142 1 1.048876
urban_area 1.319776 1 1.148815
forest 1.585130 1 1.259019
savanna 2.161837 1 1.470319
pavement 1.131427 1 1.063686
traffic 1.575984 2 1.120438model_6 <- lm(roadkill_count ~ taxa + cult_land + water + shrubby + urban_area + forest + savanna + pavement + traffic, data = d5km)
vif(model_6) GVIF Df GVIF^(1/(2*Df))
taxa 1.000000 3 1.000000
cult_land 2.007944 1 1.417019
water 1.117014 1 1.056889
shrubby 1.144921 1 1.070010
urban_area 1.264783 1 1.124626
forest 1.866391 1 1.366159
savanna 2.119969 1 1.456011
pavement 1.111799 1 1.054419
traffic 1.612147 2 1.126811Apparently, just removing pasture was enough to control the inflation.
Statistical Analyses - GLMM
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. For further details, refer to scripts taxa_and_season and Effect_of_Pavement_Taxon_Analyses_Adjustment.
GLMM - Negative binomial
Ensure that the variables are correctly interpreted as either numeric or factors:
format_data_GLMM <- 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"))
df$taxa <- factor(df$taxa, levels = c("amp", "rep", "bird", "mam"))
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 <- format_data_GLMM(d1km)
d3km <- format_data_GLMM(d3km)
d5km <- format_data_GLMM(d5km)Combining taxon and pavement into a single variable
Since this issue with NA values occurred, Professor Andre Guaraldo suggested combining the information from the ‘taxon’ and ‘pavement’ columns into a single column
taxa_pav <- function(df) {
df %>%
mutate(taxa_pav = factor(paste(taxa, pavement, sep = ""),
levels = c("mam0", "rep0", "bird0", "amp0", "rep1", "mam1", "amp1", "bird1")))
}
d1km <- taxa_pav(d1km)
d3km <- taxa_pav(d3km)
d5km <- taxa_pav(d5km)1km
taxa_1km_m0 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + forest + savanna + water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m0) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
forest + savanna + water + shrubby + road_length + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7359.3 7473.1 -3660.6 7321.3 2941
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2076 0.4556
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.191293 0.300995 -7.280 3.33e-13 ***
taxa_pavrep0 0.308599 0.382912 0.806 0.420
taxa_pavbird0 0.425592 0.375311 1.134 0.257
taxa_pavamp0 0.449274 0.372497 1.206 0.228
taxa_pavrep1 1.605705 0.307168 5.227 1.72e-07 ***
taxa_pavmam1 1.822519 0.306460 5.947 2.73e-09 ***
taxa_pavamp1 1.818086 0.306337 5.935 2.94e-09 ***
taxa_pavbird1 2.027948 0.305892 6.630 3.37e-11 ***
trafficmedium 0.107724 0.086402 1.247 0.212
traffichigh 1.499333 0.130137 11.521 < 2e-16 ***
cult_land -0.072015 0.045859 -1.570 0.116
urban_area -0.315825 0.044839 -7.044 1.87e-12 ***
forest -0.007665 0.037911 -0.202 0.840
savanna -0.186257 0.047122 -3.953 7.73e-05 ***
water 0.033585 0.033888 0.991 0.322
shrubby -0.072128 0.045512 -1.585 0.113
road_length 0.055487 0.038816 1.429 0.153
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Overdispersion
(chat <- deviance(taxa_1km_m0) / df.residual(taxa_1km_m0))[1] 0.88756323km
taxa_3km_m0 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + forest + savanna + water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m0) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
forest + savanna + water + shrubby + road_length + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7334.1 7447.8 -3648.0 7296.1 2917
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2173 0.4661
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.17809 0.30278 -7.194 6.31e-13 ***
taxa_pavrep0 0.25896 0.38809 0.667 0.5046
taxa_pavbird0 0.41982 0.37644 1.115 0.2648
taxa_pavamp0 0.47453 0.37109 1.279 0.2010
taxa_pavrep1 1.59580 0.30848 5.173 2.30e-07 ***
taxa_pavmam1 1.81544 0.30779 5.898 3.67e-09 ***
taxa_pavamp1 1.81145 0.30769 5.887 3.93e-09 ***
taxa_pavbird1 2.02947 0.30716 6.607 3.92e-11 ***
trafficmedium 0.09937 0.08753 1.135 0.2563
traffichigh 1.53917 0.14145 10.882 < 2e-16 ***
cult_land -0.01421 0.05401 -0.263 0.7924
urban_area -0.28666 0.04628 -6.195 5.84e-10 ***
forest 0.07568 0.04521 1.674 0.0941 .
savanna -0.09117 0.05500 -1.658 0.0974 .
water 0.05349 0.03651 1.465 0.1428
shrubby -0.04014 0.04066 -0.987 0.3235
road_length 0.07614 0.03554 2.143 0.0321 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Overdispersion
(chat <- deviance(taxa_3km_m0) / df.residual(taxa_3km_m0))[1] 0.88753195km
taxa_5km_m0 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + forest + savanna + water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m0) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
forest + savanna + water + shrubby + road_length + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7281.2 7394.7 -3621.6 7243.2 2893
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2342 0.4839
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.20663 0.30442 -7.249 4.21e-13 ***
taxa_pavrep0 0.30726 0.38515 0.798 0.425015
taxa_pavbird0 0.36629 0.38105 0.961 0.336420
taxa_pavamp0 0.55647 0.36726 1.515 0.129726
taxa_pavrep1 1.64623 0.30982 5.314 1.07e-07 ***
taxa_pavmam1 1.85217 0.30914 5.991 2.08e-09 ***
taxa_pavamp1 1.83939 0.30903 5.952 2.65e-09 ***
taxa_pavbird1 2.07588 0.30854 6.728 1.72e-11 ***
trafficmedium 0.08654 0.08933 0.969 0.332674
traffichigh 1.55534 0.15118 10.288 < 2e-16 ***
cult_land 0.01807 0.05302 0.341 0.733176
urban_area -0.26524 0.04562 -5.814 6.11e-09 ***
forest 0.09949 0.05034 1.976 0.048131 *
savanna -0.06372 0.05509 -1.157 0.247436
water 0.06339 0.03903 1.624 0.104318
shrubby -0.04523 0.04094 -1.105 0.269261
road_length 0.12073 0.03612 3.342 0.000831 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Overdispersion
(chat <- deviance(taxa_5km_m0) / df.residual(taxa_5km_m0))[1] 0.8837043Model selection
The variables were removed based on the highest P-value indicated in the model summary
1km
- without forest
taxa_1km_m1 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + savanna + water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m1) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
savanna + water + shrubby + road_length + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7357.3 7465.2 -3660.7 7321.3 2942
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2076 0.4556
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.18923 0.30081 -7.278 3.40e-13 ***
taxa_pavrep0 0.30891 0.38294 0.807 0.420
taxa_pavbird0 0.42578 0.37533 1.134 0.257
taxa_pavamp0 0.44917 0.37252 1.206 0.228
taxa_pavrep1 1.60514 0.30718 5.225 1.74e-07 ***
taxa_pavmam1 1.82202 0.30648 5.945 2.76e-09 ***
taxa_pavamp1 1.81777 0.30636 5.933 2.97e-09 ***
taxa_pavbird1 2.02749 0.30591 6.628 3.41e-11 ***
trafficmedium 0.10470 0.08510 1.230 0.219
traffichigh 1.49859 0.13008 11.521 < 2e-16 ***
cult_land -0.07050 0.04524 -1.558 0.119
urban_area -0.31504 0.04466 -7.054 1.74e-12 ***
savanna -0.18442 0.04624 -3.988 6.66e-05 ***
water 0.03368 0.03388 0.994 0.320
shrubby -0.07124 0.04528 -1.573 0.116
road_length 0.05542 0.03879 1.429 0.153
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
taxa_1km_m2 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + savanna + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m2) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
savanna + shrubby + road_length + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7356.3 7458.2 -3661.1 7322.3 2943
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2092 0.4574
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.19415 0.30086 -7.293 3.04e-13 ***
taxa_pavrep0 0.30906 0.38296 0.807 0.420
taxa_pavbird0 0.42606 0.37535 1.135 0.256
taxa_pavamp0 0.44897 0.37254 1.205 0.228
taxa_pavrep1 1.61597 0.30711 5.262 1.43e-07 ***
taxa_pavmam1 1.83264 0.30641 5.981 2.22e-09 ***
taxa_pavamp1 1.82874 0.30628 5.971 2.36e-09 ***
taxa_pavbird1 2.03797 0.30585 6.663 2.68e-11 ***
trafficmedium 0.09300 0.08440 1.102 0.270
traffichigh 1.51766 0.12895 11.769 < 2e-16 ***
cult_land -0.07428 0.04516 -1.645 0.100
urban_area -0.31501 0.04475 -7.039 1.94e-12 ***
savanna -0.18802 0.04619 -4.071 4.69e-05 ***
shrubby -0.06883 0.04524 -1.522 0.128
road_length 0.05404 0.03884 1.392 0.164
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
taxa_1km_m3 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + savanna + shrubby + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m3) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
savanna + shrubby + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7356.2 7452.1 -3662.1 7324.2 2944
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2091 0.4573
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.17295 0.30017 -7.239 4.52e-13 ***
taxa_pavrep0 0.31089 0.38277 0.812 0.4167
taxa_pavbird0 0.42756 0.37516 1.140 0.2544
taxa_pavamp0 0.45724 0.37227 1.228 0.2194
taxa_pavrep1 1.59277 0.30632 5.200 2.00e-07 ***
taxa_pavmam1 1.80893 0.30561 5.919 3.24e-09 ***
taxa_pavamp1 1.80613 0.30553 5.912 3.39e-09 ***
taxa_pavbird1 2.01485 0.30507 6.605 3.99e-11 ***
trafficmedium 0.09137 0.08435 1.083 0.2787
traffichigh 1.52840 0.12872 11.874 < 2e-16 ***
cult_land -0.07600 0.04513 -1.684 0.0922 .
urban_area -0.31140 0.04456 -6.988 2.78e-12 ***
savanna -0.19069 0.04613 -4.133 3.58e-05 ***
shrubby -0.06218 0.04540 -1.370 0.1708
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
taxa_1km_m4 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + cult_land + urban_area + savanna + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m4) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + cult_land + urban_area +
savanna + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7356.1 7446.0 -3663.1 7326.1 2945
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.213 0.4615
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.20109 0.29996 -7.338 2.17e-13 ***
taxa_pavrep0 0.30840 0.38283 0.806 0.4205
taxa_pavbird0 0.42715 0.37520 1.138 0.2549
taxa_pavamp0 0.45289 0.37231 1.216 0.2238
taxa_pavrep1 1.62127 0.30594 5.299 1.16e-07 ***
taxa_pavmam1 1.83704 0.30524 6.018 1.76e-09 ***
taxa_pavamp1 1.83402 0.30516 6.010 1.85e-09 ***
taxa_pavbird1 2.04324 0.30469 6.706 2.00e-11 ***
trafficmedium 0.09471 0.08459 1.120 0.2629
traffichigh 1.54789 0.12861 12.036 < 2e-16 ***
cult_land -0.08186 0.04515 -1.813 0.0698 .
urban_area -0.31275 0.04475 -6.988 2.78e-12 ***
savanna -0.20113 0.04576 -4.395 1.11e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without cultivated land
taxa_1km_m5 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + savanna + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m5) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + traffic + urban_area + savanna +
(1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7357.4 7441.3 -3664.7 7329.4 2946
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2186 0.4675
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.86
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.19501 0.30027 -7.310 2.67e-13 ***
taxa_pavrep0 0.30554 0.38285 0.798 0.4248
taxa_pavbird0 0.42498 0.37518 1.133 0.2573
taxa_pavamp0 0.45554 0.37228 1.224 0.2211
taxa_pavrep1 1.58369 0.30534 5.187 2.14e-07 ***
taxa_pavmam1 1.79813 0.30460 5.903 3.56e-09 ***
taxa_pavamp1 1.79781 0.30463 5.902 3.60e-09 ***
taxa_pavbird1 2.00542 0.30409 6.595 4.26e-11 ***
trafficmedium 0.15026 0.07933 1.894 0.0582 .
traffichigh 1.51556 0.12831 11.812 < 2e-16 ***
urban_area -0.29724 0.04412 -6.737 1.62e-11 ***
savanna -0.15478 0.03817 -4.055 5.02e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
taxa_1km_m6 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m6) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7371.5 7449.4 -3672.8 7345.5 2947
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2388 0.4887
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.85
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.22541 0.30147 -7.382 1.56e-13 ***
taxa_pavrep0 0.30545 0.38315 0.797 0.4253
taxa_pavbird0 0.42075 0.37541 1.121 0.2624
taxa_pavamp0 0.46301 0.37248 1.243 0.2138
taxa_pavrep1 1.59462 0.30633 5.206 1.93e-07 ***
taxa_pavmam1 1.81125 0.30559 5.927 3.09e-09 ***
taxa_pavamp1 1.80556 0.30564 5.908 3.47e-09 ***
taxa_pavbird1 2.01855 0.30509 6.616 3.68e-11 ***
trafficmedium 0.18181 0.08060 2.256 0.0241 *
traffichigh 1.56252 0.13138 11.893 < 2e-16 ***
urban_area -0.26064 0.04394 -5.932 2.99e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
taxa_1km_m7 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m7) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + traffic + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7406.9 7478.8 -3691.4 7382.9 2948
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2778 0.5271
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.85
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1476 0.3030 -7.087 1.37e-12 ***
taxa_pavrep0 0.3044 0.3837 0.793 0.4276
taxa_pavbird0 0.4201 0.3759 1.118 0.2637
taxa_pavamp0 0.4597 0.3730 1.232 0.2178
taxa_pavrep1 1.5524 0.3084 5.034 4.79e-07 ***
taxa_pavmam1 1.7702 0.3076 5.754 8.70e-09 ***
taxa_pavamp1 1.7510 0.3077 5.691 1.26e-08 ***
taxa_pavbird1 1.9710 0.3071 6.418 1.38e-10 ***
trafficmedium 0.1487 0.0836 1.778 0.0754 .
traffichigh 1.3597 0.1328 10.235 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
taxa_1km_m8 <- glmmTMB(roadkill_count ~ taxa_pav + (1 | polygon),
family = nbinom2,
data = d1km)
summary(taxa_1km_m8) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + (1 | polygon)
Data: d1km
AIC BIC logLik deviance df.resid
7491.2 7551.1 -3735.6 7471.2 2950
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.48 0.6928
Number of obs: 2960, groups: polygon, 370
Dispersion parameter for nbinom2 family (): 1.88
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1540 0.3103 -6.942 3.86e-12 ***
taxa_pavrep0 0.3007 0.3854 0.780 0.435
taxa_pavbird0 0.4194 0.3774 1.111 0.267
taxa_pavamp0 0.4516 0.3747 1.205 0.228
taxa_pavrep1 1.7558 0.3172 5.535 3.11e-08 ***
taxa_pavmam1 1.9754 0.3166 6.240 4.37e-10 ***
taxa_pavamp1 1.9551 0.3166 6.176 6.57e-10 ***
taxa_pavbird1 2.1589 0.3162 6.827 8.66e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13km
- without cultivated land
taxa_3km_m1 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna + water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m1) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna +
water + shrubby + road_length + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7332.2 7439.9 -3648.1 7296.2 2918
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2176 0.4665
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.17200 0.30192 -7.194 6.29e-13 ***
taxa_pavrep0 0.25866 0.38809 0.666 0.5051
taxa_pavbird0 0.41937 0.37644 1.114 0.2653
taxa_pavamp0 0.47457 0.37110 1.279 0.2010
taxa_pavrep1 1.58785 0.30700 5.172 2.31e-07 ***
taxa_pavmam1 1.80729 0.30622 5.902 3.59e-09 ***
taxa_pavamp1 1.80374 0.30630 5.889 3.89e-09 ***
taxa_pavbird1 2.02157 0.30569 6.613 3.76e-11 ***
trafficmedium 0.10293 0.08650 1.190 0.2341
traffichigh 1.52651 0.13309 11.470 < 2e-16 ***
urban_area -0.28349 0.04469 -6.344 2.24e-10 ***
forest 0.08080 0.04083 1.979 0.0478 *
savanna -0.08106 0.03938 -2.059 0.0395 *
water 0.05434 0.03638 1.493 0.1353
shrubby -0.04049 0.04065 -0.996 0.3192
road_length 0.07576 0.03552 2.133 0.0329 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
taxa_3km_m2 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna + water + road_length + (1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m2) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna +
water + road_length + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7331.2 7432.9 -3648.6 7297.2 2919
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2192 0.4682
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.18218 0.30204 -7.225 5.01e-13 ***
taxa_pavrep0 0.25716 0.38823 0.662 0.5077
taxa_pavbird0 0.41893 0.37656 1.113 0.2659
taxa_pavamp0 0.47307 0.37124 1.274 0.2026
taxa_pavrep1 1.59568 0.30709 5.196 2.03e-07 ***
taxa_pavmam1 1.81420 0.30634 5.922 3.18e-09 ***
taxa_pavamp1 1.81125 0.30639 5.912 3.39e-09 ***
taxa_pavbird1 2.02886 0.30580 6.635 3.25e-11 ***
trafficmedium 0.10671 0.08654 1.233 0.2176
traffichigh 1.54358 0.13229 11.668 < 2e-16 ***
urban_area -0.28488 0.04477 -6.364 1.97e-10 ***
forest 0.08894 0.04009 2.218 0.0265 *
savanna -0.08405 0.03934 -2.137 0.0326 *
water 0.05054 0.03628 1.393 0.1637
road_length 0.07318 0.03546 2.063 0.0391 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
taxa_3km_m3 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + water + road_length + (1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m3) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + water +
road_length + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7333.7 7429.4 -3650.8 7301.7 2920
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2257 0.4751
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.16541 0.30227 -7.164 7.85e-13 ***
taxa_pavrep0 0.25774 0.38828 0.664 0.50682
taxa_pavbird0 0.41710 0.37660 1.108 0.26806
taxa_pavamp0 0.47514 0.37130 1.280 0.20066
taxa_pavrep1 1.58150 0.30729 5.147 2.65e-07 ***
taxa_pavmam1 1.80124 0.30655 5.876 4.21e-09 ***
taxa_pavamp1 1.79524 0.30659 5.856 4.75e-09 ***
taxa_pavbird1 2.01584 0.30601 6.587 4.47e-11 ***
trafficmedium 0.09762 0.08696 1.123 0.26160
traffichigh 1.54244 0.13347 11.557 < 2e-16 ***
urban_area -0.26743 0.04430 -6.037 1.57e-09 ***
forest 0.11141 0.03899 2.858 0.00427 **
water 0.04913 0.03658 1.343 0.17918
road_length 0.07217 0.03572 2.020 0.04334 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
taxa_3km_m4 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + road_length + (1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m4) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + road_length +
(1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7333.5 7423.2 -3651.7 7303.5 2921
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2272 0.4766
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.16732 0.30230 -7.169 7.53e-13 ***
taxa_pavrep0 0.25819 0.38826 0.665 0.50606
taxa_pavbird0 0.41714 0.37659 1.108 0.26800
taxa_pavamp0 0.47577 0.37128 1.281 0.20004
taxa_pavrep1 1.59388 0.30722 5.188 2.12e-07 ***
taxa_pavmam1 1.81335 0.30648 5.917 3.29e-09 ***
taxa_pavamp1 1.80863 0.30648 5.901 3.61e-09 ***
taxa_pavbird1 2.02784 0.30595 6.628 3.40e-11 ***
trafficmedium 0.08165 0.08623 0.947 0.34371
traffichigh 1.53669 0.13364 11.499 < 2e-16 ***
urban_area -0.26208 0.04417 -5.933 2.97e-09 ***
forest 0.10816 0.03897 2.775 0.00552 **
road_length 0.07242 0.03578 2.024 0.04297 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
taxa_3km_m5 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest +
(1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m5) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7335.5 7419.3 -3653.8 7307.5 2922
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2305 0.4801
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.14903 0.30201 -7.116 1.11e-12 ***
taxa_pavrep0 0.26055 0.38800 0.672 0.50189
taxa_pavbird0 0.41970 0.37633 1.115 0.26474
taxa_pavamp0 0.48575 0.37089 1.310 0.19030
taxa_pavrep1 1.57262 0.30695 5.123 3.00e-07 ***
taxa_pavmam1 1.79087 0.30619 5.849 4.95e-09 ***
taxa_pavamp1 1.78751 0.30622 5.837 5.30e-09 ***
taxa_pavbird1 2.00600 0.30566 6.563 5.28e-11 ***
trafficmedium 0.08222 0.08651 0.950 0.34187
traffichigh 1.54710 0.13415 11.533 < 2e-16 ***
urban_area -0.25512 0.04409 -5.786 7.20e-09 ***
forest 0.11191 0.03904 2.866 0.00415 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without forest
taxa_3km_m6 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area +
(1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m6) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7341.6 7419.4 -3657.8 7315.6 2923
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2426 0.4925
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.81
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.20144 0.30237 -7.281 3.32e-13 ***
taxa_pavrep0 0.25770 0.38811 0.664 0.5067
taxa_pavbird0 0.41787 0.37640 1.110 0.2669
taxa_pavamp0 0.49337 0.37089 1.330 0.1834
taxa_pavrep1 1.57416 0.30745 5.120 3.06e-07 ***
taxa_pavmam1 1.79072 0.30671 5.839 5.27e-09 ***
taxa_pavamp1 1.78695 0.30674 5.826 5.69e-09 ***
taxa_pavbird1 2.00676 0.30618 6.554 5.60e-11 ***
trafficmedium 0.17442 0.08121 2.148 0.0317 *
traffichigh 1.58356 0.13574 11.666 < 2e-16 ***
urban_area -0.25892 0.04464 -5.800 6.63e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
taxa_3km_m7 <- glmmTMB(roadkill_count ~ taxa_pav + traffic +
(1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m7) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + traffic + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7375.0 7446.8 -3675.5 7351.0 2924
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2828 0.5318
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.8
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.12167 0.30396 -6.980 2.95e-12 ***
taxa_pavrep0 0.25761 0.38862 0.663 0.507
taxa_pavbird0 0.41735 0.37690 1.107 0.268
taxa_pavamp0 0.48851 0.37149 1.315 0.189
taxa_pavrep1 1.53690 0.30958 4.964 6.89e-07 ***
taxa_pavmam1 1.75497 0.30885 5.682 1.33e-08 ***
taxa_pavamp1 1.73912 0.30887 5.631 1.80e-08 ***
taxa_pavbird1 1.96417 0.30831 6.371 1.88e-10 ***
trafficmedium 0.13619 0.08428 1.616 0.106
traffichigh 1.31522 0.13414 9.805 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
taxa_3km_m8 <- glmmTMB(roadkill_count ~ taxa_pav+
(1 | polygon),
family = nbinom2,
data = d3km)
summary(taxa_3km_m8) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + (1 | polygon)
Data: d3km
AIC BIC logLik deviance df.resid
7453.5 7513.4 -3716.8 7433.5 2926
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.4718 0.6868
Number of obs: 2936, groups: polygon, 367
Dispersion parameter for nbinom2 family (): 1.83
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1321 0.3110 -6.856 7.10e-12 ***
taxa_pavrep0 0.2571 0.3902 0.659 0.510
taxa_pavbird0 0.4159 0.3785 1.099 0.272
taxa_pavamp0 0.4747 0.3734 1.271 0.204
taxa_pavrep1 1.7346 0.3179 5.456 4.87e-08 ***
taxa_pavmam1 1.9548 0.3173 6.161 7.23e-10 ***
taxa_pavamp1 1.9389 0.3173 6.111 9.91e-10 ***
taxa_pavbird1 2.1480 0.3169 6.777 1.22e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15km
- without cultivated land
taxa_5km_m1 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna
+ water + shrubby + road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m1) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna +
water + shrubby + road_length + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7279.3 7386.8 -3621.6 7243.3 2894
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2339 0.4836
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.21371 0.30372 -7.289 3.13e-13 ***
taxa_pavrep0 0.30784 0.38515 0.799 0.424139
taxa_pavbird0 0.36658 0.38106 0.962 0.336047
taxa_pavamp0 0.55621 0.36727 1.514 0.129914
taxa_pavrep1 1.65517 0.30874 5.361 8.28e-08 ***
taxa_pavmam1 1.86152 0.30796 6.045 1.50e-09 ***
taxa_pavamp1 1.84792 0.30805 5.999 1.99e-09 ***
taxa_pavbird1 2.08479 0.30747 6.780 1.20e-11 ***
trafficmedium 0.08247 0.08851 0.932 0.351494
traffichigh 1.57441 0.14041 11.213 < 2e-16 ***
urban_area -0.26812 0.04483 -5.981 2.22e-09 ***
forest 0.09108 0.04388 2.076 0.037901 *
savanna -0.07609 0.04142 -1.837 0.066190 .
water 0.06199 0.03880 1.597 0.110158
shrubby -0.04461 0.04089 -1.091 0.275237
road_length 0.12105 0.03610 3.353 0.000799 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without shrubby
taxa_5km_m2 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna
+ water + road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m2) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna +
water + road_length + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7278.5 7380.1 -3622.2 7244.5 2895
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2353 0.4851
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.22438 0.30380 -7.322 2.45e-13 ***
taxa_pavrep0 0.30574 0.38529 0.794 0.427474
taxa_pavbird0 0.36678 0.38117 0.962 0.335934
taxa_pavamp0 0.55606 0.36739 1.514 0.130140
taxa_pavrep1 1.66681 0.30874 5.399 6.71e-08 ***
taxa_pavmam1 1.87225 0.30799 6.079 1.21e-09 ***
taxa_pavamp1 1.85898 0.30807 6.034 1.60e-09 ***
taxa_pavbird1 2.09571 0.30750 6.815 9.40e-12 ***
trafficmedium 0.08123 0.08864 0.916 0.359474
traffichigh 1.59245 0.13969 11.400 < 2e-16 ***
urban_area -0.26744 0.04490 -5.957 2.57e-09 ***
forest 0.10133 0.04293 2.360 0.018252 *
savanna -0.07771 0.04147 -1.874 0.060926 .
water 0.05341 0.03812 1.401 0.161206
road_length 0.12109 0.03613 3.352 0.000803 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without water
taxa_5km_m3 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna
+ road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m3) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + savanna +
road_length + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7278.4 7374.0 -3623.2 7246.4 2896
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2368 0.4866
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.22828 0.30383 -7.334 2.24e-13 ***
taxa_pavrep0 0.30621 0.38529 0.795 0.426760
taxa_pavbird0 0.36697 0.38117 0.963 0.335672
taxa_pavamp0 0.55643 0.36739 1.515 0.129887
taxa_pavrep1 1.68381 0.30858 5.457 4.85e-08 ***
taxa_pavmam1 1.88935 0.30782 6.138 8.36e-10 ***
taxa_pavamp1 1.87685 0.30787 6.096 1.09e-09 ***
taxa_pavbird1 2.11297 0.30733 6.875 6.18e-12 ***
trafficmedium 0.06290 0.08776 0.717 0.473537
traffichigh 1.57604 0.13940 11.306 < 2e-16 ***
urban_area -0.25992 0.04455 -5.835 5.38e-09 ***
forest 0.09515 0.04276 2.225 0.026067 *
savanna -0.07610 0.04149 -1.834 0.066621 .
road_length 0.12050 0.03618 3.331 0.000867 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without savanna
taxa_5km_m4 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area + forest
+ road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m4) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + forest + road_length +
(1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7279.7 7369.4 -3624.9 7249.7 2897
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2417 0.4917
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.21384 0.30404 -7.281 3.30e-13 ***
taxa_pavrep0 0.30778 0.38538 0.799 0.42450
taxa_pavbird0 0.36780 0.38123 0.965 0.33467
taxa_pavamp0 0.55979 0.36746 1.523 0.12766
taxa_pavrep1 1.67157 0.30874 5.414 6.16e-08 ***
taxa_pavmam1 1.87756 0.30799 6.096 1.09e-09 ***
taxa_pavamp1 1.86297 0.30802 6.048 1.46e-09 ***
taxa_pavbird1 2.10175 0.30750 6.835 8.21e-12 ***
trafficmedium 0.05831 0.08814 0.662 0.50827
traffichigh 1.55248 0.13968 11.115 < 2e-16 ***
urban_area -0.24579 0.04404 -5.581 2.39e-08 ***
forest 0.12323 0.04016 3.069 0.00215 **
road_length 0.11600 0.03628 3.197 0.00139 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without forest
taxa_5km_m5 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area
+ road_length + (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m5) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + road_length +
(1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7287.0 7370.6 -3629.5 7259.0 2898
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2568 0.5067
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.27586 0.30465 -7.470 8.00e-14 ***
taxa_pavrep0 0.30384 0.38567 0.788 0.430802
taxa_pavbird0 0.36540 0.38141 0.958 0.338042
taxa_pavamp0 0.56591 0.36757 1.540 0.123657
taxa_pavrep1 1.67450 0.30948 5.411 6.28e-08 ***
taxa_pavmam1 1.87894 0.30874 6.086 1.16e-09 ***
taxa_pavamp1 1.86422 0.30878 6.037 1.57e-09 ***
taxa_pavbird1 2.10488 0.30826 6.828 8.60e-12 ***
trafficmedium 0.16283 0.08245 1.975 0.048274 *
traffichigh 1.60046 0.14160 11.302 < 2e-16 ***
urban_area -0.24931 0.04473 -5.574 2.49e-08 ***
road_length 0.12621 0.03681 3.429 0.000606 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without road length
taxa_5km_m6 <- glmmTMB(roadkill_count ~ taxa_pav + traffic + urban_area
+ (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m6) Family: nbinom2 ( log )
Formula:
roadkill_count ~ taxa_pav + traffic + urban_area + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7296.6 7374.3 -3635.3 7270.6 2899
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.2692 0.5188
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.78
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.24648 0.30452 -7.377 1.62e-13 ***
taxa_pavrep0 0.30588 0.38520 0.794 0.427
taxa_pavbird0 0.36662 0.38093 0.962 0.336
taxa_pavamp0 0.58262 0.36686 1.588 0.112
taxa_pavrep1 1.63582 0.30932 5.288 1.23e-07 ***
taxa_pavmam1 1.83815 0.30857 5.957 2.57e-09 ***
taxa_pavamp1 1.82599 0.30864 5.916 3.29e-09 ***
taxa_pavbird1 2.06542 0.30809 6.704 2.03e-11 ***
trafficmedium 0.17209 0.08338 2.064 0.039 *
traffichigh 1.61533 0.14365 11.245 < 2e-16 ***
urban_area -0.23622 0.04495 -5.256 1.48e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without urban area
taxa_5km_m7 <- glmmTMB(roadkill_count ~ taxa_pav + traffic
+ (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m7) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + traffic + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7323.4 7395.1 -3649.7 7299.4 2900
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.3037 0.5511
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.77
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1685 0.3057 -7.093 1.31e-12 ***
taxa_pavrep0 0.3047 0.3856 0.790 0.4294
taxa_pavbird0 0.3654 0.3813 0.958 0.3380
taxa_pavamp0 0.5775 0.3673 1.572 0.1159
taxa_pavrep1 1.5893 0.3110 5.110 3.22e-07 ***
taxa_pavmam1 1.7924 0.3103 5.777 7.60e-09 ***
taxa_pavamp1 1.7706 0.3103 5.706 1.16e-08 ***
taxa_pavbird1 2.0163 0.3098 6.509 7.57e-11 ***
trafficmedium 0.1455 0.0860 1.692 0.0906 .
traffichigh 1.3630 0.1411 9.659 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- without traffic
taxa_5km_m8 <- glmmTMB(roadkill_count ~ taxa_pav
+ (1 | polygon),
family = nbinom2,
data = d5km)
summary(taxa_5km_m8) Family: nbinom2 ( log )
Formula: roadkill_count ~ taxa_pav + (1 | polygon)
Data: d5km
AIC BIC logLik deviance df.resid
7398.1 7457.9 -3689.1 7378.1 2902
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
polygon (Intercept) 0.5 0.7071
Number of obs: 2912, groups: polygon, 364
Dispersion parameter for nbinom2 family (): 1.8
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1490 0.3130 -6.867 6.56e-12 ***
taxa_pavrep0 0.3012 0.3874 0.778 0.437
taxa_pavbird0 0.3624 0.3829 0.946 0.344
taxa_pavamp0 0.5639 0.3693 1.527 0.127
taxa_pavrep1 1.7562 0.3200 5.488 4.07e-08 ***
taxa_pavmam1 1.9564 0.3194 6.125 9.06e-10 ***
taxa_pavamp1 1.9383 0.3194 6.068 1.29e-09 ***
taxa_pavbird1 2.1692 0.3190 6.799 1.05e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICc Table
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))
n <- nobs(model)
aic <- AIC(model)
aicc <- aic + (2 * k * (k + 1)) / (n - k - 1)
return(aicc)
}Lists of all models
# 1km scale models
models_1km <- list(
taxa_1km_m0, taxa_1km_m1, taxa_1km_m2,
taxa_1km_m3, taxa_1km_m4, taxa_1km_m5,
taxa_1km_m6, taxa_1km_m7, taxa_1km_m8
)
# 3km scale models
models_3km <- list(
taxa_3km_m0, taxa_3km_m1, taxa_3km_m2,
taxa_3km_m3, taxa_3km_m4, taxa_3km_m5,
taxa_3km_m6, taxa_3km_m7, taxa_3km_m8
)
# 5km scale models
models_5km <- list(
taxa_5km_m0, taxa_5km_m1, taxa_5km_m2,
taxa_5km_m3, taxa_5km_m4, taxa_5km_m5,
taxa_5km_m6, taxa_5km_m7, taxa_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) Model AICc delta_AICc Akaike_weight
5 1km_m4 7356.147 0.00000000 2.385743e-01
4 1km_m3 7356.214 0.06664722 2.307551e-01
3 1km_m2 7356.297 0.15028873 2.213038e-01
2 1km_m1 7357.317 1.16965508 1.329341e-01
6 1km_m5 7357.417 1.27031657 1.264090e-01
1 1km_m0 7359.276 3.12872951 4.991472e-02
7 1km_m6 7371.528 15.38069995 1.090805e-04
8 1km_m7 7406.864 50.71665642 2.315483e-12
9 1km_m8 7491.182 135.03470150 1.135562e-30print(results_3km) Model AICc delta_AICc Akaike_weight
3 3km_m2 7331.166 0.000000 3.914402e-01
2 3km_m1 7332.168 1.002153 2.371650e-01
5 3km_m4 7333.466 2.299873 1.239523e-01
4 3km_m3 7333.683 2.516606 1.112222e-01
1 3km_m0 7334.099 2.933031 9.031624e-02
6 3km_m5 7335.547 4.380692 4.379357e-02
7 3km_m6 7341.612 10.445853 2.110461e-03
8 3km_m7 7375.034 43.867463 1.166721e-10
9 3km_m8 7453.523 122.357253 1.054692e-27print(results_5km) Model AICc delta_AICc Akaike_weight
4 5km_m3 7278.415 0.0000000 2.941957e-01
3 5km_m2 7278.474 0.0589771 2.856469e-01
2 5km_m1 7279.279 0.8641103 1.909840e-01
5 5km_m4 7279.754 1.3388944 1.506257e-01
1 5km_m0 7281.163 2.7478560 7.446409e-02
6 5km_m5 7286.986 8.5706596 4.050803e-03
7 5km_m6 7296.618 18.2023801 3.281252e-05
8 5km_m7 7323.366 44.9504963 5.102230e-11
9 5km_m8 7398.140 119.7243442 2.956822e-27Nakagawa’s R²
Calculation R²
The marginal R² considers only the variance of the fixed effects (without the random effects), while the conditional R² takes both the fixed and random effects into account (i.e., the total model).
models_all <- c(models_1km, models_3km, models_5km)
# Junta os nomes dos modelos em uma única lista
model_names_all <- c(model_names_1km, model_names_3km, model_names_5km)
# Aplica r2_nakagawa a todos os modelos e nomeia o resultado
r2_results <- setNames(lapply(models_all, r2_nakagawa), model_names_all)
# Imprime os resultados
cat("\nPseudo R² of the models:\n")
Pseudo R² of the models:print(r2_results)$`1km_m0`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.311
$`1km_m1`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.311
$`1km_m2`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.310
$`1km_m3`
# R2 for Mixed Models
Conditional R2: 0.433
Marginal R2: 0.308
$`1km_m4`
# R2 for Mixed Models
Conditional R2: 0.434
Marginal R2: 0.306
$`1km_m5`
# R2 for Mixed Models
Conditional R2: 0.434
Marginal R2: 0.303
$`1km_m6`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.292
$`1km_m7`
# R2 for Mixed Models
Conditional R2: 0.432
Marginal R2: 0.265
$`1km_m8`
# R2 for Mixed Models
Conditional R2: 0.462
Marginal R2: 0.188
$`3km_m0`
# R2 for Mixed Models
Conditional R2: 0.428
Marginal R2: 0.298
$`3km_m1`
# R2 for Mixed Models
Conditional R2: 0.428
Marginal R2: 0.298
$`3km_m2`
# R2 for Mixed Models
Conditional R2: 0.429
Marginal R2: 0.298
$`3km_m3`
# R2 for Mixed Models
Conditional R2: 0.430
Marginal R2: 0.294
$`3km_m4`
# R2 for Mixed Models
Conditional R2: 0.429
Marginal R2: 0.293
$`3km_m5`
# R2 for Mixed Models
Conditional R2: 0.428
Marginal R2: 0.289
$`3km_m6`
# R2 for Mixed Models
Conditional R2: 0.429
Marginal R2: 0.284
$`3km_m7`
# R2 for Mixed Models
Conditional R2: 0.428
Marginal R2: 0.258
$`3km_m8`
# R2 for Mixed Models
Conditional R2: 0.456
Marginal R2: 0.186
$`5km_m0`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.297
$`5km_m1`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.297
$`5km_m2`
# R2 for Mixed Models
Conditional R2: 0.436
Marginal R2: 0.297
$`5km_m3`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.295
$`5km_m4`
# R2 for Mixed Models
Conditional R2: 0.436
Marginal R2: 0.293
$`5km_m5`
# R2 for Mixed Models
Conditional R2: 0.439
Marginal R2: 0.287
$`5km_m6`
# R2 for Mixed Models
Conditional R2: 0.436
Marginal R2: 0.277
$`5km_m7`
# R2 for Mixed Models
Conditional R2: 0.435
Marginal R2: 0.255
$`5km_m8`
# R2 for Mixed Models
Conditional R2: 0.464
Marginal R2: 0.183General table
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)
r2_vals <- lapply(model_list, r2_nakagawa)
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),
R2_marginal = sapply(r2_vals, function(x) x$R2_marginal),
R2_conditional = sapply(r2_vals, function(x) x$R2_conditional),
stringsAsFactors = FALSE
) |>
dplyr::arrange(AICc)
}
# Renomeia modelos
names(models_1km) <- paste0("1km_m", 0:8)
names(models_3km) <- paste0("3km_m", 0:8)
names(models_5km) <- paste0("5km_m", 0:8)
# Gera a tabela final
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)
# Marca o melhor modelo
final_table$Best <- ifelse(final_table$Delta_AICc == 0, "★", "")
# Exporta para Excel
writexl::write_xlsx(final_table, "GLMM_final_model_selection.xlsx")Advanced diagnostics of the model residuals
simulationOutput <- simulateResiduals(fittedModel = taxa_1km_m4, plot = TRUE)
simulationOutput <- simulateResiduals(fittedModel = taxa_3km_m2, plot = TRUE)
simulationOutput <- simulateResiduals(fittedModel = taxa_5km_m3, plot = TRUE)
Since the outliers in the datasets are real data (and not errors) primarily related to seasonality, these data will be retained even if the model is not able to represent them very well.
Summary
# Function to extract estimates and p-values, adding an asterisk for significance
extract_est_p <- function(model_summary) {
coef_table <- model_summary$coefficients$cond # Accessing coefficients
# Create the table with estimates and p-values (rounded)
estimates <- round(coef_table[, "Estimate"], 3)
p_values <- round(coef_table[, "Pr(>|z|)"], 3)
# Add an asterisk to significant p-values
p_values <- ifelse(p_values < 0.05, paste0(p_values, "*"), as.character(p_values))
result <- data.frame(
Term = rownames(coef_table),
Estimate = estimates,
P_value = p_values,
row.names = NULL
)
return(result)
}
# Extract data for the three models
df_1km <- extract_est_p(summary(taxa_1km_m4))
df_3km <- extract_est_p(summary(taxa_3km_m2))
df_5km <- extract_est_p(summary(taxa_5km_m3))
# Comparison across scales
GLMM_scale_comparison <- df_1km %>%
rename(Estimate_1km = Estimate, P_1km = P_value) %>%
full_join(df_3km %>% rename(Estimate_3km = Estimate, P_3km = P_value), by = "Term") %>%
full_join(df_5km %>% rename(Estimate_5km = Estimate, P_5km = P_value), by = "Term")
# Display the comparison table with asterisks for significant p-values
GLMM_scale_comparison Term Estimate_1km P_1km Estimate_3km P_3km Estimate_5km P_5km
1 (Intercept) -2.201 0* -2.182 0* -2.228 0*
2 taxa_pavrep0 0.308 0.42 0.257 0.508 0.306 0.427
3 taxa_pavbird0 0.427 0.255 0.419 0.266 0.367 0.336
4 taxa_pavamp0 0.453 0.224 0.473 0.203 0.556 0.13
5 taxa_pavrep1 1.621 0* 1.596 0* 1.684 0*
6 taxa_pavmam1 1.837 0* 1.814 0* 1.889 0*
7 taxa_pavamp1 1.834 0* 1.811 0* 1.877 0*
8 taxa_pavbird1 2.043 0* 2.029 0* 2.113 0*
9 trafficmedium 0.095 0.263 0.107 0.218 0.063 0.474
10 traffichigh 1.548 0* 1.544 0* 1.576 0*
11 cult_land -0.082 0.07 NA <NA> NA <NA>
12 urban_area -0.313 0* -0.285 0* -0.260 0*
13 savanna -0.201 0* -0.084 0.033* -0.076 0.067
14 forest NA <NA> 0.089 0.027* 0.095 0.026*
15 water NA <NA> 0.051 0.164 NA <NA>
16 road_length NA <NA> 0.073 0.039* 0.121 0.001*write_xlsx(GLMM_scale_comparison, "GLMM_summary_comparison.xlsx")Summary Plot
Function to extract model coefficients
extract_glmmTMB <- function(model, label) {
sum <- summary(model)
coefs_cond <- sum$coefficients$cond
if (is.null(coefs_cond)) stop("Componente condicional não encontrado")
coefs <- as.data.frame(coefs_cond)
coefs$term <- rownames(coefs)
colnames(coefs) <- c("estimate", "std_error", "z_value", "p_value", "term")
coefs$model <- label
return(coefs[, c("term", "estimate", "std_error", "z_value", "p_value", "model")])
}Extraction of coefficients and organization into a single data frame
coef_1km <- extract_glmmTMB(taxa_1km_m4, "1 km")
coef_3km <- extract_glmmTMB(taxa_3km_m2, "3 km")
coef_5km <- extract_glmmTMB(taxa_5km_m3, "5 km")
coef_data_all <- bind_rows(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",
str_detect(term, "taxa_pavrep0") ~ "Reptile-unpaved",
str_detect(term, "taxa_pavrep1") ~ "Reptile-paved",
str_detect(term, "taxa_pavbird0") ~ "Bird-unpaved",
str_detect(term, "taxa_pavbird1") ~ "Bird-paved",
str_detect(term, "taxa_pavamp0") ~ "Amphibia-unpaved",
str_detect(term, "taxa_pavamp1") ~ "Amphibia-paved",
str_detect(term, "taxa_pavmam0") ~ "Mammal-unpaved",
str_detect(term, "taxa_pavmam1") ~ "Mammal-paved",
str_detect(term, "water") ~ "Water"
)
) %>%
filter(term != "(Intercept)")coef_data_all$p_value_num <- as.numeric(gsub("< 0.001", "0.00099", coef_data_all$p_value))
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
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") coef_data_all <- coef_data_all %>%
mutate(term_num = as.numeric(factor(term)))
coef_data_all <- coef_data_all %>%
filter(term != "(Intercept)")
coef_data_all$term <- factor(coef_data_all$term,
levels = unique(coef_data_all$term[order(coef_data_all$estimate)]))
Summary_coeff_GLMM<- 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_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_GLMM
#ggsave("GLMM_comparison.png", plot = Summary_coeff_GLMM, height = 8, width = 10, dpi = 300)Plots
Relative rate: pavement
d1km <- d1km %>%
mutate(
pavement_label = ifelse(pavement == 0, "Unpaved", "Paved"),
taxa_pav = paste(taxa, pavement, sep = "") # cria "mam0", "rep1", etc.
) %>%
mutate(
taxa_pav = recode_factor(taxa_pav,
"mam0" = "mammal_unpaved",
"rep0" = "reptile_unpaved",
"amp0" = "amphibia_unpaved",
"bird0" = "bird_unpaved",
"mam1" = "mammal_paved",
"rep1" = "reptile_paved",
"amp1" = "amphibia_paved",
"bird1" = "bird_paved"
)
)d1km <- d1km %>%
mutate(
pavement_label = factor(pavement_label,
levels = c("Unpaved", "Paved")),
taxa = factor(taxa,
levels = c("amp", "rep", "bird", "mam"))
)
taxa_rate_pavement_effect<-ggplot(d1km, aes(x = pavement_label, y = relative_rate, color = taxa)) +
stat_summary(fun = mean, geom = "point", size = 3,
position = position_dodge(width = 0.6)) +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.2,
position = position_dodge(width = 0.6)) +
scale_color_manual(values = c(
"mam" = "#D55E00", # laranja escuro
"rep" = "#0072B2", # azul
"amp" = "#009E73", # verde
"bird" = "#CC79A7" # rosa
),
labels = c(
"mam" = "Mammals",
"rep" = "Reptiles",
"amp" = "Amphibians",
"bird" = "Birds"
)
) +
labs(x = "Pavement", y = "WVC rate (individuals km⁻¹)", color = "Taxa") +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black")
) +
ylim(0, NA)Relative rate: traffic
d1km <- d1km %>%
mutate(
traffic = factor(traffic,
levels = c("low", "medium", "high")),
taxa = factor(taxa,
levels = c("amp", "rep", "bird", "mam"))
)
taxa_traffic_effect <- ggplot(
d1km,
aes(x = traffic,
y = relative_rate,
color = taxa,
group = taxa)
) +
stat_summary(
fun = mean,
geom = "point",
size = 3,
position = position_dodge(width = 0.6)
) +
stat_summary(
fun.data = mean_cl_boot,
geom = "errorbar",
width = 0.2,
position = position_dodge(width = 0.6)
) +
scale_color_manual(
values = c(
"mam" = "#D55E00",
"rep" = "#0072B2",
"amp" = "#009E73",
"bird" = "#CC79A7"
),
labels = c(
"mam" = "Mammals",
"rep" = "Reptiles",
"amp" = "Amphibians",
"bird" = "Birds"
)
) +
labs(
x = "Traffic",
y = "WVC rate (individuals km⁻¹)",
color = "Taxa"
) + scale_x_discrete(labels = c("low" = "Low", "medium" = "Medium", "high" = "High"))+
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.line = element_line(color = "black")
) +
ylim(0, NA)