How to iteratively calculate weighted standard error to report alongside a weighted mean

Question

I have a group of individuals for which I would like to report a mean and a weighted error. The data that I observe on a daily basis are two independent $iid$ random variables with unknown distributions that can be quantified as 'number of individuals per day' which i call $W$ and some measure of 'revenue' called $R$ that represents the sum of the revenue brought by all the individuals of that specific day.

Each day, those two variables are updated as they are observed. Using this data, I can compute the thing that I am intereseted in reporting which is the 'revenue per individual per day' called $X$ by computing $x_i = r_i/w_i$ on each day ($i.e.$ for each day $i$, I divide the sum of revenue for all the individuals by the number of individuals). This leaves me with the following dataframe:

W	R	X
12	100	8.33
15	110	7.33
65	740	11.38

Using this data, I would like to finally report my KPI which is the mean revenue per individual per day over this data collection period. Since not all days are equally represented with the same number of individuals bringing revenue ($e.g.$ day $3$ in the data above has a much higher individual count of $65$ compared to day $1$ which contains only $12$), I would like to represent my KPI as a weighted mean where the individual weights are the number of individuals $w_i$ such that days with higher individual count are better represented in the final mean I report like so:

$$ \bar{x}_w = \frac{\sum_i w_i \times x_i}{\sum_i w_i} $$

Then, I can calculate the weighted sample standard deviation like so (according to this note):

$$ s_w = \sqrt{\frac{ \sum_iw_i (x_i - \bar{x}_w)^2}{ \sum_i w_i -1}} $$

Finally, I can compute the standard error of the weighted mean which is defined as $se = s_w/ \sqrt{n}$.

First question

What is $n$ in the standard error calculation in my case? Is it the number of observations which is $3$ in the example dataframe that I have outlined above? Or is it the sum of the weights used to compute the weighted standard deviation $\sum_i w_i$ which is equal to $92$? In case we use the sum of the weights, the error may get quite small which is something I do not typically observe since my data has high variability. In case it is the number of observations, the error will systematically decrease when I gather more data which also seems unusual. So how may i correctly estimate this empirically with the data I have? Is there any documentation anyone can point me to for reference?

Second question

My second question is somewhat trickier. In a production environment, when I report the KPI alongside its standard error, I do not have access to the full dataframe. I only have access to the first two data points to start the calculation procedure of the weighted mean and weighted standard error deviation (and the standard error).

Given that the current observation count is $n$, I would like to iteratively calculate the weighted mean and weighted standard deviation and the standard error of the weighted mean in computer code knowing only the previous weighted mean $\bar{x}_{w,n}$, previous weighted standard deviation $s_{w,n}$, the cumulative sum of previous weights $\sum_i^n w_i$, and of course the new row of data points that contains the new observation $x_{n+1}$ and its new weight $w_{n+1}$ where $n+1$ indicates the new row.

Thank you for the help!

Answer 1

Interestingly, there is no single equation for a weighted standard error. Multiple versions have been proposed in the literature though. See for example:

Donald F. Gatz and Luther Smith (1995). "The Standard Error of a Weighted Mean Concentration - I: Bootstrapping Vs Other Methods". In: Atmospheric Environment 29.11, pp. 1185-1193

I implemented some of these method in R to use as an (unexported) function in the adjustedCurves R-package I developed, here is the code:

weighted.se <- function(x, w, se_method, na.rm=FALSE) {

  if (na.rm) {
    miss_ind <- !is.na(x)
    w <- w[miss_ind]
    x <- x[miss_ind]
  }

  n <- length(x)
  mean_Xw <- stats::weighted.mean(x=x, w=w, na.rm=na.rm)

  ## Miller (1977)
  if (se_method=="miller") {
    se <- 1/n * (1/sum(w)) * sum(w * (x - mean_Xw)^2)
  ## Galloway et al. (1984)
  } else if (se_method=="galloway") {
    se <- (n/(sum(w)^2)) * ((n*sum(w^2 * x^2) - sum(w*x)^2) / (n*(n-1)))
  ## Cochrane (1977)
  } else if (se_method=="cochrane") {
    mean_W <- mean(w)
    se <- (n/((n-1)*sum(w)^2))*(sum((w*x - mean_W*mean_Xw)^2)
                                - 2*mean_Xw*sum((w-mean_W)*(w*x-mean_W*mean_Xw))
                                + mean_Xw^2*sum((w-mean_W)^2))
  ## As implemented in Hmisc
  } else if (se_method=="Hmisc") {
    se <- (sum(w * (x - mean_Xw)^2) / (sum(w) - 1)) / n
  }
  return(sqrt(se))
}

where x is your vector of interest, w is a vector of weights with equal length, se_method specifies which method to use and na.rm specifies whether to remove missing values before performing calculations.

I understand that this doesn't fully answer your questions, but it might still be helpful to you.

Answer 2

I have kind of figured out a (almost) correct answer to my question so I will post it here and leave room for others to weigh in to improve it.

Answer to the first question

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{s_w}{\sqrt{\sum_i^n w_i}} $$

where the $s_w$ is the unbiased estimator of the standard deviation of the random variable $X$ and $\sum_i^n w_i$ is the sum of the individual weights that contribute to the unbiased estimation of $X$. The following link is a statistical note that compares how it is computed in SPSS vs WinCross and SPSS uses the sum of weights as the denominator (which happens to be almost the same as the sample size $n$ in their example). So in the example I provided in my question, the sum of weights is $\sum_i^n w_i = 92$.

Answer to the second question

I came up with the following formulas for recursive computation of the weighted mean, weighted standard deviation and the standard error on the weighted mean:

Given that the current known data points are $n$ and the next data point that triggers the update is denoted as $n+1$, we can express the weighted stats like so:

Recursive weighted mean:

$$ \bar{x}_{w,n+1} = \frac{(\sum_{i=1}^{n} w_i) \bar{x}_{w,n} + w_{n+1} \times x_{n+1}}{\sum_{i=1}^{n} w_i + w_{n+1}} $$

Recusrive weighted standard deviation

$$ s_{w,n+1} = \sqrt{\frac{(\sum_{i=1}^{n} w_i) \times (s_{w,n}^2 + [\bar{x}_{w,n} - \bar{x}_{w,n+1}]^2) + w_{n+1} (x_{n+1} - \bar{x}_{w,n+1})^2}{\sum_{i=1}^{n} w_i + w_{n+1}}} $$

Standard error of the weighted mean

$$ se_w = \frac{s_{w,n+1}}{\sqrt{\sum_{i=1}^{n} w_i + w_{n+1}}} $$

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, their unbiased estimator of the standard error with degres of freedom parameter set to $1$ is the formula that I wrote above. This answers my first question as to what should I take as a denominator when computing the unbiased estimation of the standard error on my weighted mean.

Using Python I was able to verify the implementation of the above estimators using recursive definitions vs statsmodels's weighted stats function knowing the full history of data like so:

import numpy as np
from statsmodels.stats.weightstats import DescrStatsW

def update_weighted_mean_se(current_sum_weights, current_weighted_avg, current_weighted_std, new_weight, new_x):
    ''' 
    Update the weighted statistics (mean, weighted standard deviation and weighted standard error) given the previous 
    sum of weights, previous weighted mean, previous weighted standard deviation, new weight, and new x value. 
    '''
    # new weighted mean and weighted standard deviation recursively 
    new_sum_weights = current_sum_weights + new_weight
    new_weighted_avg = (current_sum_weights*current_weighted_avg + new_weight*new_x) / new_sum_weights
    new_weighted_std = np.sqrt((current_sum_weights*(current_weighted_std**2 + (current_weighted_avg-new_weighted_avg)**2) + new_weight*(new_x-new_weighted_avg)**2)/new_sum_weights)

    # new standard error on the weighted mean
    se_w = new_weighted_std / np.sqrt(new_sum_weights)
    return new_weighted_avg, new_weighted_std, se_w

# 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)

# add new values and compute weighted stats iteratively
new_x_array = np.array([20, 30])
new_weights_array = np.array([200, 150])
for new_x, new_w in zip(new_x_array, new_weights_array):
    avg_w, std_w, se_w = update_weighted_mean_se(sum_w, avg_w, std_w, new_w, new_x)
    sum_w+=new_w

# verify new weighted stats using the formula (with ddof=1) 
weighted_stats = DescrStatsW(np.concatenate([x, new_x_array]), weights=np.concatenate([w, new_weights_array]), ddof=1)

print('iterative weighted avg = %0.5f' %avg_w)
print('iterative weighted std = %0.5f' %std_w)
print('iterative 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:
iterative weighted avg = 17.12570
iterative weighted std = 6.88164
iterative weighted se = 0.23521
statsmodels weighted avg = 17.12570
statsmodels weighted std = 6.88539
statsmodels weighted se = 0.23534

My implementation yields the correct weighted average but the standard deviation (and by extention the standard error) are only accurate up to $2$ or $3$ decimal points. This means that my implementation of the standard deviation is not exactly the same as statsmodel's and there is room for improvement. I wonder if this is just a matter of numerical precision.