让您的统计学家满意：Stata 与 R 学生的 t 检验

Question

第 1 章：按性别划分的平均年龄

我经常与对统计输出有非常具体要求的流行病学家和统计学家合作，但我经常无法在 R 中重现完全相同的内容（我们的流行病学家在 Stata 中工作）。

让我们从一个简单的例子开始，学生 t 检验。我们感兴趣的是首次诊断的平均年龄和置信区间的差异。

1) 在 R 中创建一些示例数据

set.seed(41)

cohort <- data.frame(
          id = seq(1,100),
          gender = sample(c(rep(1,33), rep(2,67)),100),
          age    = sample(seq(0,50),100, replace=TRUE)
          )

# save to import into Stata
# write.csv(cohort, "cohort.csv", row.names = FALSE)

b) 导入数据并在 Stata 中运行 t-test

import delimited "cohort.csv"
ttest age, by(gender)

我们想要的是平均值的绝对差 = 3.67 年，组合置信区间 = 95% CI：24.59 - 30.57

b) 在 R 中运行 t 检验

t.test(age~gender, data=cohort)

t.test(cohort$age[cohort$gender == 1])

t.test(cohort$age[cohort$gender == 2])

t.test(cohort$age)

肯定有另一种方法，而不是在 R 中运行 4 次 t 检验！

Answer 1

您可以尝试将所有内容放在一个函数中，并使用一些 tidyverse 魔法。当然，可以根据需要编辑输出。 boom 的 tidy 将用于很好的输出。

foo <- function(df, x, y){
  require(tidyverse)
  require(broom)
  a1 <- df %>% 
    select(ep=!!x, gr=!!y) %>% 
    mutate(gr=as.character(gr)) %>% 
    bind_rows(mutate(., gr="ALL")) %>% 
    split(.$gr) %>% 
    map(~tidy(t.test(.$ep))) %>% 
    bind_rows(.,.id = "gr") %>% 
    mutate_if(is.factor, as.character)
  tidy(t.test(as.formula(paste(x," ~ ",y)), data=df)) %>% 
    mutate_if(is.factor, as.character) %>% 
    mutate(gr="vs") %>% 
    select(gr, estimate, statistic, p.value,parameter, conf.low, conf.high, method, alternative) %>% 
    bind_rows(a1, .)}


foo(cohort, "age", "gender")
   gr  estimate statistic      p.value parameter  conf.low conf.high                  method alternative
1   1 25.121212  9.545737 6.982763e-11  32.00000  19.76068 30.481745       One Sample t-test   two.sided
2   2 28.791045 15.699854 5.700541e-24  66.00000  25.12966 32.452428       One Sample t-test   two.sided
3 ALL 27.580000 18.301678 1.543834e-33  99.00000  24.58985 30.570147       One Sample t-test   two.sided
4  vs -3.669833 -1.144108 2.568817e-01  63.37702 -10.07895  2.739284 Welch Two Sample t-test   two.sided

我建议从头开始使用这个

foo <- function(df){
 a1 <- broom::tidy(t.test(age~gender, data=df))
 a2 <- broom::tidy(t.test(df$age))
 a3 <- broom::tidy(t.test(df$age[df$gender == 1]))
 a4 <- broom::tidy(t.test(df$age[df$gender == 2]))
 list(rbind(a2, a3, a4), a1)
}

foo(cohort)
[[1]]
  estimate statistic      p.value parameter conf.low conf.high            method alternative
1 27.58000 18.301678 1.543834e-33        99 24.58985  30.57015 One Sample t-test   two.sided
2 25.12121  9.545737 6.982763e-11        32 19.76068  30.48174 One Sample t-test   two.sided
3 28.79104 15.699854 5.700541e-24        66 25.12966  32.45243 One Sample t-test   two.sided

[[2]]
   estimate estimate1 estimate2 statistic   p.value parameter  conf.low conf.high                  method alternative
1 -3.669833  25.12121  28.79104 -1.144108 0.2568817  63.37702 -10.07895  2.739284 Welch Two Sample t-test   two.sided

Answer 2

您可以制作自己的功能：

tlimits <- function(data, group){
  error <- qt(0.975, df = length(data)-1)*sd(data)/(sqrt(length(data)))
  mean <- mean(data)
  means <- tapply(data, group, mean)
  c(abs(means[1] - means[2]), mean - error, mean + error)
}

tlimits(cohort$age, cohort$gender)
        1                     
 3.669833 24.589853 30.570147 

Answer 3

我们想要的是平均值的绝对差 = 3.67 年，组合置信区间 = 95% CI：24.59 - 30.57

请注意，R 的 t.test 进行 t 检验，而您需要均值差和“组合置信区间”（即忽略分组变量的均值 CI）。所以你不想要一个 t 检验，而是别的东西。

您可以使用例如：

diff(with(cohort, tapply(age, gender, mean)))
# 3.669833 
# no point in using something more complicated e.g., t-test or lm

... 和 CI 使用，例如：

confint(lm(age~1, data=cohort))
#                2.5 %   97.5 %
# (Intercept) 24.58985 30.57015

显然，如果您经常需要，您可以轻松地将这两个步骤组合成一个功能。

doit <- function(a,b) c(diff= diff(tapply(a,b,mean)), CI=confint(lm(a~1)))
with(cohort, doit(age,gender))
#   diff.2       CI1       CI2 
# 3.669833 24.589853 30.570147