## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, echo = FALSE, message = FALSE-------------------------------------
library(svytest)
library(survey)

## -----------------------------------------------------------------------------
# Construct survey design
des <- svydesign(ids = ~1, weights = ~FINLWT21, data = svytestCE)

# Fit weighted regression
fit <- svyglm(TOTEXPCQ ~ ROOMSQ + BATHRMQ + BEDROOMQ + FAM_SIZE + AGE, design = des)

## -----------------------------------------------------------------------------
# Run all diagnostic tests
all_results <- run_all_diagnostic_tests(fit, alpha = 0.05)
print(all_results)

## -----------------------------------------------------------------------------
# Run DC test with equal residual variance
res_equal <- diff_in_coef_test(fit, var_equal = TRUE)
print(res_equal)

# Run DC test with heteroskedasticity-robust variance (HC3)
res_robust <- diff_in_coef_test(fit, var_equal = FALSE, robust_type = "HC3")
summary(res_robust)

## -----------------------------------------------------------------------------
results <- wa_test(fit, type = "DD")
print(results)

## -----------------------------------------------------------------------------
results <- wa_test(fit, type = "PS1")
print(results)

## -----------------------------------------------------------------------------
results <- wa_test(fit, type = "PS2")
print(results)

## -----------------------------------------------------------------------------
results <- wa_test(fit, type = "WF")
print(results)

## -----------------------------------------------------------------------------
linear_results <- estim_eq_test(fit, q_method = "linear")
print(linear_results)

log_results <- estim_eq_test(fit, q_method = "log")
print(log_results)

## -----------------------------------------------------------------------------
perm_mean_results <- perm_test(fit, stat = "pred_mean", B = 1000, engine = "R")
print(perm_mean_results)

library(Rcpp)
perm_mahal_results <- perm_test(fit, stat = "coef_mahal", B = 1000, engine = "C++")
print(perm_mahal_results)

## ----echo=FALSE, results='asis', message = FALSE------------------------------
library(knitr)
library(dplyr)

reject_table <- readRDS(system.file("extdata/st1_reject_table.rds", package = "svytest"))
reject_table_perm <- readRDS(system.file("extdata/perm_reject_table.rds", package = "svytest")) 

reject_table |> 
  left_join(reject_table_perm, by = c("n", "sigma", "delta", "alpha")) |>
  select(n, sigma, delta, alpha,
         DD, PS1, PS1q, PS2, PS2q, HP, PS3, pred_mean, coef_mahal) |> 
  kable(digits = 3,
        caption = "Empirical rejection rates (Study 1 design). Columns correspond to diagnostic tests; rows to simulation scenarios.")