6

線形混合効果モデルは、伝統的に次のように定式化されます。Ri = Xi × β + Zi × bi + εi ここで、β は推定された固定効果を表し、Z は変量効果を表します。したがって、X は古典的な計画行列です。R を使用して、nlme パッケージから lme を使用してモデルを適合させた後、これら 2 つの行列を抽出できるようにしたいと考えています。たとえば、nlme パッケージにもあるデータセット「Rails」には、ランダムに選択された 6 つの鉄道レールでの超音波移動時間の 3 つの個別の測定値が含まれています。次のように、レールごとにインターセプト固定効果とランダム効果を備えた単純なモデルを適合させることができます。

library(nlme)
lmemodel<-lme(travel ~ 1, random = ~ 1 | Rail, data=Rail)

X 設計マトリックスは、単一の 18x1 マトリックス (6 レール * 3 測定値) であり、次の方法で簡単に抽出できます。

 model.matrix(lmemodel, data=Rail)
   (Intercept)
1            1
2            1
3            1
4            1
5            1
6            1
7            1
8            1
9            1
10           1
11           1
12           1
13           1
14           1
15           1
16           1
17           1
18           1
attr(,"assign")
[1] 0

私がやりたいのは、変量効果の設計行列 Z を抽出することです。lme4 パッケージを使用して同じモデルに適合する場合、これは次の方法で実行できます。

library(lme4)
lmermodel<-lmer(travel ~ 1 + (1|Rail),data=Rail) 
t(lmermodel@Zt)  ##takes the transpose of lmermodel@Zt
lmermodel@X  ## extracts the X matrix

ただし、このマトリックスを lme 適合モデルから抽出する方法については途方に暮れています。

4

2 に答える 2

2

私が見る限り、Zマトリックスはlmeオブジェクトのどこにも保存されていません。最善の策はmodelStruct$reStructコンポーネントにあります（names(modelfit); str(modelfit); sapply(modelfit,class)探索するなどを試してください）が、私が知る限りそこにはありません。実際、行列が実際に明示的に構築されることは決してない可能性があることをlme.default示唆しています。Z内部的lmeには、代わりにグループ化構造で動作するようです。もちろんできます

Z <- model.matrix(~Rail-1,data=Rail)

しかし、それはおそらくあなたが考えていたものではありません...

于 2013-10-31T16:32:42.637 に答える