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 つの交互作用 (教育 ~ 所得 + 若者 + 都市 + 西部地域:若者) を持つモデルは、BIC の点で最良のようです。

coef(MOD1.subset,4)

問題は、数式を手動で記述せずに、そのモデルから ML オブジェクトを取得するにはどうすればよいかということです。

summaryHH投稿する前に、やなどの regsubsets オブジェクトに対していくつかの興味深い機能を備えたパッケージ HH を見つけました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. 、で名前を分割し:ます:
  2. それぞれの側で、データフレーム内の変数と部分的に一致するかどうかを確認しますdf
  3. df一致する場合は、必要に応じて係数に強制変換した後、一致しないビットが変数のレベルと一致することを確認します。
  4. レベルと一致する場合は、インタラクションのその側を変数名に置き換えます。
  5. 一致するものがない場合は、他に部分的に一致するものがあるかどうかを確認し、一致しない場合は失敗します

等々。それは[so]の投稿に関係するプログラミングのかなりの部分ですが、あなたが挑戦に挑戦しているなら、上記はあなたに何かを始めるための何かを与えるはずです。

于 2012-10-25T08:48:18.760 に答える
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 の場合、モデルに対して何もしません。分散は説明されていません。必要に応じて、「best」変数の最後のインデックスを for ループの「i」に置き換えることで、これをループできます。

警告: 列は、トレーニング/クロス検証/テスト データに合わせて整列する必要があります。トレーニング セットの最初のインデックスが性別で、交差検証セットの最初のインデックスが年齢である場合、間違った列をゼロにすることになります。

補足: 調整された R^2 値に関して最適なモデルを選択したことがわかります。自由に変更してください。

私が助けてくれることを願っています。乾杯。

于 2016-11-22T04:05:05.953 に答える
0

この質問を元に戻して申し訳ありませんが、私はこれに対する答えを自分で探していました。

使用される特定の基準 (AIC、BIC など) は、regsubsets の結果には影響しません。これは、関数が同じサイズのモデルに対してのみ比較され、AIC が BIC と異なるのはモデル サイズに割り当てられた「ペナルティ」のみであるためです。ただし、異なるサイズのモデルを比較することに関心がある場合は、BIC ではなく AIC を使用することをお勧めします。

regsubsets に AIC をプロットする機能があるとは思いません。ただし、AIC は次を使用して簡単に計算できます。

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

ここで、n はサンプル数、p はモデルのパラメーター数です ( aic と bic の定義については、 http: //stat.wharton.upenn.edu/~khyuns/stat431/ModelSelection.pdf の最後のページを参照してください) ）。

regsubsets をだまして新しい aic 値をプロットしようとしましたが、失敗しました。ただし、aic 値の行列を見て、おそらく 'order(aic)' を使用して順序付けするだけで、異なるサイズのモデルを簡単に比較できます。

于 2015-04-09T20:40:16.527 に答える