Mixed Effect Model - Roadkill hotspot v. coldspot

Question

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.

Answer 1

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.

Answer 2

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.