R 通过函数计算总体平均值的置信区间

Question

我打算创建一个脚本，可以生成统计信息，根据类别进行分组（在下面的示例中；SDP 和 M346）。我需要 n、平均值、标准差 (sd)、标准误差 (se) 和置信区间。除了置信区间之外，我已经准备好了所有代码。我无法使用

aggregate

函数将置信区间的计算嵌入到组中。

tab.sep.file

我试图打印按列

SDP

和

M346

分组的平均值/sd/se/置信区间。

aggregate(tab.sep.file$Hic8, FUN=mean, by=list(SDP=tab.sep.file$SDP, trait=tab.sep.file$M346))

通过将函数切换为长度、标准差并创建标准误差公式，我成功地一次性获得了分组数据的统计数据。 上面脚本的输出如下：

计算整个数据（不分组）的置信区间：

mean <- mean(tab.sep.file$Hic8)
n <- length(tab.sep.file$Hic8)
standard_deviation <- sd(tab.sep.file$Hic8)
standard_error <- standard_deviation / sqrt(n)
alpha = 0.05
degrees_of_freedom = n - 1
t_score = qt(p=alpha/2, df=degrees_of_freedom, lower.tail=F)
margin_error <- t_score * standard_error
 
# Calculating lower bound and upper bound
lower_bound <- mean_value - margin_error
upper_bound <- mean_value + margin_error
 
# Print the confidence interval
print(c(lower_bound,upper_bound))

embedded from [https://www.geeksforgeeks.org/how-to-find-confidence-intervals-in-r/][3] 

但是我应该如何像其他方法一样使用

aggregate

函数根据分组数据一次性计算置信区间？

我的主要问题是，在分组数据中，肯定会有多个

n

，而不是一个

n

。

我也尝试过使用

dplyr

工具，但输出不正确，似乎有些不对劲。

> tab.sep.file %>%
+ group_by(SDP,M346) %>%
+ summarise(mean = mean(tab.sep.file$Hic8),
+ sd = sd(tab.sep.file$Hic8),
+ n = n()) %>%
+ mutate(se = sd / sqrt(n),
+ lower.ci = mean - qt(1 - (0.05 / 2 ), n -1) * se,
+ upper.ci = mean + qt(1 - (0.05 / 2 ), n -1) * se)

Answer 1

您需要在

aggregate

中指定分组。例如通过公式。试试这个

> aggf <- \(x) {
+   n <- length(x)
+   xmn <- mean(x)
+   xsd <- sd(x)
+   se <- xsd/sqrt(n)
+   ci <- xmn + se*(-qt(p=.05/2, df=n - 1)) %*% cbind(-1, 1)
+   c(mean=xmn, sd=xsd, se=se, ci=ci)
+ }
> aggregate(Hic8 ~ sdp, dat, aggf)
   sdp   Hic8.mean     Hic8.sd     Hic8.se    Hic8.ci1    Hic8.ci2
1 SDP1  0.64358665  0.20498872  0.08368630  0.42846418  0.85870913
2 SDP2  0.32132755  0.16655345  0.09615968 -0.09241416  0.73506925
3 SDP3  0.46302347  0.37685995  0.21758019 -0.47314854  1.39919549

---

数据

> dput(dat)
structure(list(sdp = c("SDP1", "SDP1", "SDP1", "SDP2", "SDP2", 
"SDP2", "SDP1", "SDP1", "SDP1", "SDP3", "SDP3", "SDP3"), Hic8 = c(0.854164426913485, 
0.382015778450295, 0.799342160811648, 0.144684031140059, 0.475511683383957, 
0.343786930665374, 0.629821318900213, 0.412683063885197, 0.783493172377348, 
0.534960983321071, 0.798729537520558, 0.0553799034096301)), class = "data.frame", row.names = c(NA, 
-12L))