jranke / mkin
1
utils::globalVariables("ds")
2

3
#' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
4
#'
5
#' @param x An object of class [mixed.mmkin], [nlme.mmkin]
6
#' @param i A numeric index to select datasets for which to plot the individual predictions,
7
#'   in case plots get too large
8
#' @inheritParams plot.mkinfit
9
#' @param standardized Should the residuals be standardized? Only takes effect if
10
#'   `resplot = "time"`.
11
#' @param pred_over Named list of alternative predictions as obtained
12
#'   from [mkinpredict] with a compatible [mkinmod].
13
#' @param rel.height.legend The relative height of the legend shown on top
14
#' @param rel.height.bottom The relative height of the bottom plot row
15
#' @param ymax Vector of maximum y axis values
16
#' @param ncol.legend Number of columns to use in the legend
17
#' @param nrow.legend Number of rows to use in the legend
18
#' @param resplot Should the residuals plotted against time or against
19
#'   predicted values?
20
#' @param col_ds Colors used for plotting the observed data and the
21
#'   corresponding model prediction lines for the different datasets.
22
#' @param pch_ds Symbols to be used for plotting the data.
23
#' @param lty_ds Line types to be used for the model predictions.
24
#' @importFrom stats coefficients
25
#' @return The function is called for its side effect.
26
#' @author Johannes Ranke
27
#' @examples
28
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
29
#'  function(x) x$data[c("name", "time", "value")])
30
#' names(ds) <- paste0("ds ", 6:10)
31
#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
32
#'   A1 = mkinsub("SFO"), quiet = TRUE)
33
#' \dontrun{
34
#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
35
#' plot(f[, 3:4], standardized = TRUE)
36
#'
37
#' # For this fit we need to increase pnlsMaxiter, and we increase the
38
#' # tolerance in order to speed up the fit for this example evaluation
39
#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
40
#' plot(f_nlme)
41
#'
42
#' }
43
#' @export
44
plot.mixed.mmkin <- function(x,
45
  i = 1:ncol(x$mmkin),
46
  obs_vars = names(x$mkinmod$map),
47
  standardized = TRUE,
48
  xlab = "Time",
49
  xlim = range(x$data$time),
50
  resplot = c("predicted", "time"),
51
  pred_over = NULL,
52
  ymax = "auto", maxabs = "auto",
53
  ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
54
  nrow.legend = ceiling((length(i) + 1) / ncol.legend),
55
  rel.height.legend = 0.02 + 0.07 * nrow.legend,
56
  rel.height.bottom = 1.1,
57
  pch_ds = 1:length(i),
58
  col_ds = pch_ds + 1,
59
  lty_ds = col_ds,
60
  frame = TRUE, ...
61
)
62
{
63
  # Prepare parameters and data
64 1
  fit_1 <- x$mmkin[[1]]
65 1
  ds_names <- colnames(x$mmkin)
66

67 1
  backtransform = TRUE
68

69 1
  if (identical(class(x), "mixed.mmkin")) {
70 1
    degparms_pop <- mean_degparms(x$mmkin)
71

72 1
    degparms_tmp <- parms(x$mmkin, transformed = TRUE)
73 1
    degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
74 1
    residual_type = ifelse(standardized, "standardized", "residual")
75 1
    residuals <- x$data[[residual_type]]
76
  }
77

78 1
  if (inherits(x, "nlme.mmkin")) {
79 0
    degparms_i <- coefficients(x)
80 0
    degparms_pop <- nlme::fixef(x)
81 0
    residuals <- residuals(x,
82 0
      type = ifelse(standardized, "pearson", "response"))
83
  }
84

85 1
  degparms_fixed <- fit_1$fixed$value
86 1
  names(degparms_fixed) <- rownames(fit_1$fixed)
87 1
  degparms_all <- cbind(as.matrix(degparms_i),
88 1
    matrix(rep(degparms_fixed, nrow(degparms_i)),
89 1
      ncol = length(degparms_fixed),
90 1
      nrow = nrow(degparms_i), byrow = TRUE))
91 1
  degparms_all_names <- c(names(degparms_i), names(degparms_fixed))
92 1
  colnames(degparms_all) <- degparms_all_names
93

94 1
  degparms_all_pop <- c(degparms_pop, degparms_fixed)
95

96 1
  odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
97 1
  odeparms_names <- setdiff(degparms_all_names, odeini_names)
98

99 1
  observed <- cbind(x$data[c("ds", "name", "time", "value")],
100 1
    residual = residuals)
101

102 1
  solution_type = fit_1$solution_type
103

104 1
  outtimes <- sort(unique(c(x$data$time,
105 1
    seq(xlim[1], xlim[2], length.out = 50))))
106

107 1
  pred_ds <- purrr::map_dfr(i, function(ds_i)   {
108 1
    odeparms_trans <- degparms_all[ds_i, odeparms_names]
109 1
    names(odeparms_trans) <- odeparms_names # needed if only one odeparm
110 1
    if (backtransform) {
111 1
      odeparms <- backtransform_odeparms(odeparms_trans,
112 1
        x$mkinmod,
113 1
        transform_rates = fit_1$transform_rates,
114 1
        transform_fractions = fit_1$transform_fractions)
115
    } else {
116 0
      odeparms <- odeparms_trans
117
    }
118

119 1
    odeini <- degparms_all[ds_i, odeini_names]
120 1
    names(odeini) <- gsub("_0", "", odeini_names)
121

122 1
    out <- mkinpredict(x$mkinmod, odeparms, odeini,
123 1
      outtimes, solution_type = solution_type,
124 1
      atol = fit_1$atol, rtol = fit_1$rtol)
125 1
    return(cbind(as.data.frame(out), ds = ds_names[ds_i]))
126
  })
127

128 1
  odeparms_pop_trans <- degparms_all_pop[odeparms_names]
129

130 1
  if (backtransform) {
131 1
    odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
132 1
      x$mkinmod,
133 1
      transform_rates = fit_1$transform_rates,
134 1
      transform_fractions = fit_1$transform_fractions)
135
  } else {
136 0
    odeparms_pop <- odeparms_pop_trans
137
  }
138

139 1
  odeini_pop <- degparms_all_pop[odeini_names]
140 1
  names(odeini_pop) <- gsub("_0", "", odeini_names)
141

142 1
  pred_pop <- as.data.frame(
143 1
    mkinpredict(x$mkinmod, odeparms_pop, odeini_pop,
144 1
      outtimes, solution_type = solution_type,
145 1
      atol = fit_1$atol, rtol = fit_1$rtol))
146

147
  # Start of graphical section
148 1
  oldpar <- par(no.readonly = TRUE)
149 1
  on.exit(par(oldpar, no.readonly = TRUE))
150

151 1
  n_plot_rows = length(obs_vars)
152 1
  n_plots = n_plot_rows * 2
153

154
  # Set relative plot heights, so the first plot row is the norm
155 1
  rel.heights <- if (n_plot_rows > 1) {
156 1
    c(rel.height.legend, c(rep(1, n_plot_rows - 1), rel.height.bottom))
157
  } else {
158 0
    c(rel.height.legend, 1)
159
  }
160

161 1
  layout_matrix = matrix(c(1, 1, 2:(n_plots + 1)),
162 1
    n_plot_rows + 1, 2, byrow = TRUE)
163 1
  layout(layout_matrix, heights = rel.heights)
164

165 1
  par(mar = c(0.1, 2.1, 0.1, 2.1))
166

167
  # Empty plot with legend
168 0
  if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
169 1
  else lty_over <- NULL
170 1
  n_pop <- 1 + length(lty_over)
171 1
  lty_pop <- c(1, lty_over)
172

173 1
  plot(0, type = "n", axes = FALSE, ann = FALSE)
174 1
  legend("center", bty = "n", ncol = ncol.legend,
175 1
    legend = c("Population", names(pred_over), ds_names[i]),
176 1
    lty = c(lty_pop, lty_ds),
177 1
    lwd = c(rep(2, n_pop), rep(1, length(i))),
178 1
    col = c(rep(1, n_pop), col_ds),
179 1
    pch = c(rep(NA, n_pop), pch_ds))
180

181 1
  resplot <- match.arg(resplot)
182

183
  # Loop plot rows
184 1
  for (plot_row in 1:n_plot_rows) {
185

186 1
    obs_var <- obs_vars[plot_row]
187 1
    observed_row <- subset(observed, name == obs_var)
188

189
    # Set ylim to sensible default, or use ymax
190 1
    if (identical(ymax, "auto")) {
191 1
      ylim_row = c(0,
192 1
        max(c(observed_row$value, pred_ds[[obs_var]]), na.rm = TRUE))
193
    } else {
194 0
      ylim_row = c(0, ymax[plot_row])
195
    }
196

197
    # Margins for bottom row of plots when we have more than one row
198
    # This is the only row that needs to show the x axis legend
199 1
    if (plot_row == n_plot_rows) {
200 1
      par(mar = c(5.1, 4.1, 1.1, 2.1))
201
    } else {
202 1
      par(mar = c(3.0, 4.1, 1.1, 2.1))
203
    }
204

205 1
    plot(pred_pop$time, pred_pop[[obs_var]],
206 1
      type = "l", lwd = 2, lty = lty_pop,
207 1
      xlim = xlim, ylim = ylim_row,
208 1
      xlab = xlab, ylab = paste("Residues", obs_var), frame = frame)
209

210 1
    if (!is.null(pred_over)) {
211 0
      for (i_over in seq_along(pred_over)) {
212 0
        pred_frame <- as.data.frame(pred_over[[i_over]])
213 0
        lines(pred_frame$time, pred_frame[[obs_var]],
214 0
          lwd = 2, lty = lty_over[i_over])
215
      }
216
    }
217

218 1
    for (ds_i in seq_along(i)) {
219 1
      points(subset(observed_row, ds == ds_names[ds_i], c("time", "value")),
220 1
        col = col_ds[ds_i], pch = pch_ds[ds_i])
221 1
      lines(subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)),
222 1
        col = col_ds[ds_i], lty = lty_ds[ds_i])
223
    }
224

225 1
    if (identical(maxabs, "auto")) {
226 1
      maxabs = max(abs(observed_row$residual), na.rm = TRUE)
227
    }
228

229 1
    if (identical(resplot, "time")) {
230 0
      plot(0, type = "n", xlim = xlim, xlab = "Time",
231 0
        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
232 0
        ylab = if (standardized) "Standardized residual" else "Residual")
233

234 0
      abline(h = 0, lty = 2)
235

236 0
      for (ds_i in seq_along(i)) {
237 0
        points(subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
238 0
          col = col_ds[ds_i], pch = pch_ds[ds_i])
239
      }
240
    }
241

242 1
    if (identical(resplot, "predicted")) {
243 1
      plot(0, type = "n",
244 1
        xlim = c(0, max(pred_ds[[obs_var]])),
245 1
        xlab = "Predicted",
246 1
        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
247 1
        ylab = if (standardized) "Standardized residual" else "Residual")
248

249 1
      abline(h = 0, lty = 2)
250

251 1
      for (ds_i in seq_along(i)) {
252 1
        observed_row_ds <- merge(
253 1
          subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
254 1
          subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)))
255 1
        points(observed_row_ds[c(3, 2)],
256 1
          col = col_ds[ds_i], pch = pch_ds[ds_i])
257
      }
258
    }
259
  }
260
}

Read our documentation on viewing source code .

Loading