如何从regsubsets中获取LM对象

Question

假设我们想使用收入、年轻人口、城市化程度和地区作为回归变量来模拟美国州公立学校的支出（教育）。有关更多信息：?Anscombe 模型：education ~ (income+young+urban)*region

library(car)
library(leaps)

#Loading Data
data(Anscombe)
data(state)
stateinfo <- data.frame(region=state.region,row.names=state.abb)
datamodel <- data.frame(merge(stateinfo,Anscombe,by="row.names"),row.names=1)
head(datamodel)
   region education income young urban
AK   West       372   4146 439.7   484
AL  South       112   2337 362.2   584
AR  South       134   2322 351.9   500
AZ   West       207   3027 387.5   796
CA   West       273   3968 348.4   909
CO   West       192   3340 358.1   785

#Saturated Model
MOD1 <- lm(education~(.-region)*region,datamodel)
summary(MOD1)
#anova(MOD1)

#Look for the "best" model
MOD1.subset <- regsubsets(education~(.-region)*region,datamodel,nvmax=15)
plot(MOD1.subset) 

在BIC方面，具有3个变量和1个交互作用的模型（教育程度~收入+年轻+城市+西部地区:年轻）似乎是最好的。

coef(MOD1.subset,4)

问题是，我如何在不手动编写公式的情况下从该模型中获取ML对象？

在发布之前，我发现了HH包，它有一些有趣的函数适用于regsubsets对象，例如summaryHH和lm.regsubsets。

library(HH)
summaryHH(MOD1.subset)[4,]
lm.regsubsets(MOD1.subset,4)

lm.regsubsets可以实现我想要的功能，但我认为在解析交互作用方面存在一些问题，是否有其他替代方法？

- Emer

---

当我思考“最佳模型”时，我真的感到很困惑。你究竟如何解释一个交互由二元变量和仅有一个分类变量水平组成的模型？我感到不舒服，因为我觉得你并没有完全理解这个软件包背后的推理。 - Joris Meys

Answer 1

我认为这可能不太可能实现，至少需要对系数名称进行大量处理。我已经完成了95%的工作，但在交互项上失败了：

coefs <- coef(MOD1s, 4)
nams <- names(coefs)
nams <- nams[!nams %in% "(Intercept)"]
response <- as.character(as.formula(MOD1s$call[[2]])[[2]])
form <- as.formula(paste(response, paste(nams, collapse = " + "), sep = " ~ "))
df <- get(as.character(MOD1s$call[[3]]))
mod <- lm(form, data = df)

> mod <- lm(form, data = df)
Error in eval(expr, envir, enclos) : object 'regionWest' not found

这是有道理的，并且源于系数使用的名称：

> nams
[1] "income"           "young"            "urban"           
[4] "regionWest:young"

很有可能，只要付出一些努力，您就能够做到以下几点：

1. 使用冒号 : 将任何名称拆分为 ,。
2. 对于每个部分，在数据框 df 中查看是否有与之部分匹配的变量。
3. 如果有匹配，则检查未匹配部分是否在强制转换为因子后与 df 中变量的水平相匹配。
4. 如果与水平匹配，则用变量名称替换交互的该侧。
5. 如果没有匹配，则查看是否有其他部分匹配，如果没有则失败。

等等。这对于一个 [so] 的帖子来说需要许多编程，但如果您有能力应对这个挑战，上述内容应该是一个起点。

- Gavin Simpson

Answer 2

我已经想通了，这是代码。

fit <- regsubsets(y~x,data=train)
b <- data[c(predictor columns)]
best <- order(summary(fit)$adjr2,decreasing=T)[1]
a <- as.integer(summary(fit)$which[best,][1:ncol(b)+1]
new_data <- data.frame(t( t(b) * a))
fit_lm <- lm(y~x,data=new_data)

你需要将列乘以布尔值。如果输入始终为0，则对模型不会产生任何影响，也不会解释任何方差。如果想要循环执行此操作，可以通过使用for循环中的“i”替换“best”变量中的最后一个索引来实现。

警告：您的列必须与您的训练/交叉验证/测试数据相对应。如果训练集的第一个索引是性别，而交叉验证集的第一个索引是年龄，则会将错误的列清零。

另外：您可以看到我选择了最佳模型，具体取决于调整后的R ^ 2值。如有需要，请随意更改。

希望能对您有所帮助。祝好！

- Jon Joh

Answer 3

抱歉又提出这个问题，但我自己也在寻找答案。

具体使用的标准（例如 AIC，BIC）不影响 regsubsets 的结果，因为该函数仅与相同大小的模型进行比较，并且 AIC 仅通过分配给模型大小的“惩罚”与 BIC 不同。但是，如果您有兴趣比较不同大小的模型，则可能更喜欢使用 AIC 而不是 BIC。

我认为 regsubsets 没有绘制 AIC 的能力。 但是，可以使用以下方法轻松计算 AIC：

aic <-summary(leaps2)$bic + (2 - log(n))*(p+1)

n是样本数量，p是模型中的参数数量（请参见http://stat.wharton.upenn.edu/~khyuns/stat431/ModelSelection.pdf的最后一页以获取aic和bic的定义）。

我尝试欺骗regsubsets来绘制新的aic值，但没有成功。然而，您可以通过查看aic值矩阵并使用“order(aic)”对其进行排序来轻松比较不同大小的模型。

- Steve