From bd761d879a95872f82e6d8f893634a61e122a938 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 28 Oct 2019 13:23:40 +0100 Subject: Fix the cutoff for likelihood based intervals The cutoff now matches what is given by Venzon and Moolgavkar (1988). Also, confidence intervals closely match intervals obtained with stats4::confint in the test case where an stats4::mle object is created from the likelihood function in one test case. Static documentation rebuilt by pkgdown --- DESCRIPTION | 2 +- R/confint.mkinfit.R | 20 +++++-- docs/reference/confint.mkinfit.html | 31 +++++++---- man/confint.mkinfit.Rd | 17 ++++-- test.log | 16 +++--- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/test_confidence.R | 103 +++++++++++++++--------------------- 7 files changed, 105 insertions(+), 86 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 88412fa8..544af139 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,7 +19,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, purpose. Imports: stats, graphics, methods, deSolve, R6, inline, parallel, numDeriv Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, - benchmarkme, tibble + benchmarkme, tibble, stats4 License: GPL LazyLoad: yes LazyData: yes diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R index 58a3c8db..1ed689a9 100644 --- a/R/confint.mkinfit.R +++ b/R/confint.mkinfit.R @@ -1,10 +1,17 @@ #' Confidence intervals for parameters of mkinfit objects #' +#' The default method 'profile' is based on the profile likelihood for each +#' parameter. The method uses two nested optimisations. The speed of the method +#' could likely be improved by using the method of Venzon and Moolgavkar (1988). +#' #' @param object An \code{\link{mkinfit}} object #' @param parm A vector of names of the parameters which are to be given #' confidence intervals. If missing, all parameters are considered. #' @param level The confidence level required #' @param alpha The allowed error probability, overrides 'level' if specified. +#' @param 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' #' @param method The 'profile' method searches the parameter space for the #' cutoff of the confidence intervals by means of a likelihood ratio test. #' The 'quadratic' method approximates the likelihood function at the @@ -25,6 +32,9 @@ #' @importFrom stats qnorm #' @references 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") @@ -33,7 +43,7 @@ #' } #' @export confint.mkinfit <- function(object, parm, - level = 0.95, alpha = 1 - level, + level = 0.95, alpha = 1 - level, cutoff, method = c("profile", "quadratic"), transformed = TRUE, backtransform = TRUE, distribution = c("student_t", "normal"), quiet = FALSE, ...) @@ -100,15 +110,17 @@ confint.mkinfit <- function(object, parm, } if (method == "profile") { - message("Profiling the likelihood") + if (!quiet) message("Profiling the likelihood") + lci <- uci <- rep(NA, p) names(lci) <- names(uci) <- return_pnames profile_pnames <- if(missing(parm)) names(parms(object)) else parm - # We do two-sided intervals based on the likelihood ratio - cutoff <- 0.5 * qchisq(1 - (alpha / 2), 1) + if (missing(cutoff)) { + cutoff <- 0.5 * qchisq(1 - alpha, 1) + } all_parms <- parms(object) diff --git a/docs/reference/confint.mkinfit.html b/docs/reference/confint.mkinfit.html index 64ea0d4c..ec8938c5 100644 --- a/docs/reference/confint.mkinfit.html +++ b/docs/reference/confint.mkinfit.html @@ -36,7 +36,9 @@ - + @@ -133,14 +135,16 @@
-

Confidence intervals for parameters of mkinfit objects

+

The default method 'profile' is based on the profile likelihood for each +parameter. The method uses two nested optimisations. 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, method = c("profile", "quadratic"), transformed = TRUE,
-  backtransform = TRUE, distribution = c("student_t", "normal"),
-  quiet = FALSE, ...)
+ level, cutoff, method = c("profile", "quadratic"), + transformed = TRUE, backtransform = TRUE, + distribution = c("student_t", "normal"), quiet = FALSE, ...)

Arguments

@@ -162,6 +166,12 @@ confidence intervals. If missing, all parameters are considered.

+ + + +

References

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

+ 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)
@@ -213,9 +226,9 @@ the parameter estimate

 #> k_parent_sink  0.2109541  0.4440528
 #> sigma          1.9778868  7.3681380
# \dontrun{
   confint(f, method = "profile")
#> Profiling the likelihood
#>                     2.5%      97.5%
-#> parent_0      71.3471007 93.9447024
-#> k_parent_sink  0.2030765  0.4491067
-#> sigma          2.9810656  8.8633278
# }
+#> parent_0      73.0641834 92.1392181
+#> k_parent_sink  0.2170293  0.4235348
+#> sigma          3.1307772  8.0628314
# }
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 'profile' method searches the parameter space for the @@ -204,7 +214,10 @@ the parameter estimate