simRegOrd  R Documentation 
This function simulates the power of a twosample test from a proportional odds ordinal logistic model for a continuous response variable a generalization of the Wilcoxon test. The continuous data model is normal with equal variance. Nonlinear covariate adjustment is allowed, and the user can optionally specify discrete ordinal level overrides to the continuous response. For example, if the main response is systolic blood pressure, one can add two ordinal categories higher than the highest observed blood pressure to capture heart attack or death.
simRegOrd(n, nsim=1000, delta=0, odds.ratio=1, sigma,
p=NULL, x=NULL, X=x, Eyx, alpha=0.05, pr=FALSE)
n 
combined sample size (both groups combined) 
nsim 
number of simulations to run 
delta 
difference in means to detect, for continuous portion of response variable 
odds.ratio 
odds ratio to detect for ordinal overrides of continuous portion 
sigma 
standard deviation for continuous portion of response 
p 
a vector of marginal cell probabilities which must add up to one.
The 
x 
optional covariate to adjust for  a vector of length

X 
a design matrix for the adjustment covariate 
Eyx 
a function of 
alpha 
type I error 
pr 
set to 
a list containing n, delta, sigma, power, betas, se, pvals
where
power
is the estimated power (scalar), and betas, se,
pvals
are nsim
vectors containing, respectively, the ordinal
model treatment effect estimate, standard errors, and 2tailed
pvalues. When a model fit failed, the corresponding entries in
betas, se, pvals
are NA
and power
is the proportion
of nonfailed iterations for which the treatment pvalue is significant
at the alpha
level.
Frank Harrell
Department of Biostatistics
Vanderbilt University School of Medicine
fh@fharrell.com
popower
## Not run:
## First use no ordinal highend category overrides, and compare power
## to ttest when there is no covariate
n < 100
delta < .5
sd < 1
require(pwr)
power.t.test(n = n / 2, delta=delta, sd=sd, type='two.sample') # 0.70
set.seed(1)
w < simRegOrd(n, delta=delta, sigma=sd, pr=TRUE) # 0.686
## Now do ANCOVA with a quadratic effect of a covariate
n < 100
x < rnorm(n)
w < simRegOrd(n, nsim=400, delta=delta, sigma=sd, x=x,
X=cbind(x, x^2),
Eyx=function(x) x + x^2, pr=TRUE)
w$power # 0.68
## Fit a cubic spline to some simulated pilot data and use the fitted
## function as the true equation in the power simulation
require(rms)
N < 1000
set.seed(2)
x < rnorm(N)
y < x + x^2 + rnorm(N, 0, sd=sd)
f < ols(y ~ rcs(x, 4), x=TRUE)
n < 100
j < sample(1 : N, n, replace=n > N)
x < x[j]
X < f$x[j,]
w < simRegOrd(n, nsim=400, delta=delta, sigma=sd, x=x,
X=X,
Eyx=Function(f), pr=TRUE)
w$power ## 0.70
## Finally, add discrete ordinal category overrides and high end of y
## Start with no effect of treatment on these ordinal event levels (OR=1.0)
w < simRegOrd(n, nsim=400, delta=delta, odds.ratio=1, sigma=sd,
x=x, X=X, Eyx=Function(f),
p=c(.98, .01, .01),
pr=TRUE)
w$power ## 0.61 (0.3 if p=.8 .1 .1, 0.37 for .9 .05 .05, 0.50 for .95 .025 .025)
## Now assume that odds ratio for treatment is 2.5
## First compute power for clinical endpoint portion of Y alone
or < 2.5
p < c(.9, .05, .05)
popower(p, odds.ratio=or, n=100) # 0.275
## Compute power of ttest on continuous part of Y alone
power.t.test(n = 100 / 2, delta=delta, sd=sd, type='two.sample') # 0.70
## Note this is the same as the p.o. model power from simulation above
## Solve for OR that gives the same power estimate from popower
popower(rep(.01, 100), odds.ratio=2.4, n=100) # 0.706
## Compute power for continuous Y with ordinal override
w < simRegOrd(n, nsim=400, delta=delta, odds.ratio=or, sigma=sd,
x=x, X=X, Eyx=Function(f),
p=c(.9, .05, .05),
pr=TRUE)
w$power ## 0.72
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.