From 70591022c07f0e8fb4dd67789b7c8d78af8ebc18 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 2 May 2019 13:17:05 +0200 Subject: Better initials for error model parameters - Also make it possible to specify initial values for error model parameters. - Run tests - Rebuild docs --- docs/reference/mkinfit.html | 226 +++++++++++++++++++++++--------------------- 1 file changed, 117 insertions(+), 109 deletions(-) (limited to 'docs/reference/mkinfit.html') diff --git a/docs/reference/mkinfit.html b/docs/reference/mkinfit.html index d3a826b9..7693287e 100644 --- a/docs/reference/mkinfit.html +++ b/docs/reference/mkinfit.html @@ -151,6 +151,7 @@ Per default, parameters in the kinetic models are internally transformed in 
mkinfit(mkinmod, observed,
   parms.ini = "auto",
   state.ini = "auto",
+  err.ini = "auto",
   fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
   from_max_mean = FALSE,
   solution_type = c("auto", "analytical", "eigen", "deSolve"),
@@ -209,6 +210,13 @@ Per default, parameters in the kinetic models are internally transformed in
     values for the variable with the maximum observed value, and all others to 0.
     If this variable has no time zero observations, its initial value is set to 100.

     
+    
+      err.ini
+      
A named vector of initial values for the error model parameters to be
+    optimised.  If set to "auto", initial values are set to default values.
+    Otherwise, inital values for all error model parameters must be
+    given.

+    
     
       fixed_parms
       
The names of parameters that should not be optimised but rather kept at the
@@ -353,19 +361,19 @@ Per default, parameters in the kinetic models are internally transformed in
     
# 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 
-#> R version used for fitting:       3.5.3 
-#> Date of fit:     Wed Apr 10 10:10:01 2019 
-#> Date of summary: Wed Apr 10 10:10:01 2019 
+#> R version used for fitting:       3.6.0 
+#> Date of fit:     Thu May  2 12:38:56 2019 
+#> Date of summary: Thu May  2 12:38:56 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 #> 
 #> Model predictions using solution type analytical 
 #> 
-#> Fitted with method using 221 model solutions performed in 0.508 s
+#> Fitted using 221 model solutions performed in 0.455 s
 #> 
 #> Error model:
-#> NULL
+#> Constant variance 
 #> 
 #> Starting values for parameters to be optimised:
 #>              value   type
@@ -393,10 +401,10 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Parameter correlation:
 #>             parent_0  log_alpha   log_beta      sigma
-#> parent_0   1.000e+00 -1.565e-01 -3.142e-01 -1.313e-07
-#> log_alpha -1.565e-01  1.000e+00  9.564e-01 -2.634e-07
-#> log_beta  -3.142e-01  9.564e-01  1.000e+00 -2.200e-07
-#> sigma     -1.313e-07 -2.634e-07 -2.200e-07  1.000e+00
+#> parent_0   1.000e+00 -1.565e-01 -3.142e-01 -1.317e-07
+#> log_alpha -1.565e-01  1.000e+00  9.564e-01 -2.640e-07
+#> log_beta  -3.142e-01  9.564e-01  1.000e+00 -2.205e-07
+#> sigma     -1.317e-07 -2.640e-07 -2.205e-07  1.000e+00
 #> 
 #> Backtransformed parameters:
 #> Confidence intervals for internally transformed parameters are asymmetric.
@@ -408,7 +416,7 @@ Per default, parameters in the kinetic models are internally transformed in
 #> beta        1.917   4.373 3.601e-03  1.0650  3.451
 #> sigma       1.857   4.243 4.074e-03  0.7320  2.983
 #> 
-#> Chi2 error levels in percent:
+#> FOCUS Chi2 error levels in percent:
 #>          err.min n.optim df
 #> All data   6.657       3  6
 #> parent     6.657       3  6
@@ -435,7 +443,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.653       0.000       1.653 
coef(fit)
#> NULL
endpoints(fit)
#> $ff
+#>       1.551       0.000       1.552 
coef(fit)
#> NULL
endpoints(fit)
#> $ff
 #> parent_sink   parent_m1     m1_sink 
 #>    0.485524    0.514476    1.000000 
 #> 
@@ -448,71 +456,71 @@ Per default, parameters in the kinetic models are internally transformed in
 #> m1     131.760715 437.69962
 #> 
# deSolve is slower when no C compiler (gcc) was available during model generation
 print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D,
-                           solution_type = "deSolve")))
#> Warning: Observations with value of zero were removed from the data
#> Negative log-likelihood at call 1: 18915.53
-#> Negative log-likelihood at call 2: 18915.53
-#> Negative log-likelihood at call 6: 11424.02
-#> Negative log-likelihood at call 10: 11424
-#> Negative log-likelihood at call 13: 2367.052
-#> Negative log-likelihood at call 14: 2367.05
-#> Negative log-likelihood at call 19: 1314.716
-#> Negative log-likelihood at call 22: 1314.714
-#> Negative log-likelihood at call 25: 991.8311
-#> Negative log-likelihood at call 28: 991.8305
-#> Negative log-likelihood at call 30: 893.6462
-#> Negative log-likelihood at call 33: 893.6457
-#> Negative log-likelihood at call 35: 569.4049
-#> Negative log-likelihood at call 38: 569.4047
-#> Negative log-likelihood at call 40: 565.0651
-#> Negative log-likelihood at call 41: 565.065
-#> Negative log-likelihood at call 42: 565.0637
-#> Negative log-likelihood at call 45: 428.0188
-#> Negative log-likelihood at call 46: 428.0185
-#> Negative log-likelihood at call 50: 406.732
-#> Negative log-likelihood at call 52: 406.732
-#> Negative log-likelihood at call 55: 398.9115
-#> Negative log-likelihood at call 57: 398.9113
-#> Negative log-likelihood at call 60: 394.5943
-#> Negative log-likelihood at call 62: 394.5943
-#> Negative log-likelihood at call 66: 385.26
-#> Negative log-likelihood at call 67: 385.2599
-#> Negative log-likelihood at call 69: 385.2599
-#> Negative log-likelihood at call 70: 385.2597
-#> Negative log-likelihood at call 71: 374.7604
-#> Negative log-likelihood at call 72: 374.7603
-#> Negative log-likelihood at call 76: 373.199
-#> Negative log-likelihood at call 79: 373.199
-#> Negative log-likelihood at call 80: 373.199
-#> Negative log-likelihood at call 81: 372.3772
-#> Negative log-likelihood at call 84: 372.3772
-#> Negative log-likelihood at call 86: 371.2615
-#> Negative log-likelihood at call 89: 371.2615
-#> Negative log-likelihood at call 90: 371.2615
-#> Negative log-likelihood at call 92: 371.2439
-#> Negative log-likelihood at call 93: 371.2439
-#> Negative log-likelihood at call 94: 371.2439
-#> Negative log-likelihood at call 97: 371.2198
-#> Negative log-likelihood at call 98: 371.2198
-#> Negative log-likelihood at call 102: 371.2174
-#> Negative log-likelihood at call 104: 371.2174
-#> Negative log-likelihood at call 107: 371.2147
-#> Negative log-likelihood at call 110: 371.2147
-#> Negative log-likelihood at call 111: 371.2147
-#> Negative log-likelihood at call 112: 371.2145
-#> Negative log-likelihood at call 113: 371.2145
-#> Negative log-likelihood at call 116: 371.2145
-#> Negative log-likelihood at call 119: 371.2135
-#> Negative log-likelihood at call 121: 371.2135
-#> Negative log-likelihood at call 124: 371.2135
-#> Negative log-likelihood at call 126: 371.2135
-#> Negative log-likelihood at call 127: 371.2135
-#> Negative log-likelihood at call 133: 371.2134
-#> Negative log-likelihood at call 135: 371.2134
-#> Negative log-likelihood at call 138: 371.2134
-#> Negative log-likelihood at call 142: 371.2134
+                           solution_type = "deSolve")))
#> Warning: Observations with value of zero were removed from the data
#> Sum of squared residuals at call 1: 18915.53
+#> Sum of squared residuals at call 2: 18915.53
+#> Sum of squared residuals at call 6: 11424.02
+#> Sum of squared residuals at call 10: 11424
+#> Sum of squared residuals at call 13: 2367.052
+#> Sum of squared residuals at call 14: 2367.05
+#> Sum of squared residuals at call 19: 1314.716
+#> Sum of squared residuals at call 22: 1314.714
+#> Sum of squared residuals at call 25: 991.8311
+#> Sum of squared residuals at call 28: 991.8305
+#> Sum of squared residuals at call 30: 893.6462
+#> Sum of squared residuals at call 33: 893.6457
+#> Sum of squared residuals at call 35: 569.4049
+#> Sum of squared residuals at call 38: 569.4047
+#> Sum of squared residuals at call 40: 565.0651
+#> Sum of squared residuals at call 41: 565.065
+#> Sum of squared residuals at call 42: 565.0637
+#> Sum of squared residuals at call 45: 428.0188
+#> Sum of squared residuals at call 46: 428.0185
+#> Sum of squared residuals at call 50: 406.732
+#> Sum of squared residuals at call 52: 406.732
+#> Sum of squared residuals at call 55: 398.9115
+#> Sum of squared residuals at call 57: 398.9113
+#> Sum of squared residuals at call 60: 394.5943
+#> Sum of squared residuals at call 62: 394.5943
+#> Sum of squared residuals at call 66: 385.26
+#> Sum of squared residuals at call 67: 385.2599
+#> Sum of squared residuals at call 69: 385.2599
+#> Sum of squared residuals at call 70: 385.2597
+#> Sum of squared residuals at call 71: 374.7604
+#> Sum of squared residuals at call 72: 374.7603
+#> Sum of squared residuals at call 76: 373.199
+#> Sum of squared residuals at call 79: 373.199
+#> Sum of squared residuals at call 80: 373.199
+#> Sum of squared residuals at call 81: 372.3772
+#> Sum of squared residuals at call 84: 372.3772
+#> Sum of squared residuals at call 86: 371.2615
+#> Sum of squared residuals at call 89: 371.2615
+#> Sum of squared residuals at call 90: 371.2615
+#> Sum of squared residuals at call 92: 371.2439
+#> Sum of squared residuals at call 93: 371.2439
+#> Sum of squared residuals at call 94: 371.2439
+#> Sum of squared residuals at call 97: 371.2198
+#> Sum of squared residuals at call 98: 371.2198
+#> Sum of squared residuals at call 102: 371.2174
+#> Sum of squared residuals at call 104: 371.2174
+#> Sum of squared residuals at call 107: 371.2147
+#> Sum of squared residuals at call 110: 371.2147
+#> Sum of squared residuals at call 111: 371.2147
+#> Sum of squared residuals at call 112: 371.2145
+#> Sum of squared residuals at call 113: 371.2145
+#> Sum of squared residuals at call 116: 371.2145
+#> Sum of squared residuals at call 119: 371.2135
+#> Sum of squared residuals at call 121: 371.2135
+#> Sum of squared residuals at call 124: 371.2135
+#> Sum of squared residuals at call 126: 371.2135
+#> Sum of squared residuals at call 127: 371.2135
+#> Sum of squared residuals at call 133: 371.2134
+#> Sum of squared residuals at call 135: 371.2134
+#> Sum of squared residuals at call 138: 371.2134
+#> Sum of squared residuals at call 142: 371.2134
 #> Negative log-likelihood at call 152: 97.22429
 #> Optimisation successfully terminated.
 #>        User      System verstrichen 
-#>       1.136       0.000       1.135 
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff
+#>       1.152       0.000       1.153 
coef(fit.deSolve)
#> NULL
endpoints(fit.deSolve)
#> $ff
 #> parent_sink   parent_m1     m1_sink 
 #>    0.485524    0.514476    1.000000 
 #> 
@@ -543,9 +551,9 @@ Per default, parameters in the kinetic models are internally transformed in
 
# 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 
-#> R version used for fitting:       3.5.3 
-#> Date of fit:     Wed Apr 10 10:10:17 2019 
-#> Date of summary: Wed Apr 10 10:10:17 2019 
+#> R version used for fitting:       3.6.0 
+#> Date of fit:     Thu May  2 12:39:13 2019 
+#> Date of summary: Thu May  2 12:39:13 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -553,10 +561,10 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted with method using 404 model solutions performed in 1.105 s
+#> Fitted using 404 model solutions performed in 1.178 s
 #> 
 #> Error model:
-#> NULL
+#> Constant variance 
 #> 
 #> Starting values for parameters to be optimised:
 #>                     value   type
@@ -605,7 +613,7 @@ Per default, parameters in the kinetic models are internally transformed in
 #> f_parent_to_m1  0.514500  23.070 3.104e-22  0.469100 5.596e-01
 #> sigma           3.126000   8.718 2.235e-10  2.396000 3.855e+00
 #> 
-#> Chi2 error levels in percent:
+#> FOCUS Chi2 error levels in percent:
 #>          err.min n.optim df
 #> All data   6.398       4 15
 #> parent     6.459       2  7
@@ -661,9 +669,9 @@ Per default, parameters in the kinetic models are internally transformed in
 #>   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 
-#> R version used for fitting:       3.5.3 
-#> Date of fit:     Wed Apr 10 10:10:19 2019 
-#> Date of summary: Wed Apr 10 10:10:19 2019 
+#> R version used for fitting:       3.6.0 
+#> Date of fit:     Thu May  2 12:39:14 2019 
+#> Date of summary: Thu May  2 12:39:14 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -671,10 +679,10 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted with method using 558 model solutions performed in 1.602 s
+#> Fitted using 558 model solutions performed in 1.518 s
 #> 
 #> Error model:
-#> NULL
+#> Variance unique to each observed variable 
 #> 
 #> Starting values for parameters to be optimised:
 #>                   value   type
@@ -735,7 +743,7 @@ Per default, parameters in the kinetic models are internally transformed in
 #> sigma_parent    3.401000   5.985 5.662e-07  2.244000 4.559e+00
 #> sigma_m1        2.855000   6.311 2.215e-07  1.934000 3.777e+00
 #> 
-#> Chi2 error levels in percent:
+#> FOCUS Chi2 error levels in percent:
 #>          err.min n.optim df
 #> All data   6.398       4 15
 #> parent     6.464       2  7
@@ -791,9 +799,9 @@ Per default, parameters in the kinetic models are internally transformed in
 #>   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 
-#> R version used for fitting:       3.5.3 
-#> Date of fit:     Wed Apr 10 10:10:22 2019 
-#> Date of summary: Wed Apr 10 10:10:22 2019 
+#> R version used for fitting:       3.6.0 
+#> Date of fit:     Thu May  2 12:39:18 2019 
+#> Date of summary: Thu May  2 12:39:18 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent * parent
@@ -801,10 +809,10 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Model predictions using solution type deSolve 
 #> 
-#> Fitted with method using 756 model solutions performed in 3.222 s
+#> Fitted using 844 model solutions performed in 3.501 s
 #> 
 #> Error model:
-#> NULL
+#> Two-component variance function 
 #> 
 #> Starting values for parameters to be optimised:
 #>                   value   type
@@ -812,8 +820,8 @@ Per default, parameters in the kinetic models are internally transformed in
 #> k_parent         0.1000 deparm
 #> k_m1             0.1001 deparm
 #> f_parent_to_m1   0.5000 deparm
-#> sigma_low        0.5000  error
-#> rsd_high         0.0700  error
+#> sigma_low        3.0000  error
+#> rsd_high         0.0100  error
 #> 
 #> Starting values for the transformed parameters actually optimised:
 #>                     value lower upper
@@ -821,8 +829,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_low        0.500000     0   Inf
-#> rsd_high         0.070000     0   Inf
+#> sigma_low        3.000000     0   Inf
+#> rsd_high         0.010000     0   Inf
 #> 
 #> Fixed parameter values:
 #>      value  type
@@ -858,7 +866,7 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 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
 #> 
-#> Chi2 error levels in percent:
+#> FOCUS Chi2 error levels in percent:
 #>          err.min n.optim df
 #> All data   6.475       4 15
 #> parent     6.573       2  7
@@ -876,18 +884,18 @@ Per default, parameters in the kinetic models are internally transformed in
 #> 
 #> Data:
 #>  time variable observed predicted   residual
-#>     0   parent    99.46 100.73434  -1.274340
-#>     0   parent   102.04 100.73434   1.305660
+#>     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.491200
-#>     7   parent    55.13  49.82880   5.301200
+#>     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.692316
-#>    21   parent    11.64  12.19232  -0.552316
+#>    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
@@ -900,20 +908,20 @@ Per default, parameters in the kinetic models are internally transformed in
 #>     3       m1    12.96  13.22867  -0.268669
 #>     7       m1    22.97  25.36417  -2.394167
 #>     7       m1    24.47  25.36417  -0.894167
-#>    14       m1    41.69  37.00974   4.680262
-#>    14       m1    33.21  37.00974  -3.799738
+#>    14       m1    41.69  37.00974   4.680263
+#>    14       m1    33.21  37.00974  -3.799737
 #>    21       m1    44.37  41.90133   2.468668
 #>    21       m1    46.44  41.90133   4.538668
 #>    35       m1    41.22  43.45691  -2.236914
 #>    35       m1    37.95  43.45691  -5.506914
 #>    50       m1    41.19  41.34199  -0.151986
 #>    50       m1    40.01  41.34199  -1.331986
-#>    75       m1    40.09  36.61470   3.475295
-#>    75       m1    33.85  36.61470  -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.891303
-#>   120       m1    33.31  29.04130   4.268697

+#>    75       m1    40.09  36.61471   3.475294
+#>    75       m1    33.85  36.61471  -2.764706
+#>   100       m1    31.04  32.20082  -1.160824
+#>   100       m1    33.13  32.20082   0.929176
+#>   120       m1    25.15  29.04130  -3.891305
+#>   120       m1    33.31  29.04130   4.268695