20
$\begingroup$

I am looking for a program (in R or SAS or standalone, if free or low cost) that will do power analysis for ordinal logistic regression.

$\endgroup$
1
  • $\begingroup$ I agree, It would help to know exactly how many ordinal categories are being modeled here (it looks like maybe 3 to me), what beta0 and beta2 are (thresholds?), and how the -1/2, 1/4, and 1/4 were chosen. $\endgroup$
    – Craig L-S
    Jul 16, 2021 at 0:09

3 Answers 3

36
$\begingroup$

I prefer to do power analyses beyond the basics by simulation. With precanned packages, I am never quite sure what assumptions are being made.

Simulating for power is quite straight forward (and affordable) using R.

  1. decide what you think your data should look like and how you will analyze it
  2. write a function or set of expressions that will simulate the data for a given relationship and sample size and do the analysis (a function is preferable in that you can make the sample size and parameters into arguments to make it easier to try different values). The function or code should return the p-value or other test statistic.
  3. use the replicate function to run the code from above a bunch of times (I usually start at about 100 times to get a feel for how long it takes and to get the right general area, then up it to 1,000 and sometimes 10,000 or 100,000 for the final values that I will use). The proportion of times that you rejected the null hypothesis is the power.
  4. redo the above for another set of conditions.

Here is a simple example with ordinal regression:

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 )
$\endgroup$
16
  • 7
    $\begingroup$ +1, this is a very robust, universal approach. I have often used it. I'd like to suggest another feature: You can generate data for the max $N$ you would consider, then fit the model for proportions of those data by successively fitting the 1st n of them at regular intervals up to $N$ (eg, n=100, 120, 140, 160, 180, & 200). Instead of saving a p-value from each generated dataset, you can save a row of p-values. Averaging over each column gives you a quick and dirty sense of how power is changing w/ $N$, and helps you hone in on an appropriate value quickly. $\endgroup$ Feb 8, 2012 at 2:57
  • 2
    $\begingroup$ @gung: yours comment make sense, would you mind adding your codes so that less experience people in R could also benefit from it? thanks $\endgroup$
    – user12349
    Jul 2, 2012 at 14:06
  • 1
    $\begingroup$ I am looking at this again and I have a couple questions: 1) Why is x uniform on 1:10? 2) How would you generalize it to more than 1 independent variable? $\endgroup$
    – Peter Flom
    Feb 22, 2015 at 15:48
  • 2
    $\begingroup$ @PeterFlom, x had to be something, so I chose (arbitrarily) to have it uniform between 0 and 10, it could have also been normal, gamma, etc. Best would be to choose something that is similar to what we expect the real x variables to look like. To use more than 1 predictor variable, generate them independently (or from a multivariate normal, copula, etc.) then just include them all in the eta1 piece, e.g. eta1 <- beta0 + beta1*x1 + beta2*x2 + beta3*x3. $\endgroup$
    – Greg Snow
    Feb 23, 2015 at 18:07
  • 3
    $\begingroup$ @Hack-R, the above code is for ordinal logistic regression, or proportional odds logistic regression, where there are 3 ordered levels in the response variable, e.g. low, medium, and high. So beta_0 and beta_1 together create eta1 which translates to the probability of being in the medium or high group (anything above low), then beta_2 is added to give eta2 which is the probability of being in the high group. I probably should have come up with a better name for beta_2 to convey this better. $\endgroup$
    – Greg Snow
    Sep 27, 2021 at 22:08
3
$\begingroup$

Besides Snow's excellent example, I believe you can also do a power simulation by resampling from an existing dataset which has your effect. Not quite a bootstrap, since you're not sampling-with-replacement the same n, but the same idea.

So here's an example: I ran a little self-experiment which turned in a positive point-estimate but because it was little, was not nearly statistically-significant in the ordinal logistic regression. With that point-estimate, how big a n would I need? For various possible n, I many times generated a dataset & ran the ordinal logistic regression & saw how small the p-value was:

library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
    d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
    lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
    return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
   bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
   print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}

With the output (for me):

[1] 300.0000   0.1823
[1] 330.0000   0.1925
[1] 360.0000   0.2083
[1] 390.0000   0.2143
[1] 420.0000   0.2318
[1] 450.0000   0.2462
[1] 480.000   0.258
[1] 510.0000   0.2825
[1] 540.0000   0.2855
[1] 570.0000   0.3184
[1] 600.0000   0.3175

In this case, at n=600 the power was 32%. Not very encouraging.

(If my simulation approach is wrong, please someone tell me. I'm going off a few medical papers discussing power simulation for planning clinical trials, but I'm not at all certain about my precise implementation.)

$\endgroup$
1
  • $\begingroup$ We need more answers involving nootropics as the example data :) $\endgroup$
    – Hack-R
    Sep 27, 2021 at 16:58
3
$\begingroup$

I would add one other thing to Snow's answer (and this applies to any power analysis via simulation) - pay attention to whether you are looking for a 1 or 2 tailed test. Popular programs like G*Power default to 1-tailed test, and if you are trying to see if your simulations match them (always a good idea when you are learning how to do this), you will want to check that first.

To make Snow's run a 1-tailed test, I would add a parameter called "tail" to the function inputs, and put something like this in the function itself:

 #two-tail test
  if (tail==2) fit$stats[5]

  #one-tail test
  if (tail==1){
    if (fit$coefficients[5]>0) {
          fit$stats[5]/2
    } else 1

The 1-tailed version basically checks to see that the coefficient is positive, and then cuts the p-value in half.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.