MCMCEc4photo: Markov chain Monte Carlo for C4 photosynthesis parameters
In BioCro: Simulation of Biomass Crops

Description Usage Arguments Value References Examples

This function attempts to implement Markov chain Monte Carlo methods for models with no likelihoods. In this case it is done for the von Caemmerer C4 photosynthesis model. The method implemented is based on a paper by Marjoram et al. (2003). The method is described in Miguez (2007). The chain is constructed using a Gaussian random walk. This is definitely a beta version of this function and more testing and improvements are needed. The value of this function is in the possibility of examining the empirical posterior distribution (i.e. vectors) of the Vcmax and alpha parameters. Notice that you will get different results each time you run it.

1 2	MCMCEc4photo(obsDat, niter = 30000, iCa = 380, iOa = 210, iVcmax = 60, iVpmax = 120,iVpr = 80, iJmax = 400, thresh = 0.01, scale = 1)

`obsDat`	observed assimilation data, which should be a data frame or matrix. The first column should be observed net assimilation rate (μ mol m^{-2} s^{-1}). The second column should be the observed quantum flux (μ mol m^{-2} s^{-1}). The third column should be observed temperature of the leaf (Celsius). The fourth column should be the observed relative humidity in proportion (e.g. 0.7).
`niter`	number of iterations to run the chain for (default = 30000).
`iCa`	CO2 atmospheric concentration (ppm or μbar).
`iOa`	O2 atmospheric concentration (mbar).
`iVcmax`	initial value for Vcmax (default = 60).
`iVpmax`	initial value for Vpmax (default = 120).
`iVpr`	initial value for Vpr (default = 80).
`iJmax`	initial value for Jmax (default = 400).
`thresh`	this is a threshold that determines the “convergence” of the initial burn-in period. The convergence of the whole chain can be evaluated by running the model with different starting values for Vcmax and alpha. The chain should convergence, but for this, runs with similar thresholds should be compared. This threshold reflects the fact that for any given data the model will not be able to simulate it perfectly so it represents a compromise between computability and fit.
`scale`	This scale parameter controls the size of the standard deviations which generate the moves in the chain. Decrease it to increase the acceptance rate and viceversa.

a list structure with components

`RsqBI`	This is the R^2 for the “burn-in” period. This value becomes the cut off value for the acceptance in the chain.
`CoefBI`	parameter estimates after the burn-in period. These are not optimal as in the case of the optimization routine but are starting values for the chain.
`accept1`	number of iterations for the initial burn-in period.
`resuBI`	matrix of dimensions 5 by `accept1` containing the values for Vcmax and alpha and the R^2 in each iteration of the burn-in period.
`resuMC`	matrix of dimensions 5 by `accept2` containing the values for Vcmax and alpha and the R^2 in each iteration of the chain period.
`accept2`	number of accepted samples or length of the chain.
`accept3`	number of accepted moves in the chain.

P. Marjoram, J. Molitor, V. Plagnol, S. Tavare, Markov chain monte carlo without likelihoods, PNAS 100 (26) (2003) 15324–15328.

Miguez (2007) Miscanthus x giganteus production: meta-analysis, field study and mathematical modeling. PhD Thesis. University of Illinois at Urbana-Champaign.

## Not run: 
## This is an example for the MCMCEc4photo
## evaluating the convergence of the chain
## Notice that if a parameter does not seem
## to converge this does not mean that the method
## doesn't work. Careful examination is needed
## in order to evaluate the validity of the results
data(obsNaid)
res1 <- MCMCEc4photo(obsNaid,100000,thresh=0.007)
res1

## Run it a few more times
## and test the stability of the method
res2 <- MCMCEc4photo(obsNaid,100000,thresh=0.007)
res3 <- MCMCEc4photo(obsNaid,100000,thresh=0.007)

## Now plot it
plot(res1,res2,res3)
plot(res1,res2,res3,type="density")

## End(Not run)