BACI designs are commonly used in environmental monitoring. They are relevant where you want to measure the effect of an impact (e.g. the effect on benthic ecology of dredging in an area). Observations for treatment and control areas are measured BEFORE and after the impact. This function allows you to examine the power of particular BACI designs to detect differences between the control and the treatment.
1 2  power.BACI(change, change.type="M", nt, nc, parst, parsc,
distribution, test="P", alpha=0.05, nsims=1000)

change 
AFTER treatment mean minus BEFORE treatment mean or percentage change of
AFTER treatment mean relative to BEFORE treatment mean (depending on value of

change.type 
Whether the parameter 
nt 
Vector of sample sizes for treatment group. Must be of same dimension as 
nc 
Vector of sample sizes for control group. Must be of same dimension as 
parst 
Parameters for the treatment data. If 
parsc 
Parameters for the control data. If 
distribution 
The statistical distribution for the two groups. Can be either: 
test 
The statistical test used to compare the interaction between the control and treatment
means at the BEFORE and AFTER sampling occasions. If If When 
alpha 
If the pvalue for the interaction is less than 
nsims 
Number of repeat simulations used to calculate the power. Default is 1000. 
BACI designs are relevant where you want to measure the effect of an impact (e.g. the effect on benthic ecology of dredging in an area). You take a number of samples both in (treatment) and outside (control) the affected area BEFORE the impact. Those in the control area should be as similar to the treatment area as possible in terms of benthic ecology. You then sample the areas again AFTER the impact has taken place. If there is an interaction for the 2x2 crossed design then there is an effect of the impact. That is, if the control area has changed differently to the treatment area.
This function allows you to examine the power of particular BACI designs. The distribution of the measure being used can be Normal, Poisson, Negative Binomial or Lognormal.
It is also assumed that the sample sizes before and after the impact are the same, although the sample size in the treatment area can be different to the control area. Thus, if 10 and 8 samples are taken in the treatment and control areas before the impact, then 10 and 8 samples are assumed to be taken after the impact.
For the Normal, Poisson and Negative Binomial distributions, the parameter change is simply the percentage or
additive change of the treatment mean from
the BEFORE to the AFTER sampling occasions. For the Lognormal distribution, the percentage change is relative to the
BEFORE treatment mean on the nonlog scale. The BEFORE treatment mean is estimated from the mean of the
log data (parst[1]
), the standard deviation of the log data (parst[2]
) and the proposed sampling size
nt
.
The estimator used is the one proposed by Shen (2006), which did better in terms of mean squared error than both the
sample mean on the nonlog scale and the maximum likelihood estimators. This is given by
mean.before = exp(parst[1] + (nt1)*ss / (2*(nt+4)*(nt1) + 3*ss)), where ss = (nt1)*parst[2]**2
.
The Negative Binomial distribution option (parst[3]
allows the user to specify the size parameter of the
AFTER treatment distribution. One possibility is to keep the size the same for both the BEFORE and AFTER
distributions. However, because the mean changes and because the variance V = mu+mu^2/size, this means that V
will be different for the BEFORE and AFTER 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 nt
and nc
.
power 
The estimated power for the design. 
before.mean 
The treatment mean used for the before sampling. Only really of interest if

Jon Barry (email Jon.Barry@cefas.co.uk)
Shen H, Brown LD and Zhi H (2006) Efficient estimation of lognormal means with application to pharmacokinetic data. Statistics in Medicine, 25, 3023 to 3038.
power.trend
, power.groups
, size2.samevar
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36  # Data is richness (number of species) and abundance from grab samples from
# the Dogger Bank, UK. In practice, \code{nsims} would be set to at least 1000.
rich = c(15,16,37,12,15,5,13,16,17,34,23,20,22,30,85,55,13,19,30,41,22,8,43,10,38,24,17,
23,17,17,24,33,30,18,26,18,12,50,19,21,35)
abun = c(50,91,140,21,25,8,28,37,30,90,56,50,40,83,964,180,21,60,81,138,67,17,250,63,152,
68,42,69,57,67,74,96,75,44,61,49,62,281,55,50,198)
par(mfrow=c(2,2))
hist(rich)
hist(abun)
hist(sqrt(rich))
hist(log(abun))
ssize = seq(10, 50, 10)
parsc.rich = mean(rich); parst.rich = mean(rich)
parsc.abun = rep(0,2); parst.abun = parsc.abun
parst.abun[1] = mean(log(abun)); parst.abun[2] = sd(log(abun))
parsc.abun = parst.abun
power.rich = rep(0, length(ssize))
power.abun = rep(0, length(ssize))
power.rich = power.BACI(change=35, change.type="M", nt=ssize, nc=ssize,
parst=parst.rich, parsc=parsc.rich,
distribution="Poisson", test="P", nsims=50)$power
power.abun = power.BACI(change=35, change.type="M", nt=ssize, nc=ssize,
parst=parst.abun, parsc=parsc.abun, distribution="Lognormal",
test="P", nsims=50)$power
par(mfrow=c(1,1))
plot(ssize, power.rich, ylim=c(0,1), ylab="Power", xlab="Sample size", type="l")
lines(ssize, power.abun, lty=2, col=2)
legend("bottomright", legend=c("Richness power", "Abundance power"), lty=c(1,2),
col=c(1,2))
title("BACI power plots")

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.