使用R中的ANCOVA设计获取每个对象的调整值

Question

我感兴趣的是，在控制抑郁的同时，在记忆任务上是否存在群体差异(对照组和PTSD)。我运行了一个ANCOVA模型。这是我的语法：

memory <- c(0.7085, 1, 0.375, 0.4585, 1, 0.8, 0.875, 0.35, 0.8, 0.9, 1, 0.425, 0.8335, 1, 0.4665, 0.75, 1, 0.6665, 0.8, 0.9, 1, 0.375, 0.5, 0.625, 0.5, 1, 0.7085, 0.625, 0.5835, 0.8335, 0.6, 1, 0.75, 0.4585, 0.6, 0.7, 0.6665, 0.8, 1, 0.875, 0.5, 1, 0, 0.4165)
group <- factor(c(rep("controls", 31), rep("ptsd", 13)), levels=c("controls", "ptsd"))
BDI <- c(1, 7, 0, 5, 3, 2, 4, 18, 7, 8, 0, 2, 10, 17, 1, 22, 6, 5, 4, 4, 0, 2, 4, 4, 8, 14, 4, 11, 3, 9, 11, 15, 33, 20, 31, 23, 25, 20, 21, 36, 14, 41, 22, 16)

 #ANCOVA
 Anova(lm(memory ~ BDI + group), type="III")

 #Anova Table (Type III tests)

 #Response: memory
 #            Sum Sq Df  F value    Pr(>F)    
 #(Intercept) 8.3847  1 148.3227 3.331e-15 ***
 #BDI         0.0951  1   1.6829    0.2018    
 #group       0.1136  1   2.0087    0.1640    
 #Residuals   2.3177 41                       
 #---
 #Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我看到在这项任务上没有显著的组差异，p值= 0.1640。但是，我也想获得每个组的调整后的平均值：

data         <- data.frame(group, BDI, memory)
model        <- lm(memory ~ BDI + group)
BDI.predict  <- rep(mean(BDI), 2)
data.predict <- data.frame(group = group, BDI = BDI.predict)
adjmeans     <- predict(model, data.predict)
adjmeans

我看到对照组的调整平均值是0.7636093，创伤后应激障碍组的调整平均值是0.5810471。

我的问题是如何获得每个主题的调整值？

我想做的一件事是看看每一组是否执行了随机(单样本t-检验)，并且我想使用调整后的值(在控制抑郁之后)，而不是原始值。

Answer 1

John Fox的effects包将为您提供答案：

 library (effects)
 help(pac=effects)
 mod.dat <- lm(memory ~ BDI + group,
     data=data)
 mod.dat <- lm(memory ~ BDI + group,
     data=data)
 eff.dat <- allEffects(mod.dat, xlevels=list(group=c("controls","ptsd") ))
 eff.dat
#----- output
 model: memory ~ BDI + group

 BDI effect
BDI
        0        10        20        30        40 
0.6241901 0.6975066 0.7708231 0.8441396 0.9174561 

 group effect
group
 controls      ptsd 
0.7636093 0.5810471 

这是包含print.eff函数所传递的值的eff.dat对象内容的“组”部分：

str(eff.dat) #然后剪掉了第一部分：

$ group:List of 21
  ..$ term                   : chr "group"
  ..$ formula                :Class 'formula'  language memory ~ BDI + group
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  ..$ response               : chr "memory"
  ..$ variables              :List of 1
  .. ..$ group:List of 3
  .. .. ..$ name     : chr "group"
  .. .. ..$ is.factor: logi TRUE
  .. .. ..$ levels   : chr [1:2] "controls" "ptsd"
  ..$ fit                    : num [1:2, 1] 0.764 0.581
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "1" "2"
  .. .. ..$ : NULL

因此，如果你想把它们作为一个向量，只需按名称对列表进行寻址：

> eff.dat$group$fit
       [,1]
1 0.7636093
2 0.5810471