Create an nlme model for an mmkin row object

+ + +

This functions sets up a nonlinear mixed effects model for an mmkin row +object. An mmkin row object is essentially a list of mkinfit objects that +have been obtained by fitting the same model to a list of datasets.

+ +

# S3 method for mmkin
+nlme(
+  model,
+  data = sys.frame(sys.parent()),
+  fixed,
+  random = fixed,
+  groups,
+  start,
+  correlation = NULL,
+  weights = NULL,
+  subset,
+  method = c("ML", "REML"),
+  na.action = na.fail,
+  naPattern,
+  control = list(),
+  verbose = FALSE
+)
+
+# S3 method for nlme.mmkin
+print(x, ...)
+
+# S3 method for nlme.mmkin
+update(object, ...)

+ +

Arguments

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

model	An `mmkin` row object.
data	Ignored, data are taken from the mmkin model
fixed	Ignored, all degradation parameters fitted in the +mmkin model are used as fixed parameters
random	If not specified, all fixed effects are complemented +with uncorrelated random effects
groups	See the documentation of nlme
start	If not specified, mean values of the fitted degradation +parameters taken from the mmkin object are used
correlation	See the documentation of nlme
weights	passed to nlme
subset	passed to nlme
method	passed to nlme
na.action	passed to nlme
naPattern	passed to nlme
control	passed to nlme
verbose	passed to nlme
x	An nlme.mmkin object to print
...	Update specifications passed to update.nlme
object	An nlme.mmkin object to update

+ +

Value

+ +

Upon success, a fitted nlme.mmkin object, which is an nlme object +with additional elements

Examples

ds <- lapply(experimental_data_for_UBA_2019[6:10],
+ function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
+f <- mmkin("SFO", ds, quiet = TRUE, cores = 1)
+library(nlme)
+endpoints(f[[1]])
#> $distimes
+#>            DT50     DT90
+#> parent 11.96183 39.73634
+#> 
f_nlme <- nlme(f)
+print(f_nlme)
#> Nonlinear mixed-effects model fit by maximum likelihood
+#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent_sink) 
+#>   Data: "Not shown" 
+#>   Log-likelihood: -307.5269
+#>   Fixed: list(parent_0 ~ 1, log_k_parent_sink ~ 1) 
+#>          parent_0 log_k_parent_sink 
+#>         85.540979         -3.229602 
+#> 
+#> Random effects:
+#>  Formula: list(parent_0 ~ 1, log_k_parent_sink ~ 1)
+#>  Level: ds
+#>  Structure: Diagonal
+#>         parent_0 log_k_parent_sink Residual
+#> StdDev: 1.308245          1.288586 6.304923
+#> 
+#> Number of Observations: 90
+#> Number of Groups: 5 
endpoints(f_nlme)
#> $distimes
+#>            DT50     DT90
+#> parent 17.51556 58.18543
+#> 
# \dontrun{
+  f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+  update(f_nlme_2, random = parent_0 ~ 1)
#> Nonlinear mixed-effects model fit by maximum likelihood
+#>   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent_sink) 
+#>   Data: "Not shown" 
+#>   Log-likelihood: -404.3729
+#>   Fixed: list(parent_0 ~ 1, log_k_parent_sink ~ 1) 
+#>          parent_0 log_k_parent_sink 
+#>         75.933480         -3.555983 
+#> 
+#> Random effects:
+#>  Formula: parent_0 ~ 1 | ds
+#>            parent_0 Residual
+#> StdDev: 0.002416802 21.63027
+#> 
+#> Number of Observations: 90
+#> Number of Groups: 5 
  # Test on some real data
+  ds_2 <- lapply(experimental_data_for_UBA_2019[6:10],
+   function(x) x$data[c("name", "time", "value")])
+  m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+    A1 = mkinsub("SFO"), use_of_ff = "min", quiet = TRUE)
+  m_sfo_sfo_ff <- mkinmod(parent = mkinsub("SFO", "A1"),
+    A1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE)
+  m_fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+    A1 = mkinsub("SFO"), quiet = TRUE)
+  m_dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+    A1 = mkinsub("SFO"), quiet = TRUE)
+
+  f_2 <- mmkin(list("SFO-SFO" = m_sfo_sfo,
+   "SFO-SFO-ff" = m_sfo_sfo_ff,
+   "FOMC-SFO" = m_fomc_sfo,
+   "DFOP-SFO" = m_dfop_sfo),
+    ds_2, quiet = TRUE)
+  plot(f_2["SFO-SFO", 3:4]) # Separate fits for datasets 3 and 4

+  f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ])
+  # plot(f_nlme_sfo_sfo) # not feasible with pkgdown figures
+  plot(f_nlme_sfo_sfo, 3:4) # Global mixed model: Fits for datasets 3 and 4

+  # With formation fractions
+  f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ])
+  plot(f_nlme_sfo_sfo_ff, 3:4) # chi2 different due to different df attribution

+  # For more parameters, we need to increase pnlsMaxIter and the tolerance
+  # to get convergence
+  f_nlme_fomc_sfo <- nlme(f_2["FOMC-SFO", ],
+    control = list(pnlsMaxIter = 100, tolerance = 1e-4), verbose = TRUE)
#> 
+#> **Iteration 1
+#> LME step: Loglik: -394.1603, nlminb iterations: 2
+#> reStruct  parameters:
+#>        ds1        ds2        ds3        ds4        ds5 
+#> -0.2079863  0.8563823  1.7454253  1.0917707  1.2756955 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  643.8814 
+#>  fixed effects: 94.17379  -5.473189  -0.6970234  -0.202509  2.103883  
+#>  iterations: 100 
+#> Convergence crit. (must all become <= tolerance = 0.0001):
+#>     fixed  reStruct 
+#> 0.7959873 0.1447512 
+#> 
+#> **Iteration 2
+#> LME step: Loglik: -396.3824, nlminb iterations: 7
+#> reStruct  parameters:
+#>           ds1           ds2           ds3           ds4           ds5 
+#> -1.712406e-01 -2.278541e-05  1.842120e+00  1.073975e+00  1.322924e+00 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  643.8025 
+#>  fixed effects: 94.17385  -5.473491  -0.6970406  -0.2025139  2.103871  
+#>  iterations: 100 
+#> Convergence crit. (must all become <= tolerance = 0.0001):
+#>       fixed    reStruct 
+#> 5.51758e-05 1.26861e-03 
+#> 
+#> **Iteration 3
+#> LME step: Loglik: -396.3825, nlminb iterations: 7
+#> reStruct  parameters:
+#>           ds1           ds2           ds3           ds4           ds5 
+#> -0.1712500923 -0.0001515734  1.8420972550  1.0739796967  1.3229177241 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  643.7941 
+#>  fixed effects: 94.17386  -5.473523  -0.6970424  -0.2025146  2.103869  
+#>  iterations: 100 
+#> Convergence crit. (must all become <= tolerance = 0.0001):
+#>        fixed     reStruct 
+#> 5.792621e-06 1.335434e-04 
+#> 
+#> **Iteration 4
+#> LME step: Loglik: -396.3825, nlminb iterations: 7
+#> reStruct  parameters:
+#>           ds1           ds2           ds3           ds4           ds5 
+#> -0.1712517206 -0.0001651603  1.8420950864  1.0739800294  1.3229173529 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  643.7949 
+#>  fixed effects: 94.17386  -5.473521  -0.6970423  -0.2025145  2.10387  
+#>  iterations: 100 
+#> Convergence crit. (must all become <= tolerance = 0.0001):
+#>        fixed     reStruct 
+#> 4.025781e-07 9.628656e-06 
  f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ],
+    control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE)
#> 
+#> **Iteration 1
+#> LME step: Loglik: -404.9583, nlminb iterations: 1
+#> reStruct  parameters:
+#>        ds1        ds2        ds3        ds4        ds5        ds6 
+#> -0.4114357  0.9798641  1.6990035  0.7293314  0.3354323  1.7113047 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  630.3642 
+#>  fixed effects: 93.82269  -5.455991  -0.6788957  -1.862196  -4.199671  0.0553284  
+#>  iterations: 120 
+#> Convergence crit. (must all become <= tolerance = 0.0005):
+#>     fixed  reStruct 
+#> 0.7879730 0.5822574 
+#> 
+#> **Iteration 2
+#> LME step: Loglik: -407.7755, nlminb iterations: 11
+#> reStruct  parameters:
+#>          ds1          ds2          ds3          ds4          ds5          ds6 
+#> -0.371224105  0.003056163  1.789939431  0.724671132  0.301602942  1.754200482 
+#>  Beginning PNLS step: ..  completed fit_nlme() step.
+#> PNLS step: RSS =  630.364 
+#>  fixed effects: 93.82269  -5.455991  -0.6788958  -1.862196  -4.199671  0.05532834  
+#>  iterations: 120 
+#> Convergence crit. (must all become <= tolerance = 0.0005):
+#>        fixed     reStruct 
+#> 9.814652e-07 1.059239e-05 
  plot(f_2["FOMC-SFO", 3:4])
  plot(f_nlme_fomc_sfo, 3:4)

+  plot(f_2["DFOP-SFO", 3:4])
  plot(f_nlme_dfop_sfo, 3:4)

+  anova(f_nlme_dfop_sfo, f_nlme_fomc_sfo, f_nlme_sfo_sfo)
#>                 Model df       AIC       BIC    logLik   Test   L.Ratio p-value
+#> f_nlme_dfop_sfo     1 13  843.8547  884.6201 -408.9274                         
+#> f_nlme_fomc_sfo     2 11  818.5151  853.0089 -398.2576 1 vs 2  21.33957  <.0001
+#> f_nlme_sfo_sfo      3  9 1085.1821 1113.4043 -533.5910 2 vs 3 270.66697  <.0001
  anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) # if we ignore FOMC
#>                 Model df       AIC       BIC    logLik   Test  L.Ratio p-value
+#> f_nlme_dfop_sfo     1 13  843.8547  884.6201 -408.9274                        
+#> f_nlme_sfo_sfo      2  9 1085.1821 1113.4043 -533.5910 1 vs 2 249.3274  <.0001

+  endpoints(f_nlme_sfo_sfo)
#> $ff
+#> parent_sink   parent_A1     A1_sink 
+#>   0.5912432   0.4087568   1.0000000 
+#> 
+#> $distimes
+#>            DT50     DT90
+#> parent 19.13518  63.5657
+#> A1     66.02155 219.3189
+#> 
  endpoints(f_nlme_dfop_sfo)
#> $ff
+#>   parent_A1 parent_sink 
+#>   0.2768574   0.7231426 
+#> 
+#> $distimes
+#>             DT50     DT90  DT50_k1  DT50_k2
+#> parent  11.07091 104.6320 4.462384 46.20825
+#> A1     162.30518 539.1661       NA       NA
+#> 
# }
+

+ +

Create an nlme model for an mmkin row object

Arguments

Value

See also

Examples