Bayesian parameter estimation with proportion data

Question

I am trying to do a Bayesian analysis using a model that comes from the literature in non-Bayesian form: $y = \Phi\Bigg(\frac{1}{\alpha} * log(A/\beta)\Bigg)$. Because the model uses the function $\Phi$ its outcome is in the range $0-1$. My observations cannot exceed this range either (but can be 0 or 1). However, using a Gaussian in the likelihood function means that the uncertainty bounds are sometimes outside the $0-1$ range (see Figure below). Statistically speaking, is this a problem, and if so, what would be a better solution?

$R_i \sim \mathcal{N}(\mu,\sigma)$

$\mu_i = \Phi \Bigg(\frac{1}{\alpha} * log(A/\beta) \Bigg)$

$\alpha \sim \mathcal{N}(0.35, 0.1) $

$\beta\sim \mathcal{N}(70, 10) $

$\sigma \sim \mathcal{U}(0, 0.3) $

Answer 1

We are talking about probabilistic model, so you need to define a likelihood function for $R_i$. As already mentioned by others, normal distribution is not the best choice because of it's support.

Since $R_i \in [0,1]$, as in beta regression (Ferrari and Cribari-Neto, 2004), you can choose re-parametrized beta distribution with $\alpha = \phi\mu$ and $\beta = \phi(1-\mu)$

$$ f(y) = \frac{1}{\mathrm{B}(\phi\mu,\; \phi(1-\mu))}\; y^{\phi\mu-1} (1-y)^{\phi(1-\mu)-1} $$

where $E(Y) = \mu$ and $\mathrm{Var}(Y) = \frac{\mu(1-\mu)}{1+\phi}$.

What translates to

$$\begin{align} \mu_i &= \Phi \Bigg(\frac{1}{\alpha} \log(A/\beta) \Bigg) \\ R_i &\sim \mathsf{Beta}(\mu, \phi) \end{align}$$

with some priors for $\alpha$, $\beta$, and $\phi$. The problem would be that for beta distribution, the support excludes 0 and 1, so you could either (1) use inflated beta distribution, or (2) transform your data to move the extreme values away from 0 and 1, e.g. by replacing $R_i$ with something like $(R_i \times (n−1) + 0.5) / n$.

_{Ferrari, S., & Cribari-Neto, F. (2004). Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799-815.}

Answer 2

Before you get to estimation, you need to fix your model so that it respects the support of the observable outcome variable. You say in the comments that this represents percentage damage to a building, so it is restricted to be between zero and one (i.e., one-hundred percent). In its present form, your model has a regression function that is limited to values in the support, but it allows errors to go beyond this. That is a bad model, since it allows observable values well outside their possible range, and it will give prediction intervals that extend beyond the possible range.

In view of the meaning of your output variable, I would recommend you use some form of nonlinear model where the output is restricted to the interval between zero and one. If you would like to use something like your present model, then the simplest way to correct its present problem ---while preserving the intended form of your model--- is to put the error terms inside the transformation, rather than adding it afterwards. Given explanatory variables $a_1,...,a_n$ and response variables $R_1,...,R_n$, you can model the latter as:

$$R_i = \Phi \Bigg( \frac{1}{\alpha} \cdot \log \bigg( \frac{a_i}{\beta} \bigg) + \varepsilon_i \Bigg) \quad \quad \quad \varepsilon_1,...,\varepsilon_n \sim \text{IID N}(0, \sigma^2).$$

This model form still uses the function $\Phi$ (standard notation for the CDF of the standard normal distribution) as your transformation, but now the error terms appear inside this function. This ensures that the response variable is subject to random error, but it must still remain between zero and one.

From this point you can proceed with your Bayesian analysis by imposing the relevant priors and finding the posterior distribution and predictive distribution (e.g., by MCMC methods). The above model form will ensure that the predictive distribution for the response variable respects its support, so your interval estimates should not go beyond the bounds of allowable values.

Answer 3

Since you are modelling the proportion of damage to a building, the data must be between 0 and 1, inclusive. This is troubling because the go to models for this sort of data (namely, logistic regression) usually operates with an integer number of trials from which the proportion was derived.

The other option is to use a beta likelihood, but since your data can and do achieve either 0 and 1, then you have to use a ZOIB Model.

ZOIB stands for: Zero One Inflated Beta. You can fit these very easily with brms. Also, since you specify that your mean should look like an erf function, then you should pass link='probit' to the family argument. here is a small simulation example I just cooked up.

library(brms)

#Simulate some data
set.seed(0)
x = sample(1:50,200, replace = T)

a = -12
b  = 4
eta = a+ b*log(x)
p = pnorm(eta)
y = rbeta(length(p), 20*p, 20*(1-p))

model.data = data.frame(x = x, y = y)

#The model
model = brm(y~log(x),
            data = model.data,
            family = zero_one_inflated_beta(link = 'probit'))


marginal_effects(model)

The mean of this model is a normal cdf (aka the erf function), which depends on the linear predictor of the form $\beta_0 + \beta_1 \log(x)$. This is equivalent to the linear predictor you provide, and it is easy to show with application of log rules and some algebra. Also note that the ZOIB model has many more knobs to tune/learn than what I have shown in this small example. I don't have a wealth of experience with these models myself, so I can't point you in the best direction for your application as of now.

This should maintain the overall structure of your model while solving the issue you have with the credible intervals going outside the support of the data. At least, I think. I will let much smarter members of the community let me know if I have made a gross misstep.