估计类似于 SAS 的均值的 ANCOVA 差异的置信区间答案

【问题标题】：Estimate Confidence Interval for ANCOVA difference of means similar to SAS估计类似于 SAS 的均值的 ANCOVA 差异的置信区间
【发布时间】：2016-02-27 18:33:56
【问题描述】：

我需要为一些医学数据展示一个 ANCOVA 模型。我习惯于使用 SAS proc GLM，我想知道是否可以像使用 SAS 一样获得均值差异的置信区间。我正在使用来自doBy 包的LSmatrix 来获取组均值。

model <- aov(yield ~ block + N * P + K, npk)
library(doBy)
k1 <- LSmatrix(model, effect="block")
linest(model,K=k1)

  estimate       se df   t.stat      p.value block
1   54.025 1.977116 14 27.32515 1.510775e-13     1
2   57.450 1.977116 14 29.05747 6.479499e-14     2
3   60.775 1.977116 14 30.73922 2.981292e-14     3
4   50.125 1.977116 14 25.35258 4.229573e-13     4
5   50.525 1.977116 14 25.55490 3.792520e-13     5
6   56.350 1.977116 14 28.50111 8.457748e-14     6

我想获得平均值块 1 - 块 2 和置信区间（或标准错误）的差异，但我不知道如何从这里开始。

结果可能是这样的：

    estimate     se upper.CI lower.CI
1-2   -3.425 ??????    ?????    ?????
1-3   -6.750 ??????    ?????    ?????
...

如果有人能告诉我如何计算第一个，我可以管理其余的:)

提前致谢。

【问题讨论】：

澄清一下。我的问题是在哪里可以找到计算标准的数据。均值差的误差。
如果您要查找数据，请准确说明您要查找的数据。我感觉您可能真正想要的是来自 stats.stackexchange.com 的统计帮助。所有数据都在对象npk 中。 linest() 结果是一个列表，其中包含输出中的确切数据。从我的解决方案中查看tmp 的内容或tmp <- linest(model,K=k1); str(tmp)
另外，这可能有助于cran.r-project.org/web/packages/lsmeans/vignettes/…

标签： r

【解决方案1】：

我。解决方案 1

> pacman::p_load(lsmeans, multcompView)
> 
> eindzl      <- lm(yield ~ block + N * P + K, npk)
> anova(eindzl)
Analysis of Variance Table

Response: yield
          Df Sum Sq Mean Sq F value   Pr(>F)   
block      5 343.29  68.659  4.3911 0.012954 * 
N          1 189.28 189.282 12.1055 0.003684 **
P          1   8.40   8.402  0.5373 0.475637   
K          1  95.20  95.202  6.0886 0.027114 * 
N:P        1  21.28  21.282  1.3611 0.262841   
Residuals 14 218.90  15.636                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> eindzl.rg   <- ref.grid(eindzl)
> eindzl.lsm  <- lsmeans(eindzl.rg, "block")
> cld(eindzl.lsm, alpha = .05)
 block lsmean       SE df lower.CL upper.CL .group
 4     50.125 1.977116 14 45.88451 54.36549  1    
 5     50.525 1.977116 14 46.28451 54.76549  1    
 1     54.025 1.977116 14 49.78451 58.26549  12   
 6     56.350 1.977116 14 52.10951 60.59049  12   
 2     57.450 1.977116 14 53.20951 61.69049  12   
 3     60.775 1.977116 14 56.53451 65.01549   2   

Results are averaged over the levels of: N, P, K 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 6 estimates 
significance level used: alpha = 0.05 

pairs(eindzl.lsm)
 contrast estimate       SE df t.ratio p.value
 1 - 2      -3.425 2.796064 14  -1.225  0.8180
 1 - 3      -6.750 2.796064 14  -2.414  0.2160
 1 - 4       3.900 2.796064 14   1.395  0.7295
 1 - 5       3.500 2.796064 14   1.252  0.8049
 1 - 6      -2.325 2.796064 14  -0.832  0.9564
 2 - 3      -3.325 2.796064 14  -1.189  0.8348
 2 - 4       7.325 2.796064 14   2.620  0.1560
 2 - 5       6.925 2.796064 14   2.477  0.1960
 2 - 6       1.100 2.796064 14   0.393  0.9985
 3 - 4      10.650 2.796064 14   3.809  0.0190
 3 - 5      10.250 2.796064 14   3.666  0.0248
 3 - 6       4.425 2.796064 14   1.583  0.6217
 4 - 5      -0.400 2.796064 14  -0.143  1.0000
 4 - 6      -6.225 2.796064 14  -2.226  0.2856
 5 - 6      -5.825 2.796064 14  -2.083  0.3485

Results are averaged over the levels of: N, P, K 
P value adjustment: tukey method for comparing a family of 6 estimates 
>

二。解决方案2

请随意使用您自己的公式来计算上限和下限，或者如果您需要帮助确定公式，请访问CrossValidated。此示例基于标准错误并显示您需要的代码。 CrossValidated 用于统计。

pacman::p_load(data.table,doBy)

model <- aov(yield ~ block + N * P + K, npk)
k1    <- LSmatrix(model, effect="block")
linest(model,K=k1)

tmp <- linest(model,K=k1)

tmp <- cbind(as.data.frame(tmp$coef), as.data.frame(tmp$grid))

tmp2 <- as.data.frame(matrix(ncol=4,nrow=9))
setnames(tmp2, c("group","estimate","lower bound of CI", "upper bound of CI"))

for(i in 1:5){
  n                  <- i + 1
  tmp2$estimate[i]   <- tmp$estimate[tmp$block == i] -  tmp$estimate[tmp$block == n]
  tmp2$group[i]      <- paste(i,n,sep="-")
  tmp2[i,3]          <- (tmp$estimate[tmp$block == i] - tmp$se[tmp$block == i]) - (tmp$estimate[tmp$block == n] + tmp$se[tmp$block == n])
  tmp2[i,4]          <- (tmp$estimate[tmp$block == i] + tmp$se[tmp$block == i]) - (tmp$estimate[tmp$block == n] - tmp$se[tmp$block == n])

}


for(i in 1:4){
  n                   <- i + 2
  nn                  <- 5+i    
  tmp2$estimate[nn]   <- tmp$estimate[tmp$block == i] -  tmp$estimate[tmp$block == n]
  tmp2$group[nn]      <- paste(i,n,sep="-")
  tmp2[nn,3]          <- (tmp$estimate[tmp$block == i] - tmp$se[tmp$block == i]) - (tmp$estimate[tmp$block == n] + tmp$se[tmp$block == n])
  tmp2[nn,4]          <- (tmp$estimate[tmp$block == i] + tmp$se[tmp$block == i]) - (tmp$estimate[tmp$block == n] - tmp$se[tmp$block == n])
}

tmp2

  group   estimate lower bound of CI upper bound of CI
1   1-2     -3.425         -7.379232         0.5292322
2   2-3     -3.325         -7.279232         0.6292322
3   3-4     10.650          6.695768        14.6042322
4   4-5     -0.400         -4.354232         3.5542322
5   5-6     -5.825         -9.779232        -1.8707678
6   1-3     -6.750        -10.704232        -2.7957678
7   2-4      7.325          3.370768        11.2792322
8   3-5     10.250          6.295768        14.2042322
9   4-6     -6.225        -10.179232        -2.2707678

请注意，如果您没有选择这个特定的包和功能，您将能够在没有太多自定义操作的情况下获得类似的结果。通常，您在 R 中输入的内容比在 SAS 中输入的要少。

【讨论】：

据我了解，您正在计算均值差异的 CI 边界作为均值边界的差异（lower.mean.CI[block 1] - upper.mean.CI[block 2]）虽然你没有乘以标准。误差 1.96 (qnorm)。但即使你这样做了。我不认为这是计算差异置信区间的合法方法。
@eindzl 如果您想要不同的置信区间公式，您不能根据自己的喜好进行修改吗？或者提供一个方程来做到这一点。如果你不知道它应该是什么，那么 CrossValidated (stats.stackexchange.com) 将是你要去的地方。标准误差是标准差的样本估计值。我的例子是，下限是最低差异，其中两个估计值都在平均值的 1 个标准差之内，反之亦然。 StackOverflow 不是用于统计的——这是用于编程的——CrossValidated 是用于统计的。
@eindzl 我正在添加另一种您可能更喜欢的方法
谢谢@hack-r，pairs{lsmeans} 方法正是我想要的。可能我应该去交叉验证，但我认为这更多是找到正确的功能和包的问题。
@eindzl 不客气。一切都有一个 SE 站点——如果您正在寻找一个工具（一个包），那么他们实际上希望您使用名为 Software Recommendations 的 StackExchange。