inst/doc/ParameterEstimation.R

In SoilR: Models of Soil Organic Matter Decomposition

### R code from vignette source 'ParameterEstimation.Rnw'

###################################################
### code chunk number 1: ParameterEstimation.Rnw:58-64
###################################################
library(SoilR)
library(FME)
library(MASS)
library(lattice)
BorealCO2=eCO2
names(BorealCO2)<-c("time","eCO2","eCO2sd")


###################################################
### code chunk number 2: ParameterEstimation.Rnw:71-75
###################################################
plot(BorealCO2[,1:2], xlab="Days", ylab="Evolved CO2 (mgC g-1 
     soil day-1)",ylim=c(0,50))
arrows(BorealCO2[,1],BorealCO2[,2]-BorealCO2[,3],BorealCO2[,1], 
       BorealCO2[,2]+BorealCO2[,3],code=3,angle=90,length=0.1)


###################################################
### code chunk number 3: ParameterEstimation.Rnw:90-108
###################################################
days=seq(0,35)
Ctotal=mean(c(0.04683975, 0.04703255, 0.04687287))*450 #C concentration * 450 g soil
BorealCO2=data.frame(time=BorealCO2[,1],BorealCO2[,2:3]*1e-06*Ctotal) #Changed units

eCO2func=function(pars){
  mod=TwopFeedbackModel(
    t=days,
    ks=pars[1:2],
    a21=0,
    a12=0, 
    C0=Ctotal*c(pars[3],1-pars[3]), 
    In=0,
    pass=TRUE
  )
  Rt=getReleaseFlux(mod)
  return(data.frame(time=days,eCO2=rowSums(Rt)))
}



###################################################
### code chunk number 4: ParameterEstimation.Rnw:115-119
###################################################
eCO2cost=function(pars){
  modelOutput=eCO2func(pars)
  return(modCost(model=modelOutput, obs=BorealCO2, err="eCO2sd"))
}


###################################################
### code chunk number 5: ParameterEstimation.Rnw:125-130
###################################################
inipars=c(k1=1/20,k2=1/100,gamma=0.5)

eCO2fit=modFit(f=eCO2cost,p=inipars,method="Marq",
               upper=c(Inf,Inf,1),lower=c(0,0,0))



###################################################
### code chunk number 6: ParameterEstimation.Rnw:134-135
###################################################
eCO2fit$par


###################################################
### code chunk number 7: ParameterEstimation.Rnw:140-141
###################################################
fitmod=eCO2func(eCO2fit$par)


###################################################
### code chunk number 8: ParameterEstimation.Rnw:146-150
###################################################
plot(BorealCO2[,1:2], xlab="Days", ylab="Evolved CO2 (gC day-1)",ylim=c(0,9e-04))
arrows(BorealCO2[,1],BorealCO2[,2]-BorealCO2[,3],BorealCO2[,1],
       BorealCO2[,2]+BorealCO2[,3],code=3,angle=90,length=0.1)
lines(fitmod)


###################################################
### code chunk number 9: ParameterEstimation.Rnw:158-163
###################################################
var0=eCO2fit$var_ms_unweighted

eCO2mcmc=modMCMC(f=eCO2cost, p=eCO2fit$par, niter=3000, jump=NULL,  
                 var0=var0, wvar0=NULL, updatecov=100, lower=c(0,0,0),
                 upper=c(1,1,1))


###################################################
### code chunk number 10: ParameterEstimation.Rnw:169-170
###################################################
summary(eCO2mcmc)


###################################################
### code chunk number 11: ParameterEstimation.Rnw:177-178
###################################################
pairs(eCO2mcmc)


###################################################
### code chunk number 12: ParameterEstimation.Rnw:188-194
###################################################
predRange=sensRange(func=eCO2func, parInput=eCO2mcmc$par)
plot(summary(predRange),ylim=c(0,9e-04),xlab="Days",
     ylab="Evolved CO2 (g C day-1)",main="")
points(BorealCO2)
arrows(BorealCO2[,1],BorealCO2[,2]-BorealCO2[,3],BorealCO2[,1],
       BorealCO2[,2]+BorealCO2[,3],code=3,angle=90,length=0.1)


###################################################
### code chunk number 13: ParameterEstimation.Rnw:208-210
###################################################
plot(D14C~Year,data=HarvardForest14CO2,
     ylab=expression(paste(Delta^14,"C ","(\u2030)")))


###################################################
### code chunk number 14: ParameterEstimation.Rnw:232-235
###################################################
time=C14Atm_NH$YEAR
t_start=min(time)
t_end=max(time)


###################################################
### code chunk number 15: ParameterEstimation.Rnw:240-245
###################################################
inputFluxes=BoundInFlux(
	function(t0){matrix(nrow=3,ncol=1,c(270,150,0))},
	t_start,
	t_end
)


###################################################
### code chunk number 16: ParameterEstimation.Rnw:250-251
###################################################
C0=c(390,220+390+1376,90+1800+560) 


###################################################
### code chunk number 17: ParameterEstimation.Rnw:257-276
###################################################
Fc=BoundFc(C14Atm_NH,lag=0,format="Delta14C")
Mod1<-function(ks,pass=TRUE){
  At=ConstLinDecompOp(
           matrix(nrow=3,ncol=3,byrow=TRUE,c(ks[1],0,0,
                                             ks[4],ks[2],0,
                                             ks[5],0,ks[3]))
  ) 
  mod=GeneralModel_14(
	t=time,
	A=At,
	ivList=C0,
	initialValF=ConstFc(rep(0,3),"Delta14C"),
	inputFluxes=inputFluxes,
	inputFc=Fc,
	pass=TRUE
  ) 
  R14t=getF14R(mod)
  return(data.frame(time=time,R14t=R14t))
}


###################################################
### code chunk number 18: ParameterEstimation.Rnw:280-283
###################################################
DataR14t=cbind(time=HarvardForest14CO2[,1],
               R14t=HarvardForest14CO2[,2],
               sd=sd(HarvardForest14CO2[,2]))


###################################################
### code chunk number 19: ParameterEstimation.Rnw:289-307
###################################################
#Create the cost function
R14tCost <- function(pars){
  R14t <- Mod1(pars)
  return(modCost(model=R14t,obs=DataR14t,err="sd"))
}

#Fit the model to the observed data given some initial value for the parameters
Fit <- modFit(f=R14tCost,p=c(-0.5,-0.9,-0.1,0.1,0.1))
#
# Run an MCMC using the variance and covariance results from the previous optimization
var0 <- Fit$var_ms_unweighted
cov0 <- summary(Fit)$cov.scaled 
MCMC <- modMCMC(f=R14tCost, p = Fit$par, niter = 1000, jump = NULL, var0 = var0, wvar0 = 0, 
                lower=c(-3,-3,-1,0,0),upper=c(0,0,0,1,1))

#The sensitivity range is calculated from the output of the MCMC
sR=sensRange(func=Mod1, parInput=MCMC$par)



###################################################
### code chunk number 20: ParameterEstimation.Rnw:313-314
###################################################
pairs(MCMC,nsample=500)


###################################################
### code chunk number 21: ParameterEstimation.Rnw:323-328
###################################################
par(mar=c(5,5,4,1))
plot(summary(sR),xlim=c(1950,2010),ylim=c(0,1000),xlab="Year",
     ylab=expression(paste(Delta^14,"C ","(\u2030)")),main="")
points(DataR14t,pch=20)
lines(C14Atm_NH,col=4)