从具有随机截距和二次项的 GLMM 生成空间预测(栅格)

Generating a spatial prediction (raster) from a GLMM with a random intercept and a quadratic term

如标题所示,我正在尝试生成描述相对使用概率的预测栅格。为了创建一个可重现的示例,我使用了 maxlike 包中的 MaungaWhau 数据。 MaungaWhau 是一个 list,它包含两个栅格图层以及一个包含 1000 个位置的图层。所以,有了这些数据和包...

library(maxlike)
library(lme4)
data(MaungaWhau)

我们可以制作两个栅格图层,一个栅格堆栈,以及一个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)

然后我制作了一个包含 IndID 的新数据框:每个人的唯一 ID (AAA - DDD); Used:表示点是否被使用的二进制响应变量(分别为1或0);以及 SpatialPoints object.

中每个点的 ElevPrecip
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

然后我构建了一个资源选择模型 (RSF) 并将 IndID 指定为随机效应。另请注意,我为 Elev.

添加了一个二次项
#Make model
Mod <- glmer(Used ~ Elev + I(Elev^2) + Precip + (1|IndID), family=binomial, data = df, verbos = 1) 
summary(Mod)

鉴于已经使用过的和可用的点,我对解释不感兴趣。我的第一个问题是更多的确认。 raster 包小插图指出 "The names in the Raster object should exactly match those expected by the model." 在与 I(Elev^2) 拟合的二次项的实例中,我是否正确预测 Elev 将是 'looking'?似乎是这种情况,因为在下面的 predict 代码中没有与 Elev 相关的错误。

其次,如何处理随机截距项(1|IndID)?我对边际预测感兴趣,不需要考虑个人。

用下面的代码

#Change names of layers in raster stack to match model
names(rs) <- c("Elev", "Precip")

Pred <- predict(rs, Mod)

我收到错误:

Error in eval(expr, envir, enclos) : object 'IndID' not found

是否可以在不将 IndID 协变量传递给 predict 函数的情况下为 'typical' 个体生成边际预测?换句话说,我想在制作预测表面时忽略 IndID 项和对截距的相关单独调整。

与其依赖 predict 函数,不如通过首先创建固定效应的 beta 对象来手动生成预测

betas <- fixef(Mod)

然后通过将每个栅格乘以各自的 beta 系数来生成预测。

Pred <- betas[1] + (elev * betas[2]) + (elev^2 * betas[3]) + (precip * betas[4])
plot(Pred)

然后很容易添加或删除截距并手动指定 link 函数(例如 logit)。

lme4 (merMod) 对象的 predict 函数默认进行条件预测。

要进行 marginal/unconditional 预测,您需要使用 re.form 参数。您的代码如下所示:

Pred <- predict(rs, Mod, re.form = NA)

如果您还想在响应变量的范围内进行预测,您可以使用 type 参数。有关可用参数的更多详细信息,请参阅帮助页面,?predict.merMod