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).

- -

# S3 method for 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.183   0.000   1.182 
-# 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.418   0.116   0.294 
-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.403833578 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.403833583 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.0307283041 0.0290588361
-#> f_parent_to_m1 0.0046881769 0.0027780063
-#> sigma          0.0550252516 0.0327066836
-rel_diffs_untransformed_ff
-#>                        2.5%        97.5%
-#> parent_0       0.0005408689 0.0002217232
-#> k_parent       0.0046102156 0.0023732281
-#> k_m1           0.0146740690 0.0025291820
-#> f_parent_to_m1 0.0046995211 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.596039609 106.19954892
-#> k_M1            0.037605368   0.04490762
-#> k_M2            0.008568731   0.01087676
-#> f_parent_to_M1  0.021462489   0.62023882
-#> f_parent_to_M2  0.015165617   0.37975348
-#> k1              0.273897348   0.33388101
-#> k2              0.018614554   0.02250378
-#> g               0.671943411   0.73583305
-#> sigma_low       0.251283495   0.83992077
-#> rsd_high        0.040411024   0.07662008
-confint(f_tc_2, "parent_0", method = "quadratic")
-#>              2.5%    97.5%
-#> parent_0 94.59604 106.1995
-# }
-

- -