Nothing
#' Compute (pseudo-) residuals
#'
#' @description
#' This function computes (pseudo-) residuals of an \code{\link{fHMM_model}}
#' object.
#'
#' @param x
#' An object of class \code{\link{fHMM_model}}.
#' @param verbose
#' Set to \code{TRUE} (default) to print progress messages.
#'
#' @return
#' An object of class \code{\link{fHMM_model}} with residuals included.
#'
#' @export
#'
#' @examples
#' compute_residuals(dax_model_3t)
#' res <- residuals(dax_model_3t)
#' summary(res)
#'
#' @importFrom stats pt pgamma qnorm
compute_residuals <- function(x, verbose = TRUE) {
### check input
if (!inherits(x,"fHMM_model")) {
stop("'x' must be of class 'fHMM_model'.", call. = FALSE)
}
if (!isTRUE(verbose) && !isFALSE(verbose)) {
stop("'verbose' must be either TRUE or FALSE.", call. = FALSE)
}
if (is.null(x$decoding)) {
warning(
paste(
"Cannot compute residuals without decoding.",
"Please call 'decode_states()' first."
),
immediate. = TRUE, call. = FALSE
)
return(x)
}
### function that computes the residuals
cr <- function(data, mus, sigmas, dfs, sdd_name, decoding) {
stopifnot(length(data) == length(decoding))
out <- rep(NA_real_, length(data))
for (t in seq_along(data)) {
if (sdd_name == "t") {
Fxt <- stats::pt(
q = (data[t] - mus[decoding[t]]) / sigmas[decoding[t]],
df = dfs[decoding[t]]
)
}
if (sdd_name == "gamma") {
Fxt <- stats::pgamma(
q = data[t],
shape = mus[decoding[t]]^2 / sigmas[decoding[t]]^2,
scale = sigmas[decoding[t]]^2 / mus[decoding[t]]
)
}
if (sdd_name == "lnorm") {
Fxt <- stats::plnorm(
q = data[t],
meanlog = mus[decoding[t]],
sdlog = sigmas[decoding[t]]
)
}
out[t] <- stats::qnorm(Fxt)
}
return(out)
}
### compute residuals
par <- parUncon2par(x$estimate, x$data$controls)
if (!x$data$controls$hierarchy) {
residuals <- cr(
data = x$data$data, mus = par$mus, sigmas = par$sigmas,
dfs = par$dfs, sdd_name = x$data$controls$sdds[[1]]$name,
decoding = x$decoding
)
} else {
residuals <- matrix(NA_real_, nrow = nrow(x$data$data),
ncol = ncol(x$data$data))
residuals[, 1] <- cr(
data = x$data$data[, 1], mus = par$mus,
sigmas = par$sigmas, dfs = par$dfs,
sdd_name = x$data$controls$sdds[[1]]$name,
decoding = x$decoding[, 1]
)
for (t in seq_len(nrow(residuals))) {
curr <- x$decoding[t, 1]
residuals[t, -1] <- cr(
data = x$data$data[t, -1],
mus = par$mus_star[[curr]],
sigmas = par$sigmas_star[[curr]],
dfs = par$dfs_star[[curr]],
sdd_name = x$data$controls$sdds[[2]]$name,
decoding = x$decoding[t, -1]
)
}
}
### save residuals in 'x' and return 'x'
if (verbose) {
message("Computed residuals")
}
class(residuals) <- "fHMM_residuals"
x$residuals <- residuals
return(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.