calculating asymmetric standard deviation of a beta distributed variable

Question

I have data on the relative abundance of a species within a population that varies continuously on the interval [0,1]. I have divided the data into 5 categories, and would like to take estimate the mean and some level of uncertainty of the dependent variable for each category. Here is some example data in R:

#load betareg package, grab some data.
library(betareg)
data("GasolineYield", package = "betareg")

#Divide the yield variable into 5 categories.
GasolineYield$category <- cut(GasolineYield$yield, 5, labels=F)

#calculate means of the 5 categories.
aggregate(yield ~ category, data = GasolineYield, FUN = 'mean')
  category      yield
1        1 0.02213433
2        2 0.02079663
3        3 0.02013171
4        4 0.01792205
5        5         NA

Taking the mean is straightforward. However, I want to calculate a standard deviation or some other metric of uncertainty, accounting for the fact that the variance of a beta distributed dependent variable is asymmetric. How can I do this?

Answer 1

What I would do in your situation is to calculate quantiles, e.g., 90th or 95th percentile, instead of the standard deviation. If you have enough data, e.g., several hundred observations, I'd probably use a bootstrap. If you really believe your data to be well-approximated by a beta distribution, I'd use that instead.

I'll provide an example, based upon the GasolineYield dataset, but constructed slightly differently. I group the yields into three groups based upon low, medium, and high temperature, then use method-of-moments estimators to estimate the parameters of the three beta distributions. (Note that MOM estimators can fail, especially in small sample sizes, and that there is some evidence that for small sample sizes the MLE estimator is better than the MOM estimator.)

I use the data.table package for conciseness, but it should be clear what I am doing:

library(data.table)

GasolineYield <- data.table(GasolineYield)
GasolineYield[, category := cut(temp10, 3, labels=FALSE)]

# calculate sample mean, variance for each category
estimates <- GasolineYield[, .(N = .N, xbar = mean(yield), s2 = var(yield)), category]

# calculate MOM parameter estimates
estimates[, ':='(a=xbar*(xbar*(1-xbar)/s2-1), b=(1-xbar)*(xbar*(1-xbar)/s2-1))]

# calculate 75th and 90th percentiles of the estimated Beta dist'ns
estimates[, ':='(u75 = qbeta(0.75, a, b), u90 = qbeta(0.9, a, b))]

# Display results
 > estimates[, .(category, N, xbar, u75, u90)]
   category  N      xbar       u75       u90
1:        1 17 0.2214118 0.2920784 0.3789146
2:        2 10 0.1743000 0.2392397 0.3306011
3:        3  5 0.1568000 0.1904971 0.2292984

The xbar column contains the estimate for the category, and the u75 and u90 columns contain the estimated upper 75th and 90th percentiles of the data. In your case the lower quantiles would seem to be what you're more interested in:

> estimates[, .(category, N, xbar, u25 = qbeta(0.25, a, b), u10 = qbeta(0.1, a, b))]
   category  N      xbar        u25        u10
1:        1 17 0.2214118 0.13485425 0.08562626
2:        2 10 0.1743000 0.08815858 0.04789428
3:        3  5 0.1568000 0.11752094 0.09117189

No matter how low you go (the 0.001th percentile, for example), the estimates won't go below 0:

> estimates[, .(category, N, xbar, qbeta(0.00001, a, b), qbeta(0.99999, a, b))]
   category  N      xbar           V4        V5
1:        1 17 0.2214118 0.0023495621 0.7883513
2:        2 10 0.1743000 0.0002724701 0.7840161
3:        3  5 0.1568000 0.0170810810 0.4544612

... and you can see that there's no symmetry around the point estimate xbar. This comes about because you are working directly with the quantiles, which, of course, take asymmetry into account.