fitKinrespExperiment: fitKinrespExperiment
In bgctw/twKinresp: Fitting kinetic models to microbial respiration data

Description Usage Arguments Details Value Author(s) Examples

Fit microbial form to several replicated respiration time series of one experiment.

fitKinrespExperiment(rde.e, repFits, lambda = 0.9, 
    YCO2 = 1.5, weights = NULL, start = NULL, 
    tmp.names = c("none", "x0l", "r0l", "x0l+r0l", 
        "x0l+r0l+mumaxl"), showFitErrorMsg = FALSE)

`rde.e`	dataset with columns experiment, replicate, resp and time,
`repFits`	Initial coefficients mumax, x0, and r0 for each replicate
`lambda`	Ratio of growth associated (coupled) specific respiration
`YCO2`	Ratio of assimilated carbon per respired carbon (Y/(1-Y)).
`weights`	Variance function. see details
`start`	Specifying starting values in an alternative way to
`tmp.names`	scenarios
`showFitErrorMsg`	if FALSE (standard) then error msg are

Microbial parameters across replicates

Microbial parameters can be inferred from fitting against the time series of respiration of each replicate separately (fitKinrespReplicate. However, what are then parameters of the population? The average across replicate parameters will be wrong.

Alternatively, one can first average the measurement across replicates for each measurement time. However, the uncertainty of the parameters will be wrong.

A viable solution is to fit a mixed model to all the replicate data. The micoribal parameters, then are described by a mean value of the population and a variance across replicates. Several options of which parameters vary across replicates can be tested.

Variance function: If no weights are given, the measurement errors of the single observations are assumed to be identical. If measurement errors increase with the magnitude of the observations, this can be modeled by power variance function: weights=varPower(fixed=delta), with delta being a value between 0 (constant expected absolute error) and 1 (constant expected relative error).

list with components

model: best fitting result
random: best random effects scenario
fits: fits of all random effects scenarios
aics: Akaike information criterion for random effects scenarios

Thomas Wutzler <thomas.wutzler@web.de>

# data of one example treament: measurements of several replicates of one soil
# data(respWutzler10)
rde <- subset(respWutzler10, suite == "Face" & experiment == 9 )

# constrain data to unlimited growth phase
res4 <- kinrespGrowthphaseExperiment(rde, weights = varPower(fixed = 0.5) )
rde.e <- getUnlimitedGrowthData(res4)
# fit the mixed model to all replicates
res5Scen <- fitKinrespExperiment(
  rde.e, coefKinresp(res4,rde.e), weights = varPower(fixed = 0.5)
  ,showFitErrorMsg = TRUE
  )

# get the best fit parameters of the population
coefKinresp(fixef(res5Scen$model))

# examine the random-effects scenarios:
# Here the lowest AIC suggest that activity state and initial microbial biomass
# differed between replicates, but maximum growth was the same
res5Scen$aics

# plot the fits
#windows(record = TRUE)
rde.e$fitted <- fitted(res5Scen$model)
plot( resp ~ time, data = rde.e, col = rde.e$replicate )
tmp <- by( rde.e, rde.e$replicate, function(rder){
  lines(fitted~time,data = rder, col = as.numeric(rder[1,replicate]))})

# estimated microbial coefficients
(pars <- kinrespParDist(res5Scen$model))
# plot the density of r0 and density summaries
iPar = "r0"
xGrid <- seq( pars[iPar,"cf025"]*0.8, pars[iPar,"cf975"]*1.2, length.out = 80)
#fx <- dlnorm(xGrid, mean = pars[iPar,"mu"],sd = pars[iPar,"sigma"])
fx <- dlogitnorm(xGrid, mu = pars[iPar,"mu"],sigma = pars[iPar,"sigma"])
plot( fx ~ xGrid, type = "l", xlab = iPar, ylab = "density" )
abline(v = pars[iPar,c("mle","median","mean","cf025","cf975")]
       , col = c("red","green","blue","gray","gray"))