What are some reasonable ways to estimate systematic uncertainties?

Question

I have a set of data $(x_{\rm data}, y_{\rm data}, \Delta y_{{\rm data}})$, where $\Delta y_{{\rm data}}$ is the uncertainty in the measurement. I also have some model used to explain the data. The model is non-linear:$$y_{\rm model} =\mathcal{F}(\bf{p}; x_{\rm model})~,$$ $\bf{p}$ being the model parameters. I have some prior knowledge about $\bf{p}$, so I run MCMC chain to estimate the posterior distribution of $\bf{p}$, and use that distribution to find the best-fit values and the uncertainties for $\bf{p}$. This is all good. However, I have few (less than 10) candidate models to explain the dataset, and the true model is unknown. All these models give different estimate for best-fit values and the uncertainties of $\bf{p}$. Thus, I need some good method which could be used to calculate the best-fit values and the uncertainties of $\bf{p}$, such that the model uncertainty is also included. While calculating such systematic uncertainty I also need to keep in mind that some of these models fit the data poorly compared to others (measured with some metric like $\chi^2 =\sum (y_{\rm data}- y_{\rm model})^2/\Delta y_{\rm data}^2$). Any suggestions?

P.S. For model selection, I came across information criterion like the Akaike Information Criterion, so my naive thought was choose the best model from such criterion and report the best-fit value and uncertainty of the parameters for such best-fit model.

Answer 1

Based on the comment clarifying that measurement error is known, you may consider the implementation for errors-in-variables modeling in the brms package. It seems like you're using MCMC anyway, so HMC from Stan is usually just more efficient. Assuming a generic linear regression in R, you would specify something like this:

fit <- brm(out | me(out, error) ~ x1, 
           data = dat)

Just to briefly overview that code, the real work is the me() function. This is just short for "measurement error" and has the two arguments out and error. In this case, out is the outcome variable (aka, our variable with measurement error), and error is the error that we know for that variable (on a standard deviation scale).

I would say, however, that linear regression assumptions are usually not too impacted by violating perfect measurement accuracy in the outcome variable (usually more an issue with the predictors). The residuals of the model generally cover the measurement error in the outcome. That makes some intuitive sense: the residual is error not explained by the predictors, so there's no reason that the predictors should be able to explain the variance due to measurement error.

That said, if you have reason to think that the residuals are related to certain predictors, then you can run those regressions as well. This is sometimes called a variance regression or a distributional model, but it's basically a joint estimation of a linear model for the outcome and the residual variance term. Again, in brms language, this would like like the following:

fit <- brm(bf(out ~ x1, sigma ~ x1),
           data = dat)

There are additional details about distributional models in brms here.