来自 SPSS 的难民在 R 中遇到 fit.contrast 问题

机器算法验证 r 方差分析 对比
2022-03-31 12:56:06

我正试图让自己摆脱 SPSS 并切换到 R 进行心理统计，我大部分时间都到了那里，但我遇到fit.contrast了gmodels包的真正障碍。据我所知，即使提供的示例也不能正常工作——要么我严重误解了 R 中的对比（总是有可能）。

首先让我们看一下文档附带的示例代码gmodels::fit.contrast()：

set.seed(03215)
Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
Time <- factor(sample(1:3, 1000, replace=TRUE))
y <- rnorm(1000)
data <- data.frame(y, Genotype, Time)

model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
model.tables(model, "means")

fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 , df=T)

对我来说，这会产生以下输出：

                   Estimate Std. Error   t value Pr(>|t|)   lower CI  upper CI
GenotypeKO vs WT 0.01683876  0.1095764 0.1536714   0.8779 -0.1981888 0.2318664

现在考虑后面的示例代码：

model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
              contrasts=list(Genotype=make.contrasts(cm.G),
                             Time=make.contrasts(cm.T) )
)

summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
                           Time = list( "1 vs 2" = 1,
                                        "2 vs 3" = 2 ) ) )

这为我产生了这个输出：

                                  Df Sum Sq Mean Sq F value Pr(>F)
Genotype                           1    1.2  1.1687   1.121  0.290
  Genotype: KO vs WT               1    1.2  1.1687   1.121  0.290
Time                               2    0.7  0.3677   0.353  0.703
  Time: 1 vs 2                     1    0.2  0.1784   0.171  0.679
  Time: 2 vs 3                     1    0.6  0.5571   0.534  0.465
Genotype:Time                      2    1.2  0.5760   0.552  0.576
  Genotype:Time: KO vs WT.1 vs 2   1    0.3  0.3265   0.313  0.576
  Genotype:Time: KO vs WT.2 vs 3   1    0.8  0.8256   0.792  0.374
Residuals                        994 1036.6  1.0429

我在这里完全糊涂了。基因型的两个测试的 p 值应该匹配，但它们不匹配！p 值不相同，F 值不是 t 值的平方。此外，为对比度估计 (0.0168) 计算的值不是这两种方法的差异......为什么不呢？

任何有助于理解 fit.contrast 在这里做什么的帮助将不胜感激！

2个回答

我最初的评论是：为什么你需要gmodels包中的这些功能来完成这个？这些是可以在基础 R 中直接完成的基本任务，这是我推荐的。

但是要直接回答您的问题，这里的问题fit.contrast是使用“类型 3”的效果测试，而summary使用“类型 2”的测试。这可以Anova从car包中最容易地验证，它使您可以选择所需的平方和类型：

library(car)

> Anova(model, type=2)
Anova Table (Type II tests)

Response: y
Sum Sq  Df F value Pr(>F)
Genotype         1.09   1  1.0495 0.3059
Time             0.74   2  0.3526 0.7029
Genotype:Time    1.15   2  0.5524 0.5758
Residuals     1036.64 994       

> Anova(model, type=3)
Anova Table (Type III tests)

Response: y
Sum Sq  Df F value Pr(>F)
(Intercept)      0.10   1  0.0935 0.7599
Genotype         0.02   1  0.0236 0.8779
Time             1.87   2  0.8954 0.4088
Genotype:Time    1.15   2  0.5524 0.5758
Residuals     1036.64 994

可以在这里找到这两种方法之间差异的描述： https ://stat.ethz.ch/pipermail/r-help/2006-August/111854.html

强烈建议不要使用保留字作为R对象名称。对于基因型对比，很明显来自 ANOVA 的 p 值和来自 gmodels 的对比声明是相同的并且 F = t ^ 2。

library(gmodels)
set.seed(03215)
Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
Time <- factor(sample(1:3, 1000, replace=TRUE))
y <- rnorm(1000)
dat <- data.frame(y, Genotype, Time)

fit1 <- aov( y ~ Genotype + Time + Genotype:Time, data=dat)
summary(fit1)


              Df Sum Sq Mean Sq F value Pr(>F)
Genotype        1    1.2  1.1687   1.121  0.290
Time            2    0.7  0.3677   0.353  0.703
Genotype:Time   2    1.2  0.5760   0.552  0.576
Residuals     994 1036.6  1.0429

model.tables(fit1, "means")
Tables of means
Grand mean

0.01447773 

 Genotype 
           KO        WT
     -0.02154   0.04693
rep 474.00000 526.00000

 Time 
            1         2         3
      0.03267   0.03313  -0.02539
rep 350.00000 334.00000 316.00000

 Genotype:Time 
        Time
Genotype 1      2      3     
     KO    0.02   0.02  -0.11
     rep 160.00 155.00 159.00
     WT    0.04   0.04   0.06
     rep 190.00 179.00 157.00

如 (-1) (-0.02154)+(1) (0.04693) = 0.06847

fit.contrast(fit1, "Genotype", rbind("KO vs WT"=c(1, -1)), conf=0.95, df=TRUE)

                  Estimate Std. Error  t value  Pr(>|t|)  DF   lower CI upper CI
GenotypeKO vs WT 0.06869178 0.06477589 1.060453 0.2891962 994 -0.0584214 0.195805
其它你可能感兴趣的问题
上一篇带有（太多？）分类自变量的有序 logit 下一篇做什么ββ告诉我们线性回归分析？