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

## ----load-data----------------------------------------------------------------
library(lsReg)

datafile <- system.file("extdata", "simulated_data.rds", package = "lsReg")
dat <- readRDS(datafile)

basemdl <- glm(ylin ~ x1 + x2, data = dat, family = gaussian)

## ----run-all-tests------------------------------------------------------------
testtypes <- c("lrt", "score", "robustscore", "wald", "robustwald")
zvars     <- c("z5", "z9")

results <- expand.grid(testtype = testtypes, variable = zvars,
                       stringsAsFactors = FALSE)
results$statistic <- NA_real_
results$pvalue    <- NA_real_

for (tt in testtypes) {
  mem <- lsReg(basemdl, colstoadd = 1, testtype = tt)
  for (zv in zvars) {
    xr  <- as.matrix(dat[, zv, drop = FALSE])
    addcovar(mem, xr)
    stat <- if (tt == "lrt") mem$testvalue[1] else mem$testvalue[1, 1]
    pval <- if (tt == "lrt") pchisq(stat, df = 1, lower.tail = FALSE) else
                             2 * pnorm(-abs(stat))
    idx  <- results$testtype == tt & results$variable == zv
    results$statistic[idx] <- stat
    results$pvalue[idx]    <- pval
  }
}

## ----results------------------------------------------------------------------
# Reshape for a readable side-by-side comparison
res_z5 <- results[results$variable == "z5", c("testtype", "statistic", "pvalue")]
res_z9 <- results[results$variable == "z9", c("testtype", "statistic", "pvalue")]

cat("=== z5 ===\n")
print(res_z5, digits = 5, row.names = FALSE)

cat("\n=== z9 ===\n")
print(res_z9, digits = 5, row.names = FALSE)

## ----compare-scales-----------------------------------------------------------
for (zv in zvars) {
  lrt_z   <- sqrt(results$statistic[results$testtype == "lrt"  & results$variable == zv])
  wald_z  <-      results$statistic[results$testtype == "wald" & results$variable == zv]
  score_z <-      results$statistic[results$testtype == "score"& results$variable == zv]
  cat(zv, "— sqrt(LRT):", round(lrt_z, 4),
          " Wald:", round(wald_z, 4),
          " Score:", round(score_z, 4), "\n")
}