*! Help for -powersim- version 1.0.0 26July2013

help powersim

Title

Simulation-based power analysis for linear and generalized linear models

Syntax

powersim
[using filename] options
[, options] : model_cmd

Power
b(numlist) List of effect sizes, required

null(#) Value of the null hypothesis, defaults to zero

samplesize(numlist) List of sample sizes, required

alpha(#) alpha value, defaults to alpha = 0.05

position(#) column position of b in matrix e(b) of the analysis model, required

Data generating model
family(powersim##familyname:familyname) error distribution, required

link(powersim##linkname:linkname) link function, required*

cov1(n b d # [#]) covariate specifications, covariate 1

cov2(n b d # [#]) covariate specifications, covariate 2

cov3(n b d # [#]) covariate specifications, covariate 3

block22(b1 n1 b2 n2 [ b3 ]) binary variables for 2x2 block design

cons(b) constant term

corr12(#) Pearson correlation between cov1 and cov2

corr13(#) Pearson correlation between cov1 and cov3

corr23(#) Pearson correlation between cov2 and cov3

inter1(# n*n) Interaction effect specification

inter2(# n*n) Interaction effect specification

inter3(# n*n) Interaction effect specification

Simulation
nreps(#) number of Monte Carlo replications

inside Generate predictor data inside simulation loop

dryrun Dry run

detail Detailed output for dry run

gendata create a single dataset

nobs(#) Number of observations for created dataset

seed(#) Set random-number seed

expb Interpret effect sizes b as exp(b)

nodots Suppress replication dots

silent Suppress any intermediate output (implies nodots)

addscalar(string) Fetch scalar stored in e() of analysis model

maxiter(#) Set maximum number of iterations for the analysis model

force Allow analysis model command other then regress or glm

df(string) Use Students t-distribution for significance tests with df stored in e(string) of the analysis model

Saving
dofile(filename[,replace]) Save do-file to filename

saving(filename[,replace]) Save simulation results to filename

* specification of a link function is only required if the do-file is generated via command options;
if an existing do-file is provided, the link function must be defined in that do-file.

familyname

gaussian [#] Gaussian (normal), with SD of Gaussian error #; defaults to 1

igaussian [#] inverse Gaussian, with scale parameter #; defaults to 0.25

binomial [#] Binomial, with number of trials #; defaults to 1 (i.e., Bernoulli)

poisson Poisson

nbinomial [#] negative binomial, with overdispersion #; defaults to 1

gamma [#] gamma, with scale parameter #; defaults to 1

linkname

identity identity

log log

logit logit

probit probit

cloglog cloglog

power # power

opower # odds power

nbinomial negative binomial

loglog log-log

logc log-complement

covariate_distributions

normal a# b# Gaussian, with M a# and SD b#

poisson a# Poisson, with mean a#

uniform a# b# uniform, ranging from a# to b#

binomial a# bernoulli, i.e. binomial(1, a#)

chi2 a# chi-square, with a# df

studentt a# Student's t, with a# df

beta a# b# beta, with shape parameters a# and b#

gamma a# b# gamma, with shape parameter a# and scale parameter b#

nbinomial a# b# negative binomial, with number of failures a# and success probability b#

block a# equally sized groups, with number of groups a# (2-4 groups)

Description

powersim performs power analysis Monte Carlo simulations for linear and
generalized linear models. Data will be generated according to a model that
is specified via command options. Predictor data can be generated by either
using command options or by providing an existing do-file. Stata's regress
and glm commands can be used for specifying the analysis model. See
the accompanying tutorial for a more detailed introduction.

Options

Power

b(numlist) List of the true population parameters.

null(#) Value of the null hypothesis; if the expb option is specified,
the value will be interpreted as exponentiated value.

samplesize(numlist) List of sample sizes.

alpha(#) alpha level, default is alpha = 0.05.

position(#) column position of the coefficient in matrix e(b) (see ereturn list)
of the analysis model for which power wants to be estimated.

Data generating model

family(familyname) specifies the (conditional) distribution of the outcome.

link(linkname) specifies the link function.

cov1(n b d # [#]) covariate specifications for covariate 1, where
n is the name of the covariate, b its effect size, and d its distribution, followed
by one or two parameter values, depending on the distribution
(see covariate_distributions).
Type _bp instead of an effect size value in case of the covariate for
which the power wants to be simulated.

cov2(n b d # [#]) covariate specifications as described above, covariate 2.

cov3(n b d # [#]) covariate specifications as described above, covariate 3.

block22(b1 n1 b2 n2 [b3]) creates two binary variables for a 2x2 block design
with balanced group sizes; n1 and n2 are the names of the two variables and b1 and b2 are
its coefficients; a coefficient for an interaction effect of the two variables can be specified
with b3.

constant(b) b is the coefficient for the constant term; if this option is omitted,
the constant is assumed to equal zero.

corr12(#) specifies Pearson correlation coefficient for cov1 and cov2, applicable only to
normally distributed covariates; if this option is omitted, cov1 and cov2 are assumed to be
uncorrelated (orthogonal).

corr13(#) specifies Pearson correlation coefficient for cov1 and cov3, applicable only to
normally distributed covariates; if this option is omitted, cov1 and cov3 are assumed to be
uncorrelated (orthogonal).

corr23(#) specifies Pearson correlation coefficient for cov2 and cov3, applicable only to
normally distributed covariates; if this option is omitted, cov2 and cov3 are assumed to be
uncorrelated (orthogonal).

inter1(# n*n) specifies an interaction effect; for example, to generate data
with an interaction effect of 0.5 for cov2 and cov3 with names foo and bar,
type inter1(0.5 foo*bar).

inter2(# n*n) specifies an additional interaction effect.

inter3(# n*n) specifies an additional interaction effect.

Simulation

nreps(#) number of Monte Carlo replications.

inside Generates predictor data inside the simulation loop; should be
specified if the predictor data itself is regarded as a stochastic component.

dryrun Performs a dry run in order to check command inputs before running the simulations.

detail Additionally prints the contents of the data generating do-file on the screen
if the dryrun option is specified.

gendata creates a single realization of the data.

nobs(#) Determines the number of observations for the dataset requested with
the gendata option; defaults to N = 10,000.

seed(#) Sets the initial random-number seed.

expb Effect sizes specified in option b() are interpreted as exp(b); for example, this could be
useful if the power wants to be determined in the context of a logistic regression and the user wants
to specify effect sizes as odds ratios, rather than as log odds ratios.

nodots Suppresses the replication dots during the simulations. By default, one dot character is
displayed for each successful replication. A red `x' is displayed if model_cmd returns
an error (for example, if an analysis model does not converge within the specified maximum
number of iterations).

silent Suppresses any intermediate output during simulations (implies nodots).

addscalar Specifies to fetch a scalar that is available in e() of the analysis
model from each MC replication; for example, if you would like to record the value of the
log-likelihood from each replication of a glm, the option would be specified as addscalar(ll);
the scalar will then be appended to the resulting dataset which contains results from each replication.

maxiter(#) Specifies the number of maximum iterations for the analysis model, defaults to maxiter = 200;
will be ignored for model commands other than glm.

force Allow analysis model command other then regress or glm.

df(string) Use Students t-distribution for significance tests with df
stored in e(string) of the analysis model.

Saving

dofile(filename[,replace]) Save do-file to filename; if no existing do file is provided via
using then this option is required.

saving(filename[,replace]) Save simulation results from each replication to filename;

Remarks

Stata's regress and glm commands can be used to specify an analysis model under
model_cmd after the colon of the powersim command. The outcome variable in
model_cmd must be named y. With regress,
a t-distribution is used for siginificance testing and a normal distribution is used with
glm. Note that powersim removes the current data in memory and replaces
it with simulation results; make sure to save your data before running powersim. Also note
that powersim is creating a scalar with the name _bp. In case that a scalar with
the same name already exists, it will be overwritten. If multivariate
normal data is requested, Stata matrices with the names __M, __SD, and
__C are created which would replace any existing matrices with these names.

Note that both, distributional family and link function, need to be specified
for the data generating model. The available link functions are

Link function powersim option

identity link(identity)
log link(log)
logit link(logit)
probit link(probit)
complementary log-log link(cloglog)
odds power link(opower #)
power link(power #)
negative binomial link(nbinomial)
log-log link(loglog)
log-complement link(logc)

If an existing do-file is provided, the link function
must be defined in that do-file and only the distributional family
has to be specified. Available distributional families are

Family powersim option

Gaussian(normal) family(gaussian [#])
inverse Gaussian family(igaussian [#])
Binomial family(binomial [#])
Poisson family(poisson)
negative binomial family(nbinomial [#])
gamma family(gamma [#])

As it is also stated in glm, not all family/link
combinations make sense. You may choose from the following combinations:

| id log logit probit clog pow opower nbinomial loglog logc

+

Gaussian | x x x
inv. Gau. | x x x
binomial | x x x x x x x x x
Poisson | x x x
neg. bin. | x x x x
gamma | x x x

Examples

Estimating power for various effect and sample sizes for a linear effect of
a normally ditributed predictor variable which is correlated
with another predictor variable in a linear model:

powersim , ///

b(0.1(0.1)0.3) ///

alpha(0.05) ///

pos(1) ///

sample(100(20)300) ///

nreps(500) ///

family(gaussian) ///

link(identity) ///

cov1(x1 _bp normal 0 1) ///

cov2(x2 0.3 normal 0 1) ///

corr12(0.5) ///

inside ///

saving(psim_results) ///

dofile(psim_dofile) : reg y x1 x2

The assumed model in the above call to powersim can be expressed as:

y = _bp*x1 + 0.3*x2 + e, e ~ N(0,1)

where _bp is a placeholder for the effect sizes specified in b();
x1 and x2 are drawn from a multivariate normal distribution with rho = 0.5.

Using an existing do-file (here, the one that was created with the above example):

powersim using psim_dofile, ///

b(0.1(0.1)0.3) ///

alpha(0.05) ///

pos(1) ///

sample(100(20)300) ///

nreps(500) ///

family(gaussian) ///

inside ///

saving(psim_results, replace) : reg y x1 x2

Poisson model with two correlated predictor variables and power estimated
for the interaction effect:

powersim , ///

b(0.1) ///

alpha(0.05) ///

pos(3) ///

sample(300) ///

nreps(500) ///

family(poisson) ///

link(log) ///

cov1(x1 -0.25 normal 0 1) ///

cov2(x2 0.4 normal 0 1) ///

inter1(_bp x1*x2) ///

cons(0.5) ///

corr12(0.5) ///

inside ///

dofile(psim_dofile, replace) : ///

glm y c.x1##c.x2, family(poisson) link(log)

Logistic regression model for a 2x2 balanced block design with an interaction effect,
for which the statistical power is simulated:

powersim , ///

b(0.3) ///

alpha(0.05) ///

pos(8) ///

sample(2000) ///

nreps(500) ///

family(binomial) ///

link(logit) ///

block22(0.5 x1 0.5 x2 _bp)	///

cons(0.5) ///

dofile(psim_dofile, replace) : ///

glm y i.x1##i.x2, family(binomial) link(logit)

If power is estimated for several effect and/or sample sizes, power
can be plotted as a function of sample size (by effect sizes) after
simulations are completed (see powersimplot):

powersimplot

Power can also be plotted as a function of effect size if power is simulated
for varying effect sizes (see powersimplot). The following simulation
assumes a negative binomial model with overdispersion parameter alpha = 0.5
(note that if alpha does not equal 1, ml must be specified in the family()
option of glm in the analysis model command), log link function,
and a single binomial predictor variable, distributed Bi(1,0.5), for a fixed sample
size of N = 500; effect sizes here are specified as incidence rate ratios
via usage of the expb option:

powersim , ///

b(1(.05)1.5) ///

expb ///

alpha(0.05) ///

pos(1) ///

sample(500) ///

nreps(500) ///

family(nbinomial 0.5) ///

link(log) ///

cov1(x1 _bp binomial 0.5)	///

inside ///

dofile(psim_dofile, replace) : ///

glm y x1, family(nbinomial ml) link(log)

Now plot the power curve as power vs. effect size:

powersimplot, e

Simulating power for various effect and sample sizes for
the difference between two equally sized groups in an inverse
Gaussian model with scale parameter sigma = 0.5 and canonical
link function:

powersim , ///

b(0.2(0.1)0.5) ///

alpha(0.05) ///

pos(2) ///

sample(300(100)500) ///

nreps(250) ///

family(igaussian 0.5) ///

link(power -2) ///

cov1(x1 _bp block 2) ///

cons(1) ///

dofile(psim_dofile, replace) : ///

glm y i.x1, fam(igaussian) link(power -2)

Note that the random number generator for the inverse Gaussian distribution
may become slow, depending on specified predictor and scale parameters.

Saved results

powersim saves the following in r():

Scalars

r(niter) number of MC replications (requested, per sample/effect size combination)

r(alpha) alpha level

Macros

r(model) analysis model command line

r(cmd) powersim

r(effects) effect sizes

r(samples) sample sizes

r(iseed) initial random number seed

Matrices

r(power) results matrix

Acknowledgment

Ariel Linden suggested the silent option.

Author

Joerg Luedicke

Yale University and University of Florida

United States

email: joerg.luedicke@ufl.edu

Also see

Manual: [R] glm

Manual: [R] regress

Help: powersimplot