1
# Code inspired by nlme::nlme.nlsList and R/nlme_fit.R from nlmixr
2

3
# We need to assign the degradation function created in nlme.mmkin to an
4
# environment that is always accessible, also e.g. when evaluation is done by
5
# testthat or pkgdown. Therefore parent.frame() is not good enough. The
6
# following environment will be in the mkin namespace.
7
.nlme_env <- new.env(parent = emptyenv())
8

9
#' @export
10
nlme::nlme
11

12
#' Retrieve a degradation function from the mmkin namespace
13
#'
14
#' @importFrom utils getFromNamespace
15
#' @return A function that was likely previously assigned from within
16
#'   nlme.mmkin
17
#' @export
18
get_deg_func <- function() {
19 2
  return(get("deg_func", getFromNamespace(".nlme_env", "mkin")))
20
}
21

22
#' Create an nlme model for an mmkin row object
23
#'
24
#' This functions sets up a nonlinear mixed effects model for an mmkin row
25
#' object. An mmkin row object is essentially a list of mkinfit objects that
26
#' have been obtained by fitting the same model to a list of datasets.
27
#'
28
#' @param model An [mmkin] row object.
29
#' @param data Ignored, data are taken from the mmkin model
30
#' @param fixed Ignored, all degradation parameters fitted in the
31
#'   mmkin model are used as fixed parameters
32
#' @param random If not specified, all fixed effects are complemented
33
#'   with uncorrelated random effects
34
#' @param groups See the documentation of nlme
35
#' @param start If not specified, mean values of the fitted degradation
36
#'   parameters taken from the mmkin object are used
37
#' @param correlation See the documentation of nlme
38
#' @param weights passed to nlme
39
#' @param subset passed to nlme
40
#' @param method passed to nlme
41
#' @param na.action passed to nlme
42
#' @param naPattern passed to nlme
43
#' @param control passed to nlme
44
#' @param verbose passed to nlme
45
#' @importFrom stats na.fail as.formula
46
#' @return Upon success, a fitted 'nlme.mmkin' object, which is an nlme object
47
#'   with additional elements. It also inherits from 'mixed.mmkin'.
48
#' @note As the object inherits from [nlme::nlme], there is a wealth of
49
#'   methods that will automatically work on 'nlme.mmkin' objects, such as
50
#'   [nlme::intervals()], [nlme::anova.lme()] and [nlme::coef.lme()].
51
#' @export
52
#' @seealso [nlme_function()], [plot.mixed.mmkin], [summary.nlme.mmkin]
53
#' @examples
54
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
55
#'  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
56
#' f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1)
57
#' library(nlme)
58
#' f_nlme_sfo <- nlme(f["SFO", ])
59
#' f_nlme_dfop <- nlme(f["DFOP", ])
60
#' AIC(f_nlme_sfo, f_nlme_dfop)
61
#' print(f_nlme_dfop)
62
#' plot(f_nlme_dfop)
63
#' endpoints(f_nlme_dfop)
64
#'
65
#' \dontrun{
66
#'   f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1))
67
#'   update(f_nlme_2, random = parent_0 ~ 1)
68
#'   ds_2 <- lapply(experimental_data_for_UBA_2019[6:10],
69
#'    function(x) x$data[c("name", "time", "value")])
70
#'   m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
71
#'     A1 = mkinsub("SFO"), use_of_ff = "min", quiet = TRUE)
72
#'   m_sfo_sfo_ff <- mkinmod(parent = mkinsub("SFO", "A1"),
73
#'     A1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE)
74
#'   m_dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
75
#'     A1 = mkinsub("SFO"), quiet = TRUE)
76
#'
77
#'   f_2 <- mmkin(list("SFO-SFO" = m_sfo_sfo,
78
#'    "SFO-SFO-ff" = m_sfo_sfo_ff,
79
#'    "DFOP-SFO" = m_dfop_sfo),
80
#'     ds_2, quiet = TRUE)
81
#'
82
#'   f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ])
83
#'   plot(f_nlme_sfo_sfo)
84
#'
85
#'   # With formation fractions
86
#'   f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ])
87
#'   plot(f_nlme_sfo_sfo_ff)
88
#'
89
#'   # For the following fit we need to increase pnlsMaxIter and the tolerance
90
#'   # to get convergence
91
#'   f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ],
92
#'     control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE)
93
#'
94
#'   plot(f_nlme_dfop_sfo)
95
#'
96
#'   anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo)
97
#'
98
#'   endpoints(f_nlme_sfo_sfo)
99
#'   endpoints(f_nlme_dfop_sfo)
100
#'
101
#'   if (length(findFunction("varConstProp")) > 0) { # tc error model for nlme available
102
#'     # Attempts to fit metabolite kinetics with the tc error model are possible,
103
#'     # but need tweeking of control values and sometimes do not converge
104
#'
105
#'     f_tc <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, error_model = "tc")
106
#'     f_nlme_sfo_tc <- nlme(f_tc["SFO", ])
107
#'     f_nlme_dfop_tc <- nlme(f_tc["DFOP", ])
108
#'     AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc)
109
#'     print(f_nlme_dfop_tc)
110
#'   }
111
#'
112
#'   f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo,
113
#'    "DFOP-SFO" = m_dfop_sfo),
114
#'     ds_2, quiet = TRUE, error_model = "obs")
115
#'   f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ])
116
#'   print(f_nlme_sfo_sfo_obs)
117
#'   # The same with DFOP-SFO does not converge, apparently the variances of
118
#'   # parent and A1 are too similar in this case, so that the model is
119
#'   # overparameterised
120
#'   #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100))
121
#' }
122
nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
123
  fixed, random = fixed,
124
  groups, start, correlation = NULL, weights = NULL,
125
  subset, method = c("ML", "REML"),
126
  na.action = na.fail, naPattern,
127
  control = list(), verbose= FALSE)
128
{
129 0
  if (nrow(model) > 1) stop("Only row objects allowed")
130

131 2
  thisCall <- as.list(match.call())[-1]
132

133
  # Warn in case arguments were used that are overriden
134 2
  if (any(!is.na(match(names(thisCall),
135 2
               c("fixed", "data"))))) {
136 0
    warning("'nlme.mmkin' will redefine 'fixed' and 'data'")
137
  }
138

139 2
  deg_func <- nlme_function(model)
140

141 2
  assign("deg_func", deg_func, getFromNamespace(".nlme_env", "mkin"))
142

143
  # For the formula, get the degradation function from the mkin namespace
144 2
  this_model_text <- paste0("value ~ mkin::get_deg_func()(",
145 2
    paste(names(formals(deg_func)), collapse = ", "), ")")
146 2
  this_model <- as.formula(this_model_text)
147

148 2
  thisCall[["model"]] <- this_model
149

150 2
  mean_dp_start <- mean_degparms(model)
151 2
  dp_names <- names(mean_dp_start)
152

153 2
  thisCall[["data"]] <- nlme_data(model)
154

155 2
  if (missing(start)) {
156 2
    thisCall[["start"]] <- mean_degparms(model, random = TRUE)
157
  }
158

159 2
  thisCall[["fixed"]] <- lapply(as.list(dp_names), function(el)
160 2
                                 eval(parse(text = paste(el, 1, sep = "~"))))
161

162 2
  if (missing(random)) {
163 2
    thisCall[["random"]] <- pdDiag(thisCall[["fixed"]])
164
  }
165

166 2
  error_model <- model[[1]]$err_mod
167

168 2
  if (missing(weights)) {
169 2
    thisCall[["weights"]] <- switch(error_model,
170 2
      const = NULL,
171 2
      obs = varIdent(form = ~ 1 | name),
172 2
      tc = varConstProp())
173 2
    sigma <- switch(error_model,
174 2
      tc = 1,
175 2
      NULL)
176
  }
177

178 2
  control <- thisCall[["control"]]
179 2
  if (error_model == "tc") {
180 0
    control$sigma = 1
181 0
    thisCall[["control"]] <- control
182
  }
183

184 2
  fit_time <- system.time(val <- do.call("nlme.formula", thisCall))
185 2
  val$time <- fit_time
186

187 2
  val$mkinmod <- model[[1]]$mkinmod
188 2
  val$data <- thisCall[["data"]]
189 2
  val$mmkin <- model
190 2
  val$mean_dp_start <- mean_dp_start
191 2
  val$transform_rates <- model[[1]]$transform_rates
192 2
  val$transform_fractions <- model[[1]]$transform_fractions
193 2
  val$solution_type <- model[[1]]$solution_type
194 2
  val$err_mode <- error_model
195

196 2
  val$bparms.optim <- backtransform_odeparms(val$coefficients$fixed,
197 2
    val$mkinmod,
198 2
    transform_rates = val$transform_rates,
199 2
    transform_fractions = val$transform_fractions)
200

201 2
  val$bparms.fixed <- model[[1]]$bparms.fixed
202 2
  val$date.fit <- date()
203 2
  val$nlmeversion <- as.character(utils::packageVersion("nlme"))
204 2
  val$mkinversion <- as.character(utils::packageVersion("mkin"))
205 2
  val$Rversion <- paste(R.version$major, R.version$minor, sep=".")
206 2
  class(val) <- c("nlme.mmkin", "mixed.mmkin", "nlme", "lme")
207 2
  return(val)
208
}
209

210
#' @export
211
#' @rdname nlme.mmkin
212
#' @param x An nlme.mmkin object to print
213
#' @param digits Number of digits to use for printing
214
print.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
215 0
  cat( "Kinetic nonlinear mixed-effects model fit by " )
216 0
  cat( if(x$method == "REML") "REML\n" else "maximum likelihood\n")
217 0
  cat("\nStructural model:\n")
218 0
  diffs <- x$mmkin[[1]]$mkinmod$diffs
219 0
  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
220 0
  writeLines(strwrap(nice_diffs, exdent = 11))
221 0
  cat("\nData:\n")
222 0
  cat(nrow(x$data), "observations of",
223 0
    length(unique(x$data$name)), "variable(s) grouped in",
224 0
    length(unique(x$data$ds)), "datasets\n")
225 0
  cat("\nLog-", if(x$method == "REML") "restricted-" else "",
226 0
      "likelihood: ", format(x$logLik), "\n", sep = "")
227 0
  fixF <- x$call$fixed
228 0
  cat("\nFixed effects:\n",
229 0
      deparse(
230 0
  if(inherits(fixF, "formula") || is.call(fixF) || is.name(fixF))
231 0
    x$call$fixed
232
  else
233 0
    lapply(fixF, function(el) as.name(deparse(el)))), "\n")
234 0
  print(fixef(x), digits = digits, ...)
235 0
  cat("\n")
236 0
  print(summary(x$modelStruct), sigma = x$sigma, digits = digits, ...)
237 0
  invisible(x)
238
}
239

240
#' @export
241
#' @rdname nlme.mmkin
242
#' @param object An nlme.mmkin object to update
243
#' @param ... Update specifications passed to update.nlme
244
update.nlme.mmkin <- function(object, ...) {
245 2
  res <- NextMethod()
246 2
  res$mmkin <- object$mmkin
247 2
  class(res) <- c("nlme.mmkin", "nlme", "lme")
248 2
  return(res)
249
}

Read our documentation on viewing source code .

Loading