# Mixed Effect Model - Roadkill hotspot v. coldspot

This problem concerns road mortality. We have 4 highway sections of varied length that were flagged as problem areas for road mortality with large mammals (moose, deer, elk, and even cows). We used KDE+ software to map hotspot segments of varied length using ten years of collision study data. I also mapped an equal number of intervening coldspot segments where no kills have occurred. The number of hot/cold segments per highway is as follows::

Field crews collected data on ditch depth, % fence along segment length, roadside vegetation, and adjacent land features (pasture, forested, industrial). Land data was imported into a GIS (rasterized) and we obtained spatial metrics on land cover and edge lengths within a 250 m buffer run along both sides of the highway for Developed, Shrubland, Agriculture, Coniferous, Deciduous, and Mixed forest. Traffic volume, number of dwellings within 500 m, and distance to nearest town center from the center of the segment were also calculated.

Okay - we have lots of variables to deal with here - too many in my opinion, but I'm working with a team of researchers that think we need to look at all the details. I've combed through the data (data exploration) using methods described in Zuur et al. (2009) to plot and investigate normality, heterogeneity, and multi-colinearity and have reduced the data set to 11 variables.

There are two kinds of Dependent variables:

• Hotspots data only: DVkills = (number of kills / length of hot segment) *100
• Binary: Hotspot / Coldspot segment

Zuur et al. (2009) provide an example in chapter 16 of their textbook concerning amphibian road mortality using negative binomial GAM and GAMM to address the underlying question: "Is there a relationship between amphibian roadkills and any of the dependent variables?" However, they run into issues with more coefficients than data in the model and discuss shrinkage smoothers, which is beyond my capabilities.

I found this PeerJ article to be helpful and their r-code is instructive. These authors modeled the relationship between roadkill counts and predictor variables by fitting a generalized linear mixed-effect model using the lme4 package.

However, I do not understand a couple things about their model:

"To account for spatial autocorrelation we included sampling unit, nested within the most common habitat type, as random terms in the model, including an offset to account for the length of the sampling unit." (p. 11)

I understand what spatial autocorrelation is, but I don't understand how their model accounts for this. They created 41 x 10 km sampling units along the roadway they included in their analysis. Pieces of code used in their model are as follows:

# Calculate the log of the survey area (road section length) to add it as an offset in the model
RKdata2 <- within(RKdata2, {LogArea = log(RKdata2\$LENGTH)})


The model looks as follows:

glmer.nb(roadkill ~
season +
scale(wetland) +
scale(owls) +
scale(guineafowl) +
scale(speed) +
scale(traffic) +
scale(roadWidth) +
scale(infrastructure) +
offset(log(LogArea)) +
(1|habitat/section), data = RKdata2)


My questions are as follows:

1. I looked in lme4 package to learn that offset is "can be used to specify an a priori known component to be included in the linear predictor during fitting". I have four highway segments - should I take the log of their cumulative length and used this as an offset? I don't quite understand the purpose of the offset. The lme4 package keeps a hyperlink to See model.offset, but this goes nowhere - I can't find this.

2. I have four separate highway study areas and the lengths of the hotspot and coldspot segments are varied. Should I nest these as the random component? For example: (1|Hwy/segmentL), where "segmentL" is the log of the sum length of all cold/hotspot segment lengths, and would this address spatial autocorrelation?

3. I ran a binary linear regression to compare the hotspot v. coldspot segments. The linear mixed-effect model dependent variable will be the number of kills within a segment divided by the length of the segment. There are no kills in the coldspot segments, so I have to analyze this separately. Is there a better way to do this where I include both hotspots and coldspots as a binary dependent variable and number of kills in hotspots as a second dependent variable?

I'm adding a Clevland dotplot of the # of road kills by highway to show that the # of kills varies by segment and by highway. A histogram of the # of kills is also provided.

## 2 Answers

Further to @mdewey's excellent answers, I wanted to add the following.

For your question 2, if you use a statement like (1|Hwy/segment) in your model, the implications are that (i) you have multiple higways included in your study, (ii) for each highway, you have multiple segments and (iii) for each segment, you have multiple values of your dependent variable (though some segments - but not all - might possibly have a single value for the dependent variable). In other words, both Hwy and segment would be so-called random grouping variables (i.e., categorical variable or, in R parlance, factors). In practice, you would want the highways included in your study to be representative of a larger group of highways. Similarly, you would want the segments within each highway to be representative of a larger group of segments.

However, what you suggested in your question 2, (1|Hwy/segmentL), does not seem to make sense, since:

1. segmentL is NOT a grouping variable, but a continuous variable;
2. You seemingly have a single value of the dependent variable per segment (unless I am missing out on something).

One way to think about mixed effects models is as a tool for modelling repeated values of a response variable collected across levels/categories of one or more random grouping variables. The classic example is the repeated measures setting, where your grouping variable is subject and each subject's dependent variable is repeatedly measured over time, in which case you would include something like (1|subject) in the random effects part of your model. Notice that this notation implicitly assumes there are repeated values of the response variable available for each subject.

Scenario 1

If you truly have one value of the dependent variable per segment (e.g., a single roadkill count per segment), then you can't include segment as a random grouping variable in your model. But if you think your 4 highways can be conceived of as being representative of a larger group/set of highways to which you wish to generalize the findings of your study, you could include Hwy as a random grouping factor in your model: (1|Hwy). In that case, for each highway, you would have multiple roadkill counts, i.e., one roadkill per segment.

Scenario 2

If you have repeated values of the response variable per segment (e.g., a roadkill count for the daytime and a roadkill count for the night time), then you can include segment as a random grouping variable in your model along perhaps with highway like so: (1|Hwy/segment). This notation implicitly indicates that each segment contributes 2 values of the dependent variable roadkill count to the model.

If you wanted to include an offset for Scenario 1, the offset would have to be applied to the values of the dependent variables collected at the lowest level of your hierarchy (which is the segment level). So your offset would be the log-transformed segment length.

Sinilarly, if you wanted to include an offset for Scenario 2, the offset would have to be applied to the values of the dependent variables collected at the lowest level of your hierarchy (which, again, is the segment level). So your offset would again be the log-transformed segment length.

If your dependent variable is segment-level roadkill count and you work under Scenario 1, then using an offset in your model which is equal to log(length) would enable you to model the expected roadkill count per meter of highway for a typical highway. Here, length is the segment length (with length measured in meters, say). Under Scenario 2, using the same offset of log(length) would enable you to model the expected roadkill count per meter of highway for a typical segment of a typical highway.

• Thanks Isabella - I've read and fully understood your post. I added a Clevland dotplot and histogram of the # road kill counts - dependent variable. It seems that Scenario 1 is most applicable. The hr of road kill is likely too unreliable as the data was collected by a mix of maintenance contractors that work to clear through the day. We are getting new and better data for 2017-2018 that will allow this more detailed type of analysis. – Mark Thompson Nov 26 '18 at 1:29
• You're welcome, Mark! I guess you could incorporate a coarse variable like "time of day when workers removed the road kill" in your modelling if need be (am/pm, etc.). I think you are "washing away" the temporal aspect of your data, essentially treating all segment-level values of the dependent variable for each highway as coming from the same year? If you are interested in using your model to predict what happens next year, say, you could perhaps treat year itself as a random grouping variable? – Isabella Ghement Nov 26 '18 at 1:48
• I did not separate the data by year. I will be looking into this further, but we are waiting for more data and a better data set. We used the KDE+ software to identify hotspots using location and the sum of all kills across all years. The dotplot has all years lumped together. We are doing a contract report and then planning a publication. I will do a more detailed analysis for the publication where we will sort through time of day, year, and I am also thinking of doing an analysis without KDE+ and just using linear distance of each kill location to center of nearest town. – Mark Thompson Nov 26 '18 at 2:19
• (+1) which seems to solve all the problems except the issue of the cold spots which I think may have to be modelled separately. – mdewey Nov 26 '18 at 16:36
• Thank you so much, @mdewey! I really appreciate you taking the time to review and comment on my answer. I am learning a lot both from your answers and your comment! 🤗 – Isabella Ghement Nov 26 '18 at 16:47

Q1 - the purpose of an offset here would be to account for the fact that longer stretches of highway have more chance of having a roadkill. You are modelling log roadkill so you include log length as an offset. It implies that you expect twice the length to have twice the roadkill.

Q2 - you just include the factor which identifies segments, not its length. The model then includes a random intercept for segment which accounts for correlation within segment.

Q3 - this is more tricky as I do not see any way in which you can include a count of roadkill in the model when you deliberately sampled sites with zero roadkill. You probably need to model the roadkill sites separately.

• Thanks for the very simple response. This was helpful. Q1 - offset: I created EV (explanatory variable) by dividing the length of the segment by the number of kills, which I did to account for the different lengths. Perhaps I should just use the raw roadkill counts in segments and use offset in the formula? Is there a difference in the way that offset handles this data? – Mark Thompson Nov 15 '18 at 17:34
• When I use raw counts I don't get warnings for non-integers, so it solves this problem. However, I would like to understand how offset accounts for the length and kill rate - is there a resource that could describe how this works? I like to have a complete understanding of what I put into my analysis. – Mark Thompson Nov 15 '18 at 18:07
• An offset is a variable which is included in the model with a fixed coefficient rather than an estimated one. The software will include it with a coefficient of 1 which is usually what you want. – mdewey Nov 15 '18 at 18:20
• The lme4 package states: If you do see convergence warnings, and want to trouble-shoot/double-check the results, the following steps are recommended (examples are given below): • double-check the model specification and the data for mistakes • center and scale continuous predictor variables (e.g. with scale)." The authors of the cited article used scale for their applicable variables, but I don't see where they get convergence warnings - I ran their script. Where would this come from? – Mark Thompson Nov 16 '18 at 8:30
• @mdewey: Inspired by your excellent answers, I wanted to delve a bit deeper into Mark's question 2. Would you be kind enough to read my answer and let me know if it needs any further refining? Thanks a lot! – Isabella Ghement Nov 26 '18 at 0:21