1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
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 confint.nls", {
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)
# 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:
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)
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 <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE)
# f <- fits[["SFO", "FOCUS_C"]]
# f_notrans <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE, transform_rates = FALSE)
# pf <- parms(f)
# 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_nll(parent_0 = 100, k_parent_sink = 0.3, sigma = 4.7)
# f_nll(pf[1], pf[2], pf[3])
# f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = nrow(FOCUS_2006_C))
# f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = 17L)
# logLik(f_mle)
# stats4::confint(f_mle, nobs = nrow(FOCUS_2006_C))
#
# ci_mkin_1 <- confint(f,
# method = "quadratic", backtransform = FALSE,
# transformed = TRUE)
# ci_maxLik_1 <- maxLik::confint.maxLik(f_maxLik)
#
# f_tc_2_maxLik <- maxLik::maxLik(f_tc_2$ll,
# start = f_tc_2$par)
# ci_tc_2 <- confint(f_tc_2, method = "quadratic",
# backtransform = FALSE, transformed = TRUE,
# distribution = "normal")
# ci_tc_2_maxLik <- confint(f_tc_2_maxLik)
# rel_diff_ci_tc_2 <- (ci_tc_2 - ci_tc_2_maxLik)/ci_tc_2_maxLik
# # The ilr transformed variable appear to misbehave in the numeric differentiation
# expect_equivalent(
# rel_diff_ci_tc_2[c(2, 3, 4, 6, 7, 9, 10), ], rep(0, 14),
# tolerance = 0.05)
#
# ci_tc_2[, 1]
# ci_tc_2_maxLik[, 1]
})
|
