Spatial Analysis of Mule Deer Road Carcass Reports along San Juan County’s Major Roadways
Jessie Shapiro
MSGIS | University of Utah | Fall 2020 | GEOG 6000
Introduction
As urban areas expand into wildlife ranges, it has gained importance to identify areas at risk of a collision site. Roads and vehicles impact wildlife in a numerous ways but the most significant effect is road kill. With the increase in Utah’s population, this problem is growing in scope affecting many lives and wildlife populations. Wildlife tends to cross roadways more frequently during their seasonal migrations when they are in search for better food and water. In 2019, San Juan was identified as one of the high valued hot spots in the State of Utah for wildlife vehicle collisions to occur (Cramer, 2019).
To mitigate for wildlife-vehicle collisions, state agencies collaborated to build wildlife fencing and crossing structures in the high frequency collisions sites, just south of Monticello. In 2005 and 2007, the first fencing and crossings were installed along US191. Wildlife vehicle collisions persisted between 2006-2013 so another fence was installed between milepost 66-70 in 2016. In 2019, three more wildlife crossing structures and fencing to merge the 2005 and 2016 fencing was completed. A statistical analysis should be completed to test the effectiveness of these fencing projects on collisions rates.
This report will explore the influence of seasonal phenomena by partitioning the mule deer carcass report count data into the four migration seasons. This analysis will test the effectiveness of wildlife fencing, investigate influences of road carcass reports at the county spatial extent, and finally use a hierarchical Bayesian approach to assess mule deer vehicle collisions along road segments and predict the likelihood threshold of an accident occurring.
Objective:
- Explore spatial patterns and frequencies in mule deer road carcass reports between the migration seasons
- Test the effectiveness of 2005, 2016, and 2019 wildlife fencing and wildlife crossing structures
- Run a broader spatial model to investigate variables that might lead to collisions at the county spatial extent.
Dataset
The prime data source includes 2013-2020 UDWR mule deer carcass reports and major roadways for the San Juan unit split into equal one-mile road segments. In spatial analysis of crash data, it is consistently common to segment the road networks into equal distances (Li, Linhua. 2006). Carcass reports were then assigned to each segment based on its geographic location over the 2013-2020-time period.
Brownian bridge movement models for the San Juan mule deer were computed for migration corridors, summer range and winter range. These results were then converted into binary variables with zero representing absence and one representing presence of that model intersecting a road segment. From the movement models, average migration dates were calculated and used to identify the seasonal classification for each carcass report (i.e. spring, fall, summer, winter). Mule deer habitat, wildlife fencing, and wildlife crossing variables were treated similarly with zero and one representing absence and presence of the variable intersecting a road segment.
Annual average daily traffic (AADT) a measure of traffic volumes was used as a variable provided by UDOT. AADT for each season was calculated using the seasonal factors. Speed limit will also be considered as a variable. These covariates will be used to test their effect on road carcass counts at the San Juan County extent.
Point Pattern Analysis
Density plots and frequencies tests were computed to see if there were any existing spatial patterns in the data set. From the density plots, spatial overlap is observed with significant clustering in the town of Monticello. Ripley’s K was further tested and all seasonal plots suggest significant clustering and deviations from CSR. The significance level in the Monte Carlo test was 0.02 for all seasons.
Kernel Density
fit0 <- ppm(roadkill.ppp ~ 1)
fit0## Stationary multitype Poisson process
##
## Possible marks: 'fall', 'spring', 'summer' and 'winter'
##
## Uniform intensity for each mark level: 0.00000005758216
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) -16.67005 0.01810715 -16.70554 -16.63456 *** -920.6338exp(coef(fit0))## log(lambda)
## 0.00000005758216summary(roadkill.ppp)## Marked planar point pattern: 3050 points
## Average intensity 0.0000002303287 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
##
## Multitype:
## frequency proportion intensity
## fall 150 0.04918033 0.000000011327640
## spring 65 0.02131148 0.000000004908643
## summer 2145 0.70327870 0.000000161985200
## winter 690 0.22622950 0.000000052107140
##
## Window: rectangle = [602681, 671893.4] x [4117882, 4309206] units
## (69210 x 191300 units)
## Window area = 13241900000 square unitsfit1 <- ppm(roadkill.ppp ~ polynom(x, y, 2))
fit1## Error in solve.default(M) :
## system is computationally singular: reciprocal condition number = 4.76114e-37
## Error in solve.default(M) :
## system is computationally singular: reciprocal condition number = 4.76114e-37
## Nonstationary multitype Poisson process
##
## Possible marks: 'fall', 'spring', 'summer' and 'winter'
##
## Log intensity: ~x + y + I(x^2) + I(x * y) + I(y^2)
##
## Fitted trend coefficients:
## (Intercept) x y
## -53700.481953848844569 0.018795135325848 0.022713102094499
## I(x^2) I(x * y) I(y^2)
## -0.000000009207425 -0.000000001639207 -0.000000002582775
##
## Standard errors unavailable; variance-covariance matrix is singularplot(fit1,
how = 'image',
se = FALSE,
pause = FALSE)
Ripley’s K
#Fall
fall = split(roadkill.ppp)$fall
roadkill.kest = Kest(fall)
fall.mc = envelope(fall, fun='Kest', nsim=99, verbose=FALSE)
plot(fall.mc, shade=c("hi", "lo"))
summary(fall.mc)## Pointwise critical envelopes for K(r)
## and observed value for 'fall'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Upper envelope: pointwise maximum of simulated curves
## Lower envelope: pointwise minimum of simulated curves
## Significance level of Monte Carlo test: 2/100 = 0.02
## Data: fall#Spring
spring = split(roadkill.ppp)$spring
spring.kest = Kest(spring)
spring.mc = envelope(spring, fun='Kest', nsim=99, verbose=FALSE)
plot(spring.mc, shade=c("hi", "lo"))
summary(spring.mc)## Pointwise critical envelopes for K(r)
## and observed value for 'spring'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Upper envelope: pointwise maximum of simulated curves
## Lower envelope: pointwise minimum of simulated curves
## Significance level of Monte Carlo test: 2/100 = 0.02
## Data: spring#Summer
summer = split(roadkill.ppp)$summer
summer.kest = Kest(summer)
summer.mc = envelope(summer, fun='Kest', nsim=99, verbose=FALSE)
plot(summer.mc, shade=c("hi", "lo"))
summary(summer.mc)## Pointwise critical envelopes for K(r)
## and observed value for 'summer'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Upper envelope: pointwise maximum of simulated curves
## Lower envelope: pointwise minimum of simulated curves
## Significance level of Monte Carlo test: 2/100 = 0.02
## Data: summer#Winter
winter = split(roadkill.ppp)$winter
roadkill.kest = Kest(winter)
winter.mc = envelope(winter, fun='Kest', nsim=99, verbose=FALSE)
plot(winter.mc, shade=c("hi", "lo"))
summary(winter.mc)## Pointwise critical envelopes for K(r)
## and observed value for 'winter'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Upper envelope: pointwise maximum of simulated curves
## Lower envelope: pointwise minimum of simulated curves
## Significance level of Monte Carlo test: 2/100 = 0.02
## Data: winterPoisson Distribution Model
Using the mule deer road carcass report point data, a Poisson model was fit to explore the difference in counts between the different seasons. The p-value of <2.2e-16, from the goodness of fit test indicates statistically significant difference between the seasons. The coefficients indicate that carcass reports increase in summer, winter and fall and decrease in spring
seasons.mod1<- glm(n ~ Season, family = poisson, data = seasons)
summary(seasons.mod1)##
## Call:
## glm(formula = n ~ Season, family = poisson, data = seasons)
##
## Deviance Residuals:
## [1] 0 0 0 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.01064 0.08165 61.368 < 0.0000000000000002 ***
## Seasonspring -0.83625 0.14850 -5.631 0.0000000179 ***
## Seasonsummer 2.66026 0.08446 31.499 < 0.0000000000000002 ***
## Seasonwinter 1.52606 0.09009 16.940 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3491.35886279723717962 on 3 degrees of freedom
## Residual deviance: 0.00000000000024736 on 0 degrees of freedom
## AIC: 38.748
##
## Number of Fisher Scoring iterations: 2exp(coef(seasons.mod1))## (Intercept) Seasonspring Seasonsummer Seasonwinter
## 150.0000000 0.4333333 14.3000000 4.6000000anova(seasons.mod1, test="Chisq")| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 3 | 3.49e+03 | |||
| 3 | 3.49e+03 | 0 | 2.47e-13 | 0 |
It’s important to understand how road carcass reports differ between years. When modeling the interaction of seasons and year, the model produced statistically significant results in ANOVA for seasons: year (p = 1.839934187339e-5), season (p= 2.2e-16) and year (p= 1.819e-14). This suggests that the relationship of seasons on years can vary by year.
seasons.year.mod1<- glm(Count ~ Year + Season, family = poisson, data = year)
summary(seasons.year.mod1)##
## Call:
## glm(formula = Count ~ Year + Season, family = poisson, data = year)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -9.155 -1.117 0.249 1.722 4.897
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.790430 0.124522 30.440 < 0.0000000000000002 ***
## Year -0.061097 0.007997 -7.640 0.0000000000000217 ***
## Seasonspring -0.941556 0.148531 -6.339 0.0000000002310533 ***
## Seasonsummer 2.554952 0.084524 30.227 < 0.0000000000000002 ***
## Seasonwinter 1.420748 0.090152 15.759 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3670.88 on 30 degrees of freedom
## Residual deviance: 274.41 on 26 degrees of freedom
## AIC: 456.93
##
## Number of Fisher Scoring iterations: 4seasons.year.mod2<- glm(Count ~ Season * Year, family = poisson, data = year)
summary(seasons.year.mod2)##
## Call:
## glm(formula = Count ~ Season * Year, family = poisson, data = year)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3615 -1.3192 0.1229 1.8452 3.8500
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 3.064725145064810 0.496654859883746 6.171
## Seasonspring -0.478551879486663 0.840858573750286 -0.569
## Seasonsummer 3.585356604136543 0.510236595681257 7.027
## Seasonwinter 1.344175742169117 0.539805748264784 2.490
## Year -0.000000000001704 0.040824789707465 0.000
## Seasonspring:Year -0.039627696729080 0.067909840793774 -0.584
## Seasonsummer:Year -0.086246091715364 0.041923997448833 -2.057
## Seasonwinter:Year 0.003864796828683 0.044076369011406 0.088
## Pr(>|z|)
## (Intercept) 0.00000000067974 ***
## Seasonspring 0.5693
## Seasonsummer 0.00000000000211 ***
## Seasonwinter 0.0128 *
## Year 1.0000
## Seasonspring:Year 0.5595
## Seasonsummer:Year 0.0397 *
## Seasonwinter:Year 0.9301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3670.88 on 30 degrees of freedom
## Residual deviance: 249.77 on 23 degrees of freedom
## AIC: 438.3
##
## Number of Fisher Scoring iterations: 4anova(seasons.year.mod2, test="Chisq")| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 30 | 3.67e+03 | |||
| 3 | 3.34e+03 | 27 | 333 | 0 |
| 1 | 58.7 | 26 | 274 | 1.82e-14 |
| 3 | 24.6 | 23 | 250 | 1.84e-05 |
Intervention of Wildlife Fencing and Wildlife Crossings
The influence of fencing and wildlife crossings on road carcass reports was investigated at a small spatial scale, for only fenced in road segments. Mule deer count data was split based on fence installation projects and to equal sampling years. The road segments of interest stretched for approximately 14miles of road along US191 south of Monticello. From 2013 to 2020, road carcass reports on these fenced road segments declined by 27% in the summer, 41% in the winter, 23% during spring migration and 24% during fall migration.
A negative binomial model was fit for 2013-2015 data to see if the Devil’s canyon 2005 fencing project were effective in decreasing road carcass reports. The fencing proved to be statistically significant, with a negative effect on counts for summer (p = 1.483e-2, coef = -1.3809), winter (p = 5.792e-6, coef = -2.96409), fall (p= 4.177e-5, coef = -22.1040), and spring (p = 2.953e-3, coef = -21.0585). Wildlife crossings indicated to have a negative effect on road carcass counts in summer (coef = -0.3177, p = 0.6583), winter (coef = -0.04213, p =0.94530), and spring (coef =-19.7823, p = 0.124), but it did not play a significant role in the overall model fit. Fall was also insignificant (coef = 0.7577, p = 0.2125).
Spring
##
## Call:
## glm.nb(formula = PreFSpring ~ Fence05 + Cross05, data = fencing,
## init.theta = 21071.4974, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.58112 -0.23175 -0.00004 -0.00001 1.32388
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2231 0.3162 0.706 0.480
## Fence05 -21.0585 10142.7040 -0.002 0.998
## Cross05 -19.7823 10717.1145 -0.002 0.999
##
## (Dispersion parameter for Negative Binomial(21071.5) family taken to be 1)
##
## Null deviance: 18.8657 on 13 degrees of freedom
## Residual deviance: 7.6736 on 11 degrees of freedom
## AIC: 29.893
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 21071
## Std. Err.: 1052966
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -21.893## (Intercept) Fence05 Cross05
## 1.2500000000000000 0.0000000007152009 0.0000000025623579| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 13 | 18.9 | |||
| 1 | 8.84 | 12 | 10 | 0.00295 |
| 1 | 2.36 | 11 | 7.67 | 0.125 |
Summer
##
## Call:
## glm.nb(formula = PreFSum ~ Fence05 + Cross05, data = fencing,
## init.theta = 1.252556996, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0965 -0.7803 -0.3435 0.3722 1.4386
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.0307 0.3113 12.948 < 0.0000000000000002 ***
## Fence05 -1.3809 0.5189 -2.661 0.00778 **
## Cross05 -0.3177 0.7107 -0.447 0.65483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.2526) family taken to be 1)
##
## Null deviance: 22.713 on 13 degrees of freedom
## Residual deviance: 16.585 on 11 degrees of freedom
## AIC: 133.99
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.253
## Std. Err.: 0.494
##
## 2 x log-likelihood: -125.985## (Intercept) Fence05 Cross05
## 56.2976645 0.2513473 0.7278152| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 13 | 22.7 | |||
| 1 | 5.94 | 12 | 16.8 | 0.0148 |
| 1 | 0.191 | 11 | 16.6 | 0.662 |
Winter
##
## Call:
## glm.nb(formula = PreFWint ~ Fence05 + Cross05, data = fencing,
## init.theta = 11.10762994, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5876 -0.6322 -0.1880 0.6477 1.2453
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.36280 0.20731 6.574 0.0000000000491 ***
## Fence05 -2.96409 1.02931 -2.880 0.00398 **
## Cross05 -0.04213 0.61409 -0.069 0.94530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(11.1076) family taken to be 1)
##
## Null deviance: 35.809 on 13 degrees of freedom
## Residual deviance: 15.249 on 11 degrees of freedom
## AIC: 53.881
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 11.1
## Std. Err.: 24.9
##
## 2 x log-likelihood: -45.881## (Intercept) Fence05 Cross05
## 3.90710568 0.05160733 0.95874345| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 13 | 35.8 | |||
| 1 | 20.6 | 12 | 15.3 | 5.79e-06 |
| 1 | 0.0048 | 11 | 15.2 | 0.945 |
Fall
##
## Call:
## glm.nb(formula = PreFFall ~ Fence05 + Cross05, data = fencing,
## init.theta = 39107.85051, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.93647 -0.00004 -0.00003 0.09030 0.75498
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6286 0.2582 2.435 0.0149 *
## Fence05 -22.1040 11280.4018 -0.002 0.9984
## Cross05 0.7577 0.5628 1.346 0.1782
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(39107.85) family taken to be 1)
##
## Null deviance: 24.8054 on 13 degrees of freedom
## Residual deviance: 6.4617 on 11 degrees of freedom
## AIC: 35.931
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 39108
## Std. Err.: 1323045
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -27.931## (Intercept) Fence05 Cross05
## 1.8749999999999998 0.0000000002513961 2.1333333333333333| Df | Deviance | Resid. Df | Resid. Dev | Pr(>Chi) |
|---|---|---|---|---|
| 13 | 24.8 | |||
| 1 | 16.8 | 12 | 8.02 | 4.18e-05 |
| 1 | 1.55 | 11 | 6.46 | 0.212 |
Since the 2005 fence was installed, road carcass counts stayed relatively low for these road segments. A pairwise t test was used next to compare effects of fencing projects 2016 and 2019 on the mule deer carcass count data. This test subtracts the before and after values and tests these values against zero.
Summer
#2013-2015 vs. 2016-2018
t.test(fencing$PreFSum,fencing$PostFSum, paired = TRUE)##
## Paired t-test
##
## data: fencing$PreFSum and fencing$PostFSum
## t = 2.3596, df = 13, p-value = 0.0346
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.459771 33.111657
## sample estimates:
## mean of the differences
## 17.28571The results of the summer pairwise t-test comparing 2007 fencing and 2016 fencing effects proved significant (p = 0.0346, t = 2.3596) and confidence interval [1.45, 33] with a mean difference between the pairs of 17.28.
#2013-2015 vs. 2019-2020
t.test(fencing$PreFSum,fencing$NewFSum, paired = TRUE)##
## Paired t-test
##
## data: fencing$PreFSum and fencing$NewFSum
## t = 3.6992, df = 13, p-value = 0.002674
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 13.87613 52.83815
## sample estimates:
## mean of the differences
## 33.35714When comparing the effects of 2007 fencing to the 2019 fencing installation, it was also very significant (p= 0.002675, t = 3.6992) on count data, confidence interval [13.876,52] and mean difference between the pairs of 33.
Winter
#2013-2015 vs. 2016-2018
t.test(fencing$PreFWint,fencing$PostFWint, paired = TRUE)##
## Paired t-test
##
## data: fencing$PreFWint and fencing$PostFWint
## t = -3.2284, df = 13, p-value = 0.006597
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.4575967 -0.6852605
## sample estimates:
## mean of the differences
## -2.071429When comparing the effects of 2007 fencing to 2016 fencing during the winter season, it was also significant (p = 0.0065, t = -3.2), confidence interval [-3.45, -0.68] and a mean difference of -2.
Spring
#2013-2015 vs. 2016-2018
t.test(fencing$PreFSpring,fencing$PostFSprin, paired = TRUE)##
## Paired t-test
##
## data: fencing$PreFSpring and fencing$PostFSprin
## t = 0.23421, df = 13, p-value = 0.8185
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5874454 0.7303025
## sample estimates:
## mean of the differences
## 0.07142857In the spring, fencing projects weren’t statistically significant until the 2019 fence installation. Compared to 2007 to 2019, the t-test was significant (p = 0.04454, t = 2.2234), confidence interval [0.01823416, 1.26748012] with a mean difference 0.64. When comparing 2016 fencing to 2019, the pairwise t-test was significant (t=2.2804, p=0.04009), confidence interval [0.0300661, 1.1127910] and a mean difference of 0.57.
#2013-2015 vs. 2018-2020
t.test(fencing$PreFSpring,fencing$NewFSpring, paired = TRUE)##
## Paired t-test
##
## data: fencing$PreFSpring and fencing$NewFSpring
## t = 2.2234, df = 13, p-value = 0.04454
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.01823416 1.26748012
## sample estimates:
## mean of the differences
## 0.6428571Fencing effects in fall showed to be insignificant in the t-test, which could partly be due to the low sample size.
Negative Binomial Model
While these road segments significantly declined in road carcass counts from the intervention of fencing, fencing was not included as a covariate when looking at the county spatial extent. This is because areas of fencing are highly targeted because of the high road carcass counts. Therefore, when looking at a larger spatial scale, fencing would misrepresent its effectiveness on road carcass reports because of the high concentration along the identified road segments for fencing, in comparison to the entire county where there are a lot of road segments equal to zero.
To understand what might influence road carcass reports at the county scale, AADT and Brownian Bridge modeled estimations (i.e. migration corridor, summer range and winter range) were fit to a negative binomial distribution as covariates. Speed limit was excluded because of its insignificance in the tested model. This dataset is a very good candidate for a Negative Binomial Regression since we can see the large number of observations with zero count. Negative binomial models include an extra parameter to account for over dispersion in the model. An over dispersion test was computed on a Poisson distribution to validate the decision of using a negative binomial distribution model for summer (dispersion = 28.14, p = 0.01394), winter (dispersion = 7.28, p = 0.00007), fall (dispersion =3.25, p= 0.01588) and spring (dispersion = 2.131298, p = 0.02051).
The reported negative binomial regression coefficient is interpreted as follows: for a one-unit change in the independent variable, the difference in the logs of the expected counts is expected to change by the respective regression coefficient, given that the other independent variables in the model are held constant. The negative binomial model can be interpreted in the following equation:
log(Road carcass count) = Intercept + AADT + BBMM
From the ANOVA test results, AADT, summer range, winter range and migration corridors for the respective season contributed significantly to the quality of the fit of the model. As annual average daily traffic volumes increased, it positively influenced the number of road carcass reports. For the respective season and Brownian bridge movement model, all covariates were statistically significant with a positive effect on road carcass reports.
Spring<-glm.nb(SpringCnt~ AADTSprin + MigCorr, data =roads1)
Summer<-glm.nb(SumCnt~ AADTSum + SumRg, data =roads1)
Winter<-glm.nb(WintCnt~ AADTWint + WintRg, data =roads1)
Fall<-glm.nb(FallCnt~ AADTFall + MigCorr, data =roads1)
export_summs(Spring, Summer, Winter, Fall, scale = TRUE)| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| (Intercept) | -4.73 *** | -1.71 *** | -2.50 *** | -3.99 *** |
| (0.32) | (0.18) | (0.17) | (0.26) | |
| AADTSprin | 0.85 *** | |||
| (0.16) | ||||
| MigCorr | 2.44 *** | 1.86 *** | ||
| (0.60) | (0.56) | |||
| AADTSum | 2.10 *** | |||
| (0.16) | ||||
| SumRg | 2.70 *** | |||
| (0.72) | ||||
| AADTWint | 1.70 *** | |||
| (0.13) | ||||
| WintRg | 1.50 * | |||
| (0.60) | ||||
| AADTFall | 1.64 *** | |||
| (0.15) | ||||
| N | 1640 | 1640 | 1640 | 1640 |
| AIC | 272.40 | 1270.42 | 994.79 | 478.47 |
| BIC | 294.01 | 1292.03 | 1016.40 | 500.08 |
| Pseudo R2 | 0.11 | 0.09 | 0.10 | 0.13 |
| All continuous predictors are mean-centered and scaled by 1 standard deviation. *** p < 0.001; ** p < 0.01; * p < 0.05. | ||||
Traffic data for 2018 -2020 is unavailable so count data for those years were excluded from this model. AADT seasonal traffic volumes were the highest in the summer (M= 1370, Max = 19912, Min = 0) and spring (M=1197, Max = 20045, Min = 0), and lower during the winter (M = 650, Max = 8270, Min = 0) and fall (M= 1135, Max = 17387, Min = 0) months. The high traffic volumes were highly correlated with the frequency of summer road carcass counts (cor = 0.113, p = 0.000003724) 
INLA Spatial Model
This is still in the testing phase.
total_roadkill = sum(roads1$SumCnt)
total_roadsgements = nrow(roads1)
#Standard Rate - rate per segment (rate per mile)
roads1$EI = total_roadkill/ total_roadsgements
#Setting Priors
u = 0.2/0.31
alpha = 0.01
rw1_prior <- list(theta = list(prior="pc.prec", param=c(u, alpha)))
f1 <- SumCnt ~ scale(AADTSum) + scale(SumRg) + f(idarea, model = "bym", graph = g, hyper =rw1_prior)
m.pois.bym <- inla(f1,
family = "nbinomial", data = roads1,
E = roads1$EI, control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
summary(m.pois.bym)##
## Call:
## c("inla(formula = f1, family = \"nbinomial\", data = roads1, E =
## roads1$EI, ", " control.compute = list(dic = TRUE, waic = TRUE),
## control.predictor = list(compute = TRUE))" )
## Time used:
## Pre = 3.05, Running = 28.9, Post = 0.338, Total = 32.3
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 1046.387 13007.673 -20341.801 998.881 22346.073 -20529.057
## scale(AADTSum) -3.019 0.024 -7.118 -3.019 1.079 -2.869
## scale(SumRg) 0.562 0.009 -0.923 0.562 111.915 0.535
## kld
## (Intercept) 7.383
## scale(AADTSum) 2844192.932
## scale(SumRg) 711135.571
##
## Random effects:
## Name Model
## idarea BYM model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (1/overdispersion) 75.033 1.041
## Precision for idarea (iid component) 1730.797 1799.420
## Precision for idarea (spatial component) 0.536 0.109
## 0.025quant 0.5quant
## size for the nbinomial observations (1/overdispersion) 73.630 74.825
## Precision for idarea (iid component) 122.002 1193.911
## Precision for idarea (spatial component) 0.338 0.531
## 0.975quant mode
## size for the nbinomial observations (1/overdispersion) 77.525 73.952
## Precision for idarea (iid component) 6500.523 335.800
## Precision for idarea (spatial component) 0.763 0.526
##
## Expected number of effective parameters(stdev): 133.01(4.28)
## Number of equivalent replicates : 12.33
##
## Watanabe-Akaike information criterion (WAIC) ...: 8696597407.51
## Effective number of parameters .................: 4348298239.96
##
## Marginal log-Likelihood: 933.82
## Posterior marginals for the linear predictor and
## the fitted values are computed#PLot covariates
Efxplot(list(m.pois.bym))
Discussion
The significance of covariates changed depending on the spatial scale of the model. At the county scale, there were larger frequencies of extreme observations, exceeding from the Poisson distribution. A negative binomial model was used in this case because the count data variance exceeded the sample mean (Miaou, 1994). The variable fencing only occurred where there were high counts of road carcass reports. At the county scale, this inevitably would show a positive relationship between road carcass reports and fencing due to the large portion of road segments equal to zero, misrepresenting the effectiveness of wildlife fencing.
Several covariates were originally considered for the negative binomial model at the county scale. Environment characteristics of urban, suburban and rural were going to be fit, but after a comparison of classifications, the state classifies over 80% of San Juan County as rural. Most negative binomial models also include an offset variable to account for disparate exposures in the data. It was decided against using AADT as an offset because this reported wanted to directly quantify the relationship between traffic volumes and carcass frequencies. Mule deer habitat was also considered as a covariate but caused errors in the summer negative binomial model. This was likely due to the variable to dominate the fit due to huge outliers since over 50% of San Juan county roadways fell within mule deer habitat.
Conclusions
In this report, 2005, 2016 and 2019 fencing projects proved to be effective, decreasing the number of road carcass frequencies in summer and winter. The data set also showed that carcass counts varied between mule deer migrations with a noticeable trend in higher summer count frequencies, which is suggested to be due to their seasonal movements and high summer traffic volumes. At the county spatial extent, increasing traffic volumes and presence of movement’s model had a positive influence, increasing the number of road carcass counts. The point pattern analysis also clearly suggested clustering patterns in the town of Monticello for all seasons. This report was a preliminary analysis on possible influences of mule deer carcasses reports on roadways in San Juan County and gave insight on the effectiveness of mitigation techniques. Other seasonal covariates should be considered that might affect road carcass reports at the county extent. We need to continue to be proactive in identifying where collisions occur and accommodate for these wildlife movements in San Juan County that have proven to influence road carcass counts.
References
Cramer, Patricia, E. Vasquez, A.Jones. 2019. Identification of wildlife-vehicle conflict priority hotspots in Utah summary report.
Galgamuwa, U et al., Bayesian spatial modeling to incorporate unmeasured information at road segment levels with the INLA approach: A methodological advancement of estimating crash modification factors, Journal of Traffic and Transportation Engineering (English Edition), https://doi.org/10.1016/j.jtte.2019.03.003
Gelman, Andrew, Jessica Hwang, Aki Vehatri. 2013. Understanding predictive information criteria for Bayesian models.
Li, Linhua. 2006. A GIS-Based Bayesion Approach for Analyzing Spatial-Temporal Patterns of Traffic Crashes. Texas A & M Unviersity.
Miaou, S.-P., 1994. The relationship between truck accidents and geometric design of road sections: Poisson versus negative binominal regressions. Accident Analysis and Prevention, 26 (4), 471-482.