From aaa4cab7e0c7212f91147a9789af54b97fe342ca Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 29 Nov 2022 20:23:17 +0100 Subject: Complete starting values in summary for saem.mmkin fits Also update tests to the changes in mhmkin (see NEWS) --- R/summary.saem.mmkin.R | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) (limited to 'R/summary.saem.mmkin.R') diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index 2754e9f0..7337b0f3 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -225,13 +225,22 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3) obs = "Variance unique to each observed variable", tc = "Two-component variance function"), "\n") - cat("\nMean of starting values for individual parameters:\n") + cat("\nStarting values for degradation parameters:\n") print(x$mean_dp_start, digits = digits) cat("\nFixed degradation parameter values:\n") if(length(x$fixed$value) == 0) cat("None\n") else print(x$fixed, digits = digits) + cat("\nStarting values for random effects (square root of initial entries in omega):\n") + print(sqrt(x$so@model@omega.init), digits = digits) + + cat("\nStarting values for error model parameters:\n") + errparms <- x$so@model@error.init + names(errparms) <- x$so@model@name.sigma + errparms <- errparms[x$so@model@indx.res] + print(errparms, digits = digits) + cat("\nResults:\n\n") cat("Likelihood computed by importance sampling\n") print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, -- cgit v1.2.1 From cf73354f985be17352a4d538c6f9ea69f6be9aa9 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 2 Dec 2022 13:28:26 +0100 Subject: Avoid redundant warnings in summaries --- R/summary.saem.mmkin.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'R/summary.saem.mmkin.R') diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index 7337b0f3..46ab548b 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -136,10 +136,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = } # Correlation of fixed effects (inspired by summary.nlme) - varFix <- try(vcov(object$so)[1:n_fixed, 1:n_fixed]) - if (inherits(varFix, "try-error")) { + cov_so <- try(solve(object$so@results@fim), silent = TRUE) + if (inherits(cov_so, "try-error")) { object$corFixed <- NA } else { + varFix <- cov_so[1:n_fixed, 1:n_fixed] stdFix <- sqrt(diag(varFix)) object$corFixed <- array( t(varFix/stdFix)/stdFix, -- cgit v1.2.1 From a54bd290bc3884d0000c52c1b29bc557825d9eae Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 15 Dec 2022 14:50:28 +0100 Subject: List random effects correlations in output if any Update docs --- R/summary.saem.mmkin.R | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) (limited to 'R/summary.saem.mmkin.R') diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index 46ab548b..49b02a50 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -75,10 +75,21 @@ #' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo) #' print(f_saem_dfop_sfo) #' illparms(f_saem_dfop_sfo) -#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, covariance.model = diag(c(0, 0, 1, 1, 1, 0))) +#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, +#' no_random_effect = c("parent_0", "log_k_m1")) #' illparms(f_saem_dfop_sfo_2) #' intervals(f_saem_dfop_sfo_2) #' summary(f_saem_dfop_sfo_2, data = TRUE) +#' # Add a correlation between random effects of g and k2 +#' cov_model_3 <- f_saem_dfop_sfo_2$so@model@covariance.model +#' cov_model_3["log_k2", "g_qlogis"] <- 1 +#' cov_model_3["g_qlogis", "log_k2"] <- 1 +#' f_saem_dfop_sfo_3 <- update(f_saem_dfop_sfo, +#' covariance.model = cov_model_3) +#' intervals(f_saem_dfop_sfo_3) +#' # The correlation does not improve the fit judged by AIC and BIC, although +#' # the likelihood is higher with the additional parameter +#' anova(f_saem_dfop_sfo, f_saem_dfop_sfo_2, f_saem_dfop_sfo_3) #' } #' #' @export @@ -150,7 +161,8 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = # Random effects sdnames <- intersect(rownames(conf.int), paste0("SD.", pnames)) - confint_ranef <- as.matrix(conf.int[sdnames, c("estimate", "lower", "upper")]) + corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE) + confint_ranef <- as.matrix(conf.int[c(sdnames, corrnames), c("estimate", "lower", "upper")]) colnames(confint_ranef)[1] <- "est." # Error model -- cgit v1.2.1