|
1 | 1 | #' Perform Reliability Diagnostics on Survey Regression Models |
2 | 2 | #' |
3 | 3 | #' @description |
4 | | -#' This function takes a fitted survey regression model object (e.g., from `svyglm`) |
5 | | -#' and produces a tibble with key reliability and diagnostic metrics for each |
6 | | -#' coefficient. |
| 4 | +#' This function takes a fitted survey regression model object (e.g., from `svyglm` |
| 5 | +#' or `svycoxph`) and produces a tibble with key reliability and diagnostic |
| 6 | +#' metrics for each coefficient. |
7 | 7 | #' |
8 | 8 | #' @details |
9 | 9 | #' The output provides a comprehensive overview to help assess the stability and |
|
21 | 21 | #' even if precisely estimated. It is better to rely on the standard error, |
22 | 22 | #' p-value, and confidence interval width for reliability assessment. |
23 | 23 | #' |
24 | | -#' @param fit A fitted model object, typically of class `svyglm`. |
| 24 | +#' @param fit A fitted model object from the `survey` package, such as `svyglm` or `svycoxph`. |
25 | 25 | #' @param p_threshold A numeric value (between 0 and 1) for the significance threshold. Defaults to `0.05`. |
26 | 26 | #' @param rse_threshold A numeric value for flagging high Relative Standard Error (RSE). Defaults to `30`. |
27 | 27 | #' |
|
40 | 40 | #' \item \code{is_rse_high}: A logical flag, `TRUE` if `RSE_percent` is greater than or equal to `rse_threshold`. |
41 | 41 | #' } |
42 | 42 | #' |
43 | | -#' @importFrom dplyr mutate bind_cols select |
44 | | -#' @importFrom tibble as_tibble |
45 | | -#' @importFrom stats confint setNames |
| 43 | +#' @importFrom dplyr mutate select |
| 44 | +#' @importFrom tibble tibble |
| 45 | +#' @importFrom stats confint coef vcov |
46 | 46 | #' |
47 | 47 | #' @export |
48 | 48 | #' |
|
82 | 82 | #' family = quasibinomial() |
83 | 83 | #' ) |
84 | 84 | #' |
85 | | -#' # 3. Get the reliability diagnostics table |
86 | | -#' diagnostics_table <- svyglmdiag(fit) |
| 85 | +#' # 3. Get the reliability diagnostics table using the new function |
| 86 | +#' diagnostics_table <- svydiag(fit) |
87 | 87 | #' |
88 | 88 | #' # Print the resulting table |
89 | 89 | #' print(diagnostics_table) |
|
96 | 96 | #' } |
97 | 97 | #' } |
98 | 98 |
|
99 | | -svyglmdiag <- function(fit, p_threshold = 0.05, rse_threshold = 30) { |
| 99 | +svydiag <- function(fit, p_threshold = 0.05, rse_threshold = 30) { |
100 | 100 |
|
101 | | - # --- Input validation --- |
102 | | - if (!inherits(fit, "svyglm")) { |
103 | | - warning("This function is designed for 'svyglm' objects. Results may be unexpected.") |
104 | | - } |
| 101 | + # 1. Robustly extract key model components using accessor functions |
| 102 | + s_fit <- summary(fit) |
| 103 | + estimates <- stats::coef(fit) |
| 104 | + se <- sqrt(diag(stats::vcov(fit))) |
| 105 | + conf_int <- stats::confint(fit) |
105 | 106 |
|
106 | | - # 1. Get the standard model summary and confidence intervals |
107 | | - summary_fit <- summary(fit) |
108 | | - conf_int_fit <- stats::confint(fit) |
| 107 | + # P-values are most reliably extracted from the summary coefficient table. |
| 108 | + # This assumes the p-value is the last column, which is standard for most |
| 109 | + # survey models (svyglm, svycoxph, etc.). |
| 110 | + p_vals <- s_fit$coefficients[, ncol(s_fit$coefficients)] |
109 | 111 |
|
110 | 112 | # 2. Combine these into a single, informative table |
111 | | - reliability_df <- tibble::as_tibble(summary_fit$coefficients, rownames = "Term") |
112 | | - names(reliability_df) <- c("Term", "Estimate", "SE", "t.value", "p.value") |
| 113 | + reliability_df <- tibble::tibble( |
| 114 | + Term = names(estimates), |
| 115 | + Estimate = estimates, |
| 116 | + SE = se, |
| 117 | + p.value = p_vals, |
| 118 | + CI_Lower = conf_int[, 1], |
| 119 | + CI_Upper = conf_int[, 2] |
| 120 | + ) |
113 | 121 |
|
114 | | - # 3. Add CIs, calculate metrics, and add flags |
| 122 | + # 3. Calculate derived metrics, add flags, and finalize the output |
115 | 123 | reliability_df <- reliability_df %>% |
116 | | - dplyr::bind_cols(tibble::as_tibble(conf_int_fit) %>% |
117 | | - stats::setNames(c("CI_Lower", "CI_Upper"))) %>% |
118 | 124 | dplyr::mutate( |
119 | 125 | RSE_percent = (SE / abs(Estimate)) * 100, |
120 | 126 | CI_Width = CI_Upper - CI_Lower, |
|
0 commit comments