提取R中估計協方差參數的協方差矩陣？

Question

我試圖在模式混合模型上運行一個模擬，並且需要R（在非結構化下）的「估計的漸近協方差矩陣或估計協方差參數的協方差矩陣」。
我知道這將通過SAS中的AsyCov和SPSS中的混合模型來實現。
但我不知道爲什麼asyCov（metaSEM包）的結果與SAS和SPSS輸出不一致。

這裏是我的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 

Answer 1

我們真的需要一個可重複的例子，但基於此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)。

從這裏，你必須使用增量的方法來將這些碎片和轉換一起回來......

>