From 95178837d3f91e84837628446b5fd468179af2b9 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 4 Jun 2019 15:09:28 +0200 Subject: Additional algorithm "d_c", more tests, docs The new algorithm tries direct optimization of the likelihood, as well as a three step procedure. In this way, we consistently get the model with the highest likelihood for SFO, DFOP and HS for all 12 new test datasets. --- docs/reference/mkinfit.html | 205 ++++++++++++++++++++++++++------------------ 1 file changed, 121 insertions(+), 84 deletions(-) (limited to 'docs/reference/mkinfit.html') diff --git a/docs/reference/mkinfit.html b/docs/reference/mkinfit.html index bfca44fc..8cabcb21 100644 --- a/docs/reference/mkinfit.html +++ b/docs/reference/mkinfit.html @@ -71,7 +71,7 @@ Per default, parameters in the kinetic models are internally transformed in mkin - 0.9.49.4 + 0.9.49.5 @@ -163,6 +163,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"), + reweight.tol = 1e-8, reweight.max.iter = 10, trace_parms = FALSE, ...)

Arguments

@@ -321,6 +323,41 @@ Per default, parameters in the kinetic models are internally transformed in errors follow a lognormal distribution for large values, not a normal distribution as assumed by this method.

+ + 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 can be applied.

+

The default 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 + log-likelihood.

+

The algorithm "twostep" will minimize the negative log-likelihood + after an initial unweighted leas squares optimisation step.

+

The algorithm "threestep" starts with unweighted least squares, + then optimizes only the error model using the degradation model + parameters found, and then minimizes the negative log-likelihood + with free degradation and error model parameters.

+

The algorithm "fourstep" starts with unweighted least squares, + then optimizes only the error model using the degradation model + parameters found, then optimizes the degradation model again + with fixed error model parameters, and finally minimizes the negative + log-likelihood with free degradation and error model parameters.

+

The algorithm "IRLS" starts with 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.

+ + + reweight.tol +

Tolerance for the convergence criterion calculated from the error model + parameters in IRLS fits.

+ + + reweight.max.iter +

Maximum number of iterations in IRLS fits.

+ trace_parms

Should a trace of the parameter values be listed?

@@ -360,17 +397,17 @@ Per default, parameters in the kinetic models are internally transformed in

Examples

# Use shorthand notation for parent only degradation
 fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
-summary(fit)
#> mkin version used for fitting:    0.9.49.4 
+summary(fit)
#> mkin version used for fitting:    0.9.49.5 
 #> R version used for fitting:       3.6.0 
-#> Date of fit:     Wed May  8 20:50:50 2019 
-#> Date of summary: Wed May  8 20:50:50 2019 
+#> Date of fit:     Tue Jun  4 15:01:15 2019 
+#> Date of summary: Tue Jun  4 15:01:15 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.456 s
+#> Fitted using 222 model solutions performed in 0.461 s
 #> 
 #> Error model:
 #> Constant variance 
@@ -443,7 +480,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.488       0.000       1.488 
coef(fit)
#> NULL
endpoints(fit)
#> $ff
+#>       1.521       0.000       1.526 
coef(fit)
#> NULL
endpoints(fit)
#> $ff
 #> parent_sink   parent_m1     m1_sink 
 #>    0.485524    0.514476    1.000000 
 #> 
@@ -515,7 +552,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.086       0.000       1.087 
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff
+#>       1.093       0.000       1.093 
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff
 #> parent_sink   parent_m1     m1_sink 
 #>    0.485524    0.514476    1.000000 
 #> 
@@ -545,10 +582,10 @@ Per default, parameters in the kinetic models are internally transformed in
 fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data

 
# Weighted fits, including IRLS
 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.4 
+                      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.5 
 #> R version used for fitting:       3.6.0 
-#> Date of fit:     Wed May  8 20:51:06 2019 
-#> Date of summary: Wed May  8 20:51:06 2019 
+#> Date of fit:     Tue Jun  4 15:01:31 2019 
+#> Date of summary: Tue Jun  4 15:01:31 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -556,7 +593,7 @@ 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.082 s
+#> Fitted using 421 model solutions performed in 1.096 s
 #> 
 #> Error model:
 #> Constant variance 
@@ -663,10 +700,10 @@ Per default, parameters in the kinetic models are internally transformed in
 #>   100       m1    31.04  31.98163 -9.416e-01
 #>   100       m1    33.13  31.98163  1.148e+00
 #>   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.4 
+#>   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.5 
 #> R version used for fitting:       3.6.0 
-#> Date of fit:     Wed May  8 20:51:08 2019 
-#> Date of summary: Wed May  8 20:51:08 2019 
+#> Date of fit:     Tue Jun  4 15:01:34 2019 
+#> Date of summary: Tue Jun  4 15:01:34 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -674,19 +711,19 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted using 758 model solutions performed in 1.971 s
+#> Fitted using 979 model solutions performed in 2.603 s
 #> 
 #> Error model:
 #> Variance unique to each observed variable 
 #> 
 #> Starting values for parameters to be optimised:
-#>                   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
+#>                     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
 #> 
 #> Starting values for the transformed parameters actually optimised:
 #>                     value lower upper
@@ -694,8 +731,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.000000     0   Inf
-#> sigma_m1         3.000000     0   Inf
+#> sigma_parent     3.398909     0   Inf
+#> sigma_m1         2.857157     0   Inf
 #> 
 #> Fixed parameter values:
 #>      value  type
@@ -715,14 +752,14 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 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.041368
+#> 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.03925    -0.004628
+#> sigma_m1       -0.03385     -0.06627  0.02482        0.03926    -0.004628
 #>                 sigma_m1
 #> parent_0       -0.033847
 #> log_k_parent   -0.066265
-#> log_k_m1        0.024821
-#> f_parent_ilr_1  0.039255
+#> log_k_m1        0.024823
+#> f_parent_ilr_1  0.039256
 #> sigma_parent   -0.004628
 #> sigma_m1        1.000000
 #> 
@@ -786,17 +823,17 @@ Per default, parameters in the kinetic models are internally transformed in
 #>    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.481e-03
+#>    50       m1    41.19  41.19948 -9.479e-03
 #>    50       m1    40.01  41.19948 -1.189e+00
-#>    75       m1    40.09  36.44036  3.650e+00
-#>    75       m1    33.85  36.44036 -2.590e+00
-#>   100       m1    31.04  31.98774 -9.477e-01
-#>   100       m1    33.13  31.98774  1.142e+00
-#>   120       m1    25.15  28.80430 -3.654e+00
-#>   120       m1    33.31  28.80430  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.4 
+#>    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.5 
 #> R version used for fitting:       3.6.0 
-#> Date of fit:     Wed May  8 20:51:11 2019 
-#> Date of summary: Wed May  8 20:51:11 2019 
+#> Date of fit:     Tue Jun  4 15:01:43 2019 
+#> Date of summary: Tue Jun  4 15:01:43 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -804,28 +841,28 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted using 821 model solutions performed in 3.29 s
+#> Fitted using 2289 model solutions performed in 9.499 s
 #> 
 #> Error model:
 #> Two-component variance function 
 #> 
 #> Starting values for parameters to be optimised:
-#>                   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
+#>                       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
 #> 
 #> Starting values for the transformed parameters actually optimised:
-#>                     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
+#>                        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
 #> 
 #> Fixed parameter values:
 #>      value  type
@@ -856,7 +893,7 @@ Per default, parameters in the kinetic models are internally transformed in
 #>                 Estimate  t value    Pr(>t)     Lower     Upper
 #> parent_0       1.007e+02  38.4300 1.180e-28 95.400000 1.061e+02
 #> k_parent       1.006e-01 112.8000 1.718e-43  0.098760 1.024e-01
-#> k_m1           5.167e-03  10.9500 1.172e-12  0.004290 6.223e-03
+#> k_m1           5.167e-03  10.9500 1.171e-12  0.004290 6.223e-03
 #> f_parent_to_m1 5.084e-01  26.0100 2.146e-23  0.468600 5.481e-01
 #> sigma_low      3.050e-03   0.6314 2.661e-01 -0.006786 1.289e-02
 #> rsd_high       7.928e-02   8.4170 6.418e-10  0.060100 9.847e-02
@@ -879,18 +916,18 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Data:
 #>  time variable observed predicted   residual
-#>     0   parent    99.46 100.73433  -1.274329
-#>     0   parent   102.04 100.73433   1.305671
-#>     1   parent    93.50  91.09750   2.402495
-#>     1   parent    92.50  91.09750   1.402495
-#>     3   parent    63.23  74.50140 -11.271403
-#>     3   parent    68.99  74.50140  -5.511403
-#>     7   parent    52.32  49.82880   2.491205
-#>     7   parent    55.13  49.82880   5.301205
-#>    14   parent    27.27  24.64809   2.621909
-#>    14   parent    26.64  24.64809   1.991909
-#>    21   parent    11.50  12.19231  -0.692315
-#>    21   parent    11.64  12.19231  -0.552315
+#>     0   parent    99.46 100.73434  -1.274339
+#>     0   parent   102.04 100.73434   1.305661
+#>     1   parent    93.50  91.09751   2.402486
+#>     1   parent    92.50  91.09751   1.402486
+#>     3   parent    63.23  74.50141 -11.271410
+#>     3   parent    68.99  74.50141  -5.511410
+#>     7   parent    52.32  49.82880   2.491201
+#>     7   parent    55.13  49.82880   5.301201
+#>    14   parent    27.27  24.64809   2.621908
+#>    14   parent    26.64  24.64809   1.991908
+#>    21   parent    11.50  12.19232  -0.692315
+#>    21   parent    11.64  12.19232  -0.552315
 #>    35   parent     2.85   2.98327  -0.133266
 #>    35   parent     2.91   2.98327  -0.073266
 #>    50   parent     0.69   0.66013   0.029874
@@ -899,24 +936,24 @@ Per default, parameters in the kinetic models are internally transformed in
 #>    75   parent     0.06   0.05344   0.006562
 #>     1       m1     4.84   4.88645  -0.046451
 #>     1       m1     5.64   4.88645   0.753549
-#>     3       m1    12.91  13.22867  -0.318668
-#>     3       m1    12.96  13.22867  -0.268668
-#>     7       m1    22.97  25.36416  -2.394164
-#>     7       m1    24.47  25.36416  -0.894164
-#>    14       m1    41.69  37.00974   4.680265
-#>    14       m1    33.21  37.00974  -3.799735
-#>    21       m1    44.37  41.90133   2.468670
-#>    21       m1    46.44  41.90133   4.538670
-#>    35       m1    41.22  43.45691  -2.236914
-#>    35       m1    37.95  43.45691  -5.506914
-#>    50       m1    41.19  41.34199  -0.151988
-#>    50       m1    40.01  41.34199  -1.331988
-#>    75       m1    40.09  36.61471   3.475290
-#>    75       m1    33.85  36.61471  -2.764710
-#>   100       m1    31.04  32.20083  -1.160830
-#>   100       m1    33.13  32.20083   0.929170
-#>   120       m1    25.15  29.04131  -3.891312
-#>   120       m1    33.31  29.04131   4.268688

+#>     3       m1    12.91  13.22867  -0.318669
+#>     3       m1    12.96  13.22867  -0.268669
+#>     7       m1    22.97  25.36417  -2.394166
+#>     7       m1    24.47  25.36417  -0.894166
+#>    14       m1    41.69  37.00974   4.680263
+#>    14       m1    33.21  37.00974  -3.799737
+#>    21       m1    44.37  41.90133   2.468669
+#>    21       m1    46.44  41.90133   4.538669
+#>    35       m1    41.22  43.45691  -2.236913
+#>    35       m1    37.95  43.45691  -5.506913
+#>    50       m1    41.19  41.34199  -0.151985
+#>    50       m1    40.01  41.34199  -1.331985
+#>    75       m1    40.09  36.61471   3.475295
+#>    75       m1    33.85  36.61471  -2.764705
+#>   100       m1    31.04  32.20082  -1.160823
+#>   100       m1    33.13  32.20082   0.929177
+#>   120       m1    25.15  29.04130  -3.891304
+#>   120       m1    33.31  29.04130   4.268696
-

Site built with pkgdown 1.3.0.9000.

+

Site built with pkgdown 1.3.0.

-- cgit v1.2.1