6

より複雑な実験計画の混合モデルに関する質問や投稿がいくつかあるため、このより単純なモデルは、このプロセスの他の初心者だけでなく、私にも役立つと思いました.

だから、私の質問は、sas proc 混合手順から R で反復測定 ancova を定式化したいということです。

proc mixed data=df1;
FitStatistics=akaike
class GROUP person day;
model Y = GROUP X1 / solution alpha=.1 cl;
repeated / type=cs subject=person group=GROUP;
lsmeans GROUP;
run;

以下は、R で作成されたデータを使用した SAS 出力です (以下)。

.           Effect       panel    Estimate       Error      DF    t Value    Pr > |t|     Alpha       Lower       Upper
            Intercept              -9.8693      251.04       7      -0.04      0.9697       0.1     -485.49      465.75
            panel        1         -247.17      112.86       7      -2.19      0.0647       0.1     -460.99    -33.3510
            panel        2               0           .       .        .         .             .           .           .
            X1                     20.4125     10.0228       7       2.04      0.0811       0.1      1.4235     39.4016

以下は、「nlme」パッケージを使用して R でモデルを定式化した方法ですが、同様の係数推定値が得られていません。

## create reproducible example fake panel data set:
set.seed(94); subject.id = abs(round(rnorm(10)*10000,0))

set.seed(99); sds = rnorm(10,15,5);means = 1:10*runif(10,7,13);trends = runif(10,0.5,2.5)

this = NULL; set.seed(98)
for(i in 1:10) { this = c(this,rnorm(6, mean = means[i], sd = sds[i])*trends[i]*1:6)}
set.seed(97)
that = sort(rep(rnorm(10,mean = 20, sd = 3),6))

df1 = data.frame(day = rep(1:6,10), GROUP = c(rep('TEST',30),rep('CONTROL',30)),
                 Y = this,
                 X1 = that,
                 person = sort(rep(subject.id,6)))

## use package nlme
require(nlme)

## run repeated measures mixed model using compound symmetry covariance structure:
summary(lme(Y ~ GROUP + X1, random = ~ +1 | person,
            correlation=corCompSymm(form=~day|person), na.action = na.exclude,
            data = df1,method='REML'))

さて、私が今気づいたRからの出力は、からの出力に似ていますlm()：

                Value Std.Error DF    t-value p-value
(Intercept) -626.1622  527.9890 50 -1.1859379  0.2413
GROUPTEST   -101.3647  156.2940  7 -0.6485518  0.5373
X1            47.0919   22.6698  7  2.0772934  0.0764

私は仕様に近いと信じていますが、結果を一致させるためにどの部分が欠けているのかわかりません (理由の範囲内で..)。どんな助けでも大歓迎です！

更新:以下の回答のコードを使用すると、R 出力は次のようになります。

> summary(model2)

一番下までスクロールしてパラメーター推定値を確認してください。SASと同じ。

Linear mixed-effects model fit by REML
 Data: df1 
      AIC      BIC   logLik
  776.942 793.2864 -380.471

Random effects:
 Formula: ~GROUP - 1 | person
 Structure: Diagonal
        GROUPCONTROL GROUPTEST Residual
StdDev:      184.692  14.56864 93.28885

Correlation Structure: Compound symmetry
 Formula: ~day | person 
 Parameter estimate(s):
         Rho 
-0.009929987 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | GROUP 
 Parameter estimates:
    TEST  CONTROL 
1.000000 3.068837

Fixed effects: Y ~ GROUP + X1 

                Value Std.Error DF    t-value p-value
(Intercept)   -9.8706 251.04678 50 -0.0393178  0.9688
GROUPTEST   -247.1712 112.85945  7 -2.1900795  0.0647
X1            20.4126  10.02292  7  2.0365914  0.0811
4

2 に答える 2

5

以下を試してください：

model1 <- lme(
  Y ~ GROUP + X1,
  random = ~ GROUP | person,
  correlation = corCompSymm(form = ~ day | person),
  na.action = na.exclude, data = df1, method = "REML"
)
summary(model1)

この場合、オプション with はオプション with と同様の結果をもたらすと思いrandom = ~ groupvar | subjvarます。R lmerepeated / subject = subjvar group = groupvarSAS/MIXED

編集：

SAS/混合

SAS/MIXED 共分散行列

R (修正モデル2)

model2 <- lme(
  Y ~ GROUP + X1,
  random = list(person = pdDiag(form = ~ GROUP - 1)),
  correlation = corCompSymm(form = ~ day | person),
  weights = varIdent(form = ~ 1 | GROUP),
  na.action = na.exclude, data = df1, method = "REML"
)
summary(model2)

R 共分散行列

したがって、これらの共分散構造は非常に似ていると思います (σ g1 = τ g 2 + σ 1 )。

編集2：

共変量推定 (SAS/MIXED):

Variance            person          GROUP TEST        8789.23
CS                  person          GROUP TEST         125.79
Variance            person          GROUP CONTROL       82775
CS                  person          GROUP CONTROL       33297

そう

TEST group diagonal element
  = 125.79 + 8789.23
  = 8915.02
CONTROL group diagonal element
  = 33297 + 82775
  = 116072

ここで、対角要素 = σ k1 + σ k 2 .

共変量推定 (R lme):

Random effects:
 Formula: ~GROUP - 1 | person
 Structure: Diagonal
        GROUP1TEST GROUP2CONTROL Residual
StdDev:   14.56864       184.692 93.28885

Correlation Structure: Compound symmetry
 Formula: ~day | person 
 Parameter estimate(s):
         Rho 
-0.009929987 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | GROUP 
 Parameter estimates:
   1TEST 2CONTROL 
1.000000 3.068837

そう

TEST group diagonal element
  = 14.56864^2 + (3.068837^0.5 * 93.28885 * -0.009929987) + 93.28885^2
  = 8913.432
CONTROL group diagonal element
  = 184.692^2  + (3.068837^0.5 * 93.28885 * -0.009929987) + (3.068837 * 93.28885)^2
  = 116070.5

ここで、対角要素 = τ g 2 + σ 1 + σ g 2 .

于 2012-08-06T00:33:56.477 に答える
4

ああ、これはトリッキーなものになるでしょう、そしてそれが標準のnlme関数を使ってさえ可能であるなら、Pinheiro/Batesのいくつかの真剣な研究をするつもりです。

ただし、それを行う前に、これが必要なモデルであることを絶対に確認する必要があります。おそらく、データのストーリーによりよく適合する可能性のある何かが他にあります。あるいは、Rがもっと簡単にできることで、同じくらい良いのですが、まったく同じではないかもしれません。

まず、SASでこの行を使用して行っていることについての私の見解を次に示します。

repeated / type=cs subject=person group=GROUP;

これtype=cs subject=personは、同じ人物のすべての測定値の間に相関関係を引き起こし、その相関関係はすべての日のペアで同じです。これgroup=GROUPにより、各グループの相関を異ならせることができます。

対照的に、これがあなたのRコードが何をしているのかについての私の見解です：

random = ~ +1 | person,
correlation=corCompSymm(form=~day|person)

このコードは、実際には2つの異なる方法でほぼ同じ効果を追加しています。random線は各人にランダム効果を追加し、線correlationは同じ人のすべての測定値間の相関を誘導しています。ただし、これら2つのことはほとんど同じです。相関が正の場合、どちらかを含めることでまったく同じ結果が得られます。両方を含めるとどうなるかわかりませんが、必要なのは1つだけです。とにかく、このコードはすべての個人に対して同じ相関関係を持っており、各グループが独自の相関関係を持つことを許可していません。

各グループに独自の相関関係を持たせるには、2つの異なる部分からより複雑な相関関係構造を構築する必要があると思います。私はこれをやったことはありませんが、Pinheiro/Batesがやったことを覚えていると確信しています。

代わりに、人にランダム効果を追加してから、グループごとに分散を異なるものにすることを検討してweights=varIdent(form=~1|group)ください（メモリから、構文を確認してください）。これはまったく同じではありませんが、同様の話をします。SASの話は、一部の個人の測定値が他の個人の測定値よりも相関しているということです。それが何を意味するのかを考えると、相関の高い個人の測定値は、相関の低い個人の測定値よりも近くなります。対照的に、Rの話は、個人内の測定値のばらつきが異なるということです。そのことを考えると、変動性が高く、相関が低い測定値です。ですから、彼らは似たような話をしますが、反対側からやって来ます。

これらの2つのモデルが同じものの異なるパラメーター化になる可能性さえあります（しかし私は驚かれることでしょう）。私の直感では、全体的な測定のばらつきは何らかの形で異なるでしょう。ただし、それらが同じものでなくても、パラメーター化を書き留めて、それらを確実に理解し、データのストーリーを適切に記述していることを確認することは価値があります。

于 2012-08-06T01:34:23.317 に答える