aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2022-10-19 22:59:34 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2022-10-20 11:46:35 +0200
commitcdf3848b252e8aae311bbfbf1d2921a38d904ad1 (patch)
tree19aa78676245621efe3f2db20c5e4d5eaccd15e0 /R
parent273efa0704130ff594d8b780e67827acc537b25a (diff)
First working version setting up covariate models
Diffstat (limited to 'R')
-rw-r--r--R/saem.R57
1 files changed, 53 insertions, 4 deletions
diff --git a/R/saem.R b/R/saem.R
index dbc19dec..7c8e6444 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -128,6 +128,8 @@ saem.mmkin <- function(object,
conf.level = 0.6,
solution_type = "auto",
covariance.model = "auto",
+ covariates = NULL,
+ covariate_models = NULL,
no_random_effect = NULL,
nbiter.saemix = c(300, 100),
control = list(displayProgress = FALSE, print = FALSE,
@@ -144,9 +146,11 @@ saem.mmkin <- function(object,
solution_type = solution_type,
transformations = transformations,
covariance.model = covariance.model,
+ covariates = covariates,
+ covariate_models = covariate_models,
no_random_effect = no_random_effect,
...)
- d_saemix <- saemix_data(object, verbose = verbose)
+ d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose)
fit_failed <- FALSE
FIM_failed <- NULL
@@ -268,6 +272,7 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
saemix_model <- function(object, solution_type = "auto",
transformations = c("mkin", "saemix"), degparms_start = numeric(),
covariance.model = "auto", no_random_effect = NULL,
+ covariates = NULL, covariate_models = NULL,
test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
{
if (nrow(object) > 1) stop("Only row objects allowed")
@@ -612,8 +617,8 @@ saemix_model <- function(object, solution_type = "auto",
degparms_psi0 <- degparms_optim
degparms_psi0[names(degparms_start)] <- degparms_start
- psi0_matrix <- matrix(degparms_psi0, nrow = 1)
- colnames(psi0_matrix) <- names(degparms_psi0)
+ psi0_matrix <- matrix(degparms_psi0, nrow = 1,
+ dimnames = list("(Intercept)", names(degparms_psi0)))
if (covariance.model[1] == "auto") {
covariance_diagonal <- rep(1, length(degparms_optim))
@@ -624,6 +629,41 @@ saemix_model <- function(object, solution_type = "auto",
covariance.model = diag(covariance_diagonal)
}
+ if (is.null(covariate_models)) {
+ covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel()
+ } else {
+ degparms_dependent <- sapply(covariate_models, function(x) as.character(x[[2]]))
+ covariates_in_models = unique(unlist(lapply(
+ covariate_models, function(x)
+ colnames(attr(terms(x), "factors"))
+ )))
+ covariates_not_available <- setdiff(covariates_in_models, names(covariates))
+ if (length(covariates_not_available) > 0) {
+ stop("Covariate(s) ", paste(covariates_not_available, collapse = ", "),
+ " used in the covariate models not available in the covariate data")
+ }
+ psi0_matrix <- rbind(psi0_matrix,
+ matrix(0, nrow = length(covariates), ncol = ncol(psi0_matrix),
+ dimnames = list(names(covariates), colnames(psi0_matrix))))
+ covariate.model <- matrix(0, nrow = length(covariates),
+ ncol = ncol(psi0_matrix),
+ dimnames = list(
+ covariates = names(covariates),
+ degparms = colnames(psi0_matrix)))
+ if (transformations == "saemix") {
+ stop("Covariate models with saemix transformations currently not supported")
+ }
+ parms_trans <- as.data.frame(t(sapply(object, parms, transformed = TRUE)))
+ for (covariate_model in covariate_models) {
+ covariate_name <- as.character(covariate_model[[2]])
+ model_data <- cbind(parms_trans, covariates)
+ ini_model <- lm(covariate_model, data = model_data)
+ ini_coef <- coef(ini_model)
+ psi0_matrix[names(ini_coef), covariate_name] <- ini_coef
+ covariate.model[names(ini_coef)[-1], covariate_name] <- as.numeric(as.logical(ini_coef[-1]))
+ }
+ }
+
res <- saemix::saemixModel(model_function,
psi0 = psi0_matrix,
"Mixed model generated from mmkin object",
@@ -631,6 +671,7 @@ saemix_model <- function(object, solution_type = "auto",
error.model = error.model,
verbose = verbose,
covariance.model = covariance.model,
+ covariate.model = covariate.model,
...
)
attr(res, "mean_dp_start") <- degparms_optim
@@ -641,7 +682,7 @@ saemix_model <- function(object, solution_type = "auto",
#' @importFrom rlang !!!
#' @return An [saemix::SaemixData] object.
#' @export
-saemix_data <- function(object, verbose = FALSE, ...) {
+saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) {
if (nrow(object) > 1) stop("Only row objects allowed")
ds_names <- colnames(object)
@@ -653,11 +694,19 @@ saemix_data <- function(object, verbose = FALSE, ...) {
time = ds_saemix_all$time,
value = ds_saemix_all$observed,
stringsAsFactors = FALSE)
+ if (!is.null(covariates)) {
+ name.covariates <- names(covariates)
+ covariates$ds <- rownames(covariates)
+ ds_saemix <- merge(ds_saemix, covariates, sort = FALSE)
+ } else {
+ name.covariates <- character(0)
+ }
res <- saemix::saemixData(ds_saemix,
name.group = "ds",
name.predictors = c("time", "name"),
name.response = "value",
+ name.covariates = name.covariates,
verbose = verbose,
...)
return(res)

Contact - Imprint