nlmeとlme4で異なるランダム効果を指定する方法は？

Question

nlme::lme（下部のデータ）を使用して、モデルに異なるランダム効果を指定したいとします。ランダム効果は、1）interceptおよびpositionはsubjectで変化します。 2）interceptはcomparisonで変化します。これはlme4::lmerを使用して簡単です。しかし、私はまた、（positionが時間変数である）自己相関構造をモデル化するようlmeに固執したい

lmer(rating ~ 1 + position + 
    (1 + position | subject) + 
    (1 | comparison), data=d) 

> ... 
Random effects: 
Groups  Name  Std.Dev. Corr 
comparison (Intercept) 0.31877  
subject (Intercept) 0.63289  
      position 0.06254 -1.00 
Residual    0.91458  
... 

。 lmeを使って、上記と同じようにするにはどうすればよいですか？私は以下のようにエフェクトをネストします。これは私が望むものではありません。

lme(rating ~ 1 + position, 
random = list(~ 1 + position | subject, 
       ~ 1 | comparison), data=d) 

> ... 
Random effects: 
Formula: ~1 + position | subject 
Structure: General positive-definite, Log-Cholesky parametrization 
      StdDev  Corr 
(Intercept) 0.53817955 (Intr) 
position 0.04847635 -1  

Formula: ~1 | comparison %in% subject # NESTED :(
     (Intercept)  Residual 
StdDev: 0.9707665 0.0002465237 
... 

注：;

）が存在SOおよびCV here、 here上のいくつかの同様の質問があり、 hereが、私のいずれかの答えや提案はここにカウントされませいる lmerを使用していた理解していません私がしようとする意味してきた例で使用

データ

d <- structure(list(rating = c(2, 3, 4, 3, 2, 4, 4, 3, 2, 1, 3, 2, 
2, 2, 4, 2, 4, 3, 2, 2, 3, 5, 3, 4, 4, 4, 3, 2, 3, 5, 4, 5, 2, 
3, 4, 2, 4, 4, 1, 2, 4, 5, 4, 2, 3, 4, 3, 2, 2, 2, 4, 5, 4, 4, 
5, 2, 3, 4, 3, 2), subject = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5", 
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", 
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
"28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", 
"39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", 
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", 
"61", "62", "63"), class = "factor"), position = c(1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), comparison = structure(c(1L, 
7L, 9L, 8L, 3L, 4L, 10L, 2L, 5L, 6L, 2L, 6L, 4L, 5L, 8L, 10L, 
7L, 3L, 1L, 9L, 3L, 9L, 10L, 1L, 5L, 7L, 6L, 8L, 2L, 4L, 4L, 
2L, 8L, 6L, 7L, 5L, 1L, 10L, 9L, 3L, 5L, 10L, 6L, 3L, 2L, 9L, 
4L, 1L, 8L, 7L, 6L, 5L, 2L, 10L, 4L, 3L, 8L, 9L, 7L, 1L), contrasts = structure(c(1, 
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 
0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 
0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 
0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0, 
0, 0, 0, 0, 0, 0, 0, 1, -1), .Dim = c(10L, 9L), .Dimnames = list(
    c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), NULL)), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10"), class = "factor")), .Names = c("rating", 
"subject", "position", "comparison"), row.names = c(1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 111L, 112L, 113L, 114L, 115L, 116L, 
117L, 118L, 119L, 120L, 221L, 222L, 223L, 224L, 225L, 226L, 227L, 
228L, 229L, 230L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 338L, 
339L, 340L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 448L, 449L, 
450L, 551L, 552L, 553L, 554L, 555L, 556L, 557L, 558L, 559L, 560L 
), class = "data.frame") 

Answer 1

これをしばらく把握する。はるかに多くの作業がなければ、私はlme4とまったく同じモデルを得ることができないと思っていますが、私は近づくことができます。

## source("SO36643713.dat") 
library(nlme) 
library(lme4) 

これはcomparisonためsubjectのための完全なランダム斜面項（相関傾きと切片）とあなたが望んでいたモデル、およびランダム切片である：

m1 <- lmer(rating ~ 1 + position + 
       (1 + position | subject) + 
       (1 | comparison), data=d) 

これは私がことができるものですlmeで複製する方法を理解する：独立したインターセプトとスロープ（私は特に、これらのモデルを好きではないが、彼らは人々があまりにも複雑でランダム効果モデルを簡素化するための方法としてかなり一般的に使用されている。）

m2 <- lmer(rating ~ 1 + position + 
       (1 + position || subject) + 
       (1 | comparison), data=d) 

結果：

VarCorr(m2) 
## Groups  Name  Std.Dev. 
## comparison (Intercept) 0.28115 
## subject position 0.00000 
## subject.1 (Intercept) 0.28015 
## Residual    0.93905 

のためにこの特定のデータセットでは、ランダムスロープは、とにかくゼロ分散を有すると推定される。

今度はlmeに設定しましょう。キー（???）の洞察は、pdBlocked()マトリックス内のすべての用語は同じグループ化変数の中にネストされているでなければならないということです。例えば、Pinheiro and Batesのpp。163ffの交差ランダム効果の例では、ランダム効果としてブロック、ブロック内の行、ブロック内の列があります。今、私たちが合うことができる

d$dummy <- factor(1) 

：comparisonとsubjectが両方にネストされている範囲内で、何のグループ化の要因が存在しないので、私はちょうど単一のブロックでデータセット全体が含まdummy「ファクタ」を作るつもりですモデル。

m3 <- lme(rating~1+position, 
      random=list(dummy = 
       pdBlocked(list(pdIdent(~subject-1), 
           pdIdent(~position:subject), 
           pdIdent(~comparison-1)))), 
      data=d) 

我々はランダム効果の分散共分散行列における三つのブロックがあります。subjectに1つ、position行列subject相互作用のための1、およびcomparisonための1つを。新しいpdMatクラスを定義するのに手間がかかりましたが、各勾配（position:subjectXX）を対応する切片（subjectXX）に関連付ける簡単な方法を理解できませんでした。 （pdBlocked構造でこれを設定できると思うかもしれませんが、pdBlockedオブジェクト内の複数のブロックにわたって分散の推定値を同じにするような方法はありません）。

結果はほぼ同じです、彼らは異なって報告されています。

vv <- VarCorr(m3) 
vv2 <- vv[c("subject1","position:subject1","comparison1","Residual"),] 
storage.mode(vv2) <- "numeric" 
print(vv2,digits=4) 
        Variance StdDev 
subject1   7.849e-02 2.802e-01 
position:subject1 4.681e-11 6.842e-06 
comparison1  7.905e-02 2.812e-01 
Residual   8.818e-01 9.390e-01