# Sample size calculation for univariate logistic regression

How does one calculate the sample size needed for a study in which a cohort of subjects will have a single continuous variable measured at the time of a surgery and then two years later they will be classified as functional outcome or impaired outcome.

We would like to see if that measurement could have predicted the bad outcome. At some point we may want to derive a cut point in the continuous variable above which we would try to intervene to diminish the probability of the impaired outcome.

Any ideas? Any R implementation.

• Do you expect some dropouts during follow-up? Are there any other covariates to be included in your model? – chl Sep 23 '10 at 6:34
• Let me suck a dropout rate out of my thumb - 20%. We will indeed collect many variables for instance, age, trauma score but I wanted to keep things as simple as possible for the power calculation. I have often found it useful to discuss a primary model and then secondary models that are loaded with more finesse and nuance. – Farrel Sep 25 '10 at 2:11
• Ok, but usually the expected % dropout, the number of covariates, and whether covariates are measured with errors (see e.g., j.mp/9fJkhb) enter the formula (in all case, it will increase the sample size). – chl Sep 25 '10 at 7:37

Sample size calculations for logistic regression are complex. I wont attempt to summarise it here. Reasonably accessible solutions to this problem are found in:

Hsieh FY. Sample size tables for logistic regression. Statistics in Medicine. 1989 Jul;8(7):795-802.

Hsieh FY, et al. A simple method of sample size calculation for linear and logistic regression. Statistics in Medicine. 1998 Jul 30;17(14):1623-34.

An accessible discussion of the issues with example calculations can be found in the last chapter (Section 8.5 pp 339-347) of Hosmer & Lemeshow's Applied Logistic Regression.

I usually find it easier and faster to run a simulation. Papers take a long time to read, to understand and finally come to the conclusion that they don't apply in the special case one is interested in.

Therefore, I would just pick a number of subjects, simulate the covariate you are interested in (distributed as you believe it will be), simulate good/bad outcomes based on the functional form you posit (threshold effects of the covariate? nonlinearity?) with the minimum (clinically) significant effect size you would like to detect, run the result through your analysis and see whether the effect is found at your alpha. Rerun this 10,000 times and look whether you found the effect in 80% of the simulations (or whatever other power you need). Adjust the number of subjects, repeat until you have a power you are happy with.

This has the advantage of being very general, so you are not confined to a specific functional form or a specific number or distribution of covariates. You can include dropouts, see chl's comment above, either at random or influenced by covariate or outcome. You basically code the analysis you are going to do on the final sample beforehand, which sometimes helps focus my thinking on the study design. And it is easily done in R (vectorize!).

• Do you have a worked case in R? – Farrel Sep 25 '10 at 2:07
• @Farrel - here's a very short script, which assumes [0,1]-uniformly distributed covariates, an OR of 2 between the first and third quartile of the covariate and standard normal noise, leading to power .34 for n=100. I'd play around with this to see how sensitive everything is to my assumptions: runs <- 1000; nn <- 100; set.seed(2010); detections <- replicate(n=runs,expr={covariate <- runif(nn); outcome <- runif(nn)<1/(1+exp(-2*log(2)*covariate+rnorm(nn))); summary(glm(outcome~covariate,family="binomial"))$coefficients["covariate","Pr(>|z|)"] < .05}) cat("Power:",sum(detections)/runs,"\n") – Stephan Kolassa Oct 2 '10 at 20:04 • You can attach your code as a pastie (pastebin.com) or a Gist (gist.github.com) if you feel it's more convenient, and link back to it in your comment. – chl Oct 2 '10 at 20:20 • @chl: +1, thanks a lot! Here's the gist: gist.github.com/607968 – Stephan Kolassa Oct 2 '10 at 20:41 • Great code but there is a problem. I am not as smart as you are. I need it broken down stepwise. I take it is runs the number of simulations? What is nn? Is it the number of subjects in the study? Then I see you created a distribution of covariates and made them determine a yes or a no depending on a threshold. – Farrel Oct 19 '10 at 3:03 Following on from the post by Stephan Kolassa (I can't add this as a comment), I have some alternative code for a simulation. This uses the same basic structure, but is exploded a bit more, so perhaps it is a little easier to read. It also is based on the code by Kleinman and Horton to simulate the logistic regression. nn is the number in the sample. The covariate should be continuously normally distributed, and standardized to mean 0 and sd 1. We use rnorm(nn) to generate this. We select an odds ratio and store it in odds.ratio. We also pick a number for the intercept. Choice of this number governs what proportion of the sample experience the "event" (e.g. 0.1, 0.4, 0.5). You have to play around with this number until you get the right proportion. The following code gives you a proportion of 0.1 with a sample size of 950 and an OR of 1.5: 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)  summary(proportion) confirms that the proportion is ~ 0.1 Then using the same variables, the power is calculated over 10000 runs: 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 think that this code is correct - I checked it against the examples given in Hsieh, 1998 (table 2), and it seems to agree with the three examples given there. I also tested it against the example on p 342 - 343 of Hosmer and Lemeshow, where it found a power of 0.75 (compared to 0.8 in Hosmer and Lemeshow). So it may be that in some circumstances this approach underestimates power. However, when I've run the same example in this on-line calculator, I've found that it agrees with me and not the result in Hosmer and Lemeshow.

If anyone can tell us why this is the case, I'd be interested to know.

• I do have 2 questions if you dont mind.1)Is the proportion function simply to get the intercept correct? 2) what is the logic behind using ytest (comparing prob to a random uni draw)? – B_Miner Mar 14 '11 at 19:18
• @B_Miner 1) The other way around - to get the proportion right, you need to set the intercept correctly - so adjust the intercept until you get the proportion which you are expecting. 2) The logic of ytest is that we need to get a dichotomous 0 or 1 outcome. So so we compare each sample from the uniform distribution to the probability (prob) to get our dichotomous outcome. The 'runis' doesn't have to be drawn from the random uniform distribution - a binomial or other distribution may make more sense for your data. Hope this helps (sorry for the delay in reply). – Andrew Mar 24 '11 at 13:34

a simple question about sample size is: how large a sample is needed to get a 95% confidence interval no longer than 2d for the [unknown] mean of the data distribution. another variant is: how large a sample is needed to have power 0.9 at $\theta = 1$ when testing H$_0: \theta = 0$. you don't seem to specify any criterion for choosing a sample size.

actually, it sounds as tho your study will be conducted in a sequential fashion. in that case, it may pay to make that an explicit part of the experiment. sequential sampling can often be more efficient than a fixed sample-size experiment [fewer observations needed, on average].