File: ci_wald.R

package info (click to toggle)
r-cran-parameters 0.11.0-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 3,044 kB
  • sloc: sh: 15; makefile: 2
file content (100 lines) | stat: -rw-r--r-- 3,481 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#' @rdname p_value_wald
#'
#' @param ci Confidence Interval (CI) level. Default to 0.95 (95\%).
#' @param dof Degrees of Freedom. If not specified, for \code{ci_wald()}, defaults to model's residual degrees of freedom (i.e. \code{n-k}, where \code{n} is the number of observations and \code{k} is the number of parameters). For \code{p_value_wald()}, defaults to \code{Inf}.
#'
#' @inheritParams simulate_model
#' @inheritParams standard_error
#' @inheritParams model_parameters.default
#'
#' @importFrom stats qt coef
#' @export
ci_wald <- function(model, ci = .95, dof = NULL, effects = c("fixed", "random", "all"), component = c("all", "conditional", "zi", "zero_inflated", "dispersion", "precision", "scale", "smooth_terms", "full", "marginal"), robust = FALSE, ...) {
  effects <- match.arg(effects)
  component <- match.arg(component)
  out <- lapply(ci, function(i) {
    .ci_wald(model = model, ci = i, dof = dof, effects = effects, component = component, robust = robust, method = "wald", ...)
  })
  out <- do.call(rbind, out)
  row.names(out) <- NULL
  out
}


#' @importFrom insight get_parameters n_obs
#' @importFrom stats qt complete.cases
#' @keywords internal
.ci_wald <- function(model, ci, dof, effects, component, robust = FALSE, method = "wald", se = NULL, ...) {
  if (inherits(model, "emmGrid")) {
    params <- insight::get_parameters(model, effects = effects, component = component, merge_parameters = TRUE)
  } else {
    params <- insight::get_parameters(model, effects = effects, component = component)
  }

  # check if all estimates are non-NA
  params <- .check_rank_deficiency(params, verbose = FALSE)

  method <- tolower(method)

  # if we have adjusted SE, e.g. from kenward-roger, don't recompute
  # standard errors to save time...
  if (is.null(se)) {
    stderror <- if (isTRUE(robust)) {
      standard_error_robust(model, ...)
    } else {
      switch(
        method,
        "wald" = standard_error(model, component = component),
        "kenward" = ,
        "kr" = se_kenward(model),
        "ml1" = se_ml1(model),
        "betwithin" = se_betwithin(model),
        "satterthwaite" = se_satterthwaite(model),
        standard_error(model, component = component)
      )
    }

    if (.is_empty_object(stderror)) {
      return(NULL)
    }

    # filter non-matching parameters
    if (nrow(stderror) != nrow(params)) {
      params <- stderror <- merge(stderror, params, sort = FALSE)
    }
    se <- stderror$SE
  }

  if (is.null(dof)) {
    # residual df
    dof <- degrees_of_freedom(model, method = "any")
    # make sure we have a value for degrees of freedom
    if (is.null(dof) || length(dof) == 0) {
      dof <- Inf
    } else if (length(dof) > nrow(params)) {
      # filter non-matching parameters
      dof <- dof[1:nrow(params)]
    }
  }

  alpha <- (1 + ci) / 2
  fac <- suppressWarnings(stats::qt(alpha, df = dof))
  out <- cbind(
    CI_low = params$Estimate - se * fac,
    CI_high = params$Estimate + se * fac
  )

  out <- as.data.frame(out)
  out$CI <- ci * 100
  out$Parameter <- params$Parameter

  out <- out[c("Parameter", "CI", "CI_low", "CI_high")]
  if ("Component" %in% names(params)) out$Component <- params$Component
  if ("Effects" %in% names(params) && effects != "fixed") out$Effects <- params$Effects

  if (anyNA(params$Estimate)) {
    out[stats::complete.cases(out), ]
  } else {
    out
  }
}