From 04d3370364a0a7881c4ae815c20fd37b74f5639a Mon Sep 17 00:00:00 2001 From: jranke Date: Tue, 5 Mar 2013 01:00:20 +0000 Subject: - Calculate confidence intervals for parameters based on the t distribution Thanks to the CAKE developers at Tessella for the call to qt() - Calculate asymetric confidence intervals for backtransformed parameters - Overhaul and recompile vignettes git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@74 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- vignettes/examples.Rnw | 141 +++++++++++++++++++++++++++---------------------- vignettes/examples.pdf | Bin 281843 -> 293191 bytes vignettes/mkin.pdf | Bin 162843 -> 134756 bytes 3 files changed, 79 insertions(+), 62 deletions(-) (limited to 'vignettes') diff --git a/vignettes/examples.Rnw b/vignettes/examples.Rnw index 6f3cfc9d..fa35c8f1 100644 --- a/vignettes/examples.Rnw +++ b/vignettes/examples.Rnw @@ -93,13 +93,13 @@ summary(m.L1.SFO) A plot of the fit is obtained with the plot function for mkinfit objects. -<>= +<>= plot(m.L1.SFO) @ The residual plot can be easily obtained by -<>= +<>= mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time") @ @@ -113,7 +113,8 @@ summary(m.L1.FOMC) Due to the higher number of parameters, and the lower number of degrees of freedom of the fit, the $\chi^2$ error level is actually higher for the FOMC model (3.6\%) than -for the SFO model (3.4\%). +for the SFO model (3.4\%). Additionally, the covariance matrix can not be obtained, +indicating overparameterisation of the model. The $\chi^2$ error levels reported in Appendix 3 and Appendix 7 to the FOCUS kinetics report are rounded to integer percentages and partly deviate by one percentage point @@ -167,22 +168,22 @@ models generally only implement SFO kinetics. For comparison, the FOMC model is fitted as well, and the $\chi^2$ error level is checked. -<>= -m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet=TRUE) +<>= +m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet = TRUE) +par(mfrow = c(2, 1)) plot(m.L2.FOMC) -s.m.L2.FOMC <- summary(m.L2.FOMC) -s.m.L2.FOMC$errmin +mkinresplot(m.L2.FOMC) +summary(m.L2.FOMC, data = FALSE) @ The error level at which the $\chi^2$ test passes is much lower in this case. Therefore, the FOMC model provides a better description of the data, as less experimental error has to be assumed in order to explain the data. -Fitting the four parameter DFOP model does not further reduce the -$\chi^2$ error level. +Fitting the four parameter DFOP model further reduces the $\chi^2$ error level. -<>= -m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet=TRUE) +<>= +m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet = TRUE) plot(m.L2.DFOP) @ @@ -190,17 +191,18 @@ Here, the default starting parameters for the DFOP model obviously do not lead to a reasonable solution. Therefore the fit is repeated with different starting parameters. -<>= +<>= m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, parms.ini = c(k1 = 1, k2 = 0.01, g = 0.8), quiet=TRUE) plot(m.L2.DFOP) -s.m.L2.DFOP <- summary(m.L2.DFOP) -s.m.L2.DFOP$errmin +summary(m.L2.DFOP, data = FALSE) @ -Therefore, the FOMC model is clearly the best-fit model for dataset L1 based on the -$\chi^2$ error level criterion. +Here, the DFOP model is clearly the best-fit model for dataset L2 based on the +$\chi^2$ error level criterion. However, the failure to calculate the covariance +matrix indicates that the parameter estimates correlate excessively. Therefore, +the FOMC model may be preferred for this dataset. \subsection{Laboratory Data L3} @@ -216,10 +218,10 @@ FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3) SFO model, summary and plot: -<>= +<>= m.L3.SFO <- mkinfit(SFO, FOCUS_2006_L3_mkin, quiet = TRUE) -summary(m.L3.SFO) plot(m.L3.SFO) +summary(m.L3.SFO) @ The $\chi^2$ error level of 22\% as well as the plot suggest that the model @@ -227,12 +229,10 @@ does not fit very well. The FOMC model performs better: -<>= +<>= m.L3.FOMC <- mkinfit(FOMC, FOCUS_2006_L3_mkin, quiet = TRUE) plot(m.L3.FOMC) -s.m.L3.FOMC <- summary(m.L3.FOMC) -s.m.L3.FOMC$errmin -endpoints(m.L3.FOMC) +summary(m.L3.FOMC, data = FALSE) @ The error level at which the $\chi^2$ test passes is 7\% in this case. @@ -240,14 +240,15 @@ The error level at which the $\chi^2$ test passes is 7\% in this case. Fitting the four parameter DFOP model further reduces the $\chi^2$ error level considerably: -<>= +<>= m.L3.DFOP <- mkinfit(DFOP, FOCUS_2006_L3_mkin, quiet = TRUE) plot(m.L3.DFOP) -s.m.L3.DFOP <- summary(m.L3.DFOP) -s.m.L3.DFOP$errmin +summary(m.L3.DFOP, data = FALSE) @ -Therefore, the DFOP model is the best-fit model based on the $\chi^2$ error +Here, a look to the model plot, the confidence intervals of the parameters +and the correlation matrix suggest that the paramter estimates are reliable, and +the DFOP model can be used as the best-fit model based on the $\chi^2$ error level criterion for laboratory data L3. \subsection{Laboratory Data L4} @@ -264,10 +265,10 @@ FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4) SFO model, summary and plot: -<>= +<>= m.L4.SFO <- mkinfit(SFO, FOCUS_2006_L4_mkin, quiet = TRUE) -summary(m.L4.SFO) plot(m.L4.SFO) +summary(m.L4.SFO, data = FALSE) @ The $\chi^2$ error level of 3.3\% as well as the plot suggest that the model @@ -275,11 +276,10 @@ fits very well. The FOMC model for comparison -<>= +<>= m.L4.FOMC <- mkinfit(FOMC, FOCUS_2006_L4_mkin, quiet = TRUE) plot(m.L4.FOMC) -s.m.L4.FOMC <- summary(m.L4.FOMC) -s.m.L4.FOMC$errmin +summary(m.L4.FOMC, data = FALSE) @ The error level at which the $\chi^2$ test passes is slightly lower for the FOMC @@ -315,29 +315,29 @@ Step 1 (SFO for parent only) is skipped here. We start with the model 2a, with formation and decline of metabolite Z1 and the pathway from parent directly to sink included (default in mkin). -<>= +<>= Z.2a <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), Z1 = list(type = "SFO")) m.Z.2a <- mkinfit(Z.2a, FOCUS_2006_Z_mkin, quiet = TRUE) -summary(m.Z.2a, data = FALSE) plot(m.Z.2a) +summary(m.Z.2a, data = FALSE) @ As obvious from the summary, the kinetic rate constant from parent compound Z to sink is negligible. Accordingly, the exact magnitude of the fitted parameter \texttt{log k\_Z\_sink} is ill-defined and the covariance matrix is not returned. -This suggests, in agreement to the analysis in the FOCUS kinetics report, to simplify +This suggests, in agreement with the analysis in the FOCUS kinetics report, to simplify the model by removing the pathway to sink. A similar result can be obtained when formation fractions are used in the model formulation: -<>= +<>= Z.2a.ff <- mkinmod(Z0 = list(type = "SFO", to = "Z1"), Z1 = list(type = "SFO"), use_of_ff = "max") m.Z.2a.ff <- mkinfit(Z.2a.ff, FOCUS_2006_Z_mkin, quiet = TRUE) -summary(m.Z.2a.ff, data = FALSE) plot(m.Z.2a.ff) +summary(m.Z.2a.ff, data = FALSE) @ Here, the ilr transformed formation fraction fitted in the model takes a very large value, @@ -348,18 +348,17 @@ The simplified model is obtained by setting the list component \texttt{sink} to \texttt{FALSE}. This model definition is not supported when formation fractions are used. -<>= +<>= Z.3 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), Z1 = list(type = "SFO")) m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5), quiet = TRUE) -m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, solution_type = "deSolve", - quiet = TRUE) -summary(m.Z.3, data = FALSE) +#m.Z.3 <- mkinfit(Z.3, FOCUS_2006_Z_mkin, solution_type = "deSolve") plot(m.Z.3) +summary(m.Z.3, data = FALSE) @ -The first attempt to fit the model fails, as the default solution type chosen +The first attempt to fit the model failed, as the default solution type chosen by mkinfit is based on eigenvalues, and the system defined by the starting parameters is identified as being singular to the solver. This is caused by the fact that the rate constants for both state variables are the same using the @@ -367,36 +366,36 @@ default starting paramters. Setting a different starting value for one of the parameters overcomes this. Alternatively, the \Rpackage{deSolve} based model solution can be chosen, at the cost of a bit more computing time. -<>= +<>= Z.4a <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2"), Z2 = list(type = "SFO")) m.Z.4a <- mkinfit(Z.4a, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5), quiet = TRUE) -summary(m.Z.4a, data = FALSE) plot(m.Z.4a) +summary(m.Z.4a, data = FALSE) @ As suggested in the FOCUS report, the pathway to sink was removed for metabolite Z1 as -well. While this step appears questionable on the basis of the above results, it +well in the next step. While this step appears questionable on the basis of the above results, it is followed here for the purpose of comparison. Also, in the FOCUS report, it is assumed that there is additional empirical evidence that Z1 quickly and exclusively hydrolyses to Z2. Again, in order to avoid a singular system when using default starting parameters, the starting parameter for the pathway without sink term has to be adapted. -<>= +<>= Z.5 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2", sink = FALSE), Z2 = list(type = "SFO")) m.Z.5 <- mkinfit(Z.5, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.2), quiet = TRUE) -summary(m.Z.5, data = FALSE) plot(m.Z.5) +summary(m.Z.5, data = FALSE) @ Finally, metabolite Z3 is added to the model. -<>= +<>= Z.FOCUS <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2", sink = FALSE), Z2 = list(type = "SFO", to = "Z3"), @@ -404,8 +403,8 @@ Z.FOCUS <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.2, k_Z2_Z3 = 0.3), quiet = TRUE) -summary(m.Z.FOCUS, data = FALSE) plot(m.Z.FOCUS) +summary(m.Z.FOCUS, data = FALSE) @ This is the fit corresponding to the final result chosen in Appendix 7 of the @@ -428,7 +427,7 @@ reversible binding (SFORB) model for metabolite Z3. As expected, the $\chi^2$ error level is lower for metabolite Z3 using this model and the graphical fit for Z3 is improved. However, the covariance matrix is not returned. -<>= +<>= Z.mkin.1 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2", sink = FALSE), Z2 = list(type = "SFO", to = "Z3"), @@ -436,46 +435,49 @@ Z.mkin.1 <- mkinmod(Z0 = list(type = "SFO", to = "Z1", sink = FALSE), m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, parms.ini = c(k_Z0_Z1 = 0.5, k_Z1_Z2 = 0.3, k_Z2_Z3 = 0.2), quiet = TRUE) -summary(m.Z.mkin.1, data = FALSE) plot(m.Z.mkin.1) +summary(m.Z.mkin.1, data = FALSE) @ Therefore, a further stepwise model building is performed starting from the stage of parent and one metabolite, starting from the assumption that the model fit for the parent compound can be improved by using the SFORB model. -<>= +<>= Z.mkin.2 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), Z1 = list(type = "SFO")) m.Z.mkin.2 <- mkinfit(Z.mkin.2, FOCUS_2006_Z_mkin, quiet = TRUE) -summary(m.Z.mkin.2, data = FALSE) plot(m.Z.mkin.2) +summary(m.Z.mkin.2, data = FALSE) @ -The sink is for Z1 is turned off again, for the same reasons as in the original analysis. -Then, metabolite Z2 is added. +When metabolite Z2 is added, the additional sink for Z1 is turned off again, +for the same reasons as in the original analysis. -<>= +<>= Z.mkin.3 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2"), Z2 = list(type = "SFO")) m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE) -summary(m.Z.mkin.3, data = FALSE) plot(m.Z.mkin.3) +summary(m.Z.mkin.3, data = FALSE) @ +This results in a much better representation of the behaviour of the parent +compound Z0. + Finally, Z3 is added as well. This model appears overparameterised (no covariance matrix returned) if the sink for Z1 is left in the model. -<>= +<>= Z.mkin.4 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2", sink = FALSE), Z2 = list(type = "SFO", to = "Z3"), Z3 = list(type = "SFO")) m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, parms.ini = c(k_Z1_Z2 = 0.05), quiet = TRUE) -summary(m.Z.mkin.4, data = FALSE) plot(m.Z.mkin.4) +summary(m.Z.mkin.4, data = FALSE) @ The error level of the fit, but especially of metabolite Z3, can be improved if @@ -483,20 +485,35 @@ the SFORB model is chosen for this metabolite, as this model is capable of representing the tailing of the metabolite decline phase. Using the SFORB additionally for Z1 or Z2 did not further improve the result. -Therefore, the model \texttt{Z.mkin.5} is proposed as the best-fit model -for the dataset from Appendix 7 of the FOCUS report. -<>= +<>= Z.mkin.5 <- mkinmod(Z0 = list(type = "SFORB", to = "Z1", sink = FALSE), Z1 = list(type = "SFO", to = "Z2", sink = FALSE), Z2 = list(type = "SFO", to = "Z3"), Z3 = list(type = "SFORB")) m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini = c(k_Z1_Z2 = 0.2), quiet = TRUE) -summary(m.Z.mkin.5, data = FALSE) plot(m.Z.mkin.5) +summary(m.Z.mkin.5, data = FALSE) @ +Looking at the confidence intervals of the SFORB model parameters of Z3, it is +clear that nothing can be said about a degradation rate of Z3. However, this +appears to be a feature of the data. + +<>= +par(mfrow = c(2, 2)) +mkinresplot(m.Z.mkin.5, "Z0", lpos = "bottomright") +mkinresplot(m.Z.mkin.5, "Z1", lpos = "bottomright") +mkinresplot(m.Z.mkin.5, "Z2", lpos = "bottomright") +mkinresplot(m.Z.mkin.5, "Z3", lpos = "bottomright") +@ + +As expected, the residual plots are much more random than in the case of the +all SFO model for which they were shown above. In conclusion, the model +\texttt{Z.mkin.5} is proposed as the best-fit model for the dataset from +Appendix 7 of the FOCUS report. + \bibliographystyle{plainnat} \bibliography{references} diff --git a/vignettes/examples.pdf b/vignettes/examples.pdf index 3bd0c37e..f4cfb31c 100644 Binary files a/vignettes/examples.pdf and b/vignettes/examples.pdf differ diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf index 073e3610..227ea79a 100644 Binary files a/vignettes/mkin.pdf and b/vignettes/mkin.pdf differ -- cgit v1.2.1