[R] Is my simulation to compute power of a multiple ordinal logistic regression right?

[BONACHE Adrien]

Good day,
I was performing a power analysis of articles published in a journal of management using the pwr package in R. However, it seemed to be impossible to compute power for small, medium and large Effect Size for multiple ordinal logistic regression output. I have tried using G*power, but it seemed to be just useful when we have simple logistic regression output. Thus, I have tried to simulate to calculate the power based on Greg Snow and Gung answers.    library(rms)

    tmpfun <- function(n, beta1,beta2,beta3,beta4,beta5,beta6,beta7,
    beta8,beta9,beta10,beta11,beta12,beta13,beta14) {   
    var1=log(runif(n,3,33))
    var2=rbinom(n,1,0.4502)
    x=1
    while (x>0.96) {var3=rbinom(n,1,0.9474); x=mean(var3)}
    var4=rbinom(n,1,0.2749)
    var5=rbinom(n,1,0.6199)
    var6=rbinom(n,1,0.4327)
    var7=rnorm(n,1.11,0.29)
    var8=rnorm(n,0.06,0.07)
    var9=runif(n,106,536804)
    var10=runif(n,0.08,3.7)
    var11=runif(n,0.08,0.92)
    var12=rnorm(n,96.24,113.11)
    for (i in 1:n) {if (var12[i]<=0) {var12[i]=runif(1)} else   
    {var12[i]=log(var12[i])}}
    var13=runif(n,0.05,0.88)
    var14=runif(n,-1.06,0.03)
    eta1 <- beta1*var1+beta2*var2+beta3*var3+beta4*var4+beta5*var5+
    beta6*var6
    eta2=eta1+beta7*var7+beta8*var8+beta9*var9+beta10*var10
    +beta11*var11+beta12*var12+beta13*var13+beta14*var14
    p1 <- exp(eta1)/(1+exp(eta1))
    p1=replace(p1,is.na(p1),1)
    p2 <- exp(eta2)/(1+exp(eta2))
    p2=replace(p2,is.na(p2),1)
    tmp <- runif(n)
    y <- (tmp < p1)+(tmp<p2)
    fit <- lrm(y~var1+var2+var3+var4+var5+var6)
    fit$stats[5]
    }
    out <- replicate(10000, tmpfun(n=112,0.582,-0.176,-0.413,
    -0.138,0.861,-0.942,3.481,4.542,1.059,0.322,-2.562,0.599,2.4,0.732))
    mean( out < 0.05 )
Is my simulation right? How can I get power for small, medium and large ES? Those are the questions... Thank for your time and your help.
Regards,
Adrien. 
	[[alternative HTML version deleted]]