在R中计算变量组的均值和置信区间

Question

我是R语言的新手，我试图做一些我认为应该非常简单的事情，但是在线上的代码并没有帮助到我。

data <- structure(list(Group = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3), 
Time = c(1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2), mean_PctPasses = c(68.26, 
60.2666666666667, 62.05, 66.3833333333333, 59.7333333333333, 
69.7714285714286, 57.1888888888889, 63.8875, 61.1833333333333, 
59.775, 66.2666666666667, 62.12), mean_AvgPassing = c(7.3, 
7.01111111111111, 6.35, 9.26666666666667, 6.68333333333333, 
8.78571428571429, 5.87777777777778, 8.3125, 7.63333333333333, 
7.7, 8.38333333333334, 6.89), mean_AvgRush = c(0.3, -0.3, 
3.5, 0.75, 5, 1.47142857142857, 5.71111111111111, 3.3875, 
2.74, 6.6, 4.5, 5), mean_Int = c(0.2, 0.777777777777778, 
0.25, 0.5, 1.5, 0.857142857142857, 0.777777777777778, 0.75, 
0.666666666666667, 0.75, 0.833333333333333, 1.1), mean_Rate = c(99.3, 
88.5222222222222, 80.5, 106.45, 77.2333333333333, 102.885714285714, 
76.8888888888889, 100.075, 92.1166666666667, 78.55, 98.05, 
79.56)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-12L), .Names = c("Group", "Time", "mean_PctPasses", "mean_AvgPassing", 
"mean_AvgRush", "mean_Int", "mean_Rate"))

使用这个数据集，我有两个分组变量"Group"和"Time"。我想以表格的形式获得每个组合的均值和置信区间，涉及的变量是mean_PctPasses到mean_Rate，并将结果保存在表格中。我需要它以表格的形式呈现，因为我稍后会在plotly中引用它。在SPSS中完成这个任务非常容易。

我尝试了几个函数，下面是每个函数的问题。

library(rcompanion)    
ci.mean(mean_PctPasses~Group+Time, data = data)

library(DescTools)
MeanCI(data$mean_PctPasses)

library(Rmisc)
CI(data$mean_PctPasses,    ci=0.95)

MeanCI、ci.mean和CI不允许列出多个变量，并且保存为表格（仅在控制台中显示）。

by(data = data, data$Group, FUN = stat.desc)

这不允许我根据组和时间对我的数据进行分组。下面是我希望在R中构建的图表示例（显示在SPSS中）。

SPSS Example

如果有任何帮助/协助，将不胜感激。如果需要任何澄清，让我知道，我一定会编辑我的帖子。

更新

在得到了一些很好的答案（感谢Rob和Steven）之后，我觉得我需要稍微澄清一下我的问题。 我想为每个组（而不是单独地）获取所有统计信息（从mean_PctPasses到mean_Rate）。使用Rmisc显示了一个产生我想要的一个变量的统计信息的函数示例： library(Rmisc) group.UCL(mean_PctPasses~Group+Time , data, FUN=CI) 这仅为mean_PctPasses提供了以下输出 Output Using Rmisc

但是我想要的是以下内容（我已经使用Photoshop进行了编辑） 期望输出的图像

当然，这也可以显示为另一种方向（例如下面的SPSS和SEM示例）。 SPSS中的替代方向示例

- Patrick

---

根据您的编辑，您不能使用dplyr中的summarise()函数来完成这个任务吗？似乎我在下面的第一个代码块中放置的内容正是您所要求的。或者，您是在寻找一个函数来帮您完成所有工作吗？ - Steven

只是为了记录：我是rcompanion包的作者。该包中没有ci.mean函数。但是，有一个groupwiseMean函数可以为组生成置信区间，提供了几种不同的方法。例如：groupwiseMean(PctPasses ~ Group + Time, data = data)。 - Sal Mangiafico

Answer 1

假设您只想要每个组的常规非汇总t置信区间，您可以执行以下操作：

require(dplyr)
alpha <- 0.05

data %>% 
    group_by(Group, Time) %>% 
    summarize(mean = mean(mean_PctPasses),
              lower = mean(mean_PctPasses) - qt(1- alpha/2, (n() - 1))*sd(mean_PctPasses)/sqrt(n()),
              upper = mean(mean_PctPasses) + qt(1- alpha/2, (n() - 1))*sd(mean_PctPasses)/sqrt(n()))

- IceCreamToucan

Answer 2

使用R完成这个任务也很容易。

另一种方法是使用Rmisc中的CI()函数：

library(dplyr)
library(Rmisc)
library(ggplot2)

data <- 
  data %>%
  group_by(Group, Time) %>%
  dplyr::summarise(avg_PctPasses = mean(mean_PctPasses), 
            uci_PctPasses = CI(mean_PctPasses)[1], 
            lci_PctPasses = CI(mean_PctPasses)[3]) %>%
  mutate(Time = Time %>% as.factor())

承认地说，我不是CI()调用后的“神奇数字”的粉丝。

绘制数据同样简单。

data %>%
  ggplot(aes(x = Group, y = avg_PctPasses, fill = Time)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymin = lci_PctPasses, ymax = uci_PctPasses), position = "dodge")

- Steven

---

你不必使用“魔法数字”。你可以这样做，例如 CI(mean_PctPasses)['lower']。 - IceCreamToucan

@RobJensen 好的，谢谢。我会留下它们作为纪念。 - Steven

嘿@Steven，感谢您的回复。两种方法都很好，谢谢您的回应。这种格式的问题在于需要为“数据”中的每个变量重复此操作。您是否知道一种方法可以将均值和UCI / LCI输出到一个表中？我一回到实验室就会在SPSS中上传此图像。 - Patrick

Answer 3

您可能有兴趣使用基本的R图形来复制SPSS样式。

library(DescTools)

z <- with(data, 
          aggregate(mean_PctPasses, list(Time, Group), MeanCI))
z <- xtabs(x ~ Group.1 + Group.2, z)

par(mar=c(5.1,4.1,4.1,8.1))
b <- barplot(z[,,1], beside=TRUE, ylim=c(0, 140), 
             col=c("royalblue3","limegreen"), las=1, 
             xlab="Group", ylab="Mean mean_PctPasses",
             panel.first=Bg(col="grey85", border="black"))
 
ErrBars(from=z[,,2], to=z[,,3], pos=b)
legend(x="topright", legend=c("1","2"), title="Time", bty="n", 
       fill=c("royalblue3","limegreen"), inset=c(-.2, 0), xpd=TRUE)

尽管如此，您应该考虑使用点图来显示您的数据。

col <- c("royalblue3","limegreen")
PlotDot(z[,,1], args.errbars = list(from=z[,,2], to=z[,,3], mid=z[,,1]), 
        cex.pch=1.5, pch=22, bg=col, 
        lblcolor = col, lcolor = NA, 
        panel.first=abline(v=seq(0,150,10), col="grey", lty="dotted"))

- Andri Signorell

Answer 4

也许一次只处理一个变量会更容易。有一种更简单的方法。

您需要使用library()安装/加载Hmisc包。

my_data <- data %>%
    group_by(Group, Time) %>% 
    summarise(N = n(), ci = list(enframe(Hmisc::smean.cl.normal(mean_PctPasses)))) %>% 
    unnest() %>% 
    spread(name, value)

print(my_data)

这是输出结果：

使用一个变量分组（并重复所有定量变量）进行group_by()，看起来更漂亮/整洁：

 my_data <- data %>%
   group_by(Group, Time) %>% 
   summarise(N = n(), ci = list(enframe(Hmisc::smean.cl.normal(mean_PctPasses)))) %>% 
   unnest() %>% 
   spread(name, value)

 print(my_data)

输出结果:

- Anyumba

Answer 5

更新 tidyr 1.0.0

作为对之前给出的summarise解决方案的优雅替代，很高兴看到新的tidyr 1.0.0包含了一个经常被忽视的函数：unnest_wider。 有了它，您可以将代码简化为以下内容：

data.to.plot <- data %>% 
  nest(data = -"Group") %>%
  mutate(ci = map(data, ~ MeanCI(.x$mean_PctPasses))) %>% 
  unnest_wider(ci)

这提供了

# A tibble: 3 x 5
  Group data              mean lwr.ci upr.ci
  <dbl> <list>           <dbl>  <dbl>  <dbl>
1     1 <tibble [4 × 6]>  64.2   58.3   70.1
2     2 <tibble [4 × 6]>  62.6   53.9   71.4
3     3 <tibble [4 × 6]>  62.3   57.9   66.8

你可以轻松地用编程实现这个图表：

  ggplot(aes(x = Group, y = mean)) +
  geom_bar(aes (fill = Group), stat = "identity") +
  geom_errorbar(
    aes(
      ymin = lwr.ci, ymax = upr.ci,
      width = 0.5
    ),
    size = 0.5 # line thickness
  ) + 
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() 

这将为您提供

- Agile Bean