As indicated in the title, I am trying to generate a predictive raster depicting the relative probability of use. To create a reproducible example, I have used the MaungaWhau data from the maxlike package. MaungaWhau is a list that contains two raster layers as well a layer of 1000 locations. So, with these data and packages...
library(maxlike)
library(lme4)
data(MaungaWhau)
we can make two raster layers, a raster stack, as well as a SpatialPoints object.
elev <- raster(MaungaWhau$elev, xmn=0, xmx=61, ymn=0, ymx=87)
precip <- raster(MaungaWhau$precip, xmn=0, xmx=61, ymn=0, ymx=87)
rs <- stack(elev, precip)
PointDat <- SpatialPoints(MaungaWhau$xy)
I then make a new dataframe that contains IndID: the unique ID for each individual (AAA - DDD); Used: the binary response variable indicating whether the point was used or not (1 or 0, respectively); as well as the Elev and Precip values for each point from the SpatialPoints object.
df <- data.frame(IndID = sample(c("AAA", "BBB", "CCC", "DDD"), 1000, replace = T),
Used = sample(c(0,1),1000, replace = T),
Elev = extract(elev, PointDat),
Precip = extract(precip, PointDat))
head(df); tail(df)
> head(df); tail(df)
IndID Used Elev Precip
1 DDD 0 0.3798393 0.6405494
2 DDD 1 0.8830846 1.1174869
3 AAA 0 1.9282864 0.9641432
4 DDD 0 1.5024634 0.4695881
5 BBB 1 1.3089075 -0.1341483
6 BBB 1 0.5733952 0.6246699
I then build a resource selection model (RSF) and specify IndID as a random effect. Notice also, that I included a quadratic term for Elev.
#Make model
Mod <- glmer(Used ~ Elev + I(Elev^2) + Precip + (1|IndID), family=binomial, data = df, verbos = 1)
summary(Mod)
I am not interested in interpretation given the made up used and available points. My first question is more if a confirmation. The raster package vignette states that "The names in the Raster object should exactly match those expected by the model." In the instance of the quadratic term fit with I(Elev^2) am I correct that predict will be 'looking' for Elev? This seems to be the case as there in no error associated with Elev in the predict code below.
Secondly, how do I deal with the random intercept term (1|IndID)? I am interested in marginal predictions and do not need to account for individuals.
With the following code
#Change names of layers in raster stack to match model
names(rs) <- c("Elev", "Precip")
Pred <- predict(rs, Mod)
I get the error:
Error in eval(expr, envir, enclos) : object 'IndID' not found
Is it possible to generate a marginal prediction for the 'typical' individual without passing the IndID covariate to the predict function? In other words, I want to ignore the IndID term and the associated individual adjustments to the intercept when making the prediction surface.
The predict function for lme4 (merMod) objects makes conditional predictions by default.
To make marginal/unconditional predictions, you need to make use of the re.form argument. Your code would look like:
Pred <- predict(rs, Mod, re.form = NA)
If you also wanted predictions done on the scale of the response variable, you can use the type argument. See the help page for more details on the available arguments, ?predict.merMod.