r - 相当于 R 中混合的 SAS 过程

Question

我正在尝试在 R 中转换以下 SAS 代码以获得与 SAS 相同的结果。这是 SAS 代码:

    DATA plants; 
    INPUT  sample $  treatmt $ y ; 
    cards; 

    1   trt1    6.426264755 
    1   trt1    6.95419631 
    1   trt1    6.64385619 
    1   trt2    7.348728154 
    1   trt2    6.247927513 
    1   trt2    6.491853096 
    2   trt1    2.807354922 
    2   trt1    2.584962501 
    2   trt1    3.584962501 
    2   trt2    3.906890596 
    2   trt2    3 
    2   trt2    3.459431619 
    3   trt1    2 
    3   trt1    4.321928095 
    3   trt1    3.459431619 
    3   trt2    3.807354922 
    3   trt2    3 
    3   trt2    2.807354922 
    4   trt1    0 
    4   trt1    0 
    4   trt1    0 
    4   trt2    0 
    4   trt2    0 
    4   trt2    0 
    ; 
    RUN; 

    PROC MIXED ASYCOV NOBOUND  DATA=plants ALPHA=0.05 method=ML; 
    CLASS sample treatmt; 
    MODEL  y = treatmt ; 
    RANDOM int treatmt/ subject=sample ; 
    RUN; 

我从 SAS 获得以下协方差估计:

Intercept   sample  ==> 5.5795     
Treatmt  sample ==> -0.08455      
Residual         ==> 0.3181    

我在 R 中尝试了以下操作，但得到了不同的结果。

s=as.factor(sample) 
lmer(y~ 1+treatmt+(1|treatmt:s),REML=FALSE) 

Answer 1

我不知道您是否能够获得从 SAS 到 R 的准确结果，但我能够通过处理此处概述的对比度来接近结果:

lmer for SAS PROC MIXED Users :第 6 页

When comparing estimates produced by SAS PROC MIXED and by lmer one must be careful to consider the contrasts that are used to define the effects of factors. In SAS a model with an intercept and a qualitative factor is defined in terms of the intercept and the indicator variables for all but the last level of the factor. The default behaviour in S is to use the Helmert contrasts for the factor. On a balanced factor these provide a set of orthogonal contrasts. In R the default is the “treatment” contrasts which are almost the same as the SAS parameterization except that they drop the indicator of the first level, not the last level. When in doubt, check which contrasts are being used with the contrasts function. To make comparisons easier, you may find it worthwhile to declare

options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly")) at the beginning of your session.

输出:

df <- structure(list(sample = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L), 
    treatmt = c("trt1", "trt1", "trt1", "trt2", "trt2", "trt2", 
    "trt1", "trt1", "trt1", "trt2", "trt2", "trt2", "trt1", "trt1", 
    "trt1", "trt2", "trt2", "trt2", "trt1", "trt1", "trt1", "trt2", 
    "trt2", "trt2"), y = c(6.426264755, 6.95419631, 6.64385619, 
    7.348728154, 6.247927513, 6.491853096, 2.807354922, 2.584962501, 
    3.584962501, 3.906890596, 3, 3.459431619, 2, 4.321928095, 
    3.459431619, 3.807354922, 3, 2.807354922, 0, 0, 0, 0, 0, 
    0)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-24L), .Names = c("sample", "treatmt", "y"))

当前代码:

options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))
df$sample=as.factor(df$sample) 
lmer(y~ 1+treatmt+(1|treatmt:sample),REML=FALSE, data = df) 

电流输出:

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: y ~ 1 + treatmt + (1 | treatmt:sample)
   Data: df
     AIC      BIC   logLik deviance df.resid 
 80.3564  85.0686 -36.1782  72.3564       20 
Random effects:
 Groups         Name        Std.Dev.
 treatmt:sample (Intercept) 2.344   
 Residual                   0.564   
Number of obs: 24, groups:  treatmt:sample, 8
Fixed Effects:
(Intercept)  treatmttrt1  
     3.3391      -0.1072