Power analysis for ordinal logistic regression (multiple independent variables)

I have a question with regard to a simulation based power analysis for ordinal logistic regression suggested here: https://stats.stackexchange.com/a/22410/231675 (by Greg Snow)

The suggested power analysis includes only one indepednent variable, and for the sake of my study I need to include three. I re-wrote the original code, but I am not 100% sure if it is correct. It would be great if someone could have a look at it.

Here is the original simulation (with one independent variable only):

library(rms)

tmpfun <- function(n, beta0, beta1, beta2) {
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 <- lrm(y~x)
fit$stats[5] } out <- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4)) mean( out < 0.05 )  And here is my code (with three independent variables): library(rms) library(truncnorm) tmpfun <- function(n, beta0, beta1, beta2, beta3) { x1 <- rbinom(n, 1, 0.5) #condition x2 <- rnorm(n, mean = 0, sd = 1) #GP x3 <- rnorm(n, mean = .02, sd = .64) #interaction eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3 p1 <- exp(eta1)/(1+exp(eta1)) tmp <- runif(n) y <- (tmp < p1) fit <- lrm(y~x1 + x2 + x3) fit$stats[5]
}

out <- replicate(1000, tmpfun(500, 0.231, 0.714, 0.020, 0.127))
mean( out < 0.05 )

• It may be better to merge this in with the original topic. Note that when Y has more than 50 levels the R rms orm function will be faster than lrm. If you want the most general Wald statistic approach you can manipulate the result of anova(fit) as a matrix. This will handle nonlinearities, interactions, adjusted effects, etc. – Frank Harrell Mar 12 at 11:40