组合数据的标准差答案

【问题标题】：Standard deviation of combined data组合数据的标准差
【发布时间】：2023-03-08 14:59:01
【问题描述】：

我有一个包含平均值、标准差和 n 的数据集。其中一个变量的样本量相等，而另一个变量的样本量不同。

dat <- data.frame(variable = c(rep("x", 2), rep("y", 3)), replicate = c(1,2,1,2,3),
mean = c(3.4, 2.5, 6.5, 5.7, 5.1), sd = c(1.2, 0.7, 2.4, 4.0, 3.5),
n = c(3,3,5,4,6))

我需要组合 x 和 y 变量，并试图找到一种节省代码的方法来计算组合标准差，例如使用 aggregate 函数。 The equation for combined standard deviation 如下：

对于不等的样本量 (same source)：

我的组合数据框应如下所示：

variable    mean    sd
x           2.95    sd_x
y           5.76    sd_y

如何在 R 中创建一个计算组合标准差的函数？ 或者，如果有为此设计的包，它也算作答案 =)

【问题讨论】：

听起来你真的是在要求人们编写一个为你简单地实现该公式的函数。
@joran 很抱歉，如果是这样的话。那真的不是我的意思。我整个晚上都在挣扎并想问，因为没有R解决方案。如果我问的话，我认为这会对其他人有益。我首先提出了一个更长的问题，我解释了我所做的事情，但删除了它，因为它太长且难以理解。
我想有几种方法可以做到这一点。此处介绍的方法 (arxiv.org/ftp/arxiv/papers/1007/1007.1012.pdf) 给出的值与此处的弗洛德尔答案 (stackoverflow.com/questions/9222056/…) 相同。我正在查看的方法（上图）给出的值略有不同。我不知道为什么。由于缺乏知识，我会以Rudmin的解决方案来解决。

标签： r variance propagation standard-deviation

【解决方案1】：

Rudmin (2010) 指出合并数据集的精确方差是方差的平均值加上平均值的方差。 flodel has already provided an answer and function 给出了与 Rudmin 的声明相似的值。使用Rudmin的数据集和flodel's function基于Wikipedia：

df <- data.frame(mean = c(30.66667, 31.14286, 40.33333), variance = c(8.555555, 13.26531, 1.555555), n = c(6,7,3))

grand.sd   <- function(S, M, N) {sqrt(weighted.mean(S^2 + M^2, N) -
                                      weighted.mean(M, N)^2)}

grand.sd(sqrt(df$variance), df$mean, df$n)^2 

#[1] 22.83983 = Dp variance in Rudmin (2010).

但是，与function 5.38 from Headrick (2010) 相比，此解决方案给出的值略有不同（除非某处有错误）：

dat <- data.frame(variable = c(rep("x", 2), rep("y", 3)), replicate = c(1,2,1,2,3),
mean = c(3.4, 2.5, 6.5, 5.7, 5.1), sd = c(1.2, 0.7, 2.4, 4.0, 3.5),
n = c(3,3,5,4,6))

x <- subset(dat, variable == "x")

((x$n[1]^2)*(x$sd[1]^2)+
(x$n[2]^2)*(x$sd[2]^2)-
(x$n[2])*(x$sd[1]^2) -
(x$n[2])*(x$sd[2]^2) -
(x$n[1])*(x$sd[1]^2) -
(x$n[1])*(x$sd[2]^2) +
(x$n[1])*(x$n[2])*(x$sd[1]^2) +
(x$n[1])*(x$n[2])*(x$sd[2]^2) +
(x$n[1])*(x$n[2])*(x$mean[1] - x$mean[2])^2)/
((x$n[1] + x$n[2] - 1)*(x$n[1] + x$n[2]))

#[1] 1.015

grand.sd(x$sd, x$mean, x$n)^2

#[1] 1.1675

为了回答我自己的问题，所需的data.frame 将通过以下方式获得：

library(plyr)
ddply(dat, c("variable"), function(dat) c(mean=with(dat,weighted.mean(mean, n)),  sd = with(dat, grand.sd(sd, mean, n))))   

  variable     mean       sd
1        x 2.950000 1.080509
2        y 5.726667 3.382793

【讨论】：

【解决方案2】：

使用`utilities`包中的`sample.decomp`函数

此类统计问题在utilities package 的sample.decomp 函数中自动进行。该函数可以从子组矩中计算池化样本矩，或者从其他子组矩和池化矩中计算缺失的子组矩。它适用于高达四阶的分解——即样本大小、样本均值、样本方差/标准差、样本偏度和样本峰度的分解。

如何使用该函数：这里我们将展示如何为您的数据集实现该函数。

#Input sample statistics for subgroups
SIZE <- c(3, 3, 5, 4, 6)
MEAN <- c(3.4, 2.5, 6.5, 5.7, 5.1)
SD   <- c(1.2, 0.7, 2.4, 4.0, 3.5)

#Compute sample decomposition
library(utilities)
sample.decomp(n = SIZE, sample.mean = MEAN, sample.sd = SD, include.sd = TRUE)

            n sample.mean sample.sd sample.var
1           3    3.400000  1.200000   1.440000
2           3    2.500000  0.700000   0.490000
3           5    6.500000  2.400000   5.760000
4           4    5.700000  4.000000  16.000000
5           6    5.100000  3.500000  12.250000
--pooled-- 21    4.933333  2.964428   8.787833

此输出为您提供合并样本量、样本均值和样本标准差（或等效的样本方差）。

【讨论】：

使用utilities包中的sample.decomp函数

使用`utilities`包中的`sample.decomp`函数