#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# power.README
# power.fcn
# powertst.S
# pfnc.d
# pchisqnc.d
# glhpwr.d
# sscomp.d
# This archive created: Monday 15 May 1991
export PATH; PATH=/bin:/usr/bin:$PATH
cat << \SHAR_EOF > 'power.README'
power
A suite of Splus functions to compute power and sample size for
tests of the general linear hypothesis. The functions use the
power approximations of Muller and Peterson (1984). These
functions allow you to compute the approximate power of tests of
the null hypothesis
H0: c%*%beta%*%u - theta0 = 0,
in the model
y = x%*%beta + e.
Here y(Nxp) is the matrix of responses, x(Nxq) is the design
matrix, beta(qxp) is the true parameter, and e(Nxp) is a matrix
of errors whose rows are independent draws from the multivariate
normal distribution with mean 0 and variance matrix sigma(pxp).
The functions also provide for computation of the minimum sample
size needed to achieve a given level of power.
To install, go to a likely directory and type `sh filename'
(without the quotes), where `filename' is the file in which you've
stored this shell archive. This will place 7 files in your
directory:
1. power.README (this file)
2. power.fcn (containing the functions)
3. powertst.S (some numerical examples)
4. pfnc.d (documentation for the noncentral F function)
5. pchisqnc.d (documentation for the noncentral chi^2 function)
6. glhpwr.d (documentation for the power function)
7. sscomp.d (documentation for the sample size function).
Files ending in .d are help() files, and should be moved into a
.Data/.Help directory; see prompt() in the S manual. To
accommodate different styles of use, I leave storage details to
the user.
The functions were written entirely in Splus versions 3.0 and 2.3
on a Sun SPARCstation 1 running SunOS release 4.1; they have not
been tested under any other conditions. The author is willing to
help with problems. Send email to dheitjan@biostats.hmc.psu.edu
on Internet.
REFERENCE
Muller, K.E. and Peterson, B.L. (1984). Practical methods for
computing power in testing the multivariate general linear
hypothesis. Computational Statistics and Data Analysis 2,
143--158.
Written by Daniel F. Heitjan (dheitjan@biostats.hmc.psu.edu).
COPYRIGHT NOTICE
These functions are not copyrighted. The author requests,
however, that you acknowledge the source and communicate any
improvements to him.
SHAR_EOF
cat << \SHAR_EOF > 'power.fcn'
## S functions for calculating power and sample size requirements
## when testing the general linear hypothesis.
## Daniel F. Heitjan, 25 November 1991
##
pfnc_function(q,df1,df2,lm,iprec=c(6)) {
## Daniel F. Heitjan, 15 May 1991
iprec_iprec[1]
nrw_max(c(length(q),length(df1),length(df2),length(lm)))
if (length(q)1) {
for (i in 2:ninit) xx_rbind(xx,x) }
n_ninit
pwr_glhpwr(alpha,xx,bt1,sigma,cc,u,th0,test)
if (pwr>power) {
while ((pwr>power)&(n>0)) {
xx_matrix(xx[(-1:(-nrx)),],ncol=ncol(xx))
n_n-1
opwr_pwr
pwr_glhpwr(alpha,xx,bt1,sigma,cc,u,th0,test)
}
pwr_opwr
n_n+1 }
else if (pwr 'powertst.S'
## Test the power and sample size routines in power.fcn. This
## takes about five minutes on a SPARCstation 1.
## Daniel F. Heitjan, 25 November 1991
##
## First source the file containing the functions.
source('power.fcn')
## Try to recompute a portion of Table 4.2 of Kraemer and Thiemann.
## Results differ slightly because K&T use one-sided 5% tests.
alpha_c(.05)
betavec_c(.90,.99)
x_matrix(1,1,1)
bt0vec_c(6:10)/10
sigma_matrix(1,1,1)
cc_matrix(1,1,1)
u_matrix(1,1,1)
th0_matrix(0,1,1)
tbl_matrix(0,5,2)
dimnames(tbl)_list(c('.6','.7','.8','.9','1.0'),c('.90','.99'))
init_matrix(c(31,25,18,17,1,53,41,32,24,21),5,2)
## Compute the table.
for (i in 1:5) {
for (j in 1:2) {
tbl[i,j]_sscomp(alpha,betavec[j],x,matrix(bt0vec[i],1,1),
sigma,cc,u,th0,1,init[i,j])$n }}
cat('\nPart of Table 4.2 of Kraemer and Thiemann\n')
cat('Results differ slightly because the tests differ\n')
print(tbl)
##
## Now try to recompute the ANOVA example on p. 454 of Neter and
## Wasserman (1974).
cat('\n\nANOVA power calculations from Neter&Wasserman (1974)\n')
tbl_matrix(c(.72,0,.56,0,.36,0),2,3)
x_matrix(c(1,1,rep(0,8),
0,0,1,1,1,0,0,0,0,0,
0,0,0,0,0,1,1,1,0,0,
0,0,0,0,0,0,0,0,1,1),10,4)
dimnames(tbl)_list(c('N&W','Ours'),c('Ex. 1','Ex. 2','Ex. 3'))
bt0_matrix(c(-3.5,-3,2,5),4,1)
sigma_matrix(2.5^2,1,1)
cc_t(matrix(c(1,-1,0,0,0,1,-1,0,0,0,1,-1),4,3))
u_diag(1)
th0_matrix(0,3,1)
tbl[2,1]_glhpwr(alpha,x,bt0,sigma,cc,u,th0)
sigma_matrix(3.0^2,1,1)
tbl[2,2]_glhpwr(alpha,x,bt0,sigma,cc,u,th0)
alpha_.01
sigma_matrix(2.5^2,1,1)
tbl[2,3]_glhpwr(alpha,x,bt0,sigma,cc,u,th0)
print(round(tbl,3))
##
## Now let's try a repeated-measures example from p.161 of
## Morrison (1976).
cat('\n\nA repeated-measures example from Morrison (1976)\n')
cat('Morrison gives the power as 93%; the difference is due to\n')
cat('approximation of the noncentrality parameter.\n')
tbl_matrix(c(0),1,3)
dimnames(tbl)_list(c('Power'),c('Test W','Test HLT','Test PB'))
alpha_.05
x_matrix(1,8,1)
bt0_matrix(c(6,2,0),1,3)
sigma_4*matrix(1,3,3)+6*diag(3)
cc_diag(1)
u_matrix(c(1,-1,0,0,1,-1),3,2)
th0_0*cc%*%bt0%*%u
tbl[1,]_c(glhpwr(alpha,x,bt0,sigma,cc,u,th0,1),
glhpwr(alpha,x,bt0,sigma,cc,u,th0,2),
glhpwr(alpha,x,bt0,sigma,cc,u,th0,3))
print(round(tbl,3))
##
## Here's example 5B.1, p. 217 of Johnson and Wichern (1982), who
## consider the power of Hotelling's T^2 test.
cat(
'\n\nA repeated-measures example from Johnson&Wichern (1982)\n')
cat('Again differences are due to the approximation\n')
alpha_.05
x_matrix(1,26,1)
bt0.big_matrix(1,1,4)
bt0.sml_matrix(1,1,2)
sig.11_matrix(c(1,.2,.2,1),2,2)
sig.22_matrix(c(1,.1,.1,1),2,2)
sig.12_diag(2)
cc_diag(1)
u.big_diag(4)
u.sml_diag(2)
th0.big_0*cc%*%bt0.big%*%u.big
th0.sml_0*cc%*%bt0.sml%*%u.sml
mu_c(.2,.3,.4,.5)
rho_c(.7,.8,.9)
tbl_matrix(0,4,4)
dimnames(tbl)_list(c('.2','.3','.4','.5'),
c('2v','4v,r=.7','4v,r=.8','4v,r=.9'))
## Make their table, using the W test. The results differ because
## Muller and Peterson compute the noncentrality parameter slightly
## differently -- it's a little smaller than J&W's. Note that
## J&W use an exact power, whereas M&P use an approximation.
for (i in 1:4) {
tbl[i,1]_glhpwr(alpha,x,mu[i]*bt0.sml,sig.11,cc,u.sml,th0.sml)
for (j in 1:3) {
sigma_rbind(cbind(sig.11,rho[j]*sig.12),
cbind(rho[j]*sig.12,sig.22))
tbl[i,j+1]_glhpwr(alpha,x,mu[i]*bt0.big,
sigma,cc,u.big,th0.big) } }
cat('Table 5B.1, p. 218 of J&W, using the W test\n')
print(round(tbl,3))
##
SHAR_EOF
cat << \SHAR_EOF > 'pfnc.d'
.BG
.FN pfnc
.TL
Noncentral F Distribution
.DN
Using the first terms in the series expansion 26.6.20 in
Abramowitz and Stegun, compute the probability integral of
the noncentral F distribution.
.CS
pfnc(q, df1, df2, lm, iprec=c(6))
.RA
.AG q
vector of quantiles.
.AG df1
vector of numerator degrees of freedom.
.AG df2
vector of denominator degrees of freedom.
.AG lm
vector of noncentrality parameters.
.OA
.AG iprec
a parameter governing the precision of the answer; higher values
lead to greater precision.
The default is iprec=6.
For large lm, the true answer will exceed the returned value by
no more than 1-pnorm(iprec).
.RT
probability that a noncentral F variable with numerator
degrees of freedom df1, denominator degrees of freedom df2 and
noncentrality parameter lm is less than q.
.SH REFERENCES
Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical
Functions, p. 947, Dover, 1970.
.SH NOTE
Let n be the length of the longest of q, df1, df2 and lm.
If any of q, df1, df2 and lm is of length less than n, then all
values after the first in that vector are ignored.
Only the first element of iprec is used.
.SA
pf, pchisqnc.
.WR
SHAR_EOF
cat << \SHAR_EOF > 'pchisqnc.d'
***
.BG
.FN pchisqnc
.TL
Noncentral Chi-square Distribution
.DN
Using the first terms in the series expansion 26.4.25 in
Abramowitz and Stegun, compute the probability integral of the
noncentral chi-square distribution.
.CS
pchisqnc(q, df, lm, iprec=c(6))
.RA
.AG q
vector of quantiles.
.AG df
vector of degrees of freedom.
.AG lm
vector of noncentrality parameters.
.OA
.AG iprec
a parameter governing the precision of the answer; higher values
lead to greater precision.
The default is iprec=6.
For large lm, the true answer will exceed the returned value by
no more than 1-pnorm(iprec).
.RT
probability that a noncentral chi-square variable with degrees of
freedom df and noncentrality parameter lm is less than q.
.SH REFERENCES
Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical
Functions, p. 942, Dover, 1970.
.SH NOTE
Let n be the length of the longest of q, df and lm.
If any of q, df and lm is of length less than n, then all
values after the first in that vector are ignored.
Only the first element of iprec is used.
.SA
pfnc, pchisq.
.WR
SHAR_EOF
cat << \SHAR_EOF > 'glhpwr.d'
.BG
.FN glhpwr
.TL
Power Function for Tests of the General Linear Hypothesis
.DN
Use a noncentral F approximation to compute the power of a test of
the general linear hypothesis.
.CS
glhpwr(alpha, x, bt1, sigma, cc, u, th0, test=1, tol=1e-08)
.RA
.AG alpha
type I error probability of the test.
.AG x
N-by-q design matrix.
.AG bt1
q-by-p alternative hypothesis value of the regression coefficient.
.AG sigma
p-by-p variance matrix of a single observation.
.AG cc
a-by-q between-rows contrast matrix.
.AG u
p-by-b between-columns contrast matrix.
.AG th0
a-by-b null value of cc%*%beta%*%u.
.OA
.AG test
test for which power is desired. test=1 (the default) is the Wilks
lambda test; test=2 is the Hotelling-Lawley trace test; test=3 is
the Pillai-Bartlett trace test. See Muller and Peterson (1984) for
details.
.AG tol
tolerance for computing the rank of x (using qr()). The default is
1.e-8.
.RT
the approximate power of a test of the general linear hypothesis
cc%*%beta%*%u=th0 under the alternative beta=bt1.
The model is y(Nxp)=x(Nxq)%*%beta(qxp)+e(Nxp), where x is the design,
beta is the matrix of multivariate regression parameters, and e is
the error matrix, whose rows are assumed to be independent draws
from a multivariate normal with mean 0 and p-by-p variance matrix
sigma.
.SH REFERENCE
Muller, K.E. and Peterson, B.L. (1984). Practical methods for
computing power in testing the multivariate general linear
hypothesis. Computational Statistics and Data Analysis 2, 143--158.
.SH NOTE
If the computed approximate degrees of freedom are negative, a
power of 0 is returned; this may mean that the proposed design
has no degrees of freedom for error.
.SA
pfnc, sscomp.
.WR
SHAR_EOF
cat << \SHAR_EOF > 'sscomp.d'
.BG
.FN sscomp
.TL
Sample Size for Tests of the General Linear Hypothesis
.DN
Compute the minimum sample size required to have a specified power
for rejecting the null hypothesis in a test of the general linear
hypothesis.
.CS
sscomp(alpha, power, x, bt1, sigma, cc, u, th0, test=1, ninit=1)
.RA
.AG alpha
type I error probability of the test.
.AG power
desired power of the test.
.AG x
N0-by-q design matrix.
This matrix should represent the design for a sample of size 1.
For example, if there are to be two groups, x should contain two
rows, with the first row representing the design for group 1 and
the second row the design for group 2.
The sample size returned will then be the number per group.
.AG bt1
q-by-p alternative hypothesis value of the regression coefficient.
.AG sigma
p-by-p variance matrix of a single observation.
.AG cc
a-by-q between-rows contrast matrix.
.AG u
p-by-b between-columns contrast matrix.
.AG th0
a-by-b null value of cc%*%beta%*%u.
.OA
.AG test
test to use. test=1 (the default) is the Wilks lambda test; test=2
is the Hotelling-Lawley trace test; test=3 is the Pillai-Bartlett
trace test. See Muller and Peterson (1984) for details.
.AG ninit
number of replications of x to use as the starting sample size.
The sample size algorithm is simple -- it starts at ninit and keeps
adding or subtracting 1 (that is, 1 copy of x) until the desired
power is achieved.
The default is ninit=1.
.RT
list containing two elements -- n, the smallest number of copies of
x that guarantee the desired power for the test described, and
power, the power of the design with that sample size.
.SH REFERENCE
Muller, K.E. and Peterson, B.L. (1984). Practical methods for
computing power in testing the multivariate general linear
hypothesis. Computational Statistics and Data Analysis 2, 143--158.
.SH NOTE
The power function glhpwr is an approximation, so there may be
slight differences when compared with more exact methods.
This function increases (or decreases) n one unit at a time, so
well chosen starting values can save time.
Some experimentation will often be helpful.
.SA
pfnc, glhpwr.
.WR
SHAR_EOF
exit 0
# End of shell archive