Confidence intervals for parameters of mkinfit objects

+ +

The default method 'quadratic' is based on the quadratic approximation of +the curvature of the likelihood function at the maximum likelihood parameter +estimates. +The alternative method 'profile' is based on the profile likelihood for each +parameter. The 'profile' method uses two nested optimisations and can take a +very long time, even if parallelized by specifying 'cores' on unixoid +platforms. The speed of the method could likely be improved by using the +method of Venzon and Moolgavkar (1988).

+ +

Usage

# S3 method for class 'mkinfit'
+confint(
+  object,
+  parm,
+  level = 0.95,
+  alpha = 1 - level,
+  cutoff,
+  method = c("quadratic", "profile"),
+  transformed = TRUE,
+  backtransform = TRUE,
+  cores = parallel::detectCores(),
+  rel_tol = 0.01,
+  quiet = FALSE,
+  ...
+)

+ +

Arguments

+ + +

object: An mkinfit object
parm: A vector of names of the parameters which are to be given +confidence intervals. If missing, all parameters are considered.
level: The confidence level required
alpha: The allowed error probability, overrides 'level' if specified.
cutoff: Possibility to specify an alternative cutoff for the difference +in the log-likelihoods at the confidence boundary. Specifying an explicit +cutoff value overrides arguments 'level' and 'alpha'
method: The 'quadratic' method approximates the likelihood function at +the optimised parameters using the second term of the Taylor expansion, +using a second derivative (hessian) contained in the object. +The 'profile' method searches the parameter space for the +cutoff of the confidence intervals by means of a likelihood ratio test.
transformed: If the quadratic approximation is used, should it be +applied to the likelihood based on the transformed parameters?
backtransform: If we approximate the likelihood in terms of the +transformed parameters, should we backtransform the parameters with +their confidence intervals?
cores: The number of cores to be used for multicore processing. +On Windows machines, cores > 1 is currently not supported.
rel_tol: If the method is 'profile', what should be the accuracy +of the lower and upper bounds, relative to the estimate obtained from +the quadratic method?
quiet: Should we suppress the message "Profiling the likelihood"
...: Not used

Value

A matrix with columns giving lower and upper confidence limits for +each parameter.

References

Bates DM and Watts GW (1988) Nonlinear regression analysis & its applications

Pawitan Y (2013) In all likelihood - Statistical modelling and +inference using likelihood. Clarendon Press, Oxford.

Venzon DJ and Moolgavkar SH (1988) A Method for Computing +Profile-Likelihood Based Confidence Intervals, Applied Statistics, 37, +87–94.

+ +

Examples

f <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE)
+confint(f, method = "quadratic")
+#>                2.5%      97.5%
+#> parent_0 71.8242430 93.1600766
+#> k_parent  0.2109541  0.4440528
+#> sigma     1.9778868  7.3681380
+
+# \dontrun{
+confint(f, method = "profile")
+#> Profiling the likelihood
+#>                2.5%      97.5%
+#> parent_0 73.0641834 92.1392181
+#> k_parent  0.2170293  0.4235348
+#> sigma     3.1307772  8.0628314
+
+# Set the number of cores for the profiling method for further examples
+if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+  n_cores <- parallel::detectCores() - 1
+} else {
+  n_cores <- 1
+}
+if (Sys.getenv("TRAVIS") != "") n_cores = 1
+if (Sys.info()["sysname"] == "Windows") n_cores = 1
+
+SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"),
+  use_of_ff = "min", quiet = TRUE)
+SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"),
+  use_of_ff = "max", quiet = TRUE)
+f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE)
+system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
+#>    user  system elapsed 
+#>   1.182   0.004   1.186 
+# Using more cores does not save much time here, as parent_0 takes up most of the time
+# If we additionally exclude parent_0 (the confidence of which is often of
+# minor interest), we get a nice performance improvement if we use at least 4 cores
+system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile",
+  c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
+#> Profiling the likelihood
+#>    user  system elapsed 
+#>   0.429   0.171   0.324 
+ci_profile
+#>                       2.5%        97.5%
+#> parent_0      96.456003640 1.027703e+02
+#> k_parent_sink  0.040762501 5.549764e-02
+#> k_parent_m1    0.046786482 5.500879e-02
+#> k_m1_sink      0.003892605 6.702778e-03
+#> sigma          2.535612399 3.985263e+00
+ci_quadratic_transformed <- confint(f_d_1, method = "quadratic")
+ci_quadratic_transformed
+#>                       2.5%        97.5%
+#> parent_0      96.403841640 1.027931e+02
+#> k_parent_sink  0.041033378 5.596269e-02
+#> k_parent_m1    0.046777902 5.511931e-02
+#> k_m1_sink      0.004012217 6.897547e-03
+#> sigma          2.396089689 3.854918e+00
+ci_quadratic_untransformed <- confint(f_d_1, method = "quadratic", transformed = FALSE)
+ci_quadratic_untransformed
+#>                       2.5%        97.5%
+#> parent_0      96.403841645 102.79312449
+#> k_parent_sink  0.040485331   0.05535491
+#> k_parent_m1    0.046611582   0.05494364
+#> k_m1_sink      0.003835483   0.00668582
+#> sigma          2.396089689   3.85491806
+# Against the expectation based on Bates and Watts (1988), the confidence
+# intervals based on the internal parameter transformation are less
+# congruent with the likelihood based intervals. Note the superiority of the
+# interval based on the untransformed fit for k_m1_sink
+rel_diffs_transformed <- abs((ci_quadratic_transformed - ci_profile)/ci_profile)
+rel_diffs_untransformed <- abs((ci_quadratic_untransformed - ci_profile)/ci_profile)
+rel_diffs_transformed < rel_diffs_untransformed
+#>                2.5% 97.5%
+#> parent_0      FALSE FALSE
+#> k_parent_sink  TRUE FALSE
+#> k_parent_m1    TRUE FALSE
+#> k_m1_sink     FALSE FALSE
+#> sigma         FALSE FALSE
+signif(rel_diffs_transformed, 3)
+#>                   2.5%    97.5%
+#> parent_0      0.000541 0.000222
+#> k_parent_sink 0.006650 0.008380
+#> k_parent_m1   0.000183 0.002010
+#> k_m1_sink     0.030700 0.029100
+#> sigma         0.055000 0.032700
+signif(rel_diffs_untransformed, 3)
+#>                   2.5%    97.5%
+#> parent_0      0.000541 0.000222
+#> k_parent_sink 0.006800 0.002570
+#> k_parent_m1   0.003740 0.001180
+#> k_m1_sink     0.014700 0.002530
+#> sigma         0.055000 0.032700
+
+
+# Investigate a case with formation fractions
+f_d_2 <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE)
+ci_profile_ff <- confint(f_d_2, method = "profile", cores = n_cores)
+#> Profiling the likelihood
+ci_profile_ff
+#>                        2.5%        97.5%
+#> parent_0       96.456003640 1.027703e+02
+#> k_parent        0.090911032 1.071578e-01
+#> k_m1            0.003892606 6.702775e-03
+#> f_parent_to_m1  0.471328495 5.611550e-01
+#> sigma           2.535612399 3.985263e+00
+ci_quadratic_transformed_ff <- confint(f_d_2, method = "quadratic")
+ci_quadratic_transformed_ff
+#>                        2.5%        97.5%
+#> parent_0       96.403833581 102.79311649
+#> k_parent        0.090823771   0.10725430
+#> k_m1            0.004012219   0.00689755
+#> f_parent_to_m1  0.469118824   0.55959615
+#> sigma           2.396089689   3.85491806
+ci_quadratic_untransformed_ff <- confint(f_d_2, method = "quadratic", transformed = FALSE)
+ci_quadratic_untransformed_ff
+#>                        2.5%        97.5%
+#> parent_0       96.403833586 1.027931e+02
+#> k_parent        0.090491913 1.069035e-01
+#> k_m1            0.003835485 6.685823e-03
+#> f_parent_to_m1  0.469113477 5.598387e-01
+#> sigma           2.396089689 3.854918e+00
+rel_diffs_transformed_ff <- abs((ci_quadratic_transformed_ff - ci_profile_ff)/ci_profile_ff)
+rel_diffs_untransformed_ff <- abs((ci_quadratic_untransformed_ff - ci_profile_ff)/ci_profile_ff)
+# While the confidence interval for the parent rate constant is closer to
+# the profile based interval when using the internal parameter
+# transformation, the interval for the metabolite rate constant is 'better
+# without internal parameter transformation.
+rel_diffs_transformed_ff < rel_diffs_untransformed_ff
+#>                 2.5% 97.5%
+#> parent_0       FALSE FALSE
+#> k_parent        TRUE  TRUE
+#> k_m1           FALSE FALSE
+#> f_parent_to_m1  TRUE FALSE
+#> sigma           TRUE FALSE
+rel_diffs_transformed_ff
+#>                        2.5%        97.5%
+#> parent_0       0.0005408690 0.0002217233
+#> k_parent       0.0009598532 0.0009001864
+#> k_m1           0.0307283045 0.0290588367
+#> f_parent_to_m1 0.0046881768 0.0027780062
+#> sigma          0.0550252516 0.0327066836
+rel_diffs_untransformed_ff
+#>                        2.5%        97.5%
+#> parent_0       0.0005408689 0.0002217233
+#> k_parent       0.0046102155 0.0023732280
+#> k_m1           0.0146740687 0.0025291815
+#> f_parent_to_m1 0.0046995210 0.0023457712
+#> sigma          0.0550252516 0.0327066836
+
+# The profiling for the following fit does not finish in a reasonable time,
+# therefore we use the quadratic approximation
+m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
+  M1 = mkinsub("SFO"),
+  M2 = mkinsub("SFO"),
+  use_of_ff = "max", quiet = TRUE)
+DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
+f_tc_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, error_model = "tc",
+  error_model_algorithm = "direct", quiet = TRUE)
+confint(f_tc_2, method = "quadratic")
+#>                        2.5%        97.5%
+#> parent_0       94.596181875 106.19936592
+#> k_M1            0.037605432   0.04490757
+#> k_M2            0.008568745   0.01087675
+#> f_parent_to_M1  0.021464676   0.62023880
+#> f_parent_to_M2  0.015167158   0.37975350
+#> k1              0.273897535   0.33388072
+#> k2              0.018614555   0.02250379
+#> g               0.671943738   0.73583261
+#> sigma_low       0.251283679   0.83992102
+#> rsd_high        0.040411022   0.07662008
+confint(f_tc_2, "parent_0", method = "quadratic")
+#>              2.5%    97.5%
+#> parent_0 94.59618 106.1994
+# }
+