Bayesian Methods for Ecological and Environmental Modelling
11:30 am-1 pm: Session 6b/7a Process Based Models
1-2 pm: Lunch break
2-3:30 pm: Session 6b/7a Process Based Models
Model structure based on mechanistic process understanding
Often deterministic
Good for extrapolation (eg forecasting) and “what if” type questions
Examples:
Aim to understand/predict something in nature
Build a model representing our current best understanding of the processes involved
Collect evidence/data to test/validate model against
Model does not agree with data/evidence
Errors in models and data are acknowledged but often ignored
Predictions are made without quantifying their reliability
Does it matter that we ignore uncertainty?
“…earlier on today apparently a woman rang the BBC and said she’d heard there was a hurricane on the way. Well if you are watching don’t worry, there isn’t…”
“…earlier on today apparently a woman rang the BBC and said she’d heard there was a hurricane on the way. Well if you are watching don’t worry, there isn’t…”
18 died in Britain, four in France
Devastation costs more than one billion pounds
An estimated 15 million trees were lost
Worst storm since 1703, one in 200 year storm for southern Britain
A public enquiry was held and an internal enquiry was conducted by the Met Office
“…unlike 1987, we now run multiple forecasts (called ensembles) that help us to understand the probabilities involved in a forecast and give a better estimate of the risk of high impact weather events.”
Model predictions without uncertainty could miss less likely but high impact events
Variability at scales that we are currently not interested in modelling/predicting.
Do observational samples adequately represent the site/system we are interested in - experimental design
Is measurement site (a)typical of (eg pine forest) sites in Europe?
Climate/weather data, N deposition, management etc.
For example which initial soil/tree values are consistent with each other and present day observations
All PBMs are wrong
Toy model (designed as part of a theory paper)
Included in BayesianTools (R package that we will use)
Very simplified form of more complex models
To give a concrete example of inferring the parameters of a PBM
PAR - Photosynthetically active radiation
KEXT - Beer’s law light extinction coeff
LUE - Light use efficiency
LAI - Leaf area index (LAR * Cv)
LAR - Leaf area ratio
Cv/Cs - Vegetative/Soil carbon
GAMMA - Fraction of GPP that is Autotrophic Respiration
\(C_v\), \(C_r\) and \(C_s\) : Carbon in vegetation, root and soil pools
best | lower | upper | |
---|---|---|---|
KEXT | 0.500 | 2e-01 | 1e+00 |
LAR | 1.500 | 2e-01 | 3e+00 |
LUE | 0.002 | 5e-04 | 4e-03 |
GAMMA | 0.400 | 2e-01 | 6e-01 |
tauV | 1440.000 | 5e+02 | 3e+03 |
tauS | 27370.000 | 4e+03 | 5e+04 |
tauR | 1440.000 | 5e+02 | 3e+03 |
Av | 0.500 | 2e-01 | 1e+00 |
Cv | 3.000 | 0e+00 | 4e+02 |
Cs | 15.000 | 0e+00 | 1e+03 |
Cr | 3.000 | 0e+00 | 2e+02 |
Inference of model parameters and initial conditions of carbon pools (Cv, Cs, Cr)
The model structure is fixed i.e. is not being inferred
We will not use measured data from the field, but virtual data.
referenceData <- VSEM(refPars$best[1:11], PAR)
referenceData[,1] = 1000 * referenceData[,1]
aobs = referenceData
aobs[,1] = abs(referenceData[,1])
maobs = colMeans(aobs)
obs = referenceData
obs[,1] <- referenceData[,1] + rnorm(length(referenceData[,1]), sd = maobs[1]*refPars$best[12])
obs[,2] <- referenceData[,2] + rnorm(length(referenceData[,1]), sd = maobs[2]*refPars$best[12])
obs[,3] <- referenceData[,3] + rnorm(length(referenceData[,1]), sd = maobs[3]*refPars$best[12])
obs[,4] <- referenceData[,4] + rnorm(length(referenceData[,1]), sd = maobs[4]*refPars$best[12])
PBMs often preexist in a language (R, python, FORTRAN, C) Could be costly to recode in an Bayesian package’s bespoke language
For PBMs very useful to be able to code a bespoke likelihood function as we need control over calling the PBM
Few MCMC packages allow this (eg BayesianTools, NIMBLE)
Gaussian likelihood Logarithm helps avoid numerical rounding issues
likelihood <- function(par, sum = TRUE){
x = refPars$best
x[parSel] = par
predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model
predicted[,1] = 1000 * predicted[,1] # this is just rescaling
diff <- c(predicted[,1] - obs[,1])
llValues1 <- dnorm(diff, sd = x[12]*maobs[1], log = T)
diff <- c(predicted[,2] - obs[,2])
llValues2 <- dnorm(diff, sd = x[12]*maobs[2], log = T)
diff <- c(predicted[,3] - obs[,3])
llValues3 <- dnorm(diff, sd = x[12]*maobs[3], log = T)
diff <- c(predicted[,4] - obs[,4])
llValues4 <- dnorm(diff, sd = x[12]*maobs[4], log = T)
return(sum(llValues1,llValues2,llValues3,llValues4))
}
Should reflect prior thinking about uncertainty in parameters
More than one chain helps test for convergence
bayesianSetup <- createBayesianSetup(likelihood, prior,
names = rownames(refPars)[parSel])
settings <- list(iterations = 10000, nrChains = 2)
out <- runMCMC(bayesianSetup = bayesianSetup,
sampler = "DEzs", settings = settings)
Running DEzs-MCMC, chain 1 iteration 300 of 10002 . Current logp -221813.5 -2574404 -238855.1 . Please wait!
Running DEzs-MCMC, chain 1 iteration 600 of 10002 . Current logp -192612.6 -13975.48 -192521 . Please wait!
Running DEzs-MCMC, chain 1 iteration 900 of 10002 . Current logp -13777.76 -13975.48 -23585.09 . Please wait!
Running DEzs-MCMC, chain 1 iteration 1200 of 10002 . Current logp -13284.03 -12549.67 -14109.13 . Please wait!
Running DEzs-MCMC, chain 1 iteration 1500 of 10002 . Current logp -12294.14 -12376.23 -12468.31 . Please wait!
Running DEzs-MCMC, chain 1 iteration 1800 of 10002 . Current logp -11893.12 -12247.43 -11343.65 . Please wait!
Running DEzs-MCMC, chain 1 iteration 2100 of 10002 . Current logp -10823.31 -10806.83 -10925.68 . Please wait!
Running DEzs-MCMC, chain 1 iteration 2400 of 10002 . Current logp -10756.36 -10786.47 -10852.72 . Please wait!
Running DEzs-MCMC, chain 1 iteration 2700 of 10002 . Current logp -10747.73 -10786.47 -10842.44 . Please wait!
Running DEzs-MCMC, chain 1 iteration 3000 of 10002 . Current logp -10706.31 -10738.19 -10813.63 . Please wait!
Running DEzs-MCMC, chain 1 iteration 3300 of 10002 . Current logp -10684.52 -10711.33 -10698.76 . Please wait!
Running DEzs-MCMC, chain 1 iteration 3600 of 10002 . Current logp -10686.39 -10682.21 -10681.47 . Please wait!
Running DEzs-MCMC, chain 1 iteration 3900 of 10002 . Current logp -10674.72 -10682.21 -10646.98 . Please wait!
Running DEzs-MCMC, chain 1 iteration 4200 of 10002 . Current logp -10669.24 -10642.01 -10646.98 . Please wait!
Running DEzs-MCMC, chain 1 iteration 4500 of 10002 . Current logp -10642.12 -10642.01 -10646.98 . Please wait!
Running DEzs-MCMC, chain 1 iteration 4800 of 10002 . Current logp -10597.76 -10611.22 -10641.4 . Please wait!
Running DEzs-MCMC, chain 1 iteration 5100 of 10002 . Current logp -10454.25 -10611.22 -10638.23 . Please wait!
Running DEzs-MCMC, chain 1 iteration 5400 of 10002 . Current logp -9903.171 -10583.33 -10472.17 . Please wait!
Running DEzs-MCMC, chain 1 iteration 5700 of 10002 . Current logp -9867.326 -10564.4 -10438.35 . Please wait!
Running DEzs-MCMC, chain 1 iteration 6000 of 10002 . Current logp -9809.074 -10423.17 -9873.418 . Please wait!
Running DEzs-MCMC, chain 1 iteration 6300 of 10002 . Current logp -9611.766 -9773.626 -9679.813 . Please wait!
Running DEzs-MCMC, chain 1 iteration 6600 of 10002 . Current logp -9466.296 -9638.427 -9614.455 . Please wait!
Running DEzs-MCMC, chain 1 iteration 6900 of 10002 . Current logp -9375.288 -9638.427 -9303.79 . Please wait!
Running DEzs-MCMC, chain 1 iteration 7200 of 10002 . Current logp -9239.417 -9454.801 -9299.174 . Please wait!
Running DEzs-MCMC, chain 1 iteration 7500 of 10002 . Current logp -9206.008 -9366.308 -9288.549 . Please wait!
Running DEzs-MCMC, chain 1 iteration 7800 of 10002 . Current logp -9155.129 -9246.229 -9274.466 . Please wait!
Running DEzs-MCMC, chain 1 iteration 8100 of 10002 . Current logp -9105.075 -9174.73 -9259.018 . Please wait!
Running DEzs-MCMC, chain 1 iteration 8400 of 10002 . Current logp -9037.463 -9174.73 -9157.09 . Please wait!
Running DEzs-MCMC, chain 1 iteration 8700 of 10002 . Current logp -8825.387 -9166.308 -9131.233 . Please wait!
Running DEzs-MCMC, chain 1 iteration 9000 of 10002 . Current logp -8818.227 -9151.364 -9012.713 . Please wait!
Running DEzs-MCMC, chain 1 iteration 9300 of 10002 . Current logp -8782.651 -8810.989 -8703.661 . Please wait!
Running DEzs-MCMC, chain 1 iteration 9600 of 10002 . Current logp -8691.51 -8723.941 -8501.368 . Please wait!
Running DEzs-MCMC, chain 1 iteration 9900 of 10002 . Current logp -8457.715 -8521.397 -8419.513 . Please wait!
Running DEzs-MCMC, chain 1 iteration 10002 of 10002 . Current logp -8457.715 -8521.397 -8296.083 . Please wait!
Running DEzs-MCMC, chain 2 iteration 300 of 10002 . Current logp -946789.3 -587129.7 -1298392 . Please wait!
Running DEzs-MCMC, chain 2 iteration 600 of 10002 . Current logp -435046.4 -404600.4 -471514.1 . Please wait!
Running DEzs-MCMC, chain 2 iteration 900 of 10002 . Current logp -185940.4 -243414.8 -189215.8 . Please wait!
Running DEzs-MCMC, chain 2 iteration 1200 of 10002 . Current logp -185940.4 -200537.3 -189215.8 . Please wait!
Running DEzs-MCMC, chain 2 iteration 1500 of 10002 . Current logp -179087.7 -188539.4 -184431.9 . Please wait!
Running DEzs-MCMC, chain 2 iteration 1800 of 10002 . Current logp -161545.9 -147391.3 -147859.2 . Please wait!
Running DEzs-MCMC, chain 2 iteration 2100 of 10002 . Current logp -123638.8 -131547.3 -132917.7 . Please wait!
Running DEzs-MCMC, chain 2 iteration 2400 of 10002 . Current logp -84244.48 -109220.7 -88193.19 . Please wait!
Running DEzs-MCMC, chain 2 iteration 2700 of 10002 . Current logp -78351.16 -68278.21 -63522.77 . Please wait!
Running DEzs-MCMC, chain 2 iteration 3000 of 10002 . Current logp -48634.54 -65941.92 -47602.59 . Please wait!
Running DEzs-MCMC, chain 2 iteration 3300 of 10002 . Current logp -46909.14 -59629.18 -40836.18 . Please wait!
Running DEzs-MCMC, chain 2 iteration 3600 of 10002 . Current logp -34815.83 -59629.18 -32626.59 . Please wait!
Running DEzs-MCMC, chain 2 iteration 3900 of 10002 . Current logp -22225.77 -28268.33 -25164.37 . Please wait!
Running DEzs-MCMC, chain 2 iteration 4200 of 10002 . Current logp -15434.03 -19783.08 -18383.11 . Please wait!
Running DEzs-MCMC, chain 2 iteration 4500 of 10002 . Current logp -14704.24 -19783.08 -14961.86 . Please wait!
Running DEzs-MCMC, chain 2 iteration 4800 of 10002 . Current logp -12863.96 -19783.08 -12047.53 . Please wait!
Running DEzs-MCMC, chain 2 iteration 5100 of 10002 . Current logp -12863.96 -12583.28 -11722.26 . Please wait!
Running DEzs-MCMC, chain 2 iteration 5400 of 10002 . Current logp -12863.96 -12583.28 -11335.35 . Please wait!
Running DEzs-MCMC, chain 2 iteration 5700 of 10002 . Current logp -12256.05 -12583.28 -11009.57 . Please wait!
Running DEzs-MCMC, chain 2 iteration 6000 of 10002 . Current logp -10655.78 -11364.26 -10891.59 . Please wait!
Running DEzs-MCMC, chain 2 iteration 6300 of 10002 . Current logp -10120.02 -10803 -10227.24 . Please wait!
Running DEzs-MCMC, chain 2 iteration 6600 of 10002 . Current logp -10008.7 -10148.3 -10173.9 . Please wait!
Running DEzs-MCMC, chain 2 iteration 6900 of 10002 . Current logp -9928.805 -10097.65 -10173.9 . Please wait!
Running DEzs-MCMC, chain 2 iteration 7200 of 10002 . Current logp -9928.805 -10080 -10172.49 . Please wait!
Running DEzs-MCMC, chain 2 iteration 7500 of 10002 . Current logp -9922.478 -10061.75 -9956.195 . Please wait!
Running DEzs-MCMC, chain 2 iteration 7800 of 10002 . Current logp -9922.478 -10018.28 -9956.195 . Please wait!
Running DEzs-MCMC, chain 2 iteration 8100 of 10002 . Current logp -9922.478 -10018.28 -9956.195 . Please wait!
Running DEzs-MCMC, chain 2 iteration 8400 of 10002 . Current logp -9922.478 -9998.835 -9932.861 . Please wait!
Running DEzs-MCMC, chain 2 iteration 8700 of 10002 . Current logp -9912.895 -9855.134 -9932.861 . Please wait!
Running DEzs-MCMC, chain 2 iteration 9000 of 10002 . Current logp -9912.895 -9844.444 -9921.684 . Please wait!
Running DEzs-MCMC, chain 2 iteration 9300 of 10002 . Current logp -9912.895 -9818.236 -9860.527 . Please wait!
Running DEzs-MCMC, chain 2 iteration 9600 of 10002 . Current logp -9896.801 -9816.592 -9830.889 . Please wait!
Running DEzs-MCMC, chain 2 iteration 9900 of 10002 . Current logp -9848.399 -9809.455 -9830.889 . Please wait!
Running DEzs-MCMC, chain 2 iteration 10002 of 10002 . Current logp -9848.399 -9786.997 -9830.889 . Please wait!
# # # # # # # # # # # # # # # # # # # # # # # # #
## MCMC chain summary ##
# # # # # # # # # # # # # # # # # # # # # # # # #
# MCMC sampler: DEzs
# Nr. Chains: 6
# Iterations per chain: 3334
# Rejection rate: 0.97
# Effective sample size: 67
# Runtime: 4.39 sec.
# Parameters
psf MAP 2.5% median 97.5%
KEXT 2.249 0.384 0.214 0.426 0.806
LUE 1.133 0.002 0.002 0.003 0.004
tauV 3.309 693.985 618.596 2012.347 2996.495
tauS 1.226 9090.541 4076.336 12593.051 25224.597
Av 15.516 0.904 0.205 0.792 0.999
Cv 1.585 3.897 0.102 4.954 34.952
Cs 1.282 12.381 9.125 17.395 222.179
Cr 3.367 4.424 0.368 4.407 101.340
error-sd 2.093 0.455 0.481 0.663 0.699
## DIC: 980006.3
## Convergence
Gelman Rubin multivariate psrf: 84.415
## Correlations
KEXT LUE tauV tauS Av Cv Cs Cr error-sd
KEXT 1.000 0.458 0.084 -0.241 0.194 0.403 0.457 0.446 -0.493
LUE 0.458 1.000 0.153 0.209 0.261 0.016 0.011 0.048 0.070
tauV 0.084 0.153 1.000 -0.218 -0.803 0.130 0.005 -0.033 0.330
tauS -0.241 0.209 -0.218 1.000 0.377 -0.085 0.159 0.198 0.491
Av 0.194 0.261 -0.803 0.377 1.000 -0.199 -0.026 0.066 -0.149
Cv 0.403 0.016 0.130 -0.085 -0.199 1.000 0.718 0.643 -0.360
Cs 0.457 0.011 0.005 0.159 -0.026 0.718 1.000 0.763 -0.243
Cr 0.446 0.048 -0.033 0.198 0.066 0.643 0.763 1.000 -0.173
error-sd -0.493 0.070 0.330 0.491 -0.149 -0.360 -0.243 -0.173 1.000
# # # # # # # # # # # # # # # # # # # # # # # # #
## MCMC chain summary ##
# # # # # # # # # # # # # # # # # # # # # # # # #
# MCMC sampler: DEzs
# Nr. Chains: 6
# Iterations per chain: 50000
# Rejection rate: 0.943
# Effective sample size: 719
# Runtime: 229.801 sec.
# Parameters
psf MAP 2.5% median 97.5%
KEXT 1.001 0.529 0.254 0.538 0.820
LUE 1.001 0.002 0.002 0.002 0.002
tauV 1.005 1394.867 1171.209 1425.719 2844.404
tauS 1.002 25672.714 21957.750 26787.986 49039.340
Av 1.008 0.508 0.355 0.503 0.578
Cv 1.002 3.020 2.849 3.008 3.460
Cs 1.007 14.915 14.503 14.952 16.571
Cr 1.009 3.037 2.648 3.025 3.371
error-sd 1.003 0.100 0.098 0.101 0.413
## DIC: 65755.55
## Convergence
Gelman Rubin multivariate psrf: 1.012
## Correlations
KEXT LUE tauV tauS Av Cv Cs Cr error-sd
KEXT 1.000 -0.386 -0.621 0.035 0.588 -0.035 -0.040 -0.019 -0.271
LUE -0.386 1.000 0.251 -0.112 -0.174 -0.079 -0.008 0.020 -0.248
tauV -0.621 0.251 1.000 0.507 -0.821 0.072 0.089 -0.002 0.367
tauS 0.035 -0.112 0.507 1.000 -0.371 0.014 0.045 -0.004 0.109
Av 0.588 -0.174 -0.821 -0.371 1.000 0.109 0.072 0.151 -0.263
Cv -0.035 -0.079 0.072 0.014 0.109 1.000 0.789 0.688 0.420
Cs -0.040 -0.008 0.089 0.045 0.072 0.789 1.000 0.711 0.484
Cr -0.019 0.020 -0.002 -0.004 0.151 0.688 0.711 1.000 0.444
error-sd -0.271 -0.248 0.367 0.109 -0.263 0.420 0.484 0.444 1.000
best | lower | upper | |
---|---|---|---|
KEXT | 0.500 | 2e-01 | 1e+00 |
LAR | 1.500 | 2e-01 | 3e+00 |
LUE | 0.002 | 5e-04 | 4e-03 |
GAMMA | 0.400 | 2e-01 | 6e-01 |
tauV | 1440.000 | 5e+02 | 3e+03 |
tauS | 27370.000 | 4e+03 | 5e+04 |
tauR | 1440.000 | 5e+02 | 3e+03 |
Av | 0.500 | 2e-01 | 1e+00 |
Cv | 3.000 | 0e+00 | 4e+02 |
Cs | 15.000 | 0e+00 | 1e+03 |
Cr | 3.000 | 0e+00 | 2e+02 |
error-sd | 0.100 | 1e-02 | 7e-01 |
VSEM is fast but even a “simple” model needed 150k runs
Much more complex “slower/expensive” models
Sometime good enough to recode model in a more efficient language eg C or FORTRAN
More efficient inference schemes might be useful eg R package NIMBLE
We will now move onto a practical session
The web document for this can be accessed from the main page
What is a process based model PBM?
Why considering uncertainties in PBM predictions may be important
Sources of uncertainty in the Bayesian inference of PBMs
Worked example: Inference of a simple PBM (VSEM) using BayesianTools
Practical: Why it is important to represent all sources of uncertainty in Bayesian inference?