R中混合的SAS proc的等价物

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 的确切结果，但我可以通过与contrast 打交道，如下所述：

lmer for SAS PROC MIXED Users：第 6 页

比较 SAS PROC MIXED 和 lmer one 生成的估计值时 必须小心考虑用于定义 因素的影响。在 SAS 中，具有截距和定性的模型 因子根据截距和指标定义 除了因子的最后一个级别之外的所有变量。默认 S 中的行为是使用 Helmert 对比作为因子。在一个 这些平衡因素提供了一组正交对比。在 R 中 默认是“治疗”对比，几乎与 SAS 参数化，只是它们删除了第一个指标 一级，不是最后一级。如有疑问，请检查哪些对比是 与对比功能一起使用。为了便于比较， 你可能会觉得值得声明

options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly")) 在会话开始时。

输入：

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  

Answer 2

您正在使用 SAS 选项 NOBOUND，它允许对方差进行负估计，并且您得到负估计。这对于 lmer 是不可能的，它将方差限制为正数。

我们可以尝试手动获取 SAS 结果。首先，注意等效的lmer 语法是：

lmer(y ~ 1 + treatment + (1+treatment|sample), REML=FALSE, data = dat)

让我们最大化对数似然，允许负方差：

dattxt <- "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 
"

dat <- read.table(text = dattxt)
names(dat) <- c("sample", "treatment", "y")
dat$sample <- as.factor(dat$sample)

opts <- options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))

library(lme4)
fit <- lmer(y ~ 1 + treatment + (1+treatment|sample), REML=FALSE, data = dat) 

# marginal variance matrix in function of variance components
Vfun <- function(fit, vcs){
  Z <- getME(fit, "Z")
  n <- getME(fit, "n")
  l_i <- getME(fit, "l_i")
  sigma2_a <- vcs[1]
  sigma2_b <- vcs[2]
  sigma_ab <- vcs[3]
  sigma2 <- vcs[4]
  G <- matrix(c(sigma2_a, sigma_ab, sigma_ab, sigma2_b), nrow = 2)
  R <- Diagonal(n, sigma2)
  Z %*% bdiag(rep(list(G),l_i)) %*% t(Z) + R
}


# minus log-likelihood
library(mvtnorm)
logLHD <- function(params, fit){
  X <- getME(fit, "X")
  beta <- params[1:ncol(X)]
  y <- getME(fit, "y")
  vcs <- tail(params, length(params)-ncol(X))
  V <- as.matrix(Vfun(fit, vcs))
  if(any(eigen(V)$values <= 0)){
    return(runif(1, 1e7, 1e8)) # return a high-value if V is not positive
  }
  -dmvnorm(y, c(X%*%beta), sigma = V, log = TRUE)  
}

# optimization of log-likelihood
library(dfoptim)
start <- 
  c(fixef(fit), vc$sample[1,1], vc$sample[2,2], vc$sample[1,2], sigma(fit)^2)
names(start)[3:6] <- 
  c("sample.Intercept", "sample.trt1", "covariance", "sigma2")
opt <- hjkb(start, logLHD, lower=c(-Inf,-Inf,-Inf,-Inf,-Inf,0), fit=fit)

### results 
opt$par
# (Intercept) treatmenttrt1 sample.Intercept  sample.trt1 covariance     sigma2 
# 3.33912840    -0.10721533       5.50671885  -0.16909628 0.07275635 0.31812378 

残差方差与使用 SAS 获得的方差相同。要获得其他 SAS 结果，必须对我们的结果做一些体操，我不明白为什么，但我们以这种方式获得它们：

### SAS results
opt$par[["sample.Intercept"]] + opt$par[["covariance"]]
# 5.579475
opt$par[["sample.trt1"]] / 2
# -0.08454814

请注意，负方差确实可以更好地最大化对数似然：

### remark: lmer achieves a lower log-likelihood
logLik(fit)
# 'log Lik.' -27.88947 (df=6)
-opt$value
# -26.43355

如果有人能解释所需的体操，我将不胜感激......

---

编辑

抱歉，这不是好模型。型号为：

lmer(y ~ 1 + treatment + (1|sample/treatment), REML=FALSE, data = dat)

以下是 SAS 结果：

opts <- options(contrasts = c(factor = "contr.SAS", ordered = "contr.poly"))
library(lme4)
fit <- lmer(y ~ 1+treatment+(1|sample/treatment), REML=FALSE, data = dat) 
vc <- VarCorr(fit)

Vfun <- function(fit, vcs){
  Z <- getME(fit, "Z")
  n <- getME(fit, "n")
  l_i <- getME(fit, "l_i")
  G <- Diagonal(sum(l_i), rep(vcs[1:2], l_i))
  R <- Diagonal(n, vcs[3])
  Z %*% G %*% t(Z) + R
}

library(mvtnorm)
logLHD <- function(params, fit){
  X <- getME(fit, "X")
  beta <- params[1:ncol(X)]
  y <- getME(fit, "y")
  vcs <- tail(params, length(params)-ncol(X))
  V <- as.matrix(Vfun(fit, vcs))
  if(any(eigen(V)$values <= 0)) return(runif(1, 1e7, 1e8))
  -dmvnorm(y, c(X%*%beta), sigma = V, log = TRUE)  
}

library(dfoptim)
start <- c(fixef(fit), vc[[1]], vc[[2]], sigma(fit)^2)
opt <- hjkb(start, logLHD, lower=c(-Inf,-Inf,-Inf,-Inf,0), fit=fit)
opt$par[3:5]
# -0.08454877    5.57947601    0.31812697