首页 文章

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

提问于
浏览 次
3

我正在使用R包 MuMIn 进行多模型推理，并使用函数 model.avg 来平均由一组模型估计的系数 . 为了在视觉上将数据与基于平均系数的估计关系进行比较，我想使用部分残差图，类似于 car 函数的 crPlots 函数创建的图 . 我've tried three ways and I'我不确定是否合适 . 这是一个演示 .

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

选项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)

请注意，具有平均系数的完整和黑客版本的crPlots是不同的 .
crPlots Example

这里的问题是：这是否合适？结果依赖于我在_1349369中发现的黑客行为 . 除了残差和系数之外，我还需要更改模型的部分吗？

选项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"])

该图与上面的 crPlots 产生的图不同 .
home made example

部分残差具有相似的模式，但它们的值和建模关系是不同的 . 值的差异似乎是由于crPlots使用了居中的部分残差这一事实（有关R中部分残差的讨论，请参阅此answer） . 这让我想到了第三个选择 .

选项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)

Home made example with centered partial residuals

现在我们的值与上面的 crPlots 相似，但关系仍然不同 . 差异可能与拦截有关 . 但我不确定应该使用什么而不是0 .

有关哪种方法更合适的建议？是否有更直接的方法来获得基于模型平均系数的部分残差图？

非常感谢！

r plot lm r-car mumin

1 回答

  • 3

    通过查看 crPlot.lm 源代码，看起来只有函数 residuals(model, type="partial") ， predict(model, type="terms", term=var) 和与查找变量名称相关联的函数才会在模型对象上使用 . 正如@BenBolker建议的那样，这种关系看起来似乎已经退化了 . crPlot.lm 中使用的代码是： abline(lm(partial.res[,var]~.x), lty=2, lwd=lwd, col=col.lines[1]) . 因此，我认为改变模型的系数和残差足以能够在其上使用 crPlots . 我现在也可以用自制的方式重现结果 .

    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))

    comparison of crPlots and regressed

    回复于

相关问题