Potential calibration issue

Looking at 'Sensitivity analysis' vignette, the calibration seems unsuccessful at partitioning the uncertainty between process/observation error and parameter error. In the last graph, the 90% credible interval (GPP uncertainty due to the uncertainty in the parameters) is so small that it is unvisible. All the uncertainty is absorbed by the process error (grey area). This is not necessarily an issue, but it does not hint at a very healthy model either. Issues could be:

Error in computing the likelihood (unlikely but a mistake in unpacking arrays etc could always happen)
Inappropriate priors: the use of flat priors on very limited ranges might prevent MCMC to explore enough of the parameter space. Actually, using a flat prior for a variance is rather unusual. An inverse gamma is more natural in this case (although it might not change the results).
Wrong normal assumption of the process/observation error.

Point three comes from the fact that both process error (error of the representation of the physical processes in the model itself) and the observation error (error in the estimation of GPP) are gathered in a unique normal random variable. It is common to model observation errors as a normal probability. There is however no reason for the process error to be normal, and wrongly assuming so may 'force' MCMC to overemphasize the process/observation error to meet the normality constraint. This might actually be exacerbated by the choice of flat priors for the model parameters, which cannot 'counterbalance' the strong assumption on the process error.

@fabern @stineb Any thought?

Yes, we discussed this extensively with @pepaaran who implemented the Bayesian calibration (and who is a statistician by training). This result reflects that the calibration yields a narrow band for the best estimate of the parameters. I don't think this is per se "unsucessful", but may indicate that the model is relatively rigid, or that the data itself has a large error that we just cannot fit better with the given model structure. I don't expect that a different specification of the process error distribution or the priors would change that picture. Is there a justification for using a specific (non-normal) distribution for the process error? If so, of course, we can try it out.

I agree that changing the prior will not change the observed behavior, and the rigidity of the model may indeed not be an issue. I'm just explaining further my line of thought regarding the inverse beta prior before closing the issue.

Inverse gamma is usually used as prior for a variance because it is positive (as a variance is) and is a conjugate prior with two parameters that have good interpretability. The issue with using a uniform prior is that it defines hard bounds to the variance (the PDF is 0 outside of the bounds), which assumes a fairly good a priori knowledge of err_gpp. However, by definition, err_gpp is modeling the uncertainty of our model of GPP (as well as the uncertainty in GPP partitioning) which can't really be known a priori, and an uninformative prior should be used instead (which is possible with an inverse gamma distribution). The bounds used in the vignettes were probably chosen a posteriori (after having already run a calibration using the same forcing). This is frown upon as it uses the same data twice (and is called 'double dipping').
The uniform priors for the other parameters are fine, as the bounds seem to come from the literature (although using normal priors would solve the issue of being stuck at a bound that was set too narrow).

You're absolutely right. The error distribution of err_gpp should not be normal.

We discussed it with Fabian and concluded that some work could be done to add uncertainty modelling, but this would be a large undertaking. I am closing this issue as the original question was answered.

I noticed that there was actually a bug in the code that generated the calibrated figure.
It is the all-too-familiar dplyr::summarise bug. If you overwrite variables without renaming, order of operations matters. Hence recommendation is to always rename variables! The issue was here, leading to the same q95 as q50:

rsofun/analysis/03-uncertainty-estimation.R

Lines 69 to 70 in 2d8fded

    
           gpp = quantile(gpp, 0.5, na.rm = TRUE),          # get median 
        
           gpp_q95 = quantile(gpp, 0.95, na.rm = TRUE),

The impact though remains limited (see updated plots below, which look better than the preliminary from the earlier post).
However, the parameter soilmbetao is estimated below 0.1, contrary to the manuscript draft where it was 0.5. What has changed?
Furthermore, the results below use 12000 iterations of 3 MCMC chains, of which 1500 are burnin. These are the values taken from the code (analysis/02-bayesian-calibration.R).

rerun parameter calibration with main branch and update plot in this graph
understand why soilmbetao is estimated below 0.1, while previously it was 0.5
update ms_rsofun_doc_v4.pdf with

Intermediate result of calibration with 3x12000 MCMC iterations:

Ufff. Good catch of the bug. Thanks!

It makes a lot of sense if soilmbetao is close to zero: zero GPP when the bucket is empty. It will depend a lot on the bucket size if that is prescribed.

@fabern , re the summarise() bug: does this affect the calibration itself or just the displaying of the uncertainty band in this plot?

Just the plotting of the uncertainty band*. The calibration is unaffected.
But as stated: the calibrated soilmbetao is different from the manuscript pdf v4... not sure what changed.

* If you watch (very closely) in the manuscript pdf you'll see that the green are is cropped at the yellow median line, while in the plot in the post above, the green area extends above the yellow median line.

@fabern Can we close this issue?

	gpp = quantile(gpp, 0.5, na.rm = TRUE), # get median
	gpp_q95 = quantile(gpp, 0.95, na.rm = TRUE),