aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2023-03-20 22:40:20 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2023-03-20 22:40:20 +0100
commit28286eeaca84c85d2f4c3cddc9524d56d23b9aa0 (patch)
tree85f8a3156875ad320c313079d560439d28ee3ff2
parentdd525f49b852376f24851f23c36d6c50f23dbf82 (diff)
Support covariates in parms and plot.saem.mmkin
-rw-r--r--R/endpoints.R2
-rw-r--r--R/plot.mixed.mmkin.R129
-rw-r--r--R/saem.R35
3 files changed, 113 insertions, 53 deletions
diff --git a/R/endpoints.R b/R/endpoints.R
index 0c26ffae..9ea0e598 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -38,7 +38,7 @@
#' }
#'
#' @export
-endpoints <- function(fit, ...) {
+endpoints <- function(fit, ..., covariates = mean) {
mkinmod <- fit$mkinmod
obs_vars <- names(mkinmod$spec)
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 6f613b48..68f8d577 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -8,9 +8,17 @@ utils::globalVariables("ds")
#' @inheritParams plot.mkinfit
#' @param standardized Should the residuals be standardized? Only takes effect if
#' `resplot = "time"`.
-#' @param pop_curve Per default, a population curve is drawn in case
+#' @param pop_curves Per default, one population curve is drawn in case
#' population parameters are fitted by the model, e.g. for saem objects.
-#' In case there is a covariate model, no population curve is currently shown.
+#' In case there is a covariate model, the behaviour depends on the value
+#' of 'covariates'
+#' @param covariates Data frame with covariate values for all variables in
+#' any covariate models in the object. If given, it overrides 'covariate_quantiles'
+#' @param covariate_quantiles This argument only has an effect if the fitted
+#' object has covariate models. If so, the default is to show three population
+#' curves, for the 5th percentile, the 50th percentile and the 95th percentile
+#' of the covariate values used for fitting the model.
+#' @note Covariate models are currently only supported for saem.mmkin objects.
#' @param pred_over Named list of alternative predictions as obtained
#' from [mkinpredict] with a compatible [mkinmod].
#' @param test_log_parms Passed to [mean_degparms] in the case of an
@@ -68,10 +76,12 @@ plot.mixed.mmkin <- function(x,
i = 1:ncol(x$mmkin),
obs_vars = names(x$mkinmod$map),
standardized = TRUE,
+ covariates = NULL,
+ covariate_quantiles = c(0.05, 0.5, 0.95),
xlab = "Time",
xlim = range(x$data$time),
resplot = c("predicted", "time"),
- pop_curve = "auto",
+ pop_curves = "auto",
pred_over = NULL,
test_log_parms = FALSE,
conf.level = 0.6,
@@ -94,12 +104,12 @@ plot.mixed.mmkin <- function(x,
backtransform = TRUE
if (identical(class(x), "mixed.mmkin")) {
- if (identical(pop_curve, "auto")) {
- pop_curve <- FALSE
+ if (identical(pop_curves, "auto")) {
+ pop_curves <- FALSE
} else {
- pop_curve <- TRUE
+ pop_curves <- TRUE
}
- if (pop_curve) {
+ if (pop_curves) {
degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms,
conf.level = conf.level, default_log_parms = default_log_parms)
}
@@ -111,10 +121,10 @@ plot.mixed.mmkin <- function(x,
}
if (inherits(x, "nlme.mmkin")) {
- if (identical(pop_curve, "auto")) {
- pop_curve <- TRUE
+ if (identical(pop_curves, "auto")) {
+ pop_curves <- TRUE
} else {
- pop_curve <- FALSE
+ pop_curves <- FALSE
}
degparms_i <- coefficients(x)
degparms_pop <- nlme::fixef(x)
@@ -131,19 +141,30 @@ plot.mixed.mmkin <- function(x,
residual_type = ifelse(standardized, "standardized", "residual")
residuals <- x$data[[residual_type]]
- if (identical(pop_curve, "auto")) {
- if (nrow(x$sm@covariate.model) == 0) {
- pop_curve <- TRUE
+ if (identical(pop_curves, "auto")) {
+ if (length(x$covariate_models) == 0) {
+ degparms_pop <- x$so@results@fixed.effects
+ names(degparms_pop) <- degparms_i_names
+ pop_curves <- TRUE
} else {
- pop_curve <- FALSE
+ if (is.null(covariates)) {
+ covariates = as.data.frame(
+ apply(x$covariates, 2, quantile,
+ covariate_quantiles, simplify = FALSE))
+ }
+ degparms_pop <- parms(x, covariates = covariates)
+ pop_curves <- TRUE
}
- }
- if (pop_curve) {
- degparms_pop <- x$so@results@fixed.effects
- names(degparms_pop) <- degparms_i_names
+ } else {
+ pop_curves <- FALSE
}
}
+ # Make sure degparms_pop is a matrix, population curve variants in columns, if any
+ if (is.null(dim(degparms_pop))) {
+ degparms_pop <- matrix(degparms_pop)
+ }
+
degparms_fixed <- fit_1$fixed$value
names(degparms_fixed) <- rownames(fit_1$fixed)
degparms_all <- cbind(as.matrix(degparms_i),
@@ -186,28 +207,31 @@ plot.mixed.mmkin <- function(x,
names(pred_list) <- ds_names[i]
pred_ds <- vctrs::vec_rbind(!!!pred_list, .names_to = "ds")
- if (pop_curve) {
- degparms_all_pop <- c(degparms_pop, degparms_fixed)
- odeparms_pop_trans <- degparms_all_pop[odeparms_names]
-
- if (backtransform) {
- odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
- x$mkinmod,
- transform_rates = fit_1$transform_rates,
- transform_fractions = fit_1$transform_fractions)
- } else {
- odeparms_pop <- odeparms_pop_trans
- }
+ if (pop_curves) {
+ pred_list_pop <- lapply(1:ncol(degparms_pop), function(cov_i) {
+ degparms_all_pop_i <- c(degparms_pop[, cov_i], degparms_fixed)
+ odeparms_pop_trans_i <- degparms_all_pop_i[odeparms_names]
+ names(odeparms_pop_trans_i) <- odeparms_names # needed if only one odeparm
+ if (backtransform) {
+ odeparms_pop_i <- backtransform_odeparms(odeparms_pop_trans_i,
+ x$mkinmod,
+ transform_rates = fit_1$transform_rates,
+ transform_fractions = fit_1$transform_fractions)
+ } else {
+ odeparms_pop_i <- odeparms_pop_trans_i
+ }
- odeini_pop <- degparms_all_pop[odeini_names]
- names(odeini_pop) <- gsub("_0", "", odeini_names)
+ odeini <- degparms_all_pop_i[odeini_names]
+ names(odeini) <- gsub("_0", "", odeini_names)
- pred_pop <- as.data.frame(
- mkinpredict(x$mkinmod, odeparms_pop, odeini_pop,
+ out <- mkinpredict(x$mkinmod, odeparms_pop_i, odeini,
outtimes, solution_type = solution_type,
- atol = fit_1$atol, rtol = fit_1$rtol))
+ atol = fit_1$atol, rtol = fit_1$rtol)
+ })
+ names(pred_list_pop) <- colnames(degparms_pop)
+
} else {
- pred_pop <- NULL
+ pred_list_pop <- NULL
}
# Start of graphical section
@@ -233,22 +257,26 @@ plot.mixed.mmkin <- function(x,
# Empty plot with legend
if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
else lty_over <- NULL
- if (pop_curve) {
- lty_pop <- 1
- names(lty_pop) <- "Population"
+ if (pop_curves) {
+ if (is.null(covariates)) {
+ lty_pop <- 1
+ names(lty_pop) <- "Population"
+ } else {
+ lty_pop <- 1:nrow(covariates)
+ names(lty_pop) <- rownames(covariates)
+ }
} else {
lty_pop <- NULL
}
- n_pop <- length(lty_pop) + length(lty_over)
- lty_pop <- c(lty_pop, lty_over)
+ n_pop_over <- length(lty_pop) + length(lty_over)
plot(0, type = "n", axes = FALSE, ann = FALSE)
legend("center", bty = "n", ncol = ncol.legend,
legend = c(names(lty_pop), names(pred_over), ds_names[i]),
- lty = c(lty_pop, lty_ds),
- lwd = c(rep(2, n_pop), rep(1, length(i))),
- col = c(rep(1, n_pop), col_ds),
- pch = c(rep(NA, n_pop), pch_ds))
+ lty = c(lty_pop, lty_over, lty_ds),
+ lwd = c(rep(2, n_pop_over), rep(1, length(i))),
+ col = c(rep(1, n_pop_over), col_ds),
+ pch = c(rep(NA, n_pop_over), pch_ds))
resplot <- match.arg(resplot)
@@ -293,9 +321,14 @@ plot.mixed.mmkin <- function(x,
col = col_ds[ds_i], lty = lty_ds[ds_i])
}
- if (pop_curve) {
- lines(pred_pop$time, pred_pop[[obs_var]],
- type = "l", lwd = 2, lty = lty_pop)
+ if (pop_curves) {
+ for (cov_i in seq_along(pred_list_pop)) {
+ cov_name <- names(pred_list_pop)[cov_i]
+ lines(
+ pred_list_pop[[cov_i]][, "time"],
+ pred_list_pop[[cov_i]][, obs_var],
+ type = "l", lwd = 2, lty = lty_pop[cov_i])
+ }
}
if (identical(maxabs, "auto")) {
diff --git a/R/saem.R b/R/saem.R
index bddd3bfe..865f174e 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -804,8 +804,12 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
#' @export
#' @rdname saem
#' @param ci Should a matrix with estimates and confidence interval boundaries
-#' be returned? If FALSE (default), a vector of estimates is returned.
-parms.saem.mmkin <- function(object, ci = FALSE, ...) {
+#' be returned? If FALSE (default), a vector of estimates is returned if no
+#' covariates are given, otherwise a matrix of estimates is returned, with
+#' each column corresponding to a row of the data frame holding the covariates
+#' @param covariates A data frame holding covariate values for which to
+#' return parameter values. Only has an effect if 'ci' is FALSE.
+parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, covariate_quantiles = ...) {
cov.mod <- object$sm@covariance.model
n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)])
n_parms <- length(object$sm@name.modpar) +
@@ -827,6 +831,29 @@ parms.saem.mmkin <- function(object, ci = FALSE, ...) {
names(estimate) <- rownames(conf.int)
- if (ci) return(conf.int)
- else return(estimate)
+ if (ci) {
+ return(conf.int)
+ } else {
+ if (is.null(covariates)) {
+ return(estimate)
+ } else {
+ est_for_cov <- matrix(NA,
+ nrow = length(object$sm@name.modpar), ncol = nrow(covariates),
+ dimnames = (list(object$sm@name.modpar, rownames(covariates))))
+ covmods <- object$covariate_models
+ names(covmods) <- sapply(covmods, function(x) as.character(x[[2]]))
+ for (deg_parm_name in rownames(est_for_cov)) {
+ if (deg_parm_name %in% names(covmods)) {
+ covariate <- covmods[[deg_parm_name]][[3]]
+ beta_degparm_name <- paste0("beta_", covariate,
+ "(", deg_parm_name, ")")
+ est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name] +
+ covariates[[covariate]] * estimate[beta_degparm_name]
+ } else {
+ est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name]
+ }
+ }
+ return(est_for_cov)
+ }
+ }
}

Contact - Imprint