- 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

3 files changed, 79 insertions, 62 deletions

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.

 

-<<L1_SFO_plot, fig=TRUE, echo=TRUE>>=

+<<L1_SFO_plot, fig=TRUE, echo=TRUE, height=4>>=

 plot(m.L1.SFO)

 @

 

 The residual plot can be easily obtained by

 

-<<L1_SFO_residuals, fig=TRUE, echo=TRUE>>=

+<<L1_SFO_residuals, fig=TRUE, echo=TRUE, height=4>>=

 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.

 

-<<L2_FOMC, echo=TRUE, fig=TRUE>>=

-m.L2.FOMC <- mkinfit(FOMC, FOCUS_2006_L2_mkin, quiet=TRUE)

+<<L2_FOMC, echo=TRUE, fig=TRUE, height=8>>=

+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. 

 

-<<L2_DFOP, echo=TRUE, fig=TRUE>>=

-m.L2.DFOP <- mkinfit(DFOP, FOCUS_2006_L2_mkin, quiet=TRUE)

+<<L2_DFOP, echo=TRUE, fig=TRUE, height=4>>=

+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.

 

-<<L2_DFOP_2, echo=TRUE, fig=TRUE>>=

+<<L2_DFOP_2, echo=TRUE, fig=TRUE, height=4>>=

 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:

 

-<<L3_SFO, echo=TRUE, fig=TRUE>>=

+<<L3_SFO, echo=TRUE, fig=TRUE, height=4>>=

 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:

 

-<<L3_FOMC, echo=TRUE, fig=TRUE>>=

+<<L3_FOMC, echo=TRUE, fig=TRUE, height=4>>=

 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:

 

-<<L3_DFOP, echo=TRUE, fig=TRUE>>=

+<<L3_DFOP, echo=TRUE, fig=TRUE, height=4>>=

 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:

 

-<<L4_SFO, echo=TRUE, fig=TRUE>>=

+<<L4_SFO, echo=TRUE, fig=TRUE, height=4>>=

 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

 

-<<L4_FOMC, echo=TRUE, fig=TRUE>>=

+<<L4_FOMC, echo=TRUE, fig=TRUE, height=4>>=

 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).

 

-<<FOCUS_2006_Z_fits_1, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_1, echo=TRUE, fig=TRUE, height=4>>=

 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:

 

-<<FOCUS_2006_Z_fits_2, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_2, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_3, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_3, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_4, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_4, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_5, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_5, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_6, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_6, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_7, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_7, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_8, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_8, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_9, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_9, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_10, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_10, echo=TRUE, fig=TRUE, height=4>>=

 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.

 

-<<FOCUS_2006_Z_fits_11, echo=TRUE, fig=TRUE>>=

+<<FOCUS_2006_Z_fits_11, echo=TRUE, fig=TRUE, height=4>>=

 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.

+

+<<FOCUS_2006_Z_residuals_11, fig=TRUE>>=

+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
--- a/vignettes/examples.pdf
+++ b/vignettes/examples.pdf
diff --git a/vignettes/mkin.pdf b/vignettes/mkin.pdf
index 073e3610..227ea79a 100644
--- a/vignettes/mkin.pdf
+++ b/vignettes/mkin.pdf