如何计算R中的汇总标准差？

Question

我想计算数据框中所有唯一站点的合并(实际加权)标准偏差.

这些站点的值是单一物种林分的值,我想汇集平均值和sd,以便我可以比较阔叶林分和针叶林林分.
这是具有阔叶林分的值的数据框(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检验,两个人想要比较.在这种情况下,我不打算比较这些组.我的教授建议我用这个公式得到一个加权的sd.我找不到以我需要的方式包含这个公式的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)循环是完全错误的,但这是我几个小时没有成功后想出来的.

有谁能够帮我？

谢谢,

Antra

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