The Asypow S(plus) Library for Asymptotic Power Calculations

The asypow library consists of routines written in the S language that calculate power and related quantities utilizing asymptotic methods. A paper describing these methods with examples is in preparation [1]. Two methods are available. The likelihood ratio method (LR) is described in [2]. Another general method appears recently in [3]; and we designate it the SMO method after the initials of the authors.


Introduction
The asypow library consists of routines written in the S language that calculate power and related quantities utilizing asymptotic methods. A paper describing these methods with examples is in preparation [1]. Two methods are available. The likelihood ratio method (LR) is described in [2]. Another general method appears recently in [3]; and we designate it the SMO method after the initials of the authors.
Here we outline the overall steps in the use of asymptotic sample size and power calculations.
1. The statistician specifies a complete parametric alternative hypothesis model.
The specification includes the form of the model model, values of its parameters and the design. If the model involves the comparison of a single parameter over multiple groups, then the design consists of the proportion of subjects in each group. In more complex cases, covariates are involved. If the experimenter controls the covariate values, the design consists of these values and the relative number of observations at each. In an observational study in which the covariate values are not under experimental control, a probability distribution of covariate values must be specified.
2. The designer specifies the constraints on the model parameters that transform the alternative hypothesis into the null hypothesis. Each constraint either sets a parameter value to a constant or it posits the equality of several parameters without specifying their common value. Let the number of constraints be C.
If there are equality constraints, the model must be reparameterized so that there are only constant constraints. This reparameterization is achieved by replacing each set of parameters, equal according to the null hypothesis, by their differences and by the value of one of them. Fortunately, this process can be automated.
Both methods described use the model from step 1 and the constraints from step 2 to compute the noncentrality parameter, η, of a noncentral χ 2 distribution that describes the asymptotic distribution of the likelihood ratio under the alternative hypothesis.
3. The critical value, V , is computed from the significance level, α, from 1 − α = χ 2 C (V ) where χ 2 C (.) is the central χ 2 cumulative distribution function with C degrees of freedom. Power is the probability that a random variable distribued as a noncentral χ 2 with C degrees of freedom and noncentrality parameter, η, exceeds V . If sample size is to be calculated, an iterative calculation finds the requisite noncentrality parameter.

Documentation/Help
The asypow routines, unfortunately, are inherently more difficult to use than many others in Splus. This is a result of the need to specify a comparatively large amount of information for each calculation.
To mitigate this problem, a cheatsheet and an online routine, asypow.help, are available as well as the traditional help files. The routine, asypow.help prints most of the text available in the cheatsheet as items chosen from a menu.

Likelihood Ratio (LR) Method
(This presentation follows Cox and Hinkley [2, Section 9.3].) Let y = {y 1 , . . . , y n } be a realization of independent, identically distributed random variables,{Y 1 , . . . , Y n }, having a density f (y; θ), where θ is a vector of unknown parameters. The log likelihood of the observations, y, at the parameter value θ is l(y; θ) = n s=1 log(f (y s ; θ)). (1) Letθ denote the alternative hypothesis value of θ. The expected information matrix i is defined by This expectation can be written as where the integration is over the set of possible values of Y . Note that the expectation is calculated at the assumed true value of θ. If the Y i assume discrete values, f (y|θ) is the probability that Y = y, and the above integral is replaced by a sum over all possible outcomes, y.
The presence of covariates in the model requires another level of computation to obtain the information matrix. When the experimenter chooses covariate values, the information matrix is computed separately for each value, and the results are summed. In observational studies in which covariates are random, the information matrix is integrated, component by component, over the population covariate distribution.
The null hypothesis determines the values of specific components of θ, which is partitioned into (ψ, λ), where the first set of parameter values, ψ, have been set by the null hypothesis to the value ψ 0 . The remaining parameters, λ, are not constrained by the null hypothesis.
The expected information matrix, i, is partitioned in correspondence to θ into Letθ be the maximum likelihood estimate of θ,θ 0 be the value that maximizes the likelihood subject to ψ = ψ 0 , andψ be the presumed true value of ψ. The results central to LR power and sample-size calculations are as follows: The asymptotic distribution of 2[l(θ) − l(θ 0 )] is noncentral χ 2 , χ 2 C (x|η) with noncentrality parameter and degrees of freedom, C equal to the number of elements in ψ 0 , i.e., the number of constraints.

SMO Method
This method is named with the initials of its authors, [3]. The method is easily motivated from the assumption that asymptotically, 2(l(y;θ)−l(y;θ 0 )) is distributed as χ 2 C (x|η), as per the likelihood ratio calculations but with a different value of η. The model parameters are again partitioned into sets constrained and unconstrained by the null hypothesis. Let the partitioning yieldθ 0 = (ψ 0 ,λ 0 );λ 0 is that value of λ that maximizes Eθ(l(Y |(ψ 0 , λ)).
For large n,θ →θ andθ 0 →θ 0 , so 2Eθ{l(y;θ) − l(y;θ 0 )} → 2Eθ{l(y;θ) − l(y;θ 0 )} The expectation of the non-central χ 2 is C + η, so by equating expectations, we have the approximation SMO used an expansion around the null hypothesis parameter values that produced a term involving first and second derivatives of the expected log likelihood. They found that in practice, this term nearly cancels C above, and arrive at the above formula without this term. η = 2Eθ{l(y;θ) − l(y;θ 0 )} Equations (6) and (7) above give two means of estimating the noncentrality parameter. The more conservation includes the degrees of freedom term C.
In simple cases, the value ofλ 0 is immediate, but in more complex cases it must be determined through nonlinear maximization. The LR method implicitly uses an asymptotic approximation toλ 0 , a value which usually does not maximize the log likelihood under the null hypothesis. This may be one reason that the SMO method offers superior performance in some cases.
4 Example: Single Sample Poisson PROBLEM: What sample size is required to reject with power 0.8 the null hypothesis that a Poisson mean is λ = 2.0 when λ = 3.0 using a two-sided test with 0.05 significance level?

Noncentrality Parameter: LR
The only parameter of the model is the mean, λ, so the information matrix is a scalar. The Poisson density is For a single observation, the log likelihood is l(y; λ) = y * log(λ) − λ − log(y!).
The second derivative of the log likelihood with respect to lambda is Since the number of events, y, is distributed Poisson the expectation of y is λ, the information matrix is i = 1 λ .
The expected information for a single observation is 1 3 = 0.33, so for n observations, the expected information will be 0.33n. Using (4), the noncentrality parameter for a null hypothesis value of λ 0 is (λ 0 − 3.0) 2 * 0.33n, which at λ 0 = 2.0 is 0.33n.

Noncentrality Parameter: SMO
The Poisson density is For a single observation, the log likelihood is The expected value of the log likelihood of λ a , the alternative hypothesis value of λ over λ 0 is given by the Poisson density at λ a times the log likelihood at λ 0 . NOTES: is the expected value of y, i.e., λ a .
• ∞ y=0 f (y; λ a ) is the sum of the density of y, i.e., 1. Then, The most conservative SMO methods for estimating the noncentrality parameters for n observations shown in (6) is Hence, the noncentrality parameter for a null hypothesis value of λ 0 = 2.0 and alternative hypothesis of λ = 3.0 is .43n − 1.
A less conservative methods given in (7) does not include the degrees of freedom term, C, and the noncentrality parameter is .43n.

Requisite Noncentrality Parameter
We will reject the null hypothesis if twice the LR is greater than the 0.95 quantile of the (central) χ 2 1 distribution, i.e., 3.84. The probability of exceeding 3.84 is where χ 2 d (x|ν) is the cumulative distribution function of the noncentral chi-square distribution with d degrees of freedom and noncentrality parameter ν. The value of ν making the probability of exceeding 3.84 equal to 0.8 is 7.85.

Sample Size
Equating the noncentrality parameters to 7.85, we obtain n = 24 for the LR method and n = 20 for the more conservative SMO method or n = 18 for the less conservative SMO method.

The Routines
The asypow library provides routines for calculating the sample size, power or significance level for hypothesis testing. The available asypow routines can be divided into four categories: The argument <theta.ha> is named after the parameters of the distribution, for example, p.ha for the binomial distribution and lambda.ha for the poisson distribution. Some routines have extra options. The exact calling sequence and the model parameters for each routine are available in the individual help files. Each routine returns a list containing all values used in the construction of the problem. Printing the return value produces a reconstruction of the original problem. The various components of the list can be identified using the 'names()' function.

Choose a method.
The two asymptotic methods discussed are available for calculations. To choose one, set the argument method to "lr" for the likelihood ratio method or "smo" for the SMO method.  Each model allows for more than one treatment group. If there is more than one treatment group in the model, the parameter values must be set for each group. This is done by making the argument a matrix with a row for each group. For example, if there were three groups in the binomial models and the value of p was 0.4 for the first group, 0.3 for the second group and 0.9 for the third group p.ha would be set to a column vector. The routines allow for unequal sample sizes for models with more than one group. If this is the case then set the argument group.size to be a vector of the relative sample sizes. The value of the i t h component is the relative sample size of the i t h group. If this value is specified, it should be a vector whose length is the same as the number of rows in <theta.ha>, the alternative hypothesis parameters.
For example, in the three group binomial model above if 50% of the subjects are assigned to group 3 while 25% are assigned to group 1 and the final 25% to group 2, then group.size would be > group.size <-c(.25,.25,.5).
If the model has only one group or there are equal sample sizes in each group then the argument group.size does not need to be set.

Specify constraints on the parameters that transform the alterna-
tive hypothesis into the null hypothesis.
There are two types of constraints. Constant constraints set a parameter to a specific value and equality constraints posit the equality of two or more parameters.

Setting constant constraints.
Constant constraints are defined by three values, a character description of the parameter, the group number of the parameter, and the value the parameter is being set to. These are set in the arguments: The probabilities for r-1 categories. ordinal.kgp p1,p2,. . . ,p(r-1) The probabilities for r-1 categories.
Note that for the multinomial and ordinal model with r categories, there are only r-1 parameters. This is because the probabilities of the categories must add to one, hence the final probability is redundant. For a more detailed description of the model parameters available see the individual help files.
The arguments const.constraint.group and const.constraint.param do not need to be specified if there is only one group or one parameter in the model respectively. If the const.constraint.* arguments are not specified then there are no constant constraints.
For example, if the null hypothesis in a one group binomial is to set the parameter 'p' to 0.5 then the argument values would be set in the following fashion > const.constraint.group <-1 > const.constraint.param <-"p" > const.constraint.value <-0.5 Since there is only one group and one parameter value in the model the values for const.constraint.group and const.constraint.param do not need to be set.
For more than one equality constraint, the arguments equal.constraint.param and equal.constraint.group should be a list of vectors where each element in the list defines an equality constraint. For example, in a three category ordinal model the parameters are the category probabilities "p1" and "p2". Suppose the model has two groups and the null hypothesis is that the category probabilities in the 2 groups are the same, then set > equal.constraint.param <-list(c("p1","p1"),c("p2","p2")) > equal.constraint.group <-list(c (1,2),c(1,2)).

Entering constraints at the terminal.
If all five constraint arguments: const.constraint.param const.constraint.group const.constraint.value equal.constraint.group equal.constraint.value are missing, then the routines will query the user to enter the constraints interactively. This is often easier for the user.
6.1.8 Set the significance level, power and sample size.
The <model>.kgp routines can calculate the significance level, power or sample size for a hypothesis test given the other two values. The user should enter values for two of the three arguments: significance, power or sample.size. The third will be calculated by the routine. To make multiple calculations at once the values of significance, power or sample.size can be vectors. The missing argument can be specified by a -1 or not sending argument to the function.
For example, to calculate the sample size at a power of 0.8 and significance level of 0.01, 0.05 and 0.1 set > significance <-c(0.01,0.05,0.1) > power <-0.8 6.1.9 Set the degrees of freedom option for the SMO method.
The final argument for the routines is smo.df. This is a logical variable which is only relevant if method is "smo". If smo.df is TRUE sample size is calculated via the equation η = 2Eθ{l(y;θ)−l(y;θ 0 )}−C where η is the non-centrality parameter of the Chi-Square model and C is the degrees of freedom. This is the most conservative calculation. If FALSE, sample size is calculated via the equation η = 2Eθ{l(y;θ) − l(y;θ 0 )}. By default, the more conservative calculation is made.

Univariate Design Models
Popular likelihood based models posit that the parameter of the distribution depends on a covariate, x, through design parameters, a,b,c,.... A linear combination of the covariates is widely used to model this dependence. For example, let The asypow routines allow for a linear or quadratic combination. The function which ties u to the value of the parameter is called a link function. Different link functions are appropriate depending on the model.
The following steps explain the use of the asypow routines for design models, (binomial.design, poisson.design, expsurv.design, ordinal.design). A typical calling sequence for a designed models is show below where <model> is replaced with the name of the model being used. Some routines have extra options. The exact calling sequence and the model parameters for each routine are available in the individual help files.

Choose a method.
The argument method functions the same in the design routines as the no-covariate routines to select the likelihood ratio or SMO method for the calculations.

Select a design.
The parameter design specifies the function of the covariate x that will be used. Linear indicates, u = a + bx, and quadratic indicates, u = a + bx + cx 2 . Design is "linear" for a linear combination and "quadratic" for a quadratic combination. Note: For the ordinal model with r categories a linear design specifies u[i] = a[i]+bx and a quadratic design specifies u[i] = a[i]+bx+cx 2 for i = 1, . . . , r−1. The model assumes a different intercept term for each category but the other regression terms are the same for all categories.

Select a link function.
The parameter link specifies the function which transforms the linear combination of the covariate into the distribution parameters. The available link functions depend on the model.

Routine
Links Available Link Formula binomial.design logistic The link is a parameter only in the routines binomial.design and ordinal.design, where there is more than one available.

Set the alternative hypothesis.
Define the parameter values at the alternative hypothesis. For the binomial, poisson and clinical trial with exponential survival models the parameters of the model are "a" and "b" for a linear design and "a", "b" and "c" for a quadratic design. For an r category ordinal design the parameters of the model are "a1","a2",. . . up to the number of categories minus one, plus "b" for a linear design or "b" and "c" for a quadratic design.
The value of these parameters under the alternative hypothesis are defined in the argument theta.ha. theta.ha is a vector with the values of the parameters under the null hypothesis. For example, in a linear binomial design where a is 1 and b is .3 set  Each model allows for more than one treatment group. If there is more than one treatment group in the model, the parameter values must be set for each group. This is done by making the argument a matrix with a row for each group. For example, theta.ha in a three group linear binomial model might be > theta.ha <-rbind(c (1,.3),c(2,.4),c(3,.5)).

Set the covariate values.
The argument xpoints is a vector containing the covariate values in the trial. If there is more than one group in the trial and each group has different covariate values, then xpoints is a matrix where each row represents a group. In this case, xpoints must have the same number of rows as theta.ha. 6.2.6 Set the number of observations at each covariate value.
The routines allow for unequal sample sizes at each covariate value in the design. If this is the case then define the argument natx so that natx[i,j] is the relative sample size at covariate xpoints [i,j].
For example, if a trial has 2 covariate values defined in xpoints as > xpoints <-c (1,2) and 60% of the subjects have covariate value 1 while 40% of the subjects have covariate value 2 then set natx > natx <-c(0.6,0.4) 6.2.7 Set the percentage of subjects in each group.
The routines allow for unequal sample sizes for models with more than one group. If this is the case then set the argument group.size to be a vector of the relative sample sizes. The argument group.size functions the same in the design routines as the no-covariate routines. These arguments are the same as for the no-covariate routines.
6.2.9 Set the significance level, power and sample size.
The <model>.design routines can calculate the significance level, power or sample size for a hypothesis test given the other two values. The user should enter values for two of the three arguments: significance, power or sample.size. The third will be calculated by the routine. These arguments function the same in the design routines as the no-covariate routines.
6.2.10 Set the degrees of freedom option for the SMO method.
The final argument for the routines is smo.df which determines if the degrees of freedom are used in the SMO calculations. This argument functions the same in the design routines as the no-covariate routines.

Multivariate Design Models
The asypow library provides two functions which provide a binomial design where the probability of an event is a function of k multiple covariates, x [1],. . . ,x[k] via a combination of coefficient coef [1],. . . ,coef[k].
In the multivariate logistic design (mvlogistic.design) the probability of an event is a logistic function of u where u = In the multivariate log-linear design (mvloglin.design) the probability of an event is an exponential function of u where u = log(coef [1] The usual use of these routines is for tabulated data in which case the x's will all be 0 or 1 valued indicator variables.
The following steps explain the use of the asypow routines for multivariate design models, mvlogistic.design and mvloglin.design. A typical calling sequence for a designed models is show below where <model> is replaced with the name of the model being used, either logistic or loglin for log-linear. Use of these routines is similar to using other asypow routines.

Choose a method.
This argument method functions the same in the design routines as the no-covariate routines to select the likelihood ratio or SMO method for the calculations.
6.3.2 Set the alternative hypothesis.
The parameters for multivariate designs are named "coef1", "coef2",. . . up to the number of covariates in the model. The value of these parameters under the null hypothesis are defined in the argument coef.ha.

Set the design points.
The design points are the values of the k covariates included in the trial. These are set in the arguments xpoints. xpoints is a matrix of dimension (n × k), each row of which gives values of the k covariates at one of the n design points. Note: Most models will include a constant term and the column of xpoints corresponding to this term will be identically 1.

Set the proportion of subjects at each design point.
The routines allow for unequal sample sizes at each design point in xpoints. If this is the case then the define the argument rss so that rss[i] is the relative sample size at the design point defined in row i of xpoints.
For example, if for the three design points defined in xpoints above we put 50% of the subject at design point 1 and 25% of the subjects at design points 2 and 3 rss would be > rss <-c(0.5,0.25,0.25). 6.3.5 Specify constraints on the parameters that transform the alternative hypothesis into the null hypothesis.
The constraints are defined by the arguments: • const.constraint.param • const.constraint.value • equal.constraint.param.
There are no group indicators since these models are not programmed to consider multiple groups. Setting constraints works the same as it does with other asypow routines.
6.3.6 Set the significance level, power and sample size.
The mv<model>.design routines can calculate the significance level, power or sample size for a hypothesis test given the other two values. The user should enter values for two of the three arguments: significance, power or sample.size. The third will be calculated by the routine. These arguments function the same in the design routines as the no-covariate routines.