aboutsummaryrefslogtreecommitdiff
path: root/tests/testthat/test_mixed.R
blob: 6ac93295fb84ffbf4db646a44c0b39bb0d69659a (plain) (blame)
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
context("Nonlinear mixed-effects models")

test_that("Parent fits using saemix are correctly implemented", {
  skip_if(!saemix_available)

  expect_error(saem(fits), "Only row objects")
  # Some fits were done in the setup script
  mmkin_sfo_2 <- update(mmkin_sfo_1, fixed_initials = c(parent = 100))
  expect_error(update(mmkin_sfo_1, models = c("SFOOO")), "Please supply models.*")

  sfo_saem_2 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "mkin")
  sfo_saem_3 <- expect_error(saem(mmkin_sfo_2, quiet = TRUE), "at least two parameters")
  s_sfo_s1 <- summary(sfo_saem_1)
  s_sfo_s2 <- summary(sfo_saem_2)

  sfo_nlme_1 <- expect_warning(nlme(mmkin_sfo_1), "not converge")
  s_sfo_n <- summary(sfo_nlme_1)

  # Compare with input
  expect_equal(round(s_sfo_s2$confint_ranef["SD.log_k_parent", "est."], 1), 0.3)
  # k_parent is a bit different from input 0.03 here
  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3), 0.035)
  expect_equal(round(s_sfo_s2$confint_back["k_parent", "est."], 3), 0.035)

  # But the result is pretty unanimous between methods
  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
    round(s_sfo_s2$confint_back["k_parent", "est."], 3))
  expect_equal(round(s_sfo_s1$confint_back["k_parent", "est."], 3),
    round(s_sfo_n$confint_back["k_parent", "est."], 3))

  mmkin_fomc_1 <- mmkin("FOMC", ds_fomc, quiet = TRUE, error_model = "tc", cores = n_cores)
  fomc_saem_1 <- saem(mmkin_fomc_1, quiet = TRUE)
  ci_fomc_s1 <- summary(fomc_saem_1)$confint_back

  fomc_pop <- as.numeric(fomc_pop)
  expect_true(all(ci_fomc_s1[, "lower"] < fomc_pop))
  expect_true(all(ci_fomc_s1[, "upper"] > fomc_pop))

  mmkin_fomc_2 <- update(mmkin_fomc_1, state.ini = 100, fixed_initials = "parent")
  fomc_saem_2 <- saem(mmkin_fomc_2, quiet = TRUE, transformations = "mkin")
  ci_fomc_s2 <- summary(fomc_saem_2)$confint_back

  expect_true(all(ci_fomc_s2[, "lower"] < fomc_pop[2:3]))
  expect_true(all(ci_fomc_s2[, "upper"] > fomc_pop[2:3]))

  s_dfop_s1 <- summary(dfop_saemix_1)
  s_dfop_s2 <- summary(dfop_saemix_2)
  s_dfop_n <- summary(dfop_nlme_1)

  dfop_pop <- as.numeric(dfop_pop)
  expect_true(all(s_dfop_s1$confint_back[, "lower"] < dfop_pop))
  expect_true(all(s_dfop_s1$confint_back[, "upper"] > dfop_pop))
  expect_true(all(s_dfop_s2$confint_back[, "lower"] < dfop_pop))
  expect_true(all(s_dfop_s2$confint_back[, "upper"] > dfop_pop))

  dfop_mmkin_means_trans_tested <- mean_degparms(mmkin_dfop_1, test_log_parms = TRUE)
  dfop_mmkin_means_trans <- apply(parms(mmkin_dfop_1, transformed = TRUE), 1, mean)

  dfop_mmkin_means_tested <- backtransform_odeparms(dfop_mmkin_means_trans_tested, mmkin_dfop_1$mkinmod)
  dfop_mmkin_means <- backtransform_odeparms(dfop_mmkin_means_trans, mmkin_dfop_1$mkinmod)

  # We get < 20% deviations for parent_0 and k1 by averaging the transformed parameters
  # If we average only parameters passing the t-test, the deviation for k2 is also < 20%
  rel_diff_mmkin <- (dfop_mmkin_means - dfop_pop) / dfop_pop
  rel_diff_mmkin_tested <- (dfop_mmkin_means_tested - dfop_pop) / dfop_pop
  expect_true(all(rel_diff_mmkin[c("parent_0", "k1")] < 0.20))
  expect_true(all(rel_diff_mmkin_tested[c("parent_0", "k1", "k2")] < 0.20))

  # We get < 15% deviations with transformations made in mkin
  rel_diff_1 <- (s_dfop_s1$confint_back[, "est."] - dfop_pop) / dfop_pop
  expect_true(all(rel_diff_1 < 0.15))

  # We get < 20% deviations with transformations made in saemix
  rel_diff_2 <- (s_dfop_s2$confint_back[, "est."] - dfop_pop) / dfop_pop
  expect_true(all(rel_diff_2 < 0.2))

  mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores)
  hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE)
  ci_hs_s1 <- summary(hs_saem_1)$confint_back

  hs_pop <- as.numeric(hs_pop)
  # expect_true(all(ci_hs_s1[, "lower"] < hs_pop)) # k1 is overestimated
  expect_true(all(ci_hs_s1[, "upper"] > hs_pop))

  mmkin_hs_2 <- update(mmkin_hs_1, state.ini = 100, fixed_initials = "parent")
  hs_saem_2 <- saem(mmkin_hs_2, quiet = TRUE)
  ci_hs_s2 <- summary(hs_saem_2)$confint_back

  #expect_true(all(ci_hs_s2[, "lower"] < hs_pop[2:4])) # k1 again overestimated
  expect_true(all(ci_hs_s2[, "upper"] > hs_pop[2:4]))

  # HS would likely benefit from implemenation of transformations = "saemix"
})

test_that("Print methods work", {
  expect_known_output(print(fits[, 2:3], digits = 2), "print_mmkin_parent.txt")
  expect_known_output(print(mmkin_biphasic_mixed, digits = 2), "print_mmkin_biphasic_mixed.txt")
  expect_known_output(print(nlme_biphasic, digits = 1), "print_nlme_biphasic.txt")

  skip_if(!saemix_available)
  expect_known_output(print(sfo_saem_1, digits = 1), "print_sfo_saem_1.txt")
})

test_that("nlme results are reproducible to some degree", {

  test_summary <- summary(nlme_biphasic)
  test_summary$nlmeversion <- "Dummy 0.0 for testing"
  test_summary$mkinversion <- "Dummy 0.0 for testing"
  test_summary$Rversion <- "Dummy R version for testing"
  test_summary$date.fit <- "Dummy date for testing"
  test_summary$date.summary <- "Dummy date for testing"
  test_summary$time <- c(elapsed = "test time 0")

  expect_known_output(print(test_summary, digits = 1), "summary_nlme_biphasic_s.txt")

  # k1 just fails the first test (lower bound of the ci), so we need to exclude it
  dfop_no_k1 <- c("parent_0", "k_m1", "f_parent_to_m1", "k2", "g")
  dfop_sfo_pop_no_k1 <- as.numeric(dfop_sfo_pop[dfop_no_k1])
  dfop_sfo_pop <- as.numeric(dfop_sfo_pop)

  ci_dfop_sfo_n <- summary(nlme_biphasic)$confint_back

  expect_true(all(ci_dfop_sfo_n[dfop_no_k1, "lower"] < dfop_sfo_pop_no_k1))
  expect_true(all(ci_dfop_sfo_n[, "upper"] > dfop_sfo_pop))
})

test_that("saem results are reproducible for biphasic fits", {

  skip_if(!saemix_available)
  test_summary <- summary(saem_biphasic_s)
  test_summary$saemixversion <- "Dummy 0.0 for testing"
  test_summary$mkinversion <- "Dummy 0.0 for testing"
  test_summary$Rversion <- "Dummy R version for testing"
  test_summary$date.fit <- "Dummy date for testing"
  test_summary$date.summary <- "Dummy date for testing"
  test_summary$time <- c(elapsed = "test time 0")

  expect_known_output(print(test_summary, digits = 2), "summary_saem_biphasic_s.txt")

  dfop_sfo_pop <- as.numeric(dfop_sfo_pop)
  no_k1 <- c(1, 2, 3, 5, 6)
  no_k2 <- c(1, 2, 3, 4, 6)
  no_k1_k2 <- c(1, 2, 3, 6)

  ci_dfop_sfo_s_s <- summary(saem_biphasic_s)$confint_back
  # k1 and k2 are overestimated
  expect_true(all(ci_dfop_sfo_s_s[no_k1_k2, "lower"] < dfop_sfo_pop[no_k1_k2]))
  expect_true(all(ci_dfop_sfo_s_s[, "upper"] > dfop_sfo_pop))

  # k1 and k2 are not fitted well
  ci_dfop_sfo_s_m <- summary(saem_biphasic_m)$confint_back
  expect_true(all(ci_dfop_sfo_s_m[no_k2, "lower"] < dfop_sfo_pop[no_k2]))
  expect_true(all(ci_dfop_sfo_s_m[no_k1, "upper"] > dfop_sfo_pop[no_k1]))

  # I tried to only do few iterations in routine tests as this is so slow
  # but then deSolve fails at some point (presumably at the switch between
  # the two types of iterations)
  #saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve",
  # control = list(nbiter.saemix = c(10, 5), nbiter.burn = 5), quiet = TRUE)

  skip("Fitting with saemix takes around 10 minutes when using deSolve")
  saem_biphasic_2 <- saem(mmkin_biphasic, solution_type = "deSolve", quiet = TRUE)

  # As with the analytical solution, k1 and k2 are not fitted well
  ci_dfop_sfo_s_d <- summary(saem_biphasic_2)$confint_back
  expect_true(all(ci_dfop_sfo_s_d[no_k2, "lower"] < dfop_sfo_pop[no_k2]))
  expect_true(all(ci_dfop_sfo_s_d[no_k1, "upper"] > dfop_sfo_pop[no_k1]))
})

Contact - Imprint