## Logistic regression – simulation for a power calculation…

Please note-I’ve spotted a problem with the approach taken in this post – it seems to underestimate power in certain circumstances. I’ll post again with a correction or a more full explanation when I’ve sorted it.

So, I posted an answer on cross validation regarding logistic regression. I thought I’d post it in a little more depth here, with a few illustrative figures. It’s based on the approach which Stephen Kolassa described.

Power calculations for logistic regression are discussed in some detail in Hosmer and Lemeshow (Ch 8.5). One approach with R is to simulate a dataset a few thousand times, and see how often your dataset gets the p value right. If it does 95% of the time, then you have 95% power.

In this code we use the approach which Kleinman and Horton use to simulate data for a logistic regression. We then initially calculate the overall proportion of events. To change the number of events adjust odds.ratio. The independent variable is assumed to be normally distributed with mean 0 and variance 1.

nn <- 950 runs <- 10000 intercept <- log(9) odds.ratio <- 1.5 beta <- log(odds.ratio) proportion <- replicate( n = runs, expr = { xtest <- rnorm(nn) linpred <- intercept + (xtest * beta) prob <- exp(linpred)/(1 + exp(linpred)) runis <- runif(length(xtest),0,1) ytest <- ifelse(runis < prob,1,0) prop <- length(which(ytest <= 0.5))/length(ytest) } ) summary(proportion)

This plot shows how the intercept and odds ratio affect the overall proportion of events per trial:

When you’re happy that the proportion of events is right (with some prior knowledge of the dataset), you can then fit a model and calculate a p value for that model. We use R’s inbuilt function replicate to do this 10,000 times, and count the proportion where it gets it right (i.e. p < 0.05). The proportion of the time that the simulation correctly get's the p < 0.05 is essentially the power of the logistic regression for your number of cases, odds ratio and intercept.

result <- replicate( n = runs, expr = { xtest <- rnorm(nn) linpred <- intercept + (xtest * beta) prob <- exp(linpred)/(1 + exp(linpred)) runis <- runif(length(xtest),0,1) ytest <- ifelse(runis < prob,1,0) summary(model <- glm(ytest ~ xtest, family = "binomial"))$coefficients[2,4] < .05 } ) print(sum(result)/runs)

I checked it against the examples given in Hsieh, 1999. It seemed to work pretty well calculating the power to be within ~ 1% of the power of the examples given in table II of that paper.

We can do some interesting things with R. I simulated a range of odds ratios and a range of sample sizes. The plot of these looks like this (each line represents an odd ratio):-

We can also keep the odds ratio constant, but adjust the proportion of events per trial. This looks like this (each line represents an event rate):

As ever, if anyone can spot an error or suggest a simpler way to do this then let me know. I haven’t tested my simulation against any packages which calculate power for logistic regression, but if anyone can it would be great to hear from you.

Could you please help with a similar power calculation except for a log binomial regression?

Thanks,

Rebecca

@rrwpress Could you give me more idea of what you’re trying to do (e.g. design of the study, expected distribution of the predictor variable etc?). Then I can tailor it more for your needs.

Andrew

That would be great. (:

It’s an observational study with 30%-40% expected to have health impairment, which is defined by being in the top 20% in 2 out of 3 groups, number of prescriptions, number of doctor’s visits and number of sickness days. Thank you! Rebecca

So, just so I’m clear – health impairment (yes / no) is your binomial outcome, and your predictor variables are number of prescriptions, number of doctor’s visits and number of sickness days?

Yes, Health impairment is the binomial outcome. At this stage it is defined by those three groups or variables, number of prescriptions, number of doctor’s visits, and number of sickness days (but it might change in the future). Let me explain. If a person is in the top 20% for number of prescriptions, for example then they get a 1, and if are also in the top 20% of another group like number of doctor’s visits then then also get 1. So, they need to be in the top 20% of at least 2 out of 3 groups to get a 1 for health impairment. That’s the dependent variable. Then there are at least three predictor variables, one of them a baseline variable. Something general for log binomial regression or negative binomial regression with some predictor variables would be perfect since the dependent variable is quite fluid at this time. Am I being clear?

Ok. I’m not overly familiar with log binomial regression. Give me a bit of time to sort out simulating it and I’ll post something.