1
#' Confidence intervals for parameters of mkinfit objects
2
#'
3
#' The default method 'quadratic' is based on the quadratic approximation of
4
#' the curvature of the likelihood function at the maximum likelihood parameter
5
#' estimates.
6
#' The alternative method 'profile' is based on the profile likelihood for each
7
#' parameter. The 'profile' method uses two nested optimisations and can take a
8
#' very long time, even if parallelized by specifying 'cores' on unixoid
9
#' platforms. The speed of the method could likely be improved by using the
10
#' method of Venzon and Moolgavkar (1988).
11
#'
12
#' @param object An \code{\link{mkinfit}} object
13
#' @param parm A vector of names of the parameters which are to be given
14
#'   confidence intervals. If missing, all parameters are considered.
15
#' @param level The confidence level required
16
#' @param alpha The allowed error probability, overrides 'level' if specified.
17
#' @param cutoff Possibility to specify an alternative cutoff for the difference
18
#'   in the log-likelihoods at the confidence boundary. Specifying an explicit
19
#'   cutoff value overrides arguments 'level' and 'alpha'
20
#' @param method The 'quadratic' method approximates the likelihood function at
21
#'   the optimised parameters using the second term of the Taylor expansion,
22
#'   using a second derivative (hessian) contained in the object.
23
#'   The 'profile' method searches the parameter space for the
24
#'   cutoff of the confidence intervals by means of a likelihood ratio test.
25
#' @param transformed If the quadratic approximation is used, should it be
26
#'   applied to the likelihood based on the transformed parameters?
27
#' @param backtransform If we approximate the likelihood in terms of the
28
#'   transformed parameters, should we backtransform the parameters with
29
#'   their confidence intervals?
30
#' @param rel_tol If the method is 'profile', what should be the accuracy
31
#'   of the lower and upper bounds, relative to the estimate obtained from
32
#'   the quadratic method?
33
#' @param cores The number of cores to be used for multicore processing.
34
#'   On Windows machines, cores > 1 is currently not supported.
35
#' @param quiet Should we suppress the message "Profiling the likelihood"
36
#' @return A matrix with columns giving lower and upper confidence limits for
37
#'   each parameter.
38
#' @param \dots Not used
39
#' @importFrom stats qnorm
40
#' @references
41
#'   Bates DM and Watts GW (1988) Nonlinear regression analysis & its applications
42
#'
43
#'   Pawitan Y (2013) In all likelihood - Statistical modelling and
44
#'   inference using likelihood. Clarendon Press, Oxford.
45
#'
46
#'   Venzon DJ and Moolgavkar SH (1988) A Method for Computing
47
#'   Profile-Likelihood Based Confidence Intervals, Applied Statistics, 37,
48
#'   87–94.
49
#' @examples
50
#' f <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE)
51
#' confint(f, method = "quadratic")
52
#'
53
#' \dontrun{
54
#' confint(f, method = "profile")
55
#'
56
#' # Set the number of cores for the profiling method for further examples
57
#' if (identical(Sys.getenv("NOT_CRAN"), "true")) {
58
#'   n_cores <- parallel::detectCores() - 1
59
#' } else {
60
#'  n_cores <- 1
61
#' }
62
#' if (Sys.getenv("TRAVIS") != "") n_cores = 1
63
#' if (Sys.info()["sysname"] == "Windows") n_cores = 1
64
#'
65
#' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), quiet = TRUE)
66
#' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"),
67
#'   use_of_ff = "max", quiet = TRUE)
68
#' f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE)
69
#' system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
70
#' # Using more cores does not save much time here, as parent_0 takes up most of the time
71
#' # If we additionally exclude parent_0 (the confidence of which is often of
72
#' # minor interest), we get a nice performance improvement from about 50
73
#' # seconds to about 12 seconds if we use at least four cores
74
#' system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile",
75
#'   c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
76
#' ci_profile
77
#' ci_quadratic_transformed <- confint(f_d_1, method = "quadratic")
78
#' ci_quadratic_transformed
79
#' ci_quadratic_untransformed <- confint(f_d_1, method = "quadratic", transformed = FALSE)
80
#' ci_quadratic_untransformed
81
#' # Against the expectation based on Bates and Watts (1988), the confidence
82
#' # intervals based on the internal parameter transformation are less
83
#' # congruent with the likelihood based intervals. Note the superiority of the
84
#' # interval based on the untransformed fit for k_m1_sink
85
#' rel_diffs_transformed <- abs((ci_quadratic_transformed - ci_profile)/ci_profile)
86
#' rel_diffs_untransformed <- abs((ci_quadratic_untransformed - ci_profile)/ci_profile)
87
#' rel_diffs_transformed < rel_diffs_untransformed
88
#' signif(rel_diffs_transformed, 3)
89
#' signif(rel_diffs_untransformed, 3)
90
#'
91
#'
92
#' # Investigate a case with formation fractions
93
#' f_d_2 <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE)
94
#' ci_profile_ff <- confint(f_d_2, method = "profile", cores = n_cores)
95
#' ci_profile_ff
96
#' ci_quadratic_transformed_ff <- confint(f_d_2, method = "quadratic")
97
#' ci_quadratic_transformed_ff
98
#' ci_quadratic_untransformed_ff <- confint(f_d_2, method = "quadratic", transformed = FALSE)
99
#' ci_quadratic_untransformed_ff
100
#' rel_diffs_transformed_ff <- abs((ci_quadratic_transformed_ff - ci_profile_ff)/ci_profile_ff)
101
#' rel_diffs_untransformed_ff <- abs((ci_quadratic_untransformed_ff - ci_profile_ff)/ci_profile_ff)
102
#' # While the confidence interval for the parent rate constant is closer to
103
#' # the profile based interval when using the internal parameter
104
#' # transformation, the interval for the metabolite rate constant is 'better
105
#' # without internal parameter transformation.
106
#' rel_diffs_transformed_ff < rel_diffs_untransformed_ff
107
#' rel_diffs_transformed_ff
108
#' rel_diffs_untransformed_ff
109
#'
110
#' # The profiling for the following fit does not finish in a reasonable time,
111
#' # therefore we use the quadratic approximation
112
#' m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
113
#'   M1 = mkinsub("SFO"),
114
#'   M2 = mkinsub("SFO"),
115
#'   use_of_ff = "max", quiet = TRUE)
116
#' DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
117
#' f_tc_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, error_model = "tc",
118
#'   error_model_algorithm = "direct", quiet = TRUE)
119
#' confint(f_tc_2, method = "quadratic")
120
#' confint(f_tc_2, "parent_0", method = "quadratic")
121
#' }
122
#' @export
123
confint.mkinfit <- function(object, parm,
124
  level = 0.95, alpha = 1 - level, cutoff,
125
  method = c("quadratic", "profile"),
126
  transformed = TRUE, backtransform = TRUE,
127
  cores = parallel::detectCores(), rel_tol = 0.01, quiet = FALSE, ...)
128
{
129 2
  tparms <- parms(object, transformed = TRUE)
130 2
  bparms <- parms(object, transformed = FALSE)
131 2
  tpnames <- names(tparms)
132 2
  bpnames <- names(bparms)
133

134 2
  return_pnames <- if (missing(parm)) {
135 2
    if (backtransform) bpnames else tpnames
136
  } else {
137 0
    parm
138
  }
139

140 2
  p <- length(return_pnames)
141

142 2
  method <- match.arg(method)
143

144 2
  a <- c(alpha / 2, 1 - (alpha / 2))
145

146 2
  quantiles <- qt(a, object$df.residual)
147

148 2
  covar_pnames <- if (missing(parm)) {
149 2
    if (transformed) tpnames else bpnames
150
  } else {
151 0
    parm
152
  }
153

154 2
  return_parms <- if (backtransform) bparms[return_pnames]
155 2
    else tparms[return_pnames]
156

157 2
  covar_parms <- if (transformed) tparms[covar_pnames]
158 2
    else bparms[covar_pnames]
159

160 2
  if (transformed) {
161 2
    covar <- try(solve(object$hessian), silent = TRUE)
162
  } else {
163 2
    covar <- try(solve(object$hessian_notrans), silent = TRUE)
164
  }
165

166
  # If inverting the covariance matrix failed or produced NA values
167 2
  if (!is.numeric(covar) | is.na(covar[1])) {
168 0
    ses <- lci <- uci <- rep(NA, p)
169
  } else {
170 2
    ses     <- sqrt(diag(covar))[covar_pnames]
171 2
    lci    <- covar_parms + quantiles[1] * ses
172 2
    uci    <- covar_parms + quantiles[2] * ses
173 2
    if (transformed & backtransform) {
174 2
      lci_back <- backtransform_odeparms(lci,
175 2
        object$mkinmod, object$transform_rates, object$transform_fractions)
176 2
      uci_back <- backtransform_odeparms(uci,
177 2
        object$mkinmod, object$transform_rates, object$transform_fractions)
178

179 2
      return_errparm_names <- intersect(names(object$errparms), return_pnames)
180 2
      lci <- c(lci_back, lci[return_errparm_names])
181 2
      uci <- c(uci_back, uci[return_errparm_names])
182
    }
183
  }
184 2
  ci <- cbind(lower = lci, upper = uci)
185

186 2
  if (method == "profile") {
187

188 2
    ci_quadratic <- ci
189

190 0
    if (!quiet) message("Profiling the likelihood")
191

192 2
    lci <- uci <- rep(NA, p)
193 2
    names(lci) <- names(uci) <- return_pnames
194

195 2
    profile_pnames <- if(missing(parm)) names(parms(object))
196 2
      else parm
197

198 2
    if (missing(cutoff)) {
199 2
      cutoff <- 0.5 * qchisq(1 - alpha, 1)
200
    }
201

202 2
    all_parms <- parms(object)
203

204 2
    get_ci <- function(pname) {
205 2
      pnames_free <- setdiff(names(all_parms), pname)
206 2
      profile_ll <- function(x)
207
      {
208 2
        pll_cost <- function(P) {
209 2
          parms_cost <- all_parms
210 2
          parms_cost[pnames_free] <- P[pnames_free]
211 2
          parms_cost[pname] <- x
212 2
          - object$ll(parms_cost)
213
        }
214 2
        - nlminb(all_parms[pnames_free], pll_cost)$objective
215
      }
216

217 2
      cost <- function(x) {
218 2
        (cutoff - (object$logLik - profile_ll(x)))^2
219
      }
220

221 2
      lower_quadratic <- ci_quadratic["lower"][pname]
222 2
      upper_quadratic <- ci_quadratic["upper"][pname]
223 2
      ltol <- if (!is.na(lower_quadratic)) rel_tol * lower_quadratic else .Machine$double.eps^0.25
224 2
      utol <- if (!is.na(upper_quadratic)) rel_tol * upper_quadratic else .Machine$double.eps^0.25
225 2
      lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname], tol = ltol)$minimum
226 2
      uci_pname <- optimize(cost, lower = all_parms[pname],
227 2
        upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]),
228 2
        tol = utol)$minimum
229 2
      return(c(lci_pname, uci_pname))
230
    }
231 2
    ci <- t(parallel::mcmapply(get_ci, profile_pnames, mc.cores = cores))
232
  }
233

234 2
  colnames(ci) <- paste0(
235 2
    format(100 * a, trim = TRUE, scientific = FALSE, digits = 3), "%")
236

237 2
  return(ci)
238
}

Read our documentation on viewing source code .

Loading