"Averaging" variances

Question

I need to obtain some sort of "average" among a list of variances, but have trouble coming up with a reasonable solution. There is an interesting discussion about the differences among the three Pythagorean means (arithmetic, geometric, and harmonic) in this thread; however, I still don't feel any of them would be a good candidate. Any suggestions?

P.S. Some context - These variances are sample variances from $n$ subjects, each of whom went through the same experiment design with roughly the same sample size $k$. In other words, there are $n$ sampling variances $\sigma_1^2$, $\sigma_2^2$, ..., $\sigma_n^2$, corresponding to those $n$ subjects. A meta analysis has been already performed at the population level. The reason I need to obtain some kind of "average" or "summarized" sample variance is that I want to use it to calculate an index such as ICC after the meta analysis.

P.P.S. To keep the discussion more concrete, let me explain the issue with the following example in R:

library(metafor)
dat <- get(data(dat.konstantopoulos2011))
dat$district <- as.factor(dat$district)
dat$school <- as.factor(dat$school)

In the dataset there is a variance associated with each school's performance score:

str(dat)
Classes ‘escalc’ and 'data.frame':  56 obs. of  6 variables:
 $ district: Factor w/ 11 levels "11","12","18",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ school  : Factor w/ 11 levels "1","2","3","4",..: 1 2 3 4 1 2 3 4 1 2 ...
 $ year    : int  1976 1976 1976 1976 1989 1989 1989 1989 1994 1994 ...
 $ yi      : atomic  -0.18 -0.22 0.23 -0.3 0.13 -0.26 0.19 0.32 0.45 0.38 ...
 $ vi      : num  0.118 0.118 0.144 0.144 0.014 0.014 0.015 0.024 0.023 0.043 ...

Suppose that we perform a meta analysis with a hierarchical or mixed-effects model:

$y_{ij} = a + \alpha_i + \beta_j + \epsilon_{ij}$

where $\alpha_i$ and $\beta_j$ are the random effects for the $i$th school and $j$th district, respectively, and $\epsilon_{ij}$ is the measurement error with a known Gaussian distribution $N(0,v_{ij})$. This model can be analyzed as below:

(fm <- rma.mv(yi, vi, random = list(~1 | district, ~1 | school), data=dat)) 

rendering the following variance estimates for the two variance components:

Multivariate Meta-Analysis Model (k = 56; method: REML)
Variance Components: 

            estim    sqrt  nlvls  fixed    factor
sigma^2.1  0.0814  0.2853     11     no  district
sigma^2.2  0.0010  0.0308     11     no    school

The two variances in the result, sigma^2.1 and sigma^2.2, correspond to the two random-effects variables (district and school).

I would like to compute the ICC for district, and that is why I wanted to obtain a summarized variance in the first place for those individual variances, $v_{ij}$, of the measurement term $\epsilon_{ij}$. Since the total variance is

$Var(y_{ij})= Var(\alpha_i + \beta_j + \epsilon_{ij}) = \sigma_1^2 + \sigma_2^2+v_{ij}$

my original (and simple) approach was to use just the arithmetic mean:

$\frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2+mean(v_{ij})}$

but I am not sure if the arithmetic mean, $mean(v_{ij})$, is appropriate in this context.

Answer 1

Expanding the comments that you got, answer for the question in your title is already given in How to 'sum' a standard deviation? thread, and reads as follows: to get average standard deviation, first take average of variances and then take square root of it.

At face value this approach is valid, but it ignores the hierarchical nature of your data. Similar example is discussed in chapter 5 of Bayesian Data Analysis by Andrew Gelman et al (see also here), who show that it is actually wiser to use hierarchical models that rely on pooled estimates. In your case you have $n \times k$ observations, for $n$ subjects in $k$ treatments and I guess it can be assumed that there is some kind of similarity between results obtained by each subject and between each treatment. This already suggests a hierarchical model with crossed upper-level effects for treatments and for subjects. By using such model you would account for both sources of variation.

Notice that modern formulations of ICC in fact define it in terms of mixed-effects models of the kind as described above, so employing such model solves multiple problems for you and it is often the recommended approach to meta-analysis (but notice that ICC can be misleading).

---

Regarding your edit, if your model is

$$ y_{ij} = a + \alpha_i + \beta_j + \epsilon_{ij} $$

then $\alpha_i \sim \mathcal{N}(\mu_\alpha, \sigma^2_\alpha)$, $\beta_j \sim \mathcal{N}(\mu_\beta, \sigma^2_\beta)$ and $\epsilon_{ij} \sim \mathcal{N}(0, \sigma^2_\epsilon)$, so your ICC is

$$ \mathrm{ICC}_\alpha = \frac{\sigma_\alpha^2}{\sigma_\alpha^2 + \sigma_\beta^2 + \sigma_\epsilon^2} $$

The mean of errors does not come into the equation at any point. What comes to the equation is variance of each of the random effects $\alpha,\beta$ and global "noise" $\epsilon$. The idea is to estimate the share of variance taken by $\alpha$, i.e. how much of the total variance does it account for. This is how ICC was defined by its creator Ronald A. Fisher (1966) in Statistical Methods for Research Workers:

(...) the intraclass correlation will be merely the fraction of the total variance due to that cause which observations in the same class have in common.

So the numerator in the ICC formula is the variance of the effect of interest and the denominator is the total variance. Notice that mean of variances has nothing to do with total variance (sum of variances), so unless I misunderstand something, I can't see why the mean is of your interest in here.