2017-07-02 253 views
0

我试图在模式混合模型上运行一个模拟,并且需要R(在非结构化下)的“估计的渐近协方差矩阵或估计协方差参数的协方差矩阵”。
我知道这将通过SAS中的AsyCov和SPSS中的混合模型来实现。
但我不知道为什么asyCov(metaSEM包)的结果与SAS和SPSS输出不一致。提取R中估计协方差参数的协方差矩阵?

这里是我的SAS代码:

proc Mixed data=OutcomeSort method=reml asycov covtest; 
    class Subject Rep; 
    model Y= x/s covb; 
    repeated Rep/subject=Subject type=UN r; 
run; 

和我的SAS输出:

       Covariance Parameter Estimates 

               Standard   Z 
      Cov Parm Subject Estimate  Error  Value  Pr Z 

      UN(1,1)  Subject  50.6700  5.2325  9.68  <.0001 
      UN(2,1)  Subject  38.9197  6.1316  6.35  <.0001 
      UN(2,2)  Subject  109.57  11.3113  9.69  <.0001 
      UN(3,1)  Subject  37.8759  7.1731  5.28  <.0001 
      UN(3,2)  Subject  83.7478  11.5030  7.28  <.0001 
      UN(3,3)  Subject  162.44  16.7682  9.69  <.0001 
      UN(4,1)  Subject  32.3689  8.6577  3.74  0.0002 
      UN(4,2)  Subject  85.1421  13.6659  6.23  <.0001 
      UN(4,3)  Subject  149.65  18.3647  8.15  <.0001 
      UN(4,4)  Subject  247.16  25.7446  9.60  <.0001 


        Asymptotic Covariance Matrix of Estimates 

Row Cov Parm  CovP1  CovP2  CovP3  CovP4  CovP5  CovP6  CovP7  CovP8 

    1 UN(1,1) 27.3790 20.9481 15.9859 20.3577 15.5259 15.0775 17.0100 12.8662 
    2 UN(2,1) 20.9481 37.5971 45.4151 30.4336 39.4733 33.8181 29.8290 36.7129 
    3 UN(2,2) 15.9859 45.4151 127.95 34.7757 97.8981 74.9444 36.0516 100.22 
    4 UN(3,1) 20.3577 30.4336 34.7757 51.4533 50.6228 65.5963 47.2107 45.8363 
    5 UN(3,2) 15.5259 39.4733 97.8981 50.6228 132.32 145.12 49.1089 126.34 
    6 UN(3,3) 15.0775 33.8181 74.9444 65.5963 145.12 281.17 61.4499 134.74 
    7 UN(4,1) 17.0100 29.8290 36.0516 47.2107 49.1089 61.4499 74.9564 69.2153 
    8 UN(4,2) 12.8662 36.7129 100.22 45.8363 126.34 134.74 69.2153 186.76 
    9 UN(4,3) 12.4790 31.8059 76.8862 58.5769 141.49 260.08 79.1265 182.24 
10 UN(4,4) 10.2027 29.7120 78.8488 52.3062 137.66 240.73 91.0933 231.19 

           Row  CovP9  CovP10 

           1 12.4790  10.2027 
           2 31.8059  29.7120 
           3 76.8862  78.8488 
           4 58.5769  52.3062 
           5 141.49  137.66 
           6 260.08  240.73 
           7 79.1265  91.0933 
           8 182.24  231.19 
           9 337.26  401.18 
           10 401.18  662.78 

,在这里我的R代码里面:

Std.Err.Cov.par <- matrix (c(5.2325,6.1316,7.1731,8.6577,6.1316, 11.3113,11.5030,13.6659,7.1731,11.5030,16.7682,18.3647,8.6577,13.6659,18.3647,25.7446),4)) 
Std.Err.Cov.parasyCov 
     [,1] [,2] [,3] [,4] 
[1,] 5.2325 6.1316 7.1731 8.6577 
[2,] 6.1316 11.3113 11.5030 13.6659 
[3,] 7.1731 11.5030 16.7682 18.3647 
[4,] 8.6577 13.6659 18.3647 25.7446 

asyCov(Std.Err.Cov.par,n=100,cor.analysis = F) 

     x1x1  x2x1  x3x1  x4x1  x2x2  x3x2  x4x2 
x1x1 0.5366875 0.6289067 0.7357319 0.8880043 0.7369721 0.8621534 1.040591 
x2x1 0.6289067 0.9485741 1.0209965 1.2211370 1.3595297 1.4865150 1.781083 
x3x1 0.7357319 1.0209965 1.3642389 1.5504870 1.3825724 1.8164115 2.079731 
x4x1 0.8880043 1.2211370 1.5504870 2.0549314 1.6425356 2.0644151 2.706764 
x2x2 0.7369721 1.3595297 1.3825724 1.6425356 2.5079958 2.5505029 3.030071 
x3x2 0.8621534 1.4865150 1.8164115 2.0644151 2.5505029 3.1558301 3.576670 
x4x2 1.0405908 1.7810828 2.0797308 2.7067640 3.0300707 3.5766699 4.684520 
x3x3 1.0085977 1.6174143 2.3577429 2.5822236 2.5937308 3.7809433 4.140927 
x4x3 1.2173439 1.9368496 2.7139703 3.3682746 3.0814257 4.3163971 5.362250 
x4x4 1.4692936 2.3192276 3.1166576 4.3690889 3.6608210 4.9195376 6.896459 

    x3x3  x4x3  x4x4 
x1x1 1.008598 1.217344 1.469294 
x2x1 1.617414 1.936850 2.319228 
x3x1 2.357743 2.713970 3.116658 
x4x1 2.582224 3.368275 4.369089 
x2x2 2.593731 3.081426 3.660821 
x3x2 3.780943 4.316397 4.919538 
x4x2 4.140927 5.362250 6.896459 
x3x3 5.511570 6.036326 6.611043 
x4x3 6.036326 7.536536 9.267696 
x4x4 6.611043 9.267696 12.991931 
+1

什么功能(S),您使用的是混合模式?从'lme4'?或从'nlme'或其他东西? – ekstroem

+1

发布一些示例数据和使用的代码会产生不同的结果。 – Tom

+0

我使用了nlme包。 –

回答

0

我们真的需要一个可重复的例子,但基于此SAS代码:

proc Mixed data=OutcomeSort method=reml asycov covtest; 
    class Subject Rep; 
    model Y= x/s covb; 
    * options: covb displays fixed-effect var-cov matrix 
    *   s (solution) displays fixed-effect estimates 
    repeated Rep/subject=Subject type=UN r; 
    * repeated: "R-side" (residual) effects 
    * Rep governs order? of obs (I think this is 
    * irrelevant for unstructured var-cov matrices?) 
    * Subject is grouping variable 
    * unstructured var-cov matrix (this is R's default) 
    * r: display blocks of R matrix 
run; 

我建议相应的R代码是这样的:

library(nlme) 
Orthodont$fAge <- factor(Orthodont$age) 
Orthodont$nAge <- as.numeric(Orthodont$fAge) 
fit1 <- gls(distance~Sex, 
    correlation=corSymm(form=~nAge|Subject), 
    weights=varIdent(form=~1|fAge), 
    data=Orthodont) 

然而,这给沃尔德上的相关性参数方差(corStruct*),方差2-4的比率方差为1(varStruct*),和残差(lSigma);此外,这些是天平规模上的参数差异。

从皮涅罗和贝茨2000 p。 93 Google books

lme [和gls]所使用的方法是计算置信区间时,需要考虑不同的参数化。这种参数化,我们称之为自然参数化,使用标准偏差的对数和相关性的一般化logits。对于给定的相关性参数$ -1 < \ rho < 1 $,其广义logit为$ \ log [(1+ρ)/(1-ρ)] $,其可以取实线上的任何值

这些转换的逆转换是exp(lSigma)(exp(logitRho)-1)/(exp(logitRho)+1)

从这里,你必须使用增量的方法来将这些碎片和转换一起回来......

>

+0

谢谢Dr. Bolker,这种方式非常有帮助。我可以有一个问题,只是为了让自己清楚。那么,如何在R中提取“估计协方差参数的协方差矩阵”就像在SAS中一样? –