如何计算 R 中的合并标准差？

Question

我想计算我的数据框中所有唯一站点的合并（实际加权）标准差。

这些站点的值是单一物种林分的值，我想合并平均值和标准差，以便我可以比较阔叶林和针叶林。
这是具有阔叶林值的数据框 (df)：

keybl           n   mean    sd
Vest02DenmDesp  3   58.16   6.16
Vest02DenmDesp  5   54.45   7.85
Vest02DenmDesp  3   51.34   1.71
Vest02DenmDesp  3   59.57   5.11
Vest02DenmDesp  5   62.89   10.26
Vest02DenmDesp  3   77.33   2.14
Mato10GermDesp  4   41.89   12.6
Mato10GermDesp  4   11.92   1.8
Wawa07ChinDesp  18  0.097   0.004
Chen12ChinDesp  3   41.93   1.12
Hans11SwedDesp  2   1406.2  679.46
Hans11SwedDesp  2   1156.2  464.07
Hans11SwedDesp  2   4945.3  364.58

Keybl 是网站的代码。合并 SD 的公式为：

s=sqrt((n1-1)*s1^2+(n2-1)*s2^2)/(n1+n2-2))

（抱歉我不能发图片，也没有找到可以直接进入公式的链接）

其中 2 是组数，因此会因站点而异。我知道这用于 t 检验，并且要比较两组。在这种情况下，我不打算比较这些组。我的教授建议我使用这个公式来获得加权标准差。我没有找到以我需要的方式包含此公式的 R 函数，因此我尝试构建自己的。然而，我是 R 新手，不太擅长制作函数和循环，因此希望能得到您的帮助。

这是我目前得到的：

sd=function (data) {
nc1=data[z,"nc"]
sc1=data[z, "sc"]
nc2=data[z+1, "nc"]
sc2=data[z+1, "sc"]
sd1=(nc1-1)*sc1^2 + (nc2-1)*sc2^2
sd2=sd1/(nc1+nc2-length(nc1))
sqrt(sd2)
}

splitdf=split(df, with(df, df$keybl), drop = TRUE)

for (c in 1:length(splitdf)) {
for (i in 1:length(splitdf[[i]])) {
    a = (splitdf[[i]])
    b =sd(a)
    }
}

1) 函数本身不正确，因为它给出的值略低于应有的值，我不明白为什么。莫非是z+1到达最后一行的时候没有停止？如果是这样，如何纠正？

2) 循环完全错误，但这是我在几个小时没有成功后能想到的。

有人可以帮助我吗？

谢谢，

安特拉

Answer 1

您正在尝试做的事情将受益于更通用的公式，这将使其更容易。如果您不需要通过 keybl 变量将其分成几部分，那么您就完成了。

dd <- df #df is not a good name for a data.frame variable since df has a meaning in statistics

dd$df <- dd$n-1
pooledSD <- sqrt( sum(dd$sd^2 * dd$df) / sum(dd$df) )
# note, in this case I only pre-calculated df because I'll need it more than once. The sum of squares, variance, etc. are only used once.

R 中一个重要的一般原则是尽可能多地使用向量数学。在这种微不足道的情况下，它并不重要，但为了了解如何在计算速度更重要的大型 data.frame 对象上执行此操作，请继续阅读。

# First use R's vector facilities to define the variables you need for pooling.
dd$df <- dd$n-1
dd$s2 <- dd$sd^2 # sd isn't a good name for standard deviation variable even in a data.frame just because it's a bad habit to have... it's already a function and standard deviations have a standard name
dd$ss <- dd$s2 * dd$df

现在只需使用便捷函数来拆分和计算必要的总和。注意这里每个隐式循环只执行一个函数（*apply、aggregate等都是隐式循环多次执行函数）。

ds <- aggregate(ss ~ keybl, data = dd, sum)
ds$df <- tapply(dd$df, dd$keybl, sum) #two different built in methods for split apply, we could use aggregate for both if we wanted
# divide your ss by your df and voila
ds$s2 <- ds$ss / ds$df
# and also you can easly get your sd
ds$s <- sqrt(ds$s2)

而正确答案是：

           keybl           ss df           s2          s
1 Chen12ChinDesp 2.508800e+00  2 1.254400e+00   1.120000
2 Hans11SwedDesp 8.099454e+05  3 2.699818e+05 519.597740
3 Mato10GermDesp 4.860000e+02  6 8.100000e+01   9.000000
4 Vest02DenmDesp 8.106832e+02 16 5.066770e+01   7.118125
5 Wawa07ChinDesp 2.720000e-04 17 1.600000e-05   0.004000

这看起来比其他方法（如 42- 的答案）要简洁得多，但是如果您根据实际执行的 R 命令的数量展开这些方法，这会更加简洁。对于这样一个简短的问题，任何一种方式都可以，但我想我会向您展示使用最多矢量数学的方法。它还强调了为什么可以使用那些方便的隐式循环函数，以提高表现力。如果您使用for 循环来完成相同的操作，那么将所有内容都放入循环中的诱惑会更大。这在 R 中可能是个坏主意。

Answer 2

独立假设下的池化 SD（因此可以假设协方差项为零）将为：sqrt( sum_over_groups[ (var)/sum(n)-N_groups)] )

     lapply( split(dat, dat$keybl), 
          function(dd) sqrt( sum( dd$sd^2 * (dd$n-1) )/(sum(dd$n-1)-nrow(dd)) ) )
#-------------------------
$Chen12ChinDesp
[1] 1.583919

$Hans11SwedDesp
[1] Inf

$Mato10GermDesp
[1] 11.0227

$Vest02DenmDesp
[1] 9.003795

$Wawa07ChinDesp
[1] 0.004123106