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 .