基于R的模型平均系数的部分残差图。

I'm using the R package MuMIn to do multimodel inference and the function model.avg to average the coefficients estimated by a set of models. To visually compare the data to the estimated relationships based on the averaged coefficients, I want to use partial residual plots, similar to the ones created by the crPlots function of the car package. I've tried three ways and I'm not sure whether any is appropriate. Here is a demonstration.

我用R包MuMIn做多模型推理和函数模型。平均一组模型估计的系数。为了直观地比较数据与基于平均系数的估计关系，我想使用部分残差图，类似于car package的crplot函数所创建的那些残差图。我试过三种方法，但不确定是否合适。这是一个演示。

library(MuMIn)
# Loading the data
data(Cement)
# Creating a full model with all the covariates we are interested in
fullModel <- lm(y ~ ., data = Cement, na.action=na.fail)
# Getting all possible models based on the covariates of the full model
muModel <- dredge(fullModel)
# Averaging across all models
avgModel <- model.avg(muModel)
# Getting the averaged coefficients
coefMod <- coef(avgModel)
coefMod
# (Intercept)          X1          X2          X4          X3 
# 65.71487660  1.45607957  0.61085531 -0.49776089 -0.07148454 

Option 1: Using crPlots

选项1:使用crPlots

library(car) # For crPlots
# Creating a duplicate of the fullMode
hackModel <- fullModel
# Changing the coefficents to the averaged coefficients
hackModel$coefficients <- coefMod[names(coef(fullModel))]
# Changing the residuals
hackModel$residuals <- Cement$y - predict(hackModel)
# Plot the hacked model vs the full model
layout(matrix(1:8, nrow=2, byrow=TRUE))
crPlots(hackModel, layout=NA)
crPlots(fullModel, layout=NA)

Notice that the crPlots of the full and hacked versions with the average coefficients are different.

注意，具有平均系数的完整版本和被破解版本的crplot是不同的。

The question here is: Is this appropriate? The results rely on a hack that I found in this answer. Do I need to change parts of the model other than the residuals and the coefficients?

这里的问题是:这合适吗?结果取决于我在这个答案中发现的一个技巧。除了残差和系数之外，我还需要改变模型的其他部分吗?

Option 2: Homemade plots

选项2:自制的情节

# Partial residuals: residuals(hacked model) + beta*x
# X1
# Get partial residuals
prX1 <- resid(hackModel) + coefMod["X1"]*Cement$X1
# Plot the partial residuals
plot(prX1 ~ Cement$X1)
# Add modeled relationship
abline(a=0,b=coefMod["X1"])
# X2 - X4
plot(resid(hackModel) + coefMod["X2"]*X2 ~ X2, data=Cement); abline(a=0,b=coefMod["X2"])
plot(resid(hackModel) + coefMod["X3"]*X3 ~ X3, data=Cement); abline(a=0,b=coefMod["X3"])
plot(resid(hackModel) + coefMod["X4"]*X4 ~ X4, data=Cement); abline(a=0,b=coefMod["X4"])

The plot looks different from the ones produced by the crPlots above.

这个情节看起来和上面的情节不同。

The partial residuals have similar patterns, but their values and modeled relationships are different. The difference in values appears to due to the fact that crPlots used centered partial residuals (see this answer for a discussion of partial residuals in R). This brings me to my third option.

部分残差具有相似的模式，但它们的值和建模关系不同。值的差异似乎是由于使用了中心部分残差的一个事实(在R中关于部分残差的讨论，请参见这个答案)，这使我想到了第三个选项。

Option 3: Homemade plots with centered partial residuals

选项3:用居中的部分残余物自制地块

# Get the centered partial residuals
pRes <- resid(hackModel, type='partial')
# X1
# Plot the partial residuals
plot(pRes[,"X1"] ~ Cement$X1)
# Plot the component - modeled relationship
lines(coefMod["X1"]*(X1-mean(X1))~X1, data=Cement)
# X2 - X4
plot(pRes[,"X2"] ~ Cement$X2); lines(coefMod["X2"]*(X2-mean(X2))~X2, data=Cement) 
plot(pRes[,"X3"] ~ Cement$X3); lines(coefMod["X3"]*(X3-mean(X3))~X3, data=Cement)
plot(pRes[,"X4"] ~ Cement$X4); lines(coefMod["X4"]*(X4-mean(X4))~X4, data=Cement)

Now we have similar values than the crPlots above, but the relationships are still different. The difference may be related to the intercepts. But I'm not sure what I should use instead of 0.

现在我们有了与上面crplot相似的值，但是关系仍然是不同的。差异可能与拦截有关。但是我不确定应该用什么来代替0。

Any suggestions of which method is more appropriate? Is there a more straightforward way to get the partial residual plots based on model averaged coefficients?

对哪种方法更合适有什么建议吗?是否有一种更直接的方法来得到基于模型平均系数的局部残差图?

Many thanks!

很多谢谢!

1 个解决方案

#1

From looking at the crPlot.lm source code, it looks like only the functions residuals(model, type="partial"), predict(model, type="terms", term=var), and functions associated with finding the names of the variables are used on the model object. It also looks like the relationship is regressed, as @BenBolker suggested. The code used in crPlot.lm is: abline(lm(partial.res[,var]~.x), lty=2, lwd=lwd, col=col.lines[1]). Thus, I think that changing the coefficients and residuals of the model is sufficient to be able to use crPlots on it. I can now also reproduce the results in a homemade way.

从看电影的情节。lm源代码，看起来只有函数残差(model, type="partial")、预测值(model, type="terms"， term=var)以及与查找变量名相关的函数在模型对象上使用。正如@BenBolker所言，这段关系似乎也在倒退。crPlot中使用的代码。lm是:abline(lm)(部分。res[，var]~.x)， lty=2, lwd=lwd, col=col.lines[1])。因此，我认为改变模型的系数和残差是足够的，可以使用它的cr地块。我现在还可以用自制的方式复制结果。

library(MuMIn)
# Loading the data
data(Cement)
# Creating a full model with all the covariates we are interested in
fullModel <- lm(y ~ ., data = Cement, na.action=na.fail)
# Getting all possible models based on the covariates of the full model
muModel <- dredge(fullModel)
# Averaging across all models
avgModel <- model.avg(muModel)
# Getting the averaged coefficients
coefMod <- coef(avgModel)

# Option 1 - crPlots
library(car) # For crPlots
# Creating a duplicate of the fullMode
hackModel <- fullModel
# Changing the coefficents to the averaged coefficient
hackModel$coefficients <- coefMod[names(coef(fullModel))]
# Changing the residuals
hackModel$residuals <- Cement$y - predict(hackModel)

# Plot the crPlots and the regressed homemade version 
layout(matrix(1:8, nrow=2, byrow=TRUE))
par(mar=c(3.5,3.5,0.5,0.5), mgp=c(2,1,0))
crPlots(hackModel, layout=NA, ylab="Partial Res", smooth=FALSE)

# Option 4 - Homemade centered and regressed
# Get the centered partial residuals
pRes <- resid(hackModel, type='partial')
# X1 - X4 plot partial residuals and used lm for the relationship
plot(pRes[,"X1"] ~ Cement$X1); abline(lm(pRes[,"X1"]~Cement$X1))
plot(pRes[,"X2"] ~ Cement$X2); abline(lm(pRes[,"X2"]~Cement$X2))
plot(pRes[,"X3"] ~ Cement$X3); abline(lm(pRes[,"X3"]~Cement$X3))
plot(pRes[,"X4"] ~ Cement$X4); abline(lm(pRes[,"X4"]~Cement$X4))