aboutsummaryrefslogblamecommitdiff
path: root/tests/testthat/test_confidence.R
blob: 5f76c34460d8004a7f8b4c30e7f080fa8ae74b35 (plain) (tree)





















                                                                              
                                            
 



                                                                            
 


                                                                            
 



                                                                                  
 


                                                                  
 

  
                                                                                                        































                                                                                        







                                                                     
  
 




















                                                                                    
# We set up some models and fits with nls for comparisons
SFO_trans <- function(t, parent_0, log_k_parent_sink) {
  parent_0 * exp(- exp(log_k_parent_sink) * t)
}
SFO_notrans <- function(t, parent_0, k_parent_sink) {
  parent_0 * exp(- k_parent_sink * t)
}
f_1_nls_trans <- nls(value ~ SFO_trans(time, parent_0, log_k_parent_sink),
  data = FOCUS_2006_A,
  start = list(parent_0 = 100, log_k_parent_sink = log(0.1)))
f_1_nls_notrans <- nls(value ~ SFO_notrans(time, parent_0, k_parent_sink),
  data = FOCUS_2006_A,
  start = list(parent_0 = 100, k_parent_sink = 0.1))

f_1_mkin_OLS <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE)
f_1_mkin_OLS_notrans <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE,
  transform_rates = FALSE)


f_2_mkin <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = FOCUS_2006_C)

context("Confidence intervals and p-values")

test_that("The confint method 'quadratic' is consistent with the summary", {
  expect_equivalent(
    confint(f_SFO_lin_mkin_ML, method = "quadratic"),
    summary(f_SFO_lin_mkin_ML)$bpar[, c("Lower", "Upper")])

  expect_equivalent(
    confint(f_SFO_lin_mkin_ML, method = "quadratic", backtransform = FALSE),
    summary(f_SFO_lin_mkin_ML)$par[, c("Lower", "Upper")])

  f_notrans <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE, transform_rates = FALSE)
  expect_equivalent(
    confint(f_notrans, method = "quadratic", transformed = FALSE),
    summary(f_notrans)$par[, c("Lower", "Upper")])

  expect_equivalent(
    confint(f_notrans, method = "quadratic", transformed = FALSE),
    summary(f_notrans)$bpar[, c("Lower", "Upper")])

})

test_that("Quadratic confidence intervals for rate constants are comparable to values in summary.nls", {

  # Check fitted parameter values
  expect_equivalent(
    (f_1_mkin_OLS$bparms.optim -coef(f_1_nls_notrans))/f_1_mkin_OLS$bparms.optim,
    rep(0, 2), tolerance = 1e-6)
  expect_equivalent(
    (f_1_mkin_OLS$par[1:2] - coef(f_1_nls_trans))/f_1_mkin_OLS$par[1:2],
    rep(0, 2), tolerance = 1e-6)

  # Check the standard error for the transformed parameters
  se_nls <- summary(f_1_nls_trans)$coefficients[, "Std. Error"]
  # This is of similar magnitude as the standard error obtained with the mkin
  se_mkin <- summary(f_1_mkin_OLS)$par[1:2, "Std. Error"]

  se_nls_notrans <- summary(f_1_nls_notrans)$coefficients[, "Std. Error"]
  # This is also of similar magnitude as the standard error obtained with the mkin
  se_mkin_notrans <- summary(f_1_mkin_OLS_notrans)$par[1:2, "Std. Error"]

  # The difference can partly be explained by the ratio between
  # the maximum likelihood estimate of the standard error sqrt(rss/n)
  # and the estimate used in nls sqrt(rss/rdf), i.e. by the factor
  # sqrt(n/rdf).
  # In the case of mkin, the residual degrees of freedom are only taken into
  # account in the subsequent step of generating the confidence intervals for
  # the parameters (including sigma)

  # Strangely, this only works for the rate constant to less than 1%, but
  # not for the initial estimate
  expect_equivalent(se_nls[2] / se_mkin[2], sqrt(8/5), tolerance = 0.01)
  expect_equivalent(se_nls_notrans[2] / se_mkin_notrans[2], sqrt(8/5), tolerance = 0.01)

  # Another case:
  se_mkin_2 <- summary(f_2_mkin)$par[1:4, "Std. Error"]
  se_nls_2 <- summary(f_2_nls)$coefficients[, "Std. Error"]
  # Here we the ratio of standard errors can be explained by the same
  # principle up to about 3%
  expect_equivalent(
    se_nls_2[c("lrc1", "lrc2")] / se_mkin_2[c("log_k1", "log_k2")],
    rep(sqrt(nrow(FOCUS_2006_C) / (nrow(FOCUS_2006_C) - 4)), 2),
    tolerance = 0.03)
})

test_that("Likelihood profile based confidence intervals work", {
   f <- fits[["SFO", "FOCUS_C"]]

   # negative log-likelihood for use with mle
   f_nll <- function(parent_0, k_parent_sink, sigma) {
     - f$ll(c(parent_0 = as.numeric(parent_0),
         k_parent_sink = as.numeric(k_parent_sink),
         sigma = as.numeric(sigma)))
   }
   f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = nrow(FOCUS_2006_C))

   ci_mkin_1_p_0.95 <- confint(f, method = "profile", level = 0.95, quiet = TRUE)

   # Magically, we get very similar boundaries as stats4::mle
   # (we need to capture the output to avoid printing this while testing as
   # stats4::confint uses cat() for its message, instead of message(), so
   # suppressMessage() has no effect)
   msg <- capture.output(ci_mle_1_0.95 <- stats4::confint(f_mle, level = 0.95))
   rel_diff_ci <- (ci_mle_1_0.95 - ci_mkin_1_p_0.95)/ci_mle_1_0.95
   expect_equivalent(as.numeric(rel_diff_ci), rep(0, 6), tolerance = 1e-4)
})

Contact - Imprint