From c73b2f30ec836c949885784ab576e814eb8070a9 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 9 Mar 2021 17:35:47 +0100 Subject: Some improvements for borderline cases - fit_with_errors for saem() - test_log_parms for mean_degparms() and saem() --- R/nlme.R | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) (limited to 'R/nlme.R') diff --git a/R/nlme.R b/R/nlme.R index 9215aab0..d235a094 100644 --- a/R/nlme.R +++ b/R/nlme.R @@ -36,7 +36,7 @@ #' nlme_f <- nlme_function(f) #' # These assignments are necessary for these objects to be #' # visible to nlme and augPred when evaluation is done by -#' # pkgdown to generated the html docs. +#' # pkgdown to generate the html docs. #' assign("nlme_f", nlme_f, globalenv()) #' assign("grouped_data", grouped_data, globalenv()) #' @@ -130,13 +130,44 @@ nlme_function <- function(object) { #' fixed and random effects, in the format required by the start argument of #' nlme for the case of a single grouping variable ds. #' @param random Should a list with fixed and random effects be returned? +#' @param test_log_parms If TRUE, log parameters are only considered in +#' the mean calculations if their untransformed counterparts (most likely +#' rate constants) pass the t-test for significant difference from zero. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. #' @export -mean_degparms <- function(object, random = FALSE) { +mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) +{ if (nrow(object) > 1) stop("Only row objects allowed") parm_mat_trans <- sapply(object, parms, transformed = TRUE) + + if (test_log_parms) { + parm_mat_dim <- dim(parm_mat_trans) + parm_mat_dimnames <- dimnames(parm_mat_trans) + + log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE) + log_parm_names <- gsub("^log_", "", log_parm_trans_names) + + t_test_back_OK <- matrix( + sapply(object, function(o) { + suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level)) + }), nrow = length(log_parm_names)) + rownames(t_test_back_OK) <- log_parm_trans_names + + parm_mat_trans_OK <- parm_mat_trans + for (trans_parm in log_parm_trans_names) { + parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ], + parm_mat_trans[trans_parm, ], NA) + } + } else { + parm_mat_trans_OK <- parm_mat_trans + } + mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] - fixed <- apply(degparm_mat_trans, 1, mean) + degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE] + + fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE) if (random) { random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) # If we only have one parameter, apply returns a vector so we get a single row -- cgit v1.2.1