RパッケージMuMInを使用してマルチモデル推論を行い、関数を使用model.avgして一連のモデルによって推定された係数を平均化しています。平均係数に基づいて推定された関係とデータを視覚的に比較するには、パッケージのcrPlots関数によって作成されたものと同様の部分残差プロットを使用したいと考えています。car3つの方法を試しましたが、どれが適切かわかりません。これがデモンストレーションです。
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 に注意してください。
ここでの質問は次のとおりです。これは適切ですか? 結果は、この回答で見つけたハックに依存しています。残差と係数以外のモデルの部分を変更する必要がありますか?
オプション 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。
部分残差のパターンは似ていますが、それらの値とモデル化された関係は異なります。値の違いは、crPlots が中央の部分残差を使用したために表示されます (R の部分残差の説明については、この回答を参照してください)。これにより、3 番目のオプションにたどり着きます。
オプション 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)
上記と同様の値になりましたcrPlotsが、関係はまだ異なります。違いは切片に関連している可能性があります。しかし、0 の代わりに何を使用すればよいかわかりません。
どの方法がより適切であるかの提案はありますか? モデル平均係数に基づいて部分残差プロットを取得するより簡単な方法はありますか?
どうもありがとう！