General sdreport function.

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)

Value

Object of class sdreport

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.

Examples

if (FALSE) { # \dontrun{
runExample("linreg_parallel", thisR = TRUE) ## Non-random effect example
sdreport(obj) } # }

runExample("simple", thisR = TRUE)          ## Random effect example
#> Running example simple 
#> 
#> > require(TMB)
#> 
#> > dyn.load(dynlib("simple"))
#> 
#> > set.seed(123)
#> 
#> > y <- rep(1900:2010, each = 2)
#> 
#> > year <- factor(y)
#> 
#> > quarter <- factor(rep(1:4, length.out = length(year)))
#> 
#> > period <- factor((y > mean(y)) + 1)
#> 
#> > B <- model.matrix(~year + quarter - 1)
#> 
#> > A <- model.matrix(~period - 1)
#> 
#> > B <- as(B, "TsparseMatrix")
#> 
#> > A <- as(A, "TsparseMatrix")
#> 
#> > u <- rnorm(ncol(B))
#> 
#> > beta <- rnorm(ncol(A)) * 100
#> 
#> > eps <- rnorm(nrow(B), sd = 1)
#> 
#> > x <- as.numeric(A %*% beta + B %*% u + eps)
#> 
#> > obj <- MakeADFun(data = list(x = x, B = B, A = A), 
#> +     parameters = list(u = u * 0, beta = beta * 0, logsdu = 1, 
#> +         logsd0 = 1), random = .... [TRUNCATED] 
#> 
#> > opt <- nlminb(obj$par, obj$fn, obj$gr)
rep <- sdreport(obj)
summary(rep, "random")                      ## Only random effects
#>       Estimate Std. Error
#> u -0.480820910  0.5690955
#> u -0.716733126  0.5725030
#> u  0.657193869  0.5712858
#> u -0.217915665  0.5656616
#> u  0.415627618  0.5674558
#> u  1.037127885  0.5825155
#> u -0.035074217  0.5651449
#> u -0.226121203  0.5657118
#> u -0.509386497  0.5695582
#> u -0.566398890  0.5696368
#> u  0.495558060  0.5685250
#> u  0.298769093  0.5667365
#> u  0.355212228  0.5667787
#> u -0.870040865  0.5761263
#> u -0.939119181  0.5794885
#> u  0.764860482  0.5747977
#> u  0.519021610  0.5688762
#> u -0.728929771  0.5727655
#> u  0.206713961  0.5655963
#> u -0.431402685  0.5676512
#> u -0.564252361  0.5705169
#> u  0.004885871  0.5651013
#> u -0.594084588  0.5710766
#> u  0.107808819  0.5653671
#> u -0.397921289  0.5678945
#> u -0.972528543  0.5789366
#> u  0.736209122  0.5729243
#> u  0.002062473  0.5650987
#> u -0.055435165  0.5651904
#> u  0.150285655  0.5655734
#> u  0.343688815  0.5666624
#> u -0.219748682  0.5656726
#> u  0.571349148  0.5697202
#> u  0.113003595  0.5653893
#> u  0.346144329  0.5666868
#> u  0.699051183  0.5732600
#> u  0.390264636  0.5671578
#> u -0.110690270  0.5651981
#> u -0.410617606  0.5680647
#> u -0.064652817  0.5651094
#> u -0.591632297  0.5710295
#> u -0.681611393  0.5717720
#> u  0.279981034  0.5660939
#> u  1.288904634  0.5915390
#> u  0.437227319  0.5677253
#> u -1.090509371  0.5825501
#> u  0.225347782  0.5657070
#> u -0.207742015  0.5656022
#> u  0.653948462  0.5712226
#> u  0.105889628  0.5653591
#> u -0.192711036  0.5658356
#> u -0.787114852  0.5740793
#> u -0.014988335  0.5651127
#> u  1.271209265  0.5908473
#> u -0.220950148  0.5660413
#> u  0.421785211  0.5682186
#> u -0.964761102  0.5802521
#> u  0.802282866  0.5759681
#> u -0.112840062  0.5653996
#> u -0.494492438  0.5687318
#> u  0.254116357  0.5659156
#> u -0.495172714  0.5687417
#> u -0.597921307  0.5711543
#> u -0.557653387  0.5697122
#> u -0.154392475  0.5656067
#> u -0.319684402  0.5666591
#> u  0.173827384  0.5654428
#> u -0.313252784  0.5666009
#> u  0.426042971  0.5676030
#> u  1.545625872  0.6027581
#> u  0.177400280  0.5654596
#> u -1.881197574  0.6167924
#> u  0.255603469  0.5659263
#> u -0.093814064  0.5653875
#> u  0.678121255  0.5717240
#> u -0.026716735  0.5653169
#> u -0.283918237  0.5665978
#> u -0.256237237  0.5661412
#> u -0.133639832  0.5654966
#> u -0.037151967  0.5653186
#> u  0.042171530  0.5651036
#> u  0.279731843  0.5667837
#> u -0.515695888  0.5696688
#> u  0.333327537  0.5673407
#> u -0.087592533  0.5653002
#> u  0.655464007  0.5725508
#> u  0.585214670  0.5699792
#> u  0.819229711  0.5764038
#> u -0.057372431  0.5652073
#> u  1.186836748  0.5879069
#> u  0.689177729  0.5719491
#> u  0.609944946  0.5716226
#> u -0.318936507  0.5669554
#> u -1.010274399  0.5802594
#> u  0.721604204  0.5726304
#> u -0.939828995  0.5782188
#> u  2.174171499  0.6332481
#> u  1.405106972  0.5965310
#> u -0.841017389  0.5767380
#> u -0.564282134  0.5698222
#> u -0.836831688  0.5766270
#> u  0.824039254  0.5765290
#> u  0.244541199  0.5658487
#> u -0.469535770  0.5683781
#> u -0.535926363  0.5700152
#> u  0.065805804  0.5654506
#> u -0.098619340  0.5653412
#> u -1.085273582  0.5825898
#> u -0.022389003  0.5651355
#> u -0.026227891  0.5653169
#> u -0.816740910  0.5761013
#> u  0.472211066  0.1884481
#> u -1.530647054  0.2395324
#> u -0.139249811  0.2388424
summary(rep, "fixed", p.value = TRUE)       ## Only non-random effects
#>           Estimate Std. Error    z value Pr(>|z^2|)
#> beta   52.01370232  0.2016859 257.894548  0.0000000
#> beta   30.24058534  0.2018362 149.827373  0.0000000
#> logsdu -0.15777145  0.1201260  -1.313383  0.1890539
#> logsd0  0.03326068  0.0667396   0.498365  0.6182268
summary(rep, "report")                      ## Only report
#>              Estimate  Std. Error
#> sd0           1.03382  0.06899673
#> sum(exp(u)) 141.71087 15.59821443

## Bias correction
rep <- sdreport(obj, bias.correct = TRUE)
summary(rep, "report")                      ## Include bias correction
#>              Estimate  Std. Error Est. (bias.correct) Std. (bias.correct)
#> sd0           1.03382  0.06899673             1.03382                  NA
#> sum(exp(u)) 141.71087 15.59821443           163.51245                  NA

Usage

Arguments

Value

Details

See also

Examples