8

假设我们想使用收入、年轻人、城市和地区作为回归变量来模拟美国公立学校的支出（教育）。欲了解更多信息：?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:young）的模型在 BIC 方面似乎是最好的。

coef(MOD1.subset,4)

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

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

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

lm.regsubsets做我想要的，但我认为解析交互有一些问题，还有其他选择吗？

4

3 回答 3

3

我认为这是不可能的，至少在没有对系数名称进行大量处理的情况下是不可能的。我在那里完成了约 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. 如果没有匹配，查看是否有其他部分匹配，如果没有则失败

等等。对于 [so] 发布来说，这涉及到很多编程，但是如果您能够应对挑战，那么上面的内容应该可以为您提供一些起点。

于 2012-10-25T08:48:18.760 回答
0

很抱歉把这个问题重新提出来，但我自己也在寻找这个问题的答案。

使用的特定标准（例如 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)”对它们进行排序

于 2015-04-09T20:40:16.527 回答
0

我想到了。所以这里是代码。

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”来循环。

警告：您的列必须为您的训练/交叉验证/测试数据排队。如果您的训练集的第一个索引是 Gender，而您的交叉验证集的第一个索引是 Age，那么您会将错误的列归零。

旁注：您可以看到我选择了关于调整后的 R^2 值的最佳模型。随意改变它。

希望我有所帮助。干杯。

于 2016-11-22T04:05:05.953 回答