Estimating error when bootstrapping results are skewed

Question

I'm trying to figure out the most robust way to estimate the error on the best-fit parameters of an exponential model ($y = x^a$) for some data. I have heard that bootstrapping is a solid way to attempt this and so have fit models to 1000+ random samples of the data (generated with replacement) and saved the best-fit $a$ for each sample.

The problem is that the distribution of these bootstrapped values looks extremely skewed and is nowhere close to Gaussian. I know that with sample statistics like mean, standard error, etc the Central Limit Theorem would imply that the distribution of results approximates a Gaussian as the sample size increases. I don't think this would apply to model parameters, so I suppose it is not super surprising that the bootstrapped distribution is not.

I'm not sure how to express the uncertainty on my best-fit estimate of $a$ since typical measures like standard error wouldn't apply. Is there some sort of established method to estimate this sort of error on the peak of a skewed distribution? Thank you for any help.

Answer 1

There are three answers, depending on what your final goal is:

(1) The Bootstrapped Distribution: the simplest way to express the uncertainty in $a$ would be to simply plot the bootstrapped distribution of values of $a$ using a histogram, and let your audience draw their own conclusions about the relationship between x and y. This won't allow you to do any hypothesis testing around $a$, but if your goal is just to express the precision of your estimate then it's perfectly acceptable.

(2) Bias-corrected Accelerated Bootstrap (BCa): if you want to use bootstrapping to estimate confidence intervals "the right way" then this is your tool. It will allow you to use the bootstrapped distribution to perform hypothesis testing. The calculation here is a little uglier, so I would recommend using a package such as R's boot

(3) Transform the data: from the way you've phrased your question it sounds like you've been estimating $a$ by with something like $\hat{a} = \text{argmin}_a \sum_i |y_i - x_i^a|^2$, where $i$ is indexing over your datapoints. This isn't wrong per se, but it is inconvenient for the reasons you're encountering. If you're willing to change the loss function, it would be far simpler to simply fit the model $\log(y) = a \log(x)$. As you noted, the coefficient estimates of a linear regression are normally distributed, so by transforming your data you simplify the problem enormously, since now you have simple linear regression (with the minor twist the coefficient is know to be equal to 0, which doesn't matter for the results you need). This is the approach I would recommend, assuming it doesn't violate the constraints of your problem for some other reason.

Answer 2

A non-parametric way of expressing uncertainty in this situation could be to report percentiles. For example, you could find the 2.5th and 97.5th percentiles, producing an interval containing 95% of your bootstrapped estimates. You could then use the bootstrapped percentiles as error bars around your estimate. This would also show the skewedness of your estimates in your plots in simple way, as the error bars would be asymmetric around the original estimate.