2
$\begingroup$

I've been reading Simulation of logistic regression power analysis - designed experiments, http://sas-and-r.blogspot.com/2009/06/example-72-simulate-data-from-logistic.html, and Power analysis for ordinal logistic regression and I'm still a little lost on how to do a power analysis for my data. I want to be able to determine what N I should have if I have an interaction between a categorical variable (with 3 levels) and a continuous variable.

Simulation of logistic regression power analysis - designed experiments provides some information, but I can't figure out how to simulate the relationship between the categorical and continuous variables and outcome.

set.seed(1)

repetitions = 1000
N = 10000
n = N/8
var1 = c( .03, .03, .03, .03, .06, .06, .09, .09)
var2 = c( 0, 0, 0, 1, 0, 1, 0, 1)
rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)

var1 = rep(var1, times=n)
var2 = rep(var2, times=n)
var12 = var1**2
var1x2 = var1 *var2
var12x2 = var12*var2

significant = matrix(nrow=repetitions, ncol=7)

startT = proc.time()[3]
for(i in 1:repetitions){
responses = rbinom(n=N, size=1, prob=rates)
model = glm(responses~var1+var2+var12+var1x2+var12x2,
family=binomial(link="logit"))
significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)<br> > significant[i,6] = sum(significant[i,1:5])<br> > modelDev = model$null.deviance-model$deviance
significant[i,7] = (1-pchisq(modelDev, 5))<.05
}
endT = proc.time()[3]
endT-startT

sum(significant[,1])/repetitions # pre-specified effect power for var1
[1] 0.042
sum(significant[,2])/repetitions # pre-specified effect power for var2
[1] 0.017
sum(significant[,3])/repetitions # pre-specified effect power for var12
[1] 0.035
sum(significant[,4])/repetitions # pre-specified effect power for var1X2
[1] 0.019
sum(significant[,5])/repetitions # pre-specified effect power for var12X2
[1] 0.022
sum(significant[,7])/repetitions # power for likelihood ratio test of model
[1] 0.168
sum(significant[,6]==5)/repetitions # all effects power
[1] 0.001
sum(significant[,6]>0)/repetitions # any effect power
[1] 0.065
sum(significant[,4]&significant[,5])/repetitions # power for interaction terms
[1] 0.017

I feel like I should be able to adapt the code from Power analysis for ordinal logistic regression and that this would be a better, more succinct option

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 )

but I'm not completely sure how to do so. I'm assuming tmpfun(100,-1/2, 1/4,1/4) is specifying the N and betas that you want, but how do I adjust tmpfun to another (categorical) variable and include an interaction term? Ultimately the equation should include 6 betas: the intercept, the beta for x, the beta for z1, the beta for z2, the interaction term between x and z1, and the interaction term between x and z2.

Finally, I can't find any reliable sources on what sorts of "effect sizes" I should be using as small or medium.

Let me know if I can provide more information!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.