Questions about a model calibration

[florianhartig]

Questions by @Nastaran-ch moved into extra issue from #204:

Dear Florian,

I am planning to use Bayesian Tools package in R to calibrate a model for my PhD research project. I have read useful and applicable examples you have prepared in the below website: https://cran.r-project.org/web/packages/BayesianTools/vignettes/BayesianTools.html

I follow the instructions in Bayesian Tools package in my codes but have some problems:

For some models, I get different MAP output for parameters in each run. I don’t get good convergence for all parameters when I plot the trace of each parameter. I used “dnorm” in “singlelikelihood”, for all models and datasets. Should I choose different distribution function for this part or normal function is good enough? Kind regards, Nastaran

`My_Data=read.csv('example.csv')

library(Metrics) library(ie2misc) library(BayesianTools)

likelihood1 <- function(param){ pred = param[1]+param2)))

singlelikelihoods = dnorm(My_Data$G, mean = pred, sd = 1/(param[4]^2), log = T) return(sum(singlelikelihoods)) } setUp1 <- createBayesianSetup(likelihood1, lower = c(-10,-10,-5,0.01), upper = c(20,20,10,30)) settings = list(nrChains = 4,iterations = 15000, message = FALSE) out1 <- runMCMC(bayesianSetup = setUp1, sampler = "Metropolis", settings = settings) #sampler: SMC,Metropolis,DEzs,DREAM,Twalk M1 = marginalLikelihood(out1) M1 print(out1) plot(out1) summary(out1) correlationPlot(out1) marginalPlot(out1, prior = TRUE) marginalLikelihood(out1) DIC(out1) Matrix_out1=getSample(out1, start = 100, end = NULL, thin = 5, whichParameters = 1:4) MAP(out1)`

[florianhartig]

Hi @Nastaran-ch,

your first 2 questions (different MAP and lack of convergence) are related - if your chains are not converged, it is logical that you will get different MAP values. Thus, the first thing you have to do is to achieve convergence. How to do this is a matter of experience, and of course, there is always the option to run the sampler for longer, but here are a few hints:

• Use DEzs instead of Metropolis
• Inspect the trace plots, first of all your samplers have to run into the relevant posterior space, all samples before should be removed as burn-in
• Then sample long enough to achieve convergence

There are many more things one can do but I would start with this as a base.

Regarding the likelihood: whether a normal likelihood is a good idea or not depends on your residuals, i.e. your residuals should be normally distributed, as in a linear regression, and if that is not the case, you should consider more complicated likelihoods. You can either check residuals based on the MAP, or (more general) via posterior predictive simulations, e.g. with the DHARMa package, see https://cran.r-project.org/web/packages/DHARMa/index.html.

[ghost]

Hi Florian,

Thanks for your prompt reply and the help. The DEzs and burn-in worked very well for convergence.

I could not find support for all distribution functions by ‘getDharmaResiduals’, so to estimate the best fit between different distribution functions, I used glm with 'gaussian' , 'Gamma' , 'poisson' , 'binomial' , 'quasi' , 'quasibinomial' , 'quasipoisson'. The results showed 'gaussian' and 'quasi' have lowest residual sum of square. Should dnorm function be changed for ‘singlelikelihoods’ or that’s fine? Is there any other suggestion to improve the calibration processes?

Thanks in advance for your time and help.

[florianhartig]

Hi @Nastaran-ch,

the goal is not to get the lowest SumSq, but that the residuals conform to your distributional assumptions. So, for example, if you take your model with normal residuals, you could

a) simply look at the residuals MAP pred - Obs. Those should be normal, or if you assume a different distribution in the likelihood, they should conform to that.

b) the more elaborate way is to create posterior predictive simulations (PPP) by either predicting from the MAP or the posterior, and then add the error that you assume in the likelihood on top of your predictions (using the fitted parameters). If you have those simulations, you can read them in with DHARMa::createDHARMa and diagnose as described in the DHARMa package