2

線形モデルのリストで AIC を使用してステップワイズ回帰を実行したいと考えています。アイデアは、線形モデルのリストを使用してから、各リスト要素に stepAIC を適用することです。失敗します。

問題を突き止めてみました。私は問題を見つけたと思います。ただし、原因がわかりません。コードを試して、3 つのケースの違いを確認してください。

require(MASS)
n<-30 
x1<-rnorm(n, mean=0, sd=1) #create rv x1 
x2<-rnorm(n, mean=1, sd=1)
x3<-rnorm(n, mean=2, sd=1)
epsilon<-rnorm(n,mean=0,sd=1) # random error variable 
dat<-as.data.frame(cbind(x1,x2,x3,epsilon)) # combine to a data frame
dat$id<-c(rep(1,10),rep(2,10),rep(3,10)) 
# y is combination from all three x and a random uniform variable
dat$y<-x1+x2+x3+epsilon 
# apply lm() only resulting in a list of models
dat.lin.model.lst<-lapply(split(dat,dat$id),function(d) lm(y~x1+x2+x3,data=d)) 
stepAIC(dat.lin.model.lst[[1]]) # FAIL!!!
# apply function stepAIC(lm())-  works
dat.lin.model.stepAIC.lst<-lapply(split(dat,dat$id),function(d) stepAIC(lm(y~x1+x2+x3,data=d))) 
# create model for particular group with id==1
k<-which(dat$id==1) # manually select records with id==1
lin.model.id1<-lm(dat$y[k]~dat$x1[k]+dat$x2[k]+dat$x3[k]) 
stepAIC(lin.model.id1) # check stepAIC - works!

stepAIC() には data.frame "dat" からの元のデータが必要であると確信しています。それは私が前に考えていたことです。(私が正しいことを願っています)しかし、元のデータフレームを渡すことができる stepAIC() にはパラメーターがありません。明らかに、リストにラップされていない単純なモデルの場合、モデルを渡すだけで十分です。(コードの最後の3行)だから私は疑問に思っています:

  • Q1: stepAIC は元のデータ "dat" (パラメーターとして渡されたモデル データだけでなく) を見つける方法をどのように知っていますか?
  • Q2: stepAIC() に、ヘルプ ページに明示的に記載されていない別のパラメーターがあることをどのように知ることができますか? (私の英語が下手すぎて見つけられないだけかもしれません)
  • Q3: そのパラメーターを stepAIC() に渡すにはどうすればよいですか?

これは、適用関数の環境のどこかにあり、データを渡す必要があります。lm() または stepAIC() のいずれか、および生データへのポインター/リンクがどこかで失われる必要があります。R の環境が何をするのかよくわかりません。私にとっては、ローカル変数をグローバル変数から分離するようなものでした。しかし、おそらくもっと複雑です。上記の問題に関して私にそれを説明できる人はいますか? 正直なところ、私はR のドキュメントをあまり読んでいません。より良い理解は私を助けるでしょう。

OLD: 複数のサブグループに分割できるデータフレーム df にデータがあります。そのために、df$id というグループ ID を作成しました。lm() は、最初のサブグループに期待される係数を返します。各サブグループの基準として AIC を個別に使用して段階的回帰を実行したいと考えています。各サブグループ (id) のモデルになる lmList {lme4} を使用します。しかし、リスト要素に stepAIC{MASS} を使用すると、エラーがスローされます。下記参照。

質問は次のとおりです。手順/構文にどのような間違いがありますか? 単一のモデルの結果が得られますが、lmList で作成されたモデルの結果は得られません。lmList() は lm() とは異なるモデルに関する情報を保存しますか?
ただし、ヘルプには次のように記載されています。 class "lmList": 共通モデルを持つクラス lm のオブジェクトのリスト。

>lme4.list.lm<-lmList(formula=Scherkraft.N~Gap.um+Standoff.um+Voidflaeche.px |df$id,data = df)
>lme4.list.lm[[1]]
Call: lm(formula = formula, data = data)
Coefficients:
(Intercept)          Gap.um     Standoff.um  Voidflaeche.px  
  62.306133       -0.009878        0.026317       -0.015048  

>stepAIC(lme4.list.lm[[1]], direction="backward") 
#stepAIC on first element on the list of linear models
Start:  AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
                 Df Sum of Sq    RSS    AIC
- Standoff.um     1      2.81 7187.3 293.14
- Gap.um          1     29.55 7214.0 293.37
<none>                        7184.4 295.12
- Voidflaeche.px  1    604.38 7788.8 297.97  

Error in terms.formula(formula, data = data) : 
'data' argument is of the wrong type

明らかに、リストでは何かが機能しません。しかし、それが何であるかはわかりません。同じモデル(少なくとも同じ係数)を作成する基本パッケージで同じことをしようとしたので。結果は以下のとおりです。

>lin.model<-lm(Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px,df[which(df$id==1),]) 
# id is in order, so should be the same subgroup as for the first list element in lmList

Coefficients:  
(Intercept)    Gap.um  Standoff.um  Voidflaeche.px  
  62.306133 -0.009878     0.026317       -0.015048  

これは、私の linear.model で stepAIC を使用して返されたものです。私の知る限り、赤池情報量基準を使用して、特定のデータが与えられた場合に適合と一般化のバランスが取れているモデルを推定できます。

>stepAIC(lin.model,direction="backward")
Start:  AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
                 Df Sum of Sq    RSS    AIC
- Standoff.um     1      2.81 7187.3 293.14  
- Gap.um          1     29.55 7214.0 293.37
<none>                        7184.4 295.12
- Voidflaeche.px  1    604.38 7788.8 297.97  

Step:  AIC=293.14
Scherkraft.N ~ Gap.um + Voidflaeche.px
                 Df Sum of Sq    RSS    AIC
- Gap.um          1     28.51 7215.8 291.38
 <none>                        7187.3 293.14
- Voidflaeche.px  1    717.63 7904.9 296.85

Step:  AIC=291.38
Scherkraft.N ~ Voidflaeche.px
                 Df Sum of Sq    RSS    AIC
<none>                        7215.8 291.38
- Voidflaeche.px  1    795.46 8011.2 295.65
Call: lm(formula = Scherkraft.N ~ Voidflaeche.px, data = df[which(df$id == 1), ])

Coefficients:
(Intercept)  Voidflaeche.px  
   71.7183         -0.0151  

モデルを使用する必要がある出力から読み取りました: Scherkraft.N ~ Voidflaeche.pxこれは最小の AIC であるためです。まあ、誰かが出力を簡単に説明できればいいのですが。段階的回帰についての私の理解 (後方消去を仮定) は、すべての回帰変数が初期モデルに含まれているということです。次に、最も重要でないものが削除されます。決定する基準は AIC です。など...どういうわけか、テーブルを正しく解釈するのに問題があります。誰かが私の解釈を確認できるといいですね。「-」(マイナス) は、削除されたリグレッサーを表します。一番上に「開始」モデルがあり、下の表では、RSS と AIC が計算されて除去の可能性が示されています。したがって、最初の表の最初の行は、モデルScherkraft.N~Gap.um+Standoff.um+Voidflaeche.pxを示しています。 - Standoff.umは AIC 293.14 になります。Standoff.um のないものを選択してください: Scherkraft.N~Gap.um+Voidflaeche.px

編集:
lmList{lme4} を dlply() に置き換えて、モデルのリストを作成しました。それでもstepAICはリストに対応していません。別のエラーがスローされます。実際、stepAICが実行する必要があるデータに問題があると思います。モデルデータだけから各ステップのAIC値をどのように計算するのか疑問に思っていました. 元のデータを使用してモデルを構築し、毎回 1 つのリグレッサーを除外します。その上で、AIC を計算して比較します。元のデータにアクセスできない場合、stepAIC はどのように機能しますか。(元のデータを stepAIC に渡すパラメーターが表示されません)。それでも、単純なモデルでは機能するのに、リストにラップされたモデルでは機能しない理由がわかりません。

>model.list.all <- dlply(df, .id, function(x) 
  {return(lm(Scherkraft.N~Gap.um+Standoff.um+Voidflaeche.px,data=x)) })
>stepAIC(model.list.all[[1]])
Start:  AIC=295.12
Scherkraft.N ~ Gap.um + Standoff.um + Voidflaeche.px
                 Df Sum of Sq    RSS    AIC
- Standoff.um     1      2.81 7187.3 293.14
- Gap.um          1     29.55 7214.0 293.37
<none>                        7184.4 295.12
- Voidflaeche.px  1    604.38 7788.8 297.97
Error in is.data.frame(data) : object 'x' not found
4

2 に答える 2

4

デバッグを非常に困難にするためにバージョン管理で何が変更されたのかはわかりませんが、1 つの解決策はdo.call、呼び出しの式を実行前に評価する を使用することです。これはd、呼び出しで保存するだけでなく、作業を行うためupdatestepAIC検索する必要があるdため、データ フレーム自体の完全な表現を保存することを意味します。

つまり、する

do.call("lm", list(y~x1+x2+x3, data=d))

それ以外の

lm(y~x1+x2+x3, data=d)

callおそらく次のように、モデルの要素を見ると、それが何をしようとしているのかがわかります。

dat.lin.model.lst <- lapply(split(dat, dat$id), function(d)
                            do.call("lm", list(y~x1+x2+x3, data=d)) )
dat.lin.model.lst[[1]]$call

また、グローバル環境でデータ フレームのリストを作成し、呼び出しを作成して各データ フレームを順番に探すこともできますupdatestepAICこれは、環境チェーンが常にグローバル環境に戻るためです。このような:

dats <- split(dat, dat$id)
dat.lin.model.list <- lapply(seq_along(dats), function(d)
            do.call("lm", list(y~x1+x2+x3, data=call("[[", quote(dats),i))) )

変更内容を確認するには、もう一度実行dat.lin.model.lst[[1]]$callします。

于 2012-04-23T14:32:33.217 に答える