r - 如何从regsubsets中获取LM对象

假设我们要使用收入，年轻人，城市和地区作为回归变量来对美国州公立学校支出(教育)进行建模。有关更多信息:?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)

就BIC而言，具有3个变量和1个交互作用的模型(education〜income + young + urban + RegionWest:young)似乎是最好的。

coef(MOD1.subset,4)

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

在发布之前，我发现了HH软件包，它对regsubsets对象具有一些有趣的功能，例如summaryHH和lm.regsubsets。

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

lm.regsubsets可以满足我的要求，但是我认为解析交互存在一些问题，还有其他选择吗？

最佳答案

我认为这是不可能的，至少在没有大量处理系数名称的情况下也是如此。我到达那里的方式达到了〜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"

通过一些努力，您很可能可以做到:

在:和

上用:拆分任何名称
两侧的

，查看是否与数据框中的变量部分匹配df

(如果存在匹配项)，如果需要将其强制为一个因子

，请检查不匹配的位是否匹配df中变量的级别

(如果与某个级别匹配)，则用变量名

替换交互作用的这一面

如果不存在匹配项，请查看是否存在任何其他部分匹配项，如果不存在

则失败

等等。 [so]发布涉及大量的编程，但是如果您要应对挑战，那么以上内容应该为您提供一些起点。