gasfluxes (version 0.4)

HMR.fit: HMR fit

Description

Fit the HMR model using the Golub-Pereyra algorithm for partially linear least-squares models.

Usage

HMR.fit(t, C, A = 1, V, serie = "", k = log(1.5), verbose = TRUE,
  plot = FALSE, maxiter = 100, ...)

Arguments

t

time values (usually in hours)

C

concentration values

A

area covered by the chamber

V

effective volume of the chamber

serie

id of the flux measurement

k

starting value for nls function

verbose

logical, TRUE prints message after each flux calculation

plot

logical, mainly intended for use in gasfluxes

maxiter

see nls.control

…

further parameters, currently none

Value

A list of

f0

flux estimate

f0.se

standard error of flux estimate

f0.p

p-value of flux estimate

kappa, phi

other parameters of the HMR model

AIC

Akaike information criterion

AICc

Akaike information criterion with small sample correction

RSE

residual standard error (sigma from summary.nls)

diagnostics

error or warning messages

Details

The HMR model (Pedersen et al., 2010) is \(C(t)=\phi+f_0 \frac{e^{-\kappa t}}{-\kappa V/A}\). To ensure the lower bound \(\kappa > 0\), the substitution \(\kappa = e^k\) is used. The resulting reparameterized model is then fit using nls with algorithm = "plinear". This is computationally more efficient than the manual implementation in the HMR package and results in almost identical flux values. Flux standard errors and p-values differ strongly from those reported by the HMR package <= version 0.3.1, but are equal to those reported by later versions.

The default starting value \(k = log(\kappa)\) assumes that time is in hours. If you use a different time unit, you should adjust it accordingly.

There have been demands to return the initial concentration as predicted by the model as this is useful for checking plausibility. However, this can be easily calculated from the parameters and the equation of the model by setting t = 0, i.e., \(C_0=\phi+f_0\).

Note that nls is used internally and thus this function should not be used with artificial "zero-residual" data.

References

Pedersen, A.R., Petersen, S.O., Schelde, K., 2010. A comprehensive approach to soil-atmosphere trace-gas flux estimation with static chambers. European Journal of Soil Science 61(6), 888-902.

Examples

# NOT RUN {
#a single fit
t <- c(0, 1/3, 2/3, 1)
C <- c(320, 341, 352, 359)
print(fit <- HMR.fit(t, C, 1, 0.3, "a"))
plot(C ~ t)
curve({fit$phi + fit$f0 * exp(-fit$kappa * x)/(-fit$kappa*0.3)}, 
      from = 0, to = 1, add = TRUE)
#compare with fitting function from HMR package 0.3.1
  gasfluxes:::.HMR.fit1(t, C,
                  1, 0.3, "a",
                  ngrid = 1000, LR.always = FALSE, FollowHMR = TRUE, 
                  JPG = FALSE, PS = FALSE, PHMR = FALSE, npred = 500, 
                  xtxt = "", ytxt = "", pcttxt = "",
                  MSE.zero = 10 * max(.Machine$double.eps, .Machine$double.neg.eps), 
                  bracketing.tol = 1e-07, bracketing.maxiter = 1000)[2:4]

# }
# NOT RUN {
#a dataset of 1329 chamber N2O flux measurements
data(fluxMeas)
fluxMeas[, n := length(time), by=serie]
print(fluxMeas)
fluxes <- fluxMeas[n > 3, HMR.fit(time, C, A, V, serie), by=serie]
print(fluxes)
plot(f0.se ~ f0, data = fluxes)
#one very large f0.se value (and several infinite ones not shown in the plot)
fluxes[is.finite(f0.se),][which.max(f0.se),]
plot(C~time, data=fluxMeas[serie=="ID940",])
print(tmp <- fluxes[is.finite(f0.se),][which.max(f0.se),])
curve({tmp[, phi] + tmp[, f0] * exp(-tmp[, kappa] * x)/
      (-tmp[, kappa]*fluxMeas[serie=="ID940", V[1]]/
      fluxMeas[serie=="ID940",A[1]])}, 
      from = 0, to = 1, add = TRUE)
plot(f0.se ~ f0, data = fluxes[f0.se < 1e4,], pch = 16)
boxplot(fluxes[f0.se < 1e4, sqrt(f0.se)])
# }
# NOT RUN {
# }