jranke / mkin
1
#' Summary method for class "mkinfit"
2
#'
3
#' Lists model equations, initial parameter values, optimised parameters with
4
#' some uncertainty statistics, the chi2 error levels calculated according to
5
#' FOCUS guidance (2006) as defined therein, formation fractions, DT50 values
6
#' and optionally the data, consisting of observed, predicted and residual
7
#' values.
8
#'
9
#' @param object an object of class \code{\link{mkinfit}}.
10
#' @param x an object of class \code{summary.mkinfit}.
11
#' @param data logical, indicating whether the data should be included in the
12
#'   summary.
13
#' @param distimes logical, indicating whether DT50 and DT90 values should be
14
#'   included.
15
#' @param alpha error level for confidence interval estimation from t
16
#'   distribution
17
#' @param digits Number of digits to use for printing
18
#' @param \dots optional arguments passed to methods like \code{print}.
19
#' @importFrom stats qt pt cov2cor
20
#' @return The summary function returns a list with components, among others
21
#'   \item{version, Rversion}{The mkin and R versions used}
22
#'   \item{date.fit, date.summary}{The dates where the fit and the summary were
23
#'     produced}
24
#'   \item{diffs}{The differential equations used in the model}
25
#'   \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
26
#'   \item{bpar}{Optimised and backtransformed
27
#'     parameters}
28
#'   \item{data}{The data (see Description above).}
29
#'   \item{start}{The starting values and bounds, if applicable, for optimised
30
#'     parameters.}
31
#'   \item{fixed}{The values of fixed parameters.}
32
#'   \item{errmin }{The chi2 error levels for
33
#'     each observed variable.}
34
#'   \item{bparms.ode}{All backtransformed ODE
35
#'     parameters, for use as starting parameters for related models.}
36
#'   \item{errparms}{Error model parameters.}
37
#'   \item{ff}{The estimated formation fractions derived from the fitted
38
#'      model.}
39
#'   \item{distimes}{The DT50 and DT90 values for each observed variable.}
40
#'   \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
41
#'   The print method is called for its side effect, i.e. printing the summary.
42
#' @author Johannes Ranke
43
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
44
#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
45
#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
46
#'   EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
47
#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
48
#' @examples
49
#'
50
#'   summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE))
51
#'
52
#' @export
53
summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
54 1
  param  <- object$par
55 1
  pnames <- names(param)
56 1
  bpnames <- names(object$bparms.optim)
57 1
  epnames <- names(object$errparms)
58 1
  p      <- length(param)
59 1
  mod_vars <- names(object$mkinmod$diffs)
60 1
  covar  <- try(solve(object$hessian), silent = TRUE)
61 1
  covar_notrans  <- try(solve(object$hessian_notrans), silent = TRUE)
62 1
  rdf <- object$df.residual
63

64 1
  if (!is.numeric(covar) | is.na(covar[1])) {
65 0
    covar <- NULL
66 0
    se <- lci <- uci <- rep(NA, p)
67
  } else {
68 1
    rownames(covar) <- colnames(covar) <- pnames
69 1
    se     <- sqrt(diag(covar))
70 1
    lci    <- param + qt(alpha/2, rdf) * se
71 1
    uci    <- param + qt(1-alpha/2, rdf) * se
72
  }
73

74 1
  beparms.optim <- c(object$bparms.optim, object$par[epnames])
75 1
  if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) {
76 1
    covar_notrans <- NULL
77 1
    se_notrans <- tval <- pval <- rep(NA, p)
78
  } else {
79 1
    rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)
80 1
    se_notrans <- sqrt(diag(covar_notrans))
81 1
    tval  <- beparms.optim / se_notrans
82 1
    pval  <- pt(abs(tval), rdf, lower.tail = FALSE)
83
  }
84

85 1
  names(se) <- pnames
86

87 1
  param <- cbind(param, se, lci, uci)
88 1
  dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
89

90 1
  bparam <- cbind(Estimate = beparms.optim, se_notrans,
91 1
                  "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)
92

93
  # Transform boundaries of CI for one parameter at a time,
94
  # with the exception of sets of formation fractions (single fractions are OK).
95 1
  f_names_skip <- character(0)
96 1
  for (box in mod_vars) { # Figure out sets of fractions to skip
97 1
    f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
98 1
    n_paths <- length(f_names)
99 1
    if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
100
  }
101

102 1
  for (pname in pnames) {
103 1
    if (!pname %in% f_names_skip) {
104 1
      par.lower <- param[pname, "Lower"]
105 1
      par.upper <- param[pname, "Upper"]
106 1
      names(par.lower) <- names(par.upper) <- pname
107 1
      bpl <- backtransform_odeparms(par.lower, object$mkinmod,
108 1
                                            object$transform_rates,
109 1
                                            object$transform_fractions)
110 1
      bpu <- backtransform_odeparms(par.upper, object$mkinmod,
111 1
                                            object$transform_rates,
112 1
                                            object$transform_fractions)
113 1
      bparam[names(bpl), "Lower"] <- bpl
114 1
      bparam[names(bpu), "Upper"] <- bpu
115
    }
116
  }
117 1
  bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]
118

119 1
  ans <- list(
120 1
    version = as.character(utils::packageVersion("mkin")),
121 1
    Rversion = paste(R.version$major, R.version$minor, sep="."),
122 1
    date.fit = object$date,
123 1
    date.summary = date(),
124 1
    solution_type = object$solution_type,
125 1
    warnings = object$summary_warnings,
126 1
    use_of_ff = object$mkinmod$use_of_ff,
127 1
    error_model_algorithm = object$error_model_algorithm,
128 1
    df = c(p, rdf),
129 1
    covar = covar,
130 1
    covar_notrans = covar_notrans,
131 1
    err_mod = object$err_mod,
132 1
    niter = object$iterations,
133 1
    calls = object$calls,
134 1
    time = object$time,
135 1
    par = param,
136 1
    bpar = bparam)
137

138 1
  if (!is.null(object$version)) {
139 1
    ans$fit_version <- object$version
140 1
    ans$fit_Rversion <- object$Rversion
141 1
    if (ans$fit_version >= "0.9.49.6") {
142 1
      ans$AIC = AIC(object)
143 1
      ans$BIC = BIC(object)
144 1
      ans$logLik = logLik(object)
145
    }
146
  }
147

148 1
  ans$diffs <- object$mkinmod$diffs
149 1
  if(data) ans$data <- object$data
150 1
  ans$start <- object$start
151 1
  ans$start_transformed <- object$start_transformed
152

153 1
  ans$fixed <- object$fixed
154

155 1
  ans$errmin <- mkinerrmin(object, alpha = 0.05)
156

157 1
  if (object$calls > 0) {
158 1
    if (!is.null(ans$covar)){
159 1
      Corr <- cov2cor(ans$covar)
160 1
      rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
161 1
      ans$Corr <- Corr
162
    } else {
163 0
      warning("Could not calculate correlation; no covariance matrix")
164
    }
165
  }
166

167 1
  ans$bparms.ode <- object$bparms.ode
168 1
  ans$shapiro.p <- object$shapiro.p
169 1
  ep <- endpoints(object)
170 1
  if (length(ep$ff) != 0)
171 1
    ans$ff <- ep$ff
172 1
  if (distimes) ans$distimes <- ep$distimes
173 1
  if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
174 0
  if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
175 1
  class(ans) <- c("summary.mkinfit", "summary.modFit")
176 1
  return(ans)
177
}
178

179
#' @rdname summary.mkinfit
180
#' @export
181
print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
182 1
  if (is.null(x$fit_version)) {
183 0
    cat("mkin version:   ", x$version, "\n")
184 0
    cat("R version:      ", x$Rversion, "\n")
185
  } else {
186 1
    cat("mkin version used for fitting:   ", x$fit_version, "\n")
187 1
    cat("R version used for fitting:      ", x$fit_Rversion, "\n")
188
  }
189

190 1
  cat("Date of fit:    ", x$date.fit, "\n")
191 1
  cat("Date of summary:", x$date.summary, "\n")
192

193 1
  cat("\nEquations:\n")
194 1
  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
195 1
  writeLines(strwrap(nice_diffs, exdent = 11))
196 1
  df  <- x$df
197 1
  rdf <- df[2]
198

199 1
  cat("\nModel predictions using solution type", x$solution_type, "\n")
200

201 1
  cat("\nFitted using", x$calls, "model solutions performed in", x$time[["elapsed"]],  "s\n")
202

203 1
  if (!is.null(x$err_mod)) {
204 1
    cat("\nError model: ")
205 1
    cat(switch(x$err_mod,
206 1
               const = "Constant variance",
207 1
               obs = "Variance unique to each observed variable",
208 1
               tc = "Two-component variance function"), "\n")
209

210 1
    cat("\nError model algorithm:", x$error_model_algorithm, "\n")
211 0
    if (!is.null(x$d_3_message)) cat(x$d_3_message, "\n")
212
  }
213

214 1
  cat("\nStarting values for parameters to be optimised:\n")
215 1
  print(x$start)
216

217 1
  cat("\nStarting values for the transformed parameters actually optimised:\n")
218 1
  print(x$start_transformed)
219

220 1
  cat("\nFixed parameter values:\n")
221 1
  if(length(x$fixed$value) == 0) cat("None\n")
222 1
  else print(x$fixed)
223

224
  # We used to only have one warning - show this for summarising old objects
225 0
   if (!is.null(x[["warning"]])) cat("\n\nWarning:", x$warning, "\n\n")
226

227 1
  if (length(x$warnings) > 0) {
228 0
    cat("\n\nWarning(s):", "\n")
229 0
    cat(x$warnings, sep = "\n")
230
  }
231

232 1
  if (!is.null(x$AIC)) {
233 1
    cat("\nResults:\n\n")
234 1
    print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
235 1
      row.names = " "))
236
  }
237

238 1
  cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
239 1
  print(signif(x$par, digits = digits))
240

241 1
  if (x$calls > 0) {
242 1
    cat("\nParameter correlation:\n")
243 1
    if (!is.null(x$covar)){
244 1
      print(x$Corr, digits = digits, ...)
245
    } else {
246 0
      cat("No covariance matrix")
247
    }
248
  }
249

250 1
  cat("\nBacktransformed parameters:\n")
251 1
  cat("Confidence intervals for internally transformed parameters are asymmetric.\n")
252 1
  if ((x$version) < "0.9-36") {
253 0
    cat("To get the usual (questionable) t-test, upgrade mkin and repeat the fit.\n")
254 0
    print(signif(x$bpar, digits = digits))
255
  } else {
256 1
    cat("t-test (unrealistically) based on the assumption of normal distribution\n")
257 1
    cat("for estimators of untransformed parameters.\n")
258 1
    print(signif(x$bpar[, c(1, 3, 4, 5, 6)], digits = digits))
259
  }
260

261 1
  cat("\nFOCUS Chi2 error levels in percent:\n")
262 1
  x$errmin$err.min <- 100 * x$errmin$err.min
263 1
  print(x$errmin, digits=digits,...)
264

265 1
  printSFORB <- !is.null(x$SFORB)
266 1
  if(printSFORB){
267 1
    cat("\nEstimated Eigenvalues of SFORB model(s):\n")
268 1
    print(x$SFORB, digits=digits,...)
269
  }
270

271 1
  printff <- !is.null(x$ff)
272 1
  if(printff){
273 1
    cat("\nResulting formation fractions:\n")
274 1
    print(data.frame(ff = x$ff), digits=digits,...)
275
  }
276

277 1
  printdistimes <- !is.null(x$distimes)
278 1
  if(printdistimes){
279 1
    cat("\nEstimated disappearance times:\n")
280 1
    print(x$distimes, digits=digits,...)
281
  }
282

283 1
  printdata <- !is.null(x$data)
284 1
  if (printdata){
285 1
    cat("\nData:\n")
286 1
    print(format(x$data, digits = digits, ...), row.names = FALSE)
287
  }
288

289 1
  invisible(x)
290
}

Read our documentation on viewing source code .

Loading