Ordinal Logistic Regression Power Analysis (original) (raw)

# adopted from Greg Snow's answer on CV:
#  https://stats.stackexchange.com/questions/22406/power-analysis-for-ordinal-logistic-regression/22410#22410

Step 1: Define design and NA condition to solve

library(SimDesign)

Design <- createDesign(n = NA,
                       beta0 = -1/2, beta1 = 1/4, beta2 = 1/4)
Design

# A tibble: 1 × 4
      n beta0 beta1 beta2
  <dbl> <dbl> <dbl> <dbl>
1    NA  -0.5  0.25  0.25

Step 2 — Define analyse and summarise functions

Analyse <- function(condition, dat, fixed_objects = NULL) {
    Attach(condition)
    x <- runif(n, 0, 10)
    eta1 <- beta0 + beta1*x
    eta2 <- eta1 + beta2
    p1 <- exp(eta1)/(1+exp(eta1))
    p2 <- exp(eta2)/(1+exp(eta2))
    tmp <- runif(n)
    y <- (tmp < p1) + (tmp < p2)
    fit <- rms::lrm(y~x)
    fit$stats[5]
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    ret <- ERR(results, alpha = .05)  ## empirical rejection rate given alpha
    ret
}

Step 3 — Optimize over the rows in design

# find n associated with f(n) = 1 - B = .90 power
# terminate when prediction CI between [.895, .905]
solved <- SimSolve(design=Design, b=.9, interval=c(30,300), parallel=TRUE,
                   analyse=Analyse, summarise=Summarise,
                   predCI.tol=.01, maxiter = 200, verbose=FALSE)
solved

# A tibble: 1 × 4
       n beta0 beta1 beta2
   <dbl> <dbl> <dbl> <dbl>
1 96.631  -0.5  0.25  0.25

summary(solved)  # note that prediction CI is within [.895, .905]

$root
[1] 96.63137

$predCI.root
  CI_2.5  CI_97.5 
95.29939 97.81995 

$b
[1] 0.9

$predCI.b
[1] 0.8951701 0.9046310

$terminated_early
[1] TRUE

$time
[1] 33.62s

$iterations
[1] 80

$total.replications
[1] 22900

$tab
           y   x reps
5  0.8837736  92 2650
6  0.8909910  93 1110
7  0.8993399  94 3030
8  0.8964602  95 2260
9  0.8859873  96 1570
10 0.8960606  97 3300
11 0.9071429  98 2100
12 0.9082667  99 3750
13 0.9186667 100 1500
15 0.9093750 102  320

plot(solved, type = 'history')