# calculating asymmetric standard deviation of a beta distributed variable

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?

• The distribution is asymmetric; the variance is not - it's just related to spread, and is agnostic with respect to skew. Perhaps you are interested in calculating something like an upper 95% quantile? – jbowman Aug 5 '17 at 15:23
• @jbowman I just want to deal with the fact that the mean of a bunch of values between 0.8 and 1 may be 0.99, but a straight forward calculation of the standard deviation may return a value of greater than 0.01, resulting in a range of uncertainty that incldues values greater than 1, which are not possible. – colin Aug 5 '17 at 15:28
• Are you looking for a description of the data or estimates of the accuracy of the mean? – jbowman Aug 5 '17 at 16:08
• @jbowman is an estimate of the accuracy of the mean not a description of the data? – colin Aug 5 '17 at 16:17
• Not really, For example, I can have a Normal distribution with mean 0, standard deviation 1, and with a sample size of 100 the s.d. of the sample mean is 0.1, whereas with a sample size of 10,000 the s.d. of the sample mean is 0.01. Neither, independent of the sample size, tells me anything about the s.d. of the underlying distribution. – jbowman Aug 5 '17 at 16:19

## 1 Answer

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.

• Thanks for the detailed response. However, I really need an estimate of the uncertainty around the mean. This will only ever provide information regarding the distribution of the data with respect to the median. – colin Aug 5 '17 at 22:42
• It provides estimation about how likely a data point is to be above whatever level you pick. It has nothing to do with the median. As @maarten buis noted above, the requirements of your question seem inconsistent; you want an estimate of the std. deviation, but you don't want it large enough so that the mean +/- your estimate is outside the interval $(0,1)$... doesn't make sense, statistically, at least not in the context of a beta distribution. You appear, as Maarten pointed out, to really be wanting a low and high quantile, not an estimate of the sd. – jbowman Aug 6 '17 at 5:06