Skip to main content
Added the online resource, and also stated the difference in power with the different methods.
Source Link
Andrew
  • 6.3k
  • 5
  • 26
  • 23

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'd love any feedback on I also tested it against the code, and examplesexample on p 342 - 343 of Hosmer and Lemeshow, where it breaks downfound 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 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'd love any feedback on the code, and examples of where it breaks down.

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.

Changed the reference to "above" as answers can move around, and explained what I meant by "explode"
Source Link
Andrew
  • 6.3k
  • 5
  • 26
  • 23

Following on from the post by Stephan Kolassa above (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.

Following on from the post by Stephan Kolassa above (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. It also is based on the code by Kleinman and Horton to simulate the logistic regression.

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.

Source Link
Andrew
  • 6.3k
  • 5
  • 26
  • 23

Following on from the post by Stephan Kolassa above (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. 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'd love any feedback on the code, and examples of where it breaks down.