0
$\begingroup$

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 )
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Frank Harrell Mar 12 at 11:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.