sdreport: General sdreport function.
In TMB: Template Model Builder: A General Random Effect Tool Inspired by 'ADMB'
sdreport R Documentation
General sdreport function.
Description
After optimization of an AD model, sdreport is used to
calculate standard deviations of all model parameters, including
non linear functions of random effects and parameters specified
through the ADREPORT() macro from the user template.
Usage
sdreport(
obj,
par.fixed = NULL,
hessian.fixed = NULL,
getJointPrecision = FALSE,
bias.correct = FALSE,
bias.correct.control = list(sd = FALSE, split = NULL, nsplit = NULL),
ignore.parm.uncertainty = FALSE,
getReportCovariance = TRUE,
skip.delta.method = FALSE
)
Arguments
obj
Object returned by MakeADFun
par.fixed
Optional. Parameter estimate (will be known to obj when an optimization has been carried out).
hessian.fixed
Optional. Hessian wrt. parameters (will be calculated from obj if missing).
getJointPrecision
Optional. Return full joint precision matrix of random effects and parameters?
bias.correct
logical indicating if bias correction should be applied
bias.correct.control
a list of bias correction options; currently sd, split and nsplit are used - see details.
ignore.parm.uncertainty
Optional. Ignore estimation variance of parameters?
getReportCovariance
Get full covariance matrix of ADREPORTed variables?
skip.delta.method
Skip the delta method? (FALSE by default)
Details
First, the Hessian wrt. the parameter vector (\theta) is
calculated. The parameter covariance matrix is approximated by
V(\hat\theta)=-\nabla^2 l(\hat\theta)^{-1}
where l
denotes the log likelihood function (i.e. -obj$fn). If
ignore.parm.uncertainty=TRUE then the Hessian calculation
is omitted and a zero-matrix is used in place of
V(\hat\theta).
For non-random effect models the standard delta-method is used to
calculate the covariance matrix of transformed parameters. Let
\phi(\theta) denote some non-linear function of
\theta. Then
V(\phi(\hat\theta))\approx \nabla\phi
V(\hat\theta) \nabla\phi'
The covariance matrix of reported variables
V(\phi(\hat\theta)) is returned by default. This can cause
high memory usage if many variables are ADREPORTed. Use
getReportCovariance=FALSE to only return standard errors.
In case standard deviations are not required one can completely skip
the delta method using skip.delta.method=TRUE.
For random effect models a generalized delta-method is used. First
the joint covariance of random effect and parameter estimation error is approximated
by
V \left( \begin{array}{cc} \hat u - u \cr \hat\theta - \theta \end{array} \right) \approx
\left( \begin{array}{cc} H_{uu}^{-1} & 0 \cr 0 & 0 \end{array} \right) +
J V(\hat\theta) J'
where H_{uu} denotes random effect block of the full joint
Hessian of obj$env$f and J denotes the Jacobian of
\left( \begin{array}{cc}\hat u(\theta) \cr \theta \end{array} \right) wrt. \theta.
Here, the first term represents the expected conditional variance
of the estimation error given the data and the second term represents the variance
of the conditional mean of the estimation error given the data.
Now the delta method can be applied on a general non-linear
function \phi(u,\theta) of random effects u and
parameters \theta:
V\left(\phi(\hat u,\hat\theta) - \phi(u,\theta) \right)\approx \nabla\phi V \left( \begin{array}{cc}
\hat u - u \cr \hat\theta - \theta \end{array} \right) \nabla\phi'
The full joint covariance is not returned by default, because it
may require large amounts of memory. It may be obtained by
specifying getJointPrecision=TRUE, in which case V
\left( \begin{array}{cc} \hat u - u \cr \hat\theta - \theta \end{array} \right) ^{-1} will be part of the
output. This matrix must be manually inverted using
solve(jointPrecision) in order to get the joint covariance
matrix. Note, that the parameter order will follow the original
order (i.e. obj$env$par).
Using \phi(\hat u,\theta) as estimator of
\phi(u,\theta) may result in substantial bias. This may be
the case if either \phi is non-linear or if the distribution
of u given x (data) is sufficiently non-symmetric. A
generic correction is enabled with bias.correct=TRUE. It is
based on the identity
E_{\theta}[\phi(u,\theta)|x] =
\partial_\varepsilon\left(\log \int \exp(-f(u,\theta) +
\varepsilon \phi(u,\theta))\:du\right)_{|\varepsilon=0}
stating that the conditional expectation can be written as a
marginal likelihood gradient wrt. a nuisance parameter
\varepsilon.
The marginal likelihood is replaced by its Laplace approximation.
If bias.correct.control$sd=TRUE the variance of the
estimator is calculated using
V_{\theta}[\phi(u,\theta)|x] =
\partial_\varepsilon^2\left(\log \int \exp(-f(u,\theta) +
\varepsilon \phi(u,\theta))\:du\right)_{|\varepsilon=0}
A further correction is added to this variance to account for the
effect of replacing \theta by the MLE \hat\theta
(unless ignore.parm.uncertainty=TRUE).
Bias correction can be be performed in chunks in order to reduce
memory usage or in order to only bias correct a subset of
variables. First option is to pass a list of indices as
bias.correct.control$split. E.g. a list
list(1:2,3:4) calculates the first four ADREPORTed
variables in two chunks.
The internal function obj$env$ADreportIndex()
gives an overview of the possible indices of ADREPORTed variables.
Second option is to pass the number of
chunks as bias.correct.control$nsplit in which case all
ADREPORTed variables are bias corrected in the specified number of
chunks.
Also note that skip.delta.method may be necessary when bias
correcting a large number of variables.
Value
Object of class sdreport
See Also
summary.sdreport, print.sdreport, as.list.sdreport
Examples
## Not run:
runExample("linreg_parallel", thisR = TRUE) ## Non-random effect example
sdreport(obj)
## End(Not run)
runExample("simple", thisR = TRUE) ## Random effect example
rep <- sdreport(obj)
summary(rep, "random") ## Only random effects
summary(rep, "fixed", p.value = TRUE) ## Only non-random effects
summary(rep, "report") ## Only report
## Bias correction
rep <- sdreport(obj, bias.correct = TRUE)
summary(rep, "report") ## Include bias correction
TMB documentation built on Sept. 11, 2024, 7:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| sdreport | R Documentation |
General sdreport function.
Description
After optimization of an AD model, sdreport is used to
calculate standard deviations of all model parameters, including
non linear functions of random effects and parameters specified
through the ADREPORT() macro from the user template.
Usage
sdreport(
obj,
par.fixed = NULL,
hessian.fixed = NULL,
getJointPrecision = FALSE,
bias.correct = FALSE,
bias.correct.control = list(sd = FALSE, split = NULL, nsplit = NULL),
ignore.parm.uncertainty = FALSE,
getReportCovariance = TRUE,
skip.delta.method = FALSE
)
Arguments
obj |
Object returned by |
par.fixed |
Optional. Parameter estimate (will be known to |
hessian.fixed |
Optional. Hessian wrt. parameters (will be calculated from |
getJointPrecision |
Optional. Return full joint precision matrix of random effects and parameters? |
bias.correct |
logical indicating if bias correction should be applied |
bias.correct.control |
a |
ignore.parm.uncertainty |
Optional. Ignore estimation variance of parameters? |
getReportCovariance |
Get full covariance matrix of ADREPORTed variables? |
skip.delta.method |
Skip the delta method? ( |
Details
First, the Hessian wrt. the parameter vector (\theta) is
calculated. The parameter covariance matrix is approximated by
V(\hat\theta)=-\nabla^2 l(\hat\theta)^{-1}
where l
denotes the log likelihood function (i.e. -obj$fn). If
ignore.parm.uncertainty=TRUE then the Hessian calculation
is omitted and a zero-matrix is used in place of
V(\hat\theta).
For non-random effect models the standard delta-method is used to
calculate the covariance matrix of transformed parameters. Let
\phi(\theta) denote some non-linear function of
\theta. Then
V(\phi(\hat\theta))\approx \nabla\phi
V(\hat\theta) \nabla\phi'
The covariance matrix of reported variables
V(\phi(\hat\theta)) is returned by default. This can cause
high memory usage if many variables are ADREPORTed. Use
getReportCovariance=FALSE to only return standard errors.
In case standard deviations are not required one can completely skip
the delta method using skip.delta.method=TRUE.
For random effect models a generalized delta-method is used. First the joint covariance of random effect and parameter estimation error is approximated by
V \left( \begin{array}{cc} \hat u - u \cr \hat\theta - \theta \end{array} \right) \approx
\left( \begin{array}{cc} H_{uu}^{-1} & 0 \cr 0 & 0 \end{array} \right) +
J V(\hat\theta) J'
where H_{uu} denotes random effect block of the full joint
Hessian of obj$env$f and J denotes the Jacobian of
\left( \begin{array}{cc}\hat u(\theta) \cr \theta \end{array} \right) wrt. \theta.
Here, the first term represents the expected conditional variance
of the estimation error given the data and the second term represents the variance
of the conditional mean of the estimation error given the data.
Now the delta method can be applied on a general non-linear
function \phi(u,\theta) of random effects u and
parameters \theta:
V\left(\phi(\hat u,\hat\theta) - \phi(u,\theta) \right)\approx \nabla\phi V \left( \begin{array}{cc}
\hat u - u \cr \hat\theta - \theta \end{array} \right) \nabla\phi'
The full joint covariance is not returned by default, because it
may require large amounts of memory. It may be obtained by
specifying getJointPrecision=TRUE, in which case V
\left( \begin{array}{cc} \hat u - u \cr \hat\theta - \theta \end{array} \right) ^{-1} will be part of the
output. This matrix must be manually inverted using
solve(jointPrecision) in order to get the joint covariance
matrix. Note, that the parameter order will follow the original
order (i.e. obj$env$par).
Using \phi(\hat u,\theta) as estimator of
\phi(u,\theta) may result in substantial bias. This may be
the case if either \phi is non-linear or if the distribution
of u given x (data) is sufficiently non-symmetric. A
generic correction is enabled with bias.correct=TRUE. It is
based on the identity
E_{\theta}[\phi(u,\theta)|x] =
\partial_\varepsilon\left(\log \int \exp(-f(u,\theta) +
\varepsilon \phi(u,\theta))\:du\right)_{|\varepsilon=0}
stating that the conditional expectation can be written as a
marginal likelihood gradient wrt. a nuisance parameter
\varepsilon.
The marginal likelihood is replaced by its Laplace approximation.
If bias.correct.control$sd=TRUE the variance of the
estimator is calculated using
V_{\theta}[\phi(u,\theta)|x] =
\partial_\varepsilon^2\left(\log \int \exp(-f(u,\theta) +
\varepsilon \phi(u,\theta))\:du\right)_{|\varepsilon=0}
A further correction is added to this variance to account for the
effect of replacing \theta by the MLE \hat\theta
(unless ignore.parm.uncertainty=TRUE).
Bias correction can be be performed in chunks in order to reduce
memory usage or in order to only bias correct a subset of
variables. First option is to pass a list of indices as
bias.correct.control$split. E.g. a list
list(1:2,3:4) calculates the first four ADREPORTed
variables in two chunks.
The internal function obj$env$ADreportIndex()
gives an overview of the possible indices of ADREPORTed variables.
Second option is to pass the number of
chunks as bias.correct.control$nsplit in which case all
ADREPORTed variables are bias corrected in the specified number of
chunks.
Also note that skip.delta.method may be necessary when bias
correcting a large number of variables.
Value
Object of class sdreport
See Also
summary.sdreport, print.sdreport, as.list.sdreport
Examples
## Not run:
runExample("linreg_parallel", thisR = TRUE) ## Non-random effect example
sdreport(obj)
## End(Not run)
runExample("simple", thisR = TRUE) ## Random effect example
rep <- sdreport(obj)
summary(rep, "random") ## Only random effects
summary(rep, "fixed", p.value = TRUE) ## Only non-random effects
summary(rep, "report") ## Only report
## Bias correction
rep <- sdreport(obj, bias.correct = TRUE)
summary(rep, "report") ## Include bias correction
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.