通过分类变量和连续变量的交互可视化 GLMM 预测

Question

我正在 R 中使用 GLMM，其中混合了连续变量和分类变量，并具有一些交互作用。我使用 MuMIn 中的 dredge 和 model.avg 函数来获取每个变量的效果估计。我的问题是如何最好地绘制结果。我想制作一个图表，显示一个变量（森林）对我的数据的影响，其中趋势线反映森林参数估计，但我无法弄清楚如何将分类变量和交互变量保持在“平均值”，以便趋势线仅反映森林的影响。

这是模型和绘图设置：

#load packages and document
cuckoo<-read.table("http://www.acsu.buffalo.edu/~ciaranwi/home_range.txt", 
header=T,sep="\t")
require(lme4)
require(MuMIn)
as.factor (cuckoo$ID)
as.factor (cuckoo$Sex)
as.factor(cuckoo$MS_bin)
options(na.action = "na.fail")

# create global model and fit
fm<- lmer(log(KD_95)~ MS_bin + Forest + NDVI + Sex + Precip + MS_bin*Forest 
+ MS_bin*NDVI  + MS_bin*Sex + MS_bin*Precip + Argos + Sample + (1|ID), data 
= cuckoo, REML = FALSE)

# dredge but always include argos and sample
KD95<-dredge(fm,fixed=c("Argos","Sample"))

# model averaging 
avgmod<-model.avg(KD95, fit=TRUE)
summary(avgmod)

#plot data
plot(cuckoo$Forest, (log(cuckoo$KD_95)),
 xlab = "Mean percentage of forest cover",
 ylab = expression(paste(plain("Log of Kernel density estimate, 95%    
utilisation, km"^{2}))),
 pch = c(15,17)[as.numeric(cuckoo$MS_bin)],  
 main = "Forest cover",
 col="black", 
 ylim=c(14,23))
legend(80,22, c("Breeding","Nonbreeding"), pch=c(15, 17),  cex=0.7)

然后我就陷入了如何包含趋势线的困境。到目前为止我已经：

#parameter estimates from model.avg
argos_est<- -1.6
MS_est<- -1.77
samp_est<-0.01
forest_est<--0.02
sex_est<-0.0653
precip_est<-0.0004
ndvi_est<--0.00003
model_intercept<-22.7

#calculate mean values for parameters
argos_mean<-mean(cuckoo$Argos)
samp_mean<-mean(cuckoo$Sample)
forest_mean<-mean(cuckoo$Forest)
ndvi_mean<-mean(cuckoo$NDVI)
precip_mean<-mean(cuckoo$Precip)

#calculate the intercept and add trend line
intercept<-(model_intercept + (forest_est*cuckoo$Forest) +    
(argos_est*argos_mean) + (samp_est * samp_mean) + (ndvi_est*ndvi_mean) +  
(precip_est*precip_mean) )

abline(intercept, forest_est)

但这没有考虑交互作用或分类变量，并且截距看起来太高了。有什么想法吗？

Answer 1

在流程方面，您可以利用 R 在模型对象中存储大量有关模型的信息并具有从模型对象中获取信息的函数这一事实，使您的编码变得更加容易。例如，

coef(avgmod)

将为您提供模型系数，而

predict(avgmod)

将为您提供用于拟合模型的数据框中每个观测值的模型预测。

为了可视化对我们感兴趣的数据值的特定组合的预测，请创建一个新的数据框，其中包含我们想要保持不变的变量的平均值，以及我们想要改变的变量的一系列值（例如

 Forest

）。

expand.grid

使用下面列出的值的所有组合创建一个数据框。

pred.data = expand.grid(Argos=mean(cuckoo$Argos), Sample=mean(cuckoo$Sample), 
                        Precip=mean(cuckoo$Precip), NDVI=mean(cuckoo$NDVI), 
                        Sex="M", Forest=seq(0,100,10), MS_bin=unique(cuckoo$MS_bin), 
                        ID=unique(cuckoo$ID))

现在我们使用

predict

函数将 log(KD_95) 的预测添加到该数据框。

predict

负责计算您输入的任何数据的模型预测（假设您给它一个包含模型中所有变量的数据框）。

pred.data$lgKD_95_pred = predict(avgmod, newdata=pred.data)

现在我们绘制结果。

geom_point

绘制点，如原始图中所示，然后

geom_line

添加

MS_bin

每个级别的预测（以及 Sex="M"）。

library(ggplot2)

ggplot() +
  geom_point(data=cuckoo, aes(Forest, log(KD_95), shape=factor(MS_bin), 
             colour=factor(MS_bin), size=3)) +
  geom_line(data=pred.data, aes(Forest, lKD_95_pred, colour=factor(MS_bin)))

结果如下：

更新： 要绘制男性和女性的回归线，只需在

pred.data

中包含 Sex="F" 并在图中添加

Sex

作为美感。在下面的示例中，我在绘制点时使用不同的形状来标记

Sex

，并使用不同的线型来标记回归线的

Sex

。

pred.data = expand.grid(Argos=mean(cuckoo$Argos), Sample=mean(cuckoo$Sample), 
                        Precip=mean(cuckoo$Precip), NDVI=mean(cuckoo$NDVI), 
                        Sex=unique(cuckoo$Sex), Forest=seq(0,100,10), MS_bin=unique(cuckoo$MS_bin), 
                        ID=unique(cuckoo$ID))

pred.data$lgKD_95_pred = predict(avgmod, newdata=pred.data)

ggplot() +
  geom_point(data=cuckoo, aes(Forest, log(KD_95), shape=Sex, 
                              colour=factor(MS_bin)), size=3) +
  geom_line(data=pred.data, aes(Forest, lgKD_95_pred, linetype=Sex, 
                                colour=factor(MS_bin))) 

Answer 2

我希望我没有错过重点，但如果你想要一个线性趋势，你实际上不必手动计算所有内容，而是获取你绘制的内容并拟合 y~x 线性回归模型，如下所示：

model = lm(log(cuckoo$KD_95)~cuckoo$Forest)

model

# Call:
#   lm(formula = log(cuckoo$KD_95) ~ cuckoo$Forest)
# 
# Coefficients:
#   (Intercept)  cuckoo$Forest  
#      17.13698       -0.01461 

abline(17.13698 ,  -0.01461, col="red")

红线使用回归拟合的截距和斜率。黑线是您的手动过程。