Power analysis for a fitted ordinal regression model

[jarbet]

I know post-hoc power analysis using your observed data is controversial, but nevertheless I want to try estimating power for the effect of a continuous predictor in an ordinal regression model for my dataset. I used an LLM to review the RESI vignette and manuscript, and came up with the following approach. Please let me know if this looks correct!

# Load required packages
library(MASS)   # for polr() ordinal regression
library(RESI)   # for RESI functions
#> Registered S3 method overwritten by 'clubSandwich':
#>   method    from    
#>   bread.mlm sandwich

# Prepare data
data(mtcars)
mtcars$gear <- ordered(mtcars$gear)  # treat gear as ordinal

# Fit ordinal regression model
model_polr <- polr(gear ~ hp + wt, data = mtcars, Hess = TRUE)
summary(model_polr)
#> Call:
#> polr(formula = gear ~ hp + wt, data = mtcars, Hess = TRUE)
#> 
#> Coefficients:
#>       Value Std. Error t value
#> hp  0.02124   0.008842   2.403
#> wt -3.30999   0.987097  -3.353
#> 
#> Intercepts:
#>     Value   Std. Error t value
#> 3|4 -7.6948  2.3863    -3.2247
#> 4|5 -4.4913  1.8843    -2.3836
#> 
#> Residual Deviance: 43.46576 
#> AIC: 51.46576


# Extract coefficients and standard errors
coefs <- summary(model_polr)$coefficients
wald_stats <- coefs[, "t value"]
#df <- nrow(mtcars) - length(coefs[, "t value"])  # approximate degrees of freedom

# Residual degrees of freedom
df.resid <- model_polr$df.residual;

### Use RESI's t2S function to compute signed effect sizes

# Compute RESI for each predictor
resi_hp <- t2S(t = wald_stats["hp"], rdf = df.resid, n = nrow(mtcars))
resi_wt <- t2S(t = wald_stats["wt"], rdf = df.resid, n = nrow(mtcars))

cat("RESI for hp:", resi_hp, "\n")
#> RESI for hp: 0.4132453
cat("RESI for wt:", resi_wt, "\n")
#> RESI for wt: -0.5767323


# Power calculation for hp
alpha <- 0.05
df_test <- 1 # number of parameters being tested
lambda <- nrow(mtcars) * resi_hp^2

power_hp <- 1 - pchisq(qchisq(1 - alpha, df_test), df_test, ncp = lambda)
cat("Estimated power for hp:", round(power_hp, 3), "\n")
#> Estimated power for hp: 0.647

^{Created on 2025-08-21 with reprex v2.1.1}