SE of weighted mean

Question

$X$ is a random variable with unknown distribution. A number of experiments are conducted to estimate $X$. Each experiment has a different reliability measure in estimating $X$. These $n$ experiments resulted in following sample set $\{x_1, x_2, x_3, ... , x_n\}$ with corresponding non-zero weights being $\{w_1, w_2, w_3, ... , w_n\}$. The higher weight corresponds to higher reliability. Note ${\sum_{i=1}^n{w_i}}$ can be greater than $1$.

The best unbiased estimator of true value of $X$ is the weighted mean of sample,
$\hat{X} = \bar{x}_w$,
where, $\bar{x}_w = \frac{\sum_{i=1}^n{w_ix_i}}{\sum_{i=1}^n{w_i}}$

The estimator for variance of $X$ from its true mean is,
$\hat{\sigma^2} = \bar{\sigma^2}_w$,
where, $\bar{\sigma^2}_w = \frac{\sum_{i=1}^n {w_i(x_i-\bar{x}_w)^2}}{\sum_{i=1}^n{w_i}}$,

What would be the best estimate of Standard Error of the sampling distribution of $\bar{x}_w$. Would it be $\frac{\bar{\sigma}_w}{\sqrt{n}}$. If yes, can someone help derive/explain it.

Answer 1

Well, under i.i.d. assumptions for the $x_i$ the variance of the weighted sample mean relative to the population variance would have a factor $\left(\sum_iw_i^2\right)/\left(\sum_iw_i\right)^2$, which reduces to $1/n$ if all of the weights are equal. And there would be a comparable substitution for the "$1-1/n$" factor on the unbiased sample variance as well. See here for example.

---

For explanation, we have $\bar{x}=\frac{1}{W}\sum_iw_ix_i$, where $W=\sum_iw_i$.

So then using $\mathbb{V}[\,]$ to denote variance, we have \begin{align} \mathbb{V}[\bar{x}] &= \mathbb{V}\left[\tfrac{1}{W}\sum_iw_ix_i\right] \\ &=\sum_i\mathbb{V}\left[\frac{\,w_i\,}{W}\,x_i\right] \\ &=\sum_i\left[\left(\frac{\,w_i}{W}\right)^2\,\mathbb{V}[x_i]\right] \\ &=\mathbb{V}[x]\frac{\sum_iw_i^2}{W^2} \\ \end{align}

Answer 2

I know this is an old one but I keep encountring the same question over and over. Apparently, there is no consensus as to the definition of the standard error of the weighted mean. Even different statistical softwares use different definitions. However, the most coherent answer that I keep seeing is this for an unbiased estimation of the standard error on a weighted mean:

$$ se= \frac{\sigma_w}{\sqrt{\sum_i^n w_i}} $$

where the $\sigma_w$ is the unbiased estimator of the standard deviation of your random variable $X$ and $\sum_i^n w_i$ is the sum of the individual weights that contribute to your unbiased estimation of $X$. The unbiased estimator of the standard deviation of your random variable with degres of freedom $=1$ is the following:

$$ \sigma_w = \sqrt{\frac{\sum_i^n w_i x_i}{\sum_i^n w_i - 1}} $$

Here's a link to a note that compares how it is computed in SPSS vs WinCross.

Python's statsmodels implemented a class that computes all sorts of weighted statistics including the standard deviation and standard error (method under the name std_mean here in their source code. As we can see from their implementation, they use either a biased estimation of the standard error if the degres of freedom is equal to $0$ like so:

$$ se= \frac{\sigma_w}{\sqrt{\sum_i^n w_i-1}} $$

or an unbiased estimator of the standard error (which is your case) if the degres of freedom parameter is given which activates a condition to apply a degres of freedom correction to the standard deviation first like so:

$$ \sigma_w \leftarrow \sigma_w \times \sqrt{\frac{\sum_i^n w_i- ddof}{\sum_i^n w_i}} $$

For $ddof=1$ and if you plugin the new value of the corrected $\sigma_w$ in the biased estimation of the standard error, then you get the formula for the unbiased estimation of the standard error of the weighted mean $se = \sigma_w/\sqrt{\sum_i^n w_i}$

Here's how to numerically verify your estimators using manual definitions vs statsmodels's implementation if you use python:

# make sure you install statsmodels using pip install statsmodels
import numpy as np
from statsmodels.stats.weightstats import DescrStatsW

# define the x measurements and their weights
x = np.array([10, 12, 15.2, 12.5, 11])
w = np.array([100, 120, 108, 80, 98])

# calculate the unbiased estimators of avg, std and se (with ddof=1)
sum_w = np.sum(w)
avg_w = np.sum(w * x) / sum_w
std_w = np.sqrt(np.sum(w*(x-avg_w)**2) / (sum_w-1))
se_w = std_w / np.sqrt(sum_w)

# calculate the weighted stats using scipy's formula (with ddof=1) 
weighted_stats = DescrStatsW(x, weights=w, ddof=1)

print('manual weighted avg = %0.5f' %avg_w)
print('manual weighted std = %0.5f' %std_w)
print('manual weighted se = %0.5f' %se_w)
print('statsmodels weighted avg = %0.5f' %weighted_stats.mean)
print('statsmodels weighted std = %0.5f' %weighted_stats.std)
print('statsmodels weighted se = %0.5f' %weighted_stats.std_mean)

>>> OUTPUT:
manual weighted avg = 12.17312
manual weighted std = 1.78484
manual weighted se = 0.07935
statsmodels weighted avg = 12.17312
statsmodels weighted std = 1.78484
statsmodels weighted se = 0.07935