如何从regsubsets获取LM对象

Question

让我们想象一下,我们想用收入,年轻人,城市和地区作为回归者来模拟美国州立公立学校的支出(教育).欲了解更多信息:?Anscombe 模型:教育〜(收入+年轻人+城市)*地区

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) 

具有3个变量和1个交互(教育〜收入+年轻+城市+ RegionWest:年轻)的模型似乎是BIC方面的最佳模型.

coef(MOD1.subset,4)

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

在发布之前,我发现了包HH,它为regsubsets对象提供了一些有趣的功能,例如summaryHH和lm.regsubsets.

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

lm.regsubsets做我想要的但我认为解析交互有一些问题,任何替代方案？

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. :用on分割任何名称:，
2. 对于每一侧，查看是否与数据框中的变量部分匹配df
3. df如果存在匹配，则检查未匹配的位是否与强制转换为因子后的变量级别（如果需要）匹配，
4. 如果与级别匹配，则用变量名称替换交互的那一侧，
5. 如果没有匹配，则查看是否还有其他部分匹配，如果没有则失败

等等。一个帖子涉及相当多的编程，但如果您能够应对挑战，那么上面的内容应该会给您一些开始的机会。