类SAS均值ANCOVA差的置信区间估计

Question

我需要为一些医学数据提供一个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的差值和置信区间(或std。错误)，但我不知道如何从这里继续。

结果可能是这样的：

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

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

提前谢谢。

Answer 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 
> 

II.解决方案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中键入的少得多。