2

有谁知道如何在具有因子设计的 ANCOVA 模型中进行事后测试？

我有两个向量，由 23 个基线值（协变量）和 23 个治疗后的值（因变量）组成，并且我有两个具有两个水平的因子。我创建了一个 ANCOVA 模型并计算了调整后的平均值、标准误差和置信区间。例子：

library(effects)

baseline = c(0.7672,1.846,0.6487,0.4517,0.5599,0.2255,0.5946,1.435,0.5374,0.4901,1.258,0.5445,1.078,1.142,0.5,1.044,0.7824,1.059,0.6802,0.8003,0.5547,1.003,0.9213)
after_treatment = c(0.4222,1.442,0.8436,0.5544,0.8818,0.08789,0.6291,1.23,0.4093,0.7828,-0.04061,0.8686,0.8525,0.8036,0.3758,0.8531,0.2897,0.8127,1.213,0.05276,0.7364,1.001,0.8974)

age = factor(c(rep(c("Young","Old"),11),"Young")) 
treatment = factor(c(rep("Drug",12),rep("Placebo",11)))

ANC = aov(after_treatment ~ baseline + treatment*age)

effect_treatage = effect("treatment*age",ANC)
data.frame(effect_treatage)

  treatment   age       fit        se     lower     upper
1      Drug   Old 0.8232137 0.1455190 0.5174897 1.1289377
2   Placebo   Old 0.6168641 0.1643178 0.2716452 0.9620831
3      Drug Young 0.5689036 0.1469175 0.2602413 0.8775659
4   Placebo Young 0.7603360 0.1462715 0.4530309 1.0676410

现在我想测试调整后的方法之间是否存在差异

新安慰剂：新药物
旧安慰剂：旧药物
新安慰剂：旧药物
旧安慰剂：新药物

所以我尝试了：

library(multcomp)
pH = glht(ANC, linfct = mcp(treatment*age="Tukey"))
# Error: unexpected '=' in "ph = glht(ANC_nback, linfct = mcp(treat*age="

和：

pH = TukeyHSD(ANC)
# Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep3'
# In addition: Warning message:
# In replications(paste("~", xx), data = mf) : non-factors ignored: baseline

有谁知道如何解决这个问题？

非常感谢！

PS 更多信息见

R：如何以图形方式绘制调整后的均值、SE、CI ANCOVA

4

3 回答 3

3

如果您想使用. multcomp_ lsmeans_multcomp?lsmlsmeansglht()

baseline = c(0.7672,1.846,0.6487,0.4517,0.5599,0.2255,0.5946,1.435,0.5374,0.4901,1.258,0.5445,1.078,1.142,0.5,1.044,0.7824,1.059,0.6802,0.8003,0.5547,1.003,0.9213)
after_treatment = c(0.4222,1.442,0.8436,0.5544,0.8818,0.08789,0.6291,1.23,0.4093,0.7828,-0.04061,0.8686,0.8525,0.8036,0.3758,0.8531,0.2897,0.8127,1.213,0.05276,0.7364,1.001,0.8974)
age = factor(c(rep(c("Young","Old"),11),"Young")) 
treatment = factor(c(rep("Drug",12),rep("Placebo",11)))
Treat <- data.frame(baseline, after_treatment, age, treatment)

ANC  <-  aov(after_treatment ~ baseline + treatment*age, data=Treat)

library(multcomp)
library(lsmeans)

summary(glht(ANC, linfct = lsm(pairwise ~ treatment * age)))
## Note: df set to 18
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: aov(formula = after_treatment ~ baseline + treatment * age, data = Treat)
## 
## Linear Hypotheses:
##                                  Estimate Std. Error t value Pr(>|t|)
## Drug,Old - Placebo,Old == 0       0.20635    0.21913   0.942    0.783
## Drug,Old - Drug,Young == 0        0.25431    0.20698   1.229    0.617
## Drug,Old - Placebo,Young == 0     0.06288    0.20647   0.305    0.990
## Placebo,Old - Drug,Young == 0     0.04796    0.22407   0.214    0.996
## Placebo,Old - Placebo,Young == 0 -0.14347    0.22269  -0.644    0.916
## Drug,Young - Placebo,Young == 0  -0.19143    0.20585  -0.930    0.789
## (Adjusted p values reported -- single-step method)

这消除了重新参数化的需要。lsmeans单独使用也可以达到同样的效果：

lsmeans(ANC, list(pairwise ~ treatment * age))
## $`lsmeans of treatment, age`
##  treatment age      lsmean        SE df  lower.CL  upper.CL
##  Drug      Old   0.8232137 0.1455190 18 0.5174897 1.1289377
##  Placebo   Old   0.6168641 0.1643178 18 0.2716452 0.9620831
##  Drug      Young 0.5689036 0.1469175 18 0.2602413 0.8775659
##  Placebo   Young 0.7603360 0.1462715 18 0.4530309 1.0676410
## 
## Confidence level used: 0.95 
## 
## $`pairwise differences of contrast`
##  contrast                       estimate        SE df t.ratio p.value
##  Drug,Old - Placebo,Old       0.20634956 0.2191261 18   0.942  0.7831
##  Drug,Old - Drug,Young        0.25431011 0.2069829 18   1.229  0.6175
##  Drug,Old - Placebo,Young     0.06287773 0.2064728 18   0.305  0.9899
##  Placebo,Old - Drug,Young     0.04796056 0.2240713 18   0.214  0.9964
##  Placebo,Old - Placebo,Young -0.14347183 0.2226876 18  -0.644  0.9162
##  Drug,Young - Placebo,Young  -0.19143238 0.2058455 18  -0.930  0.7893
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
于 2015-07-08T19:55:05.930 回答
2

您需要在;中使用which参数 TukeyHSD“列出应计算间隔的拟合模型中的项”。这是必需的，因为您在模型中有一个非因子变量（“基线”）。变量包含时会引起麻烦，which未指定时为默认值。

ANC = aov(after_treatment ~ baseline + treatment*age)
TukeyHSD(ANC, which = c("treatment:age"))

如果您希望使用更灵活的glht，请参阅第 3 节，第 8 页 -这里

于 2014-05-13T10:56:50.447 回答
2

这里可以重新参数化：

treatAge <- interaction(treatment, age)
ANC1 <- aov(after_treatment ~ baseline + treatAge)
#fits are equivalent:
all.equal(logLik(ANC), logLik(ANC1))
#[1] TRUE

library(multcomp)
summary(glht(ANC1, linfct = mcp(treatAge="Tukey")))

#    Simultaneous Tests for General Linear Hypotheses
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#
#Fit: aov(formula = after_treatment ~ baseline + treatAge)
#
#Linear Hypotheses:
#                                 Estimate Std. Error t value Pr(>|t|)
#Placebo.Old - Drug.Old == 0      -0.20635    0.21913  -0.942    0.783
#Drug.Young - Drug.Old == 0       -0.25431    0.20698  -1.229    0.617
#Placebo.Young - Drug.Old == 0    -0.06288    0.20647  -0.305    0.990
#Drug.Young - Placebo.Old == 0    -0.04796    0.22407  -0.214    0.996
#Placebo.Young - Placebo.Old == 0  0.14347    0.22269   0.644    0.916
#Placebo.Young - Drug.Young == 0   0.19143    0.20585   0.930    0.789
#(Adjusted p values reported -- single-step method)
于 2014-05-13T10:56:19.940 回答