-
Notifications
You must be signed in to change notification settings - Fork 9
/
metab.kalman.Rd
109 lines (85 loc) · 5.23 KB
/
metab.kalman.Rd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/metab.kalman.R
\name{metab.kalman}
\alias{metab.kalman}
\title{Metabolism calculated from parameters estimated using a Kalman filter}
\usage{
metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...)
}
\arguments{
\item{do.obs}{Vector of dissovled oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L}}
\item{do.sat}{Vector of dissolved oxygen saturation values based on water temperature. Calculate using \link{o2.at.sat}}
\item{k.gas}{Vector of kGAS values calculated from any of the gas flux models
(e.g., \link{k.cole}) and converted to kGAS using \link{k600.2.kGAS}}
\item{z.mix}{Vector of mixed-layer depths in meters. To calculate, see \link{ts.meta.depths}}
\item{irr}{Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s}}
\item{wtr}{Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature}
\item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}}
}
\value{
A data.frame with columns corresponding to components of metabolism
\describe{
\item{GPP}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}}
\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}}
\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}}
}
Use \link{attributes} to access more model output:
\item{smoothDO}{smoothed time series of oxygen concentration (\eqn{mg O[2] L^{-1}}{mg O2 / L}), from Kalman smoother}
\item{params}{parameters estimated by the Kalman filter (\eqn{c_1, c_2, Q, H}{c1, c2, Q, H})}
}
\description{
A state space model accounting for process and observation error, with the maximum likelihood of parameteres estimated using a Kalman filter.
Also provides a smoothed time series of oxygen concentration.
}
\details{
The model has four parameters, \eqn{c_1, c_2, Q, H}{c1, c2, Q, H}, and consists of equations involving the prediction of upcoming state conditional on information of the previous state (\eqn{a_{t|t-1}}{a[t|t-1]}, \eqn{P_{t|t-1}}{P[t|t-1]}), as well as updates of those predictions that are conditional upon information of the current state (\eqn{a_{t|t}}{a[t|t]}, \eqn{P_{t|t}}{P[t|t]}). \eqn{a} is the
\deqn{v=k.gas/z.mix}{v=k.gas/z.mix}
\deqn{a_t = c_1*irr_{t-1} + c_2*log_e(wtr_{t-1}) + v_{t-1}*do.sat_{t-1}}{a[t] = c1*irr[t-1] + c2*log(wtr[t-1]) + v[t-1]*do.sat[t-1]}
\deqn{beta = e^{-v}}{beta = exp(-v)}
\deqn{do.obs_t = a_t/v_{t-1} + -e^{-v_{t-1}}*a_t/v_{t-1} + beta_{t-1}*do.obs_{t-1} + epsilon_t}{do.obs[t] = a[t]/v[t-1] + -exp(-v[t-1])*a[t]/v[t-1] + beta[t-1]*do.obs[t-1] + epsilon[t]}
The above model is used during model fitting, but if gas flux is not integrated between time steps, those equations simplify to the following:
\deqn{F_{t-1} = k.gas_{t-1}*(do.sat_{t-1} - do.obs_{t-1})/z.mix_{t-1}}{F[t-1] = k.gas[t-1]*(do.sat[t-1] - do.obs[t-1])/z.mix[t-1]}
\deqn{do.obs_t=do.obs_{t-1}+c_1*irr_{t-1}+c_2*log_e(wtr_{t-1}) + F_{t-1} + epsilon_t}{do.obs[t] = do.obs[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + F[t-1] + epsilon[t]}
The parameters are fit using maximum likelihood, and the optimization (minimization of the negative log likelihood function) is performed by \code{optim} using default settings.
GPP is then calculated as \code{mean(c1*irr, na.rm=TRUE)*freq}, where \code{freq} is the number of observations per day, as estimated from the typical size between time steps. Thus, generally \code{freq==length(do.obs)}.
Similarly, R is calculated as \code{mean(c2*log(wtr), na.rm=TRUE)*freq}.
NEP is the sum of GPP and R.
}
\note{
If observation error is substantial, consider applying a Kalman filter to the water temperature time series by supplying
\code{wtr} as the output from \link{temp.kalman}
}
\examples{
library(rLakeAnalyzer)
doobs <- load.ts(system.file('extdata',
'sparkling.doobs', package="LakeMetabolizer"))
wtr <- load.ts(system.file('extdata',
'sparkling.wtr', package="LakeMetabolizer"))
wnd <- load.ts(system.file('extdata',
'sparkling.wnd', package="LakeMetabolizer"))
irr <- load.ts(system.file('extdata',
'sparkling.par', package="LakeMetabolizer"))
#Subset a day
Sys.setenv(TZ='GMT')
mod.date <- as.POSIXct('2009-07-08', 'GMT')
doobs <- doobs[trunc(doobs$datetime, 'day') == mod.date, ]
wtr <- wtr[trunc(wtr$datetime, 'day') == mod.date, ]
wnd <- wnd[trunc(wnd$datetime, 'day') == mod.date, ]
irr <- irr[trunc(irr$datetime, 'day') == mod.date, ]
k600 <- k.cole.base(wnd[,2])
k.gas <- k600.2.kGAS.base(k600, wtr[,3], 'O2')
do.sat <- o2.at.sat.base(wtr[,3], altitude=300)
metab.kalman(irr=irr[,2], z.mix=rep(1, length(k.gas)),
do.sat=do.sat, wtr=wtr[,2],
k.gas=k.gas, do.obs=doobs[,2])
}
\references{
Batt, Ryan D. and Stephen R. Carpenter. 2012. \emph{Free-water lake metabolism:
addressing noisy time series with a Kalman filter}. Limnology and Oceanography: Methods 10: 20-30. doi: 10.4319/lom.2012.10.20
}
\seealso{
\link{temp.kalman}, \link{watts.in}, \link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.mle}, \link{metab.bayesian}
}
\author{
Ryan Batt, Luke A. Winslow
}