如何计算加权标准误差并将它们绘制在条形图中？答案

【问题标题】：How can I calculate weighted standard errors and plot them in a bar plot?如何计算加权标准误差并将它们绘制在条形图中？
【发布时间】：2020-05-16 03:05:12
【问题描述】：

我有一个计数数据框。我想计算加权比例，绘制比例，并为这些加权比例绘制标准误差线。

我的数据框示例：

head(df[1:4,])
  badge year total b_1 b_2 b_3 b_4 b_5 b_6 b_7 b_8 b_9 b_10
1    15 2014    14   3   2   1   1   1   1   1   1   1    1
2    15 2015   157  13  12  11   8   6   6   6   5   5    5
3    15 2016    15   5   3   1   1   1   1   1   1   1    0
4  2581 2014    13   1   1   1   1   1   1   1   1   1    1

数据包含在给定年份中，警察在十个不同的警察节拍（b_1、b_2、...）中响应的 911 次呼叫的计数。因此，15 号警官在 2014 年总共响应了 14 个电话，其中 3 个在第 1 拍中，2 在第 2 拍中，依此类推。

基本上，我想要的是获得每个节拍内发生的呼叫的总体比例。但我希望这些比例由电话总数加权。

到目前为止，我已经能够通过将每个 b_ 列和总列中的值相加并计算比例来计算这一点。我在一个简单的条形图中绘制了这些。我一直无法弄清楚如何计算按总计加权的标准误差。

我对如何绘制数据没有偏好。我主要对获得正确的标准错误感兴趣。

这是我目前的代码：

sums_by_beat <- apply(df[, grep('b_', colnames(df2))], 2, sum)
props_by_beat <- sums_by_beat / sum(df$total)
# Bar plot of proportions by beat
barplot(props_by_beat, main='Distribution of Calls by Beat', 
        xlab="Nth Most Common Division", ylim=c(0,1), 
        names.arg=1:length(props_by_beat), ylab="Percent of Total Calls")

还有我的数据的 30 行样本：

df <- structure(list(badge = c(15, 15, 15, 2581, 2581, 2745, 2745, 
3162, 3162, 3162, 3396, 3650, 3650, 3688, 3688, 3688, 3698, 3698, 
3698, 3717, 3717, 3717, 3740, 3740, 3740, 3813, 3873, 3907, 3930, 
4007), year = c(2014, 2015, 2016, 2014, 2015, 2015, 2016, 2014, 
2015, 2016, 2016, 2014, 2015, 2014, 2015, 2016, 2014, 2015, 2016, 
2014, 2015, 2016, 2014, 2015, 2016, 2016, 2015, 2014, 2014, 2014
), total = c(14, 157, 15, 13, 29, 1, 1, 754, 1172, 1039, 14, 
1, 2, 34, 57, 146, 3, 7, 28, 593, 1036, 1303, 461, 952, 1370, 
1, 4, 41, 5, 451), b_1 = c(3, 13, 5, 1, 3, 1, 1, 33, 84, 83, 
2, 1, 2, 5, 10, 14, 2, 7, 7, 39, 72, 75, 42, 69, 81, 1, 1, 7, 
1, 36), b_2 = c(2, 12, 3, 1, 2, 0, 0, 33, 61, 52, 2, 0, 0, 3, 
6, 8, 1, 0, 2, 37, 65, 70, 29, 65, 75, 0, 1, 5, 1, 23), b_3 = c(1, 
11, 1, 1, 2, 0, 0, 32, 57, 45, 2, 0, 0, 3, 5, 8, 0, 0, 2, 34, 
62, 67, 28, 50, 73, 0, 1, 3, 1, 22), b_4 = c(1, 8, 1, 1, 2, 0, 
0, 31, 44, 39, 2, 0, 0, 3, 3, 7, 0, 0, 2, 34, 61, 67, 26, 42, 
72, 0, 1, 3, 1, 21), b_5 = c(1, 6, 1, 1, 1, 0, 0, 30, 42, 37, 
1, 0, 0, 3, 3, 7, 0, 0, 1, 33, 53, 61, 23, 42, 67, 0, 0, 2, 1, 
21), b_6 = c(1, 6, 1, 1, 1, 0, 0, 30, 40, 36, 1, 0, 0, 2, 2, 
6, 0, 0, 1, 32, 53, 61, 22, 41, 63, 0, 0, 2, 0, 21), b_7 = c(1, 
6, 1, 1, 1, 0, 0, 26, 39, 35, 1, 0, 0, 2, 2, 6, 0, 0, 1, 30, 
47, 58, 22, 39, 62, 0, 0, 2, 0, 21), b_8 = c(1, 5, 1, 1, 1, 0, 
0, 26, 39, 33, 1, 0, 0, 2, 2, 6, 0, 0, 1, 30, 47, 58, 21, 38, 
59, 0, 0, 2, 0, 19), b_9 = c(1, 5, 1, 1, 1, 0, 0, 24, 34, 33, 
1, 0, 0, 2, 2, 5, 0, 0, 1, 30, 43, 57, 20, 37, 57, 0, 0, 2, 0, 
15), b_10 = c(1, 5, 0, 1, 1, 0, 0, 23, 34, 32, 1, 0, 0, 1, 2, 
5, 0, 0, 1, 27, 40, 56, 18, 36, 55, 0, 0, 2, 0, 14)), row.names = c(NA, 
30L), class = "data.frame")

【问题讨论】：

标签： r bar-chart

【解决方案1】：

没有（据我所知）内置 R 函数来计算加权平均值的标准误差，但计算起来相当简单 - 有一些假设在您描述的情况下可能有效.
参见，例如： https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Standard_error

加权平均值的标准误差

如果用于计算加权均值的元素是来自具有相同方差 v 的总体样本，则加权样本均值的方差估计为：

var_m = v^2 * sum( wnorm^2 )   # wnorm = weights normalized to sum to 1

并且加权均值的标准误等于方差的平方根。

sem = sqrt( var_m )

因此，我们需要根据加权数据计算样本方差。

加权方差

加权总体方差（或有偏的样本方差）计算如下：

pop_v = sum( w * (x-mean)^2 ) / sum( w )

但是，如果（如您所描述的情况），我们正在处理从总体中提取的样本，而不是总体本身，我们需要进行调整以获得无偏样本方差。
如果权重代表用于计算加权均值和方差的每个元素的观测频率，则调整为：

v = pop_v * sum( w ) / ( sum( w ) -1 )

但是，这里的情况并非如此，因为权重是每个警察的 911 呼叫的总频率，而不是每个节拍的呼叫。所以在这种情况下，权重对应于每个元素的reliabilities，调整为：

v = pop_v * sum( w )^2 / ( sum( w )^2 - sum( w^2) )

weighted.var 和 weighted.sem 函数

将所有这些放在一起，我们可以定义 weighted.var 和 weighted.sem 函数，类似于基本 R 函数 weighted.mean （请注意，几个 R 包，例如“Hmisc”，已经包含更多通用函数来计算加权方差）：

weighted.var = function(x,w,type="reliability") {
    m=weighted.mean(x,w)
    if(type=="frequency"){ return( sum(w*(x-m)^2)/(sum(w)-1) ) }
    else { return( sum(w*(x-m)^2)*sum(w)/(sum(w)^2-sum(w^2)) ) }
}
weighted.sem = function(x,w,...) { return( sqrt(weighted.var(x,w,...)*sum(w^2)/sum(w)^2) ) }

应用于问题中的 911 通话数据

在问题的情况下，我们要计算加权平均值和加权 sem 的元素对应于每个警察在每个节拍中的呼叫比例。
所以（终于...）：

props = t(apply(df,1,function(row) row[-(1:3)]/row[3]))
wmean_props = apply(props,2,function(col) weighted.mean(col,w=df[,3]))
wsem_props = apply(props,2,function(col) weighted.sem(col,w=df[,3]))

【讨论】：

【解决方案2】：

您的“比例”实际上不是加权（total）观察值的平均值吗？然后我们可以简单地计算出相应的加权colMeans。

df2 <- df[, grep('b_', colnames(df))]

means.w <- colMeans(df2 / df$total)

对于误差线，我们可以使用 1 - alpha/2 的 quantiles，即对于 alpha==.05，我们使用 c(.025, .975)。分析 sds 会产生负值。

q.w <- t(apply(df2 / df$total, 2, quantile, c(.025, .975)))

现在，我们存储 barplot 返回不可见的 x 位置，

# Bar plot of proportions by beat
b <- barplot(means.w, main='Distribution of Calls by Beat', 
             xlab="Nth Most Common Division", ylim=c(0,1), 
             names.arg=1:length(means.w), ylab="Percent of Total Calls")

并使用arrows 构建误差线。

arrows(b, q.w[,1], b, q.w[,2], length=.02, angle=90, code=3)

【讨论】：