Sample size calculation for univariate logistic regression

Question

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.

Answer 1

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!).

Answer 2

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.

Answer 3

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.

Answer 4

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].

farrel: i'm adding this in reply to your comment.

to get at a sample size, one usually specifies some sort of precision criterion for an estimate [such as length of a CI] OR power at a specified alternative of a test to be carried out on the data. you seem to have mentioned both of these criteria. there is nothing wrong with that, in principle: you just have to then do two sample size calculations - one to achieve the desired estimation precision - and another to get the desired power at the stated alternative. then the larger of the two sample sizes is what is required. [btw - other than saying 80% power - you don't seem to have mentioned what test you plan to perform - or the alternative at which you want the 80% power.]

as for using sequential analysis: if subjects are enrolled in the study all at the same time, then a fixed sample size makes sense. but if the subjects are few and far between, it may take a year or two [or more] to get the required number enrolled. thus the trial could go on for three or four years [or more]. in that case, a sequential scheme offers the possibility of stopping sooner than that - if the effect[s] you are looking for become statistically significant earlier in the trial.