From c58ccd73951b2000a7a254fb36bbd9f0733db6cd Mon Sep 17 00:00:00 2001
From: Johannes Ranke 
plot(mixed(f_parent_mkin_const["SFO", ]))
 
Using biexponential decline (DFOP) results in a slightly more random scatter of the residuals:
+plot(mixed(f_parent_mkin_const["DFOP", ]))
 
The population curve (bold line) in the above plot results from taking the mean of the individual transformed parameters, i.e. of log k1 @@ -187,7 +187,7 @@ dominates the average. This is alleviated if only rate constants that pass the t-test for significant difference from zero (on the untransformed scale) are considered in the averaging:
+plot(mixed(f_parent_mkin_const["DFOP", ]), test_log_parms = TRUE)
 
While this is visually much more satisfactory, such an average procedure could introduce a bias, as not all results from the individual @@ -199,7 +199,7 @@ degradation model and the error model (see below).
predicted residues is reduced by using the two-component error model: +plot(mixed(f_parent_mkin_tc["DFOP", ]), test_log_parms = TRUE)
 
However, note that in the case of using this error model, the fits to the Flaach and BBA 2.3 datasets appear to be ill-defined, indicated by @@ -289,7 +289,7 @@ effects does not improve the fits.
The selected model (DFOP with two-component error) fitted to the data assuming no correlations between random effects is shown below.
-plot(f_parent_nlme_dfop_tc)plot(f_parent_nlme_dfop_tc)
 
-library(saemix)Loading required package: npde
-Package saemix, version 3.3, March 2024
-  please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
-
-Attaching package: 'saemix'
-The following objects are masked from 'package:npde':
-
-    kurtosis, skewness
-
-saemix_control <- saemixControl(nbiter.saemix = c(800, 300), nb.chains = 15,
+library(saemix)
+saemix_control <- saemixControl(nbiter.saemix = c(800, 300), nb.chains = 15,
     print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)
 saemix_control_moreiter <- saemixControl(nbiter.saemix = c(1600, 300), nb.chains = 15,
     print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)
@@ -324,7 +315,7 @@ Attaching package: 'saemix'
     print = FALSE, save = FALSE, save.graphs = FALSE, displayProgress = FALSE)The convergence plot for the SFO model using constant variance is shown below.
-+@@ -332,19 +323,19 @@ shown below.f_parent_saemix_sfo_const <- mkin::saem(f_parent_mkin_const["SFO", ], quiet = TRUE, control = saemix_control, transformations = "saemix") plot(f_parent_saemix_sfo_const$so, plot.type = "convergence")Obviously the selected number of iterations is sufficient to reach convergence. This can also be said for the SFO fit using the two-component error model.
-+f_parent_saemix_sfo_tc <- mkin::saem(f_parent_mkin_tc["SFO", ], quiet = TRUE, control = saemix_control, transformations = "saemix") plot(f_parent_saemix_sfo_tc$so, plot.type = "convergence")
When fitting the DFOP model with constant variance (see below), parameter convergence is not as unambiguous.
-+f_parent_saemix_dfop_const <- mkin::saem(f_parent_mkin_const["DFOP", ], quiet = TRUE, control = saemix_control, transformations = "saemix") plot(f_parent_saemix_dfop_const$so, plot.type = "convergence")-
+print(f_parent_saemix_dfop_const)Kinetic nonlinear mixed-effects model fit by SAEM Structural model: @@ -380,14 +371,14 @@ this model. also observe that the estimated variance of k2 becomes very small, while being ill-defined, as illustrated by the excessive confidence interval ofSD.k2. -+f_parent_saemix_dfop_tc <- mkin::saem(f_parent_mkin_tc["DFOP", ], quiet = TRUE, control = saemix_control, transformations = "saemix") f_parent_saemix_dfop_tc_moreiter <- mkin::saem(f_parent_mkin_tc["DFOP", ], quiet = TRUE, control = saemix_control_moreiter, transformations = "saemix") plot(f_parent_saemix_dfop_tc$so, plot.type = "convergence")-
+print(f_parent_saemix_dfop_tc)Kinetic nonlinear mixed-effects model fit by SAEM Structural model: @@ -429,7 +420,7 @@ message.The four combinations (SFO/const, SFO/tc, DFOP/const and DFOP/tc) and the version with increased iterations can be compared using the model comparison function of the saemix package:
-+AIC_parent_saemix <- saemix::compare.saemix( f_parent_saemix_sfo_const$so, f_parent_saemix_sfo_tc$so, @@ -437,7 +428,7 @@ comparison function of the saemix package: f_parent_saemix_dfop_tc$so, f_parent_saemix_dfop_tc_moreiter$so)-Likelihoods calculated by importance sampling+ @@ -451,7 +442,7 @@ DFOP tc more iterations 665.85 663.76algorithms implemented in saemix, the likelihood from Gaussian quadrature is added to the best fit, and the AIC values obtained from the three methods are compared. -+f_parent_saemix_dfop_tc$so <- saemix::llgq.saemix(f_parent_saemix_dfop_tc$so) AIC_parent_saemix_methods <- c( @@ -475,7 +466,7 @@ iterations makes a lot of difference. When using the LAPACK version coming with Debian Bullseye, the AIC based on Gaussian quadrature is almost the same as the one obtained with the other methods, also when using defaults for the fit. -+f_parent_saemix_dfop_tc_defaults <- mkin::saem(f_parent_mkin_tc["DFOP", ]) f_parent_saemix_dfop_tc_defaults$so <- saemix::llgq.saemix(f_parent_saemix_dfop_tc_defaults$so) @@ -495,7 +486,7 @@ using defaults for the fit.The following table gives the AIC values obtained with both backend packages using the same control parameters (800 iterations burn-in, 300 iterations second phase, 15 chains).
-++kable(AIC_all)AIC_all <- data.frame( check.names = FALSE, "Degradation model" = c("SFO", "SFO", "DFOP", "DFOP"), @@ -506,7 +497,7 @@ iterations second phase, 15 chains). saemix_is = sapply(list(f_parent_saemix_sfo_const$so, f_parent_saemix_sfo_tc$so, f_parent_saemix_dfop_const$so, f_parent_saemix_dfop_tc$so), AIC, method = "is") ) -kable(AIC_all)