r - 让统计学家满意 : Stata vs. 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 检验

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