From 7624a2b8398b4ad665a3b0b622488e1893a5ee7c Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 21 Oct 2019 12:11:34 +0200 Subject: Refactor mkinfit, infrastructure work mkinfit objects now include an ll() function to calculate the log-likelihood. Part of the code was refactored, hopefully making it easier to read and maintain. IRLS is currently the default algorithm for the error model "obs", for no particular reason. This may be subject to change when I get around to investigate. Slow tests are now in a separate subdirectory and will probably only be run by my own Makefile target. Formatting of test logs is improved. Roundtripping error model parameters works with a precision of 10% when we use lots of replicates in the synthetic data (see slow tests). This is not new in this commit, but as I think it is reasonable this closes #7. --- docs/reference/mkinfit.html | 172 +++++++++++++++++++++----------------------- 1 file changed, 83 insertions(+), 89 deletions(-) (limited to 'docs/reference/mkinfit.html') diff --git a/docs/reference/mkinfit.html b/docs/reference/mkinfit.html index 0767d2f8..e07d1f13 100644 --- a/docs/reference/mkinfit.html +++ b/docs/reference/mkinfit.html @@ -167,8 +167,8 @@ Per default, parameters in the kinetic models are internally transformed in quiet = FALSE, atol = 1e-8, rtol = 1e-10, n.outtimes = 100, error_model = c("const", "obs", "tc"), - error_model_algorithm = c("d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", - "OLS"), + error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", + "fourstep", "IRLS", "OLS"), reweight.tol = 1e-8, reweight.max.iter = 10, trace_parms = FALSE, ...) @@ -330,11 +330,12 @@ Per default, parameters in the kinetic models are internally transformed in error_model_algorithm -

If the error model is "const", the error model algorithm is ignored, - because no special algorithm is needed and unweighted (also known as - ordinary) least squares fitting (listed as "OLS" in the summary) can be - applied.

-

The default algorithm "d_3" will directly minimize the negative +

If "auto", the selected algorithm depends on the error model. + If the error model is "const", nonlinear least squares fitting ("OLS") is + selected. If the error model is "obs", iteratively reweighted least squares + fitting ("IRLS") is selected. If the error model is "tc", the "d_3" + algorithm is selected.

+

The algorithm "d_3" will directly minimize the negative log-likelihood and - independently - also use the three step algorithm described below. The fit with the higher likelihood is returned.

The algorithm "direct" will directly minimize the negative @@ -354,9 +355,7 @@ Per default, parameters in the kinetic models are internally transformed in unweighted least squares, and then iterates optimization of the error model parameters and subsequent optimization of the degradation model using those error model parameters, - until the error model parameters converge.

-

The algorithm "OLS" (Ordinary Least Squares) is automatically selected when - the error model is "const" and results in an unweighted least squares fit.

+ until the error model parameters converge.

reweight.tol @@ -408,33 +407,31 @@ Per default, parameters in the kinetic models are internally transformed in fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) summary(fit)
#> mkin version used for fitting: 0.9.49.6 #> R version used for fitting: 3.6.1 -#> Date of fit: Thu Sep 19 09:50:54 2019 -#> Date of summary: Thu Sep 19 09:50:54 2019 +#> Date of fit: Mon Oct 21 12:07:39 2019 +#> Date of summary: Mon Oct 21 12:07:39 2019 #> #> Equations: #> d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent #> #> Model predictions using solution type analytical #> -#> Fitted using 222 model solutions performed in 0.458 s +#> Fitted using 222 model solutions performed in 0.459 s #> #> Error model: Constant variance #> #> Error model algorithm: OLS #> #> Starting values for parameters to be optimised: -#> value type -#> parent_0 85.100000 state -#> alpha 1.000000 deparm -#> beta 10.000000 deparm -#> sigma 1.857444 error +#> value type +#> parent_0 85.1 state +#> alpha 1.0 deparm +#> beta 10.0 deparm #> #> Starting values for the transformed parameters actually optimised: #> value lower upper #> parent_0 85.100000 -Inf Inf #> log_alpha 0.000000 -Inf Inf #> log_beta 2.302585 -Inf Inf -#> sigma 1.857444 0 Inf #> #> Fixed parameter values: #> None @@ -490,7 +487,7 @@ Per default, parameters in the kinetic models are internally transformed in m1 = mkinsub("SFO"))
#> Successfully compiled differential equation model from auto-generated C code.
# Fit the model to the FOCUS example dataset D using defaults print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "eigen", quiet = TRUE)))
#> Warning: Observations with value of zero were removed from the data
#> User System verstrichen -#> 1.479 0.002 1.482
coef(fit)
#> NULL
endpoints(fit)
#> $ff +#> 1.462 0.000 1.463
coef(fit)
#> NULL
endpoints(fit)
#> $ff #> parent_sink parent_m1 m1_sink #> 0.485524 0.514476 1.000000 #> @@ -563,7 +560,7 @@ Per default, parameters in the kinetic models are internally transformed in #> Sum of squared residuals at call 126: 371.2134 #> Sum of squared residuals at call 135: 371.2134 #> Negative log-likelihood at call 145: 97.22429
#> Optimisation successfully terminated.
#> User System verstrichen -#> 1.053 0.000 1.054
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff +#> 1.04 0.00 1.04
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff #> parent_sink parent_m1 m1_sink #> 0.485524 0.514476 1.000000 #> @@ -599,8 +596,8 @@ Per default, parameters in the kinetic models are internally transformed in SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), use_of_ff = "max")
#> Successfully compiled differential equation model from auto-generated C code.
f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
summary(f.noweight)
#> mkin version used for fitting: 0.9.49.6 #> R version used for fitting: 3.6.1 -#> Date of fit: Thu Sep 19 09:51:10 2019 -#> Date of summary: Thu Sep 19 09:51:10 2019 +#> Date of fit: Mon Oct 21 12:07:54 2019 +#> Date of summary: Mon Oct 21 12:07:54 2019 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -608,19 +605,18 @@ Per default, parameters in the kinetic models are internally transformed in #> #> Model predictions using solution type deSolve #> -#> Fitted using 421 model solutions performed in 1.138 s +#> Fitted using 421 model solutions performed in 1.07 s #> #> Error model: Constant variance #> #> Error model algorithm: OLS #> #> Starting values for parameters to be optimised: -#> value type -#> parent_0 100.750000 state -#> k_parent 0.100000 deparm -#> k_m1 0.100100 deparm -#> f_parent_to_m1 0.500000 deparm -#> sigma 3.125504 error +#> value type +#> parent_0 100.7500 state +#> k_parent 0.1000 deparm +#> k_m1 0.1001 deparm +#> f_parent_to_m1 0.5000 deparm #> #> Starting values for the transformed parameters actually optimised: #> value lower upper @@ -628,7 +624,6 @@ Per default, parameters in the kinetic models are internally transformed in #> log_k_parent -2.302585 -Inf Inf #> log_k_m1 -2.301586 -Inf Inf #> f_parent_ilr_1 0.000000 -Inf Inf -#> sigma 3.125504 0 Inf #> #> Fixed parameter values: #> value type @@ -718,8 +713,8 @@ Per default, parameters in the kinetic models are internally transformed in #> 120 m1 25.15 28.78984 -3.640e+00 #> 120 m1 33.31 28.78984 4.520e+00
f.obs <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "obs", quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
summary(f.obs)
#> mkin version used for fitting: 0.9.49.6 #> R version used for fitting: 3.6.1 -#> Date of fit: Thu Sep 19 09:51:12 2019 -#> Date of summary: Thu Sep 19 09:51:12 2019 +#> Date of fit: Mon Oct 21 12:07:56 2019 +#> Date of summary: Mon Oct 21 12:07:56 2019 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -727,21 +722,20 @@ Per default, parameters in the kinetic models are internally transformed in #> #> Model predictions using solution type deSolve #> -#> Fitted using 979 model solutions performed in 2.565 s +#> Fitted using 897 model solutions performed in 2.363 s #> #> Error model: Variance unique to each observed variable #> -#> Error model algorithm: d_3 -#> Direct fitting and three-step fitting yield approximately the same likelihood +#> Error model algorithm: IRLS #> #> Starting values for parameters to be optimised: -#> value type -#> parent_0 100.750000 state -#> k_parent 0.100000 deparm -#> k_m1 0.100100 deparm -#> f_parent_to_m1 0.500000 deparm -#> sigma_parent 3.398909 error -#> sigma_m1 2.857157 error +#> value type +#> parent_0 100.7500 state +#> k_parent 0.1000 deparm +#> k_m1 0.1001 deparm +#> f_parent_to_m1 0.5000 deparm +#> sigma_parent 3.0000 error +#> sigma_m1 3.0000 error #> #> Starting values for the transformed parameters actually optimised: #> value lower upper @@ -749,8 +743,8 @@ Per default, parameters in the kinetic models are internally transformed in #> log_k_parent -2.302585 -Inf Inf #> log_k_m1 -2.301586 -Inf Inf #> f_parent_ilr_1 0.000000 -Inf Inf -#> sigma_parent 3.398909 0 Inf -#> sigma_m1 2.857157 0 Inf +#> sigma_parent 3.000000 0 Inf +#> sigma_m1 3.000000 0 Inf #> #> Fixed parameter values: #> value type @@ -767,16 +761,16 @@ Per default, parameters in the kinetic models are internally transformed in #> #> Parameter correlation: #> parent_0 log_k_parent log_k_m1 f_parent_ilr_1 sigma_parent -#> parent_0 1.00000 0.51078 -0.19133 -0.59997 0.035670 -#> log_k_parent 0.51078 1.00000 -0.37458 -0.59239 0.069833 -#> log_k_m1 -0.19133 -0.37458 1.00000 0.74398 -0.026158 -#> f_parent_ilr_1 -0.59997 -0.59239 0.74398 1.00000 -0.041369 -#> sigma_parent 0.03567 0.06983 -0.02616 -0.04137 1.000000 -#> sigma_m1 -0.03385 -0.06627 0.02482 0.03926 -0.004628 +#> parent_0 1.00000 0.51078 -0.19132 -0.59997 0.035676 +#> log_k_parent 0.51078 1.00000 -0.37457 -0.59239 0.069834 +#> log_k_m1 -0.19132 -0.37457 1.00000 0.74398 -0.026158 +#> f_parent_ilr_1 -0.59997 -0.59239 0.74398 1.00000 -0.041371 +#> sigma_parent 0.03568 0.06983 -0.02616 -0.04137 1.000000 +#> sigma_m1 -0.03385 -0.06626 0.02482 0.03926 -0.004628 #> sigma_m1 #> parent_0 -0.033847 #> log_k_parent -0.066265 -#> log_k_m1 0.024823 +#> log_k_m1 0.024824 #> f_parent_ilr_1 0.039256 #> sigma_parent -0.004628 #> sigma_m1 1.000000 @@ -790,7 +784,7 @@ Per default, parameters in the kinetic models are internally transformed in #> k_parent 0.098970 22.850 1.099e-21 0.090530 1.082e-01 #> k_m1 0.005245 8.046 1.732e-09 0.004072 6.756e-03 #> f_parent_to_m1 0.513600 23.560 4.352e-22 0.469300 5.578e-01 -#> sigma_parent 3.401000 5.985 5.662e-07 2.244000 4.559e+00 +#> sigma_parent 3.401000 5.985 5.661e-07 2.244000 4.559e+00 #> sigma_m1 2.855000 6.311 2.215e-07 1.934000 3.777e+00 #> #> FOCUS Chi2 error levels in percent: @@ -811,47 +805,47 @@ Per default, parameters in the kinetic models are internally transformed in #> #> Data: #> time variable observed predicted residual -#> 0 parent 99.46 99.65417 -1.942e-01 -#> 0 parent 102.04 99.65417 2.386e+00 -#> 1 parent 93.50 90.26332 3.237e+00 -#> 1 parent 92.50 90.26332 2.237e+00 -#> 3 parent 63.23 74.05306 -1.082e+01 -#> 3 parent 68.99 74.05306 -5.063e+00 -#> 7 parent 52.32 49.84325 2.477e+00 -#> 7 parent 55.13 49.84325 5.287e+00 +#> 0 parent 99.46 99.65420 -1.942e-01 +#> 0 parent 102.04 99.65420 2.386e+00 +#> 1 parent 93.50 90.26335 3.237e+00 +#> 1 parent 92.50 90.26335 2.237e+00 +#> 3 parent 63.23 74.05308 -1.082e+01 +#> 3 parent 68.99 74.05308 -5.063e+00 +#> 7 parent 52.32 49.84326 2.477e+00 +#> 7 parent 55.13 49.84326 5.287e+00 #> 14 parent 27.27 24.92971 2.340e+00 #> 14 parent 26.64 24.92971 1.710e+00 #> 21 parent 11.50 12.46890 -9.689e-01 #> 21 parent 11.64 12.46890 -8.289e-01 #> 35 parent 2.85 3.11925 -2.692e-01 #> 35 parent 2.91 3.11925 -2.092e-01 -#> 50 parent 0.69 0.70679 -1.679e-02 -#> 50 parent 0.63 0.70679 -7.679e-02 +#> 50 parent 0.69 0.70678 -1.678e-02 +#> 50 parent 0.63 0.70678 -7.678e-02 #> 75 parent 0.05 0.05952 -9.523e-03 #> 75 parent 0.06 0.05952 4.772e-04 #> 1 m1 4.84 4.81075 2.925e-02 #> 1 m1 5.64 4.81075 8.292e-01 -#> 3 m1 12.91 13.04196 -1.320e-01 -#> 3 m1 12.96 13.04196 -8.196e-02 -#> 7 m1 22.97 25.06847 -2.098e+00 -#> 7 m1 24.47 25.06847 -5.985e-01 +#> 3 m1 12.91 13.04197 -1.320e-01 +#> 3 m1 12.96 13.04197 -8.197e-02 +#> 7 m1 22.97 25.06848 -2.098e+00 +#> 7 m1 24.47 25.06848 -5.985e-01 #> 14 m1 41.69 36.70308 4.987e+00 #> 14 m1 33.21 36.70308 -3.493e+00 #> 21 m1 44.37 41.65115 2.719e+00 #> 21 m1 46.44 41.65115 4.789e+00 #> 35 m1 41.22 43.29465 -2.075e+00 #> 35 m1 37.95 43.29465 -5.345e+00 -#> 50 m1 41.19 41.19948 -9.479e-03 +#> 50 m1 41.19 41.19948 -9.477e-03 #> 50 m1 40.01 41.19948 -1.189e+00 #> 75 m1 40.09 36.44035 3.650e+00 #> 75 m1 33.85 36.44035 -2.590e+00 -#> 100 m1 31.04 31.98773 -9.477e-01 -#> 100 m1 33.13 31.98773 1.142e+00 -#> 120 m1 25.15 28.80429 -3.654e+00 -#> 120 m1 33.31 28.80429 4.506e+00
f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
summary(f.tc)
#> mkin version used for fitting: 0.9.49.6 +#> 100 m1 31.04 31.98772 -9.477e-01 +#> 100 m1 33.13 31.98772 1.142e+00 +#> 120 m1 25.15 28.80428 -3.654e+00 +#> 120 m1 33.31 28.80428 4.506e+00
f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
summary(f.tc)
#> mkin version used for fitting: 0.9.49.6 #> R version used for fitting: 3.6.1 -#> Date of fit: Thu Sep 19 09:51:22 2019 -#> Date of summary: Thu Sep 19 09:51:22 2019 +#> Date of fit: Mon Oct 21 12:08:06 2019 +#> Date of summary: Mon Oct 21 12:08:06 2019 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -859,7 +853,7 @@ Per default, parameters in the kinetic models are internally transformed in #> #> Model predictions using solution type deSolve #> -#> Fitted using 2289 model solutions performed in 9.24 s +#> Fitted using 2289 model solutions performed in 9.175 s #> #> Error model: Two-component variance function #> @@ -867,22 +861,22 @@ Per default, parameters in the kinetic models are internally transformed in #> Direct fitting and three-step fitting yield approximately the same likelihood #> #> Starting values for parameters to be optimised: -#> value type -#> parent_0 1.007500e+02 state -#> k_parent 1.000000e-01 deparm -#> k_m1 1.001000e-01 deparm -#> f_parent_to_m1 5.000000e-01 deparm -#> sigma_low 5.641148e-03 error -#> rsd_high 8.430766e-02 error +#> value type +#> parent_0 100.7500 state +#> k_parent 0.1000 deparm +#> k_m1 0.1001 deparm +#> f_parent_to_m1 0.5000 deparm +#> sigma_low 0.1000 error +#> rsd_high 0.1000 error #> #> Starting values for the transformed parameters actually optimised: -#> value lower upper -#> parent_0 100.750000000 -Inf Inf -#> log_k_parent -2.302585093 -Inf Inf -#> log_k_m1 -2.301585593 -Inf Inf -#> f_parent_ilr_1 0.000000000 -Inf Inf -#> sigma_low 0.005641148 0 Inf -#> rsd_high 0.084307660 0 Inf +#> value lower upper +#> parent_0 100.750000 -Inf Inf +#> log_k_parent -2.302585 -Inf Inf +#> log_k_m1 -2.301586 -Inf Inf +#> f_parent_ilr_1 0.000000 -Inf Inf +#> sigma_low 0.100000 0 Inf +#> rsd_high 0.100000 0 Inf #> #> Fixed parameter values: #> value type -- cgit v1.2.1