power.groups {emon}

Power for comparing mean of two groups

Description

Calculates the power by simulation for comparing the mean of two groups of independent observations.

Usage

power.groups(change, change.type="M", n1, n2, pars1, pars2,
           distribution, test, alternative="two", alpha=0.05, nsims=1000, nreps=999)

Arguments

change

Mean of second group minus mean of first group (i.e. mu2-mu1) or percentage change in mu1 to create mu2 (depending on value of change.type).

change.type

Whether the parameter change represents an additive ("A") or percentage ("M") change.

n1

Vector of sample sizes for group 1. Must be of same dimension as n2.

n2

Vector of sample sizes for group 2. Must be of same dimension as n1.

pars1

Parameters for the treatment data. If distribution="Normal", pars1[1] contains the mean for group 1 and pars1[2] contains the standard deviation. If distribution="Poisson", pars1[1] contains the mean for group 1. If distribution="Lognormal", pars1[1] contains the group 1 mean of the natural log data and pars1[2] contains the standard deviation of the log data. If distribution="Negbin", pars1[1] contains the mean and pars1[2] contains the group 1 size parameter.

pars2

pars2[1] is the standard deviation for group 2 if distribution="Normal". If distribution="Lognormal", pars2[1] is the standard deviation of the log data for group 2. For distribution="Negbin", pars2[1] gives the group 2 size.

distribution

The statistical distribution for the two groups. Can be either: "Normal", "Poisson", "Lognormal" or "Negbin".

test

The statistical test used to compare the group means. If test="NP" then the test will be a non-parametric randomisation test, in the spirit of Manly (1997), using the function permute.groups. If test="P", then parametric tests are made to compare the group means. If distribution="Normal", a two sample t-test is carried out. If the standard deviations (defined by pars1[2] and pars2[1] are equal, then the t-test calculates the usual pooled standard deviation. However, if the standard deviations are not equal then the default method used by t.test is adopted.

When distribution="Lognormal", natural logs are taken of the simulated data and a t-test used in a similar way as to when distribution="Lognormal". When distribution="Poisson", the difference in deviances between the null model and the model with group membership fitted as factor is compared againt a chi-squared distributiuon on 1 degree of freedom. A similar (but not quite) method is used for when distribution="Negbin". The Generalised Linear Model function glm.nb for the Negative Binomial distribution is used. The p-value for comparing the two groups is taken from the analysis of deviance table after the model with group membership is fitted as a factor. This p-value, however, uses the same value for the size parameter, as estimated from the null model, for group member deviance. This seems to be the correct thing to do as estimating separate size parameters for the two models mucks up the nesting of the models.

alternative

A character string specifying the alternative hypothesis, must be one of "two.sided" (default), "greater" or "less". You can specify just the initial letter. As an example, "less" means that the alternative is that the mean of the first group is less than the mean of the second.

alpha

The type 1 error for assessing statistical significance (default is 0.05) in the power simulations.

nsims

Number of repeat simulations to estimate power (default is 1000).

nreps

Number of repeat permutations for randomisation test (default is 999).

Details

The Negative Binomial distribution option allows the user to specify the size parameter for both groups 1 and 2. One possibility is to keep the size the same for both groupss. However, because the mean is different between the groups and because the variance V = mu+mu^2/size, this means that V will be different for the group 1 and group 2 distributions. If you want to keep the variance the same, you can use the function size2.samevar.

Several powers can be calculated per call to this function by specifying more than one values for the sample sizes n1 and n2.

Value

The power is returned. This is the proportion of the nreps simulations that returned a p-value less than the type1 error.

Author(s)

Jon Barry: [e-mail]

References

Manly BFJ (1997) Randomization, bootstrap and monte carlo methods in biology: 2nd edition. Chapman and Hall, London, 399 pages.

See Also

permute.groups, glm.nb, size2.samevar

Examples

library(MASS)

# In practice, \code{nsims} would be set to at least 1000

power.groups(change=2.5, change.type="A", n1=20, n2=20, pars1=c(10,2),
       pars2=2, test='P', distribution="Normal", nsims=50)

power.groups(change=2.5, change.type="A", n1=seq(5,25,5), n2=seq(5,25,5), pars1=c(10,2),
       pars2=2, test='P', distribution="Normal", nsims=50)

power.groups(change=25, change.type="M", n1=20, n2=20, pars1=10,
       test='P', distribution="Poisson", nsims=50)

power.groups(change=4, change.type="A", n1=20, n2=20, pars1=c(1,2),
       pars2=2, test='P', distribution="Lognormal", nsims=50)


# Keeping size constant
power.groups(change=100, change.type="M", n1=20, n2=20, pars1=c(5,2),
       pars2=2, test='P', distribution="Negbin", nsims=50)

# Keeping variance constant
s2 = size2.samevar(mu1=5, mu2=10, s1=2)   # 13.333
power.groups(change=100, change.type="M", n1=20, n2=20, pars1=c(5,2),
       pars2=s2, test='P', distribution="Negbin", nsims=50)