Sample size analysis example

Crayweed restoration project

Researchers within the Operation Crayweed Restoration Project in Sydney are restoring the locally extinct macro-algae Phyllospora comosa (“crayweed”: see Coleman et al., 2008; Campbell et al., 2014; Vergés et al., 2020 and are interested in the effect of this restoration on associated ecological communities (Marzinelli et al., 2014, 2016; Wood et al., 2019). Pilot data have already been collected, where the abundance of fish species in nine open ocean sites have been recorded in 2015. We are interested in observing if there is a change in mean fish abundance between control sites (those in Sydney without crayweed) and restored sites (similar sites where crayweed has recently been transplanted). There are plans to collect more data in the future, however there is an upper bound of approximately 24 possible spatially independent restored or control coastal bays/sites within the Sydney region and surroundings.

In this example we attempt to answer the following experimental design questions:

1. How many sites are required to likely detect \(20\%\) differences between treatments?
2. Under the maximum independent sampling design of 24 sites, what are the size of effects that are likely to be detected?

#to install latest version of ecopower on github:
#devtools::install_github("BenMaslen/ecopower")

#load libraries
library(ecopower)
library(mvabund)
library(ecoCopula)
library(ggplot2)
library(RColorBrewer)

#read in data
data("fish")

#create mvabund abundance matrix
fish_data <- fish
abund <- mvabund(fish_data[, 1:34])

#set Site.Type as a factor
fish_data$Site.Type <- factor(fish_data$Site.Type, levels = c("control", "restored", "reference"))

#fit the model
fit       <- manyglm(abund ~ Site.Type, family = "negative.binomial", data = fish_data)
fit.cord  <- cord(fit)

#check model assumptions
plot(fit)

#specify 'increaser' and 'decreaser' species:
increasers <- c("Aplodactylus.lophodon", "Atypichthys.strigatus", "Cheilodactylus.fuscus",
                "Olisthops.cyanomelas", "Pictilabrius.laticlavius")
decreasers <- c("Abudefduf.sp", "Acanthurus.nigrofuscus", "Chromis.hypsilepis",
                "Naso.unicornis", "Parma.microlepis", "Parupeneus.signatus",
                "Pempheris.compressa", "Scorpis.lineolatus", "Trachinops.taeniatus")

#specify effect size of interest
coeff.alt <- effect_alt(fit, effect_size = 1.2, increasers, decreasers, term = "Site.Type")

fish_data_sim            <- fish_data
fish_data_sim$Site.Type[fish_data_sim$Site.Type=="reference"] <- "restored"
fish_data_sim$Site.Type  <- factor(fish_data_sim$Site.Type)

powersim(fit.cord, N=100, coeff.alt, term="Site.Type", nsim=50, ncrit=49, newdata = fish_data_sim, ncores=2)
#> Time elapsed: 0 hr 0 min 20 sec
#> Power: 0.48

Produce a single power curve:

In order to answer our first experimental design question, we can estimate power across a range of sample sizes for an effect size of \(20\%\) differences in mean abundance. By plotting the results, we can observe the sample size required to reliably detect (with a conventional power target of \(80\%\)) the effect of interest.

Users could also plot multiple power curves to get an idea of the variability in the power estimates (as in Maslen et al. 2023).

Note: The below code takes ~10minutes-1hour to run depending on the speed of your computer and the number of logical processors you have. As such, this code has not been evaluated by default. Instead, the resultant figure has been added for reference. If you want to run this code and reduce computation time, you can reduce nsim or ncrit, however this will also make power estimates more variable.

First we estimate power at each sample size of interest:

#specify sample sizes to simulate under
sample_sizes    <- c(10, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300)
power_estimates <- rep(NA, length(sample_sizes))

#loop through sample sizes and estimate power at each step
for (i in 1:length(sample_sizes)){
  power_estimates[i] <- powersim(fit.cord, N = sample_sizes[i], coeff.alt, term = "Site.Type",
                                 nsim = 1000, ncrit = 4999, newdata = fish_data_sim)
}

Next, let’s save the results as a data frame and produce a plot against the conventional power target of \(80\%\):

#store the results in a data.frame
powercurve_dat <-  data.frame(sample_sizes = sample_sizes, Power = unlist(power_estimates))

#plot power curves
ggplot(powercurve_dat, aes(x = sample_sizes, y = Power)) +
  geom_line(size = 1.06) + theme_bw() + geom_hline(yintercept = 0.8, linetype = "dashed", color = "red", size=1) +
  xlab("N") + ylab("Power") + scale_x_continuous(breaks = sample_sizes) +
  scale_y_continuous(breaks = c(0.2, 0.4, 0.6, 0.8, 1.0))

This suggests that we require a sample size of at-least \(N \approx 150\) in order to likely detect \(20\%\) changes in mean abundances, using a conventional power target of \(80\%\).

Thus we have answered our first experimental design question!

Multiple effect size power curve:

In order to answer the second experimental design question, we can estimate and then plot multiple effect size curves for different effect size specifications (let’s do this for effect sizes of \(10\%, 20\%, \ldots, 100\%\) changes in mean abundance by specifying effect_size\(=1.1, 1.2, \ldots, 2\). We will also do this across sample sizes that we can sample using a balanced sampling design (from \(N=4,6,\ldots,24\)).

As we are simulating from copula models with small sample sizes, we will also increase n.samp (number of simulations for importance sampling in copula modelling) and decrease nlv to 1 (number of latent factors in the copula) as recommended in the cord function description.

Note: The below code takes ~30minutes-2hours to run depending on the speed of your computer and the number of logical processors you have. As such, this code has not been evaluated by default. Instead, the resultant figure has been added for reference. If you want to run this code and reduce computation time, you can reduce nsim or ncrit, however this will also make power estimates more variable.

First we estimate power at each sample size and effect size of interest:

#specify effect sizes and sample sizes to loop through
sample_sizes <- c(3:12)*2
effect_sizes <- c(1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2)
mult_power_estimates <- rep(NA,length(sample_sizes)*length(effect_sizes))

#loop through sample sizes and effect sizes
for (j in 1:length(effect_sizes)){
  coeff.alt <- effect_alt(fit, effect_size = effect_sizes[j], increasers, decreasers, term = "Site.Type")
  for (i in 1:length(sample_sizes)){
    try(mult_power_estimates[10*(j-1)+i] <- powersim(fit.cord, N = sample_sizes[i], coeff.alt, term = "Site.Type", nsim = 1000,
                                                     ncrit = 4999, newdata = fish_data_sim, nlv = 1, n.samp = 500))
  }
}

Next, let’s save the results as a data frame and produce a plot of the resulting power estimates:

#store the results in a data.frame
mult_powercurve_dat <- data.frame(sample_sizes = rep(sample_sizes, times = length(effect_sizes)),
                                  power_estimates = unlist(mult_power_estimates),
                                  Perc_Change = factor(rep(c("10%", "20%", "30%", "40%", "50%", "60%",
                                                             "70%", "80%", "90%", "100%"), each = 10),
                                                     levels = c("100%", "90%", "80%", "70%", "60%", "50%",
                                                              "40%", "30%", "20%", "10%")),
                                  effect_sizes = rep(effect_sizes, each = length(sample_sizes)))



#plot multiple power curve
ggplot(mult_powercurve_dat, aes(x = sample_sizes, y = power_estimates, colour = Perc_Change)) +
  geom_line(linewidth = 1.05) + theme_bw() + geom_hline(yintercept = 0.8, linetype = "dashed", color = "red", size = 1) +
  xlab("N") + ylab("Power") + scale_x_continuous(breaks = sample_sizes) +
  scale_y_continuous(breaks = c(0.2, 0.4, 0.6, 0.8, 1.0)) + scale_color_brewer(palette = "Paired", direction = -1, name = "% change")

In the above plot we can observe that under the maximum balanced sample size of \(N=24\) sites, the smallest effect size to be observed under a conventional power target of \(80\%\) is effect_size\(\approx 1.7\), or equivalently \(70 \%\) changes in the mean abundances of fish species. The pilot survey included \(N_{\text{pilot}}=9\) sites, which would not be able to reliably detect any of the simulated effect sizes.

Thus we have answered our second experimental design question!

Note that we could also have just produced a single power curve here for a sample size of \(N=24\) across our effect sizes of interest, which would also have saved a lot of time (and still answered our second experimental design question). The multiple effect size plot that we produced however will also give us insight into the size of effects we can detect across the entire range of balanced sample sizes.

Sensitivity analysis:

As mentioned in Maslen et al. 2023 it is also a good idea to check how sensitive our prior assumptions are in generating our effect size of interest and estimating power.

Let’s produce some power curves where we play around with the most abundant species and alter the overdispersion parameter, to see how sensitive our power estimates are based on these parameters.

First let’s create a new cord object and halve the overdispersion by doubling overdispersion parameter theta.

#new cord model
fit.cord_doub_theta <- cord(fit)

#change theta to be doubled
fit.cord_doub_theta$obj$theta <- fit.cord$obj$theta*2

Secondly, let’s try investigate what happens when we change the most abundant species Atypichthys.strigatus from an ‘increaser’ to a ‘no effect’ taxa, by removing it from the ‘increaser’ species list.

increasers_no_abund <- c("Aplodactylus.lophodon", "Cheilodactylus.fuscus",
                         "Olisthops.cyanomelas", "Pictilabrius.laticlavius" )

Now, let’s re specify the effect size of interest as \(20\%\) changes in mean abundance via setting effect_size\(=1.2\). We will also specify another coefficient matrix without the most abundant species in it.

coeff.alt_no_abund  <- effect_alt(fit, effect_size = 1.2, increasers_no_abund, decreasers, term = "Site.Type")
coeff.alt           <- effect_alt(fit, effect_size = 1.2, increasers, decreasers, term = "Site.Type")

Now we can estimate power with these changes.

Note: The below code takes ~20minutes-2hours to run depending on the speed of your computer and the number of logical processors you have. As such, this code has not been evaluated by default. Instead, the resultant figure has been added for reference. If you want to run this code and reduce computation time, you can reduce nsim or ncrit, however this will also make power estimates more variable.

sample_sizes <- c(10, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300)
power_estimatesa <- rep(NA,length(sample_sizes))
power_estimatesb <- rep(NA,length(sample_sizes))

for (i in 1:length(sample_sizes)){
  power_estimatesa[i] <- powersim(fit.cord, N = sample_sizes[i], coeff.alt_no_abund, term = "Site.Type",
                                  nsim = 1000, ncrit = 4999, newdata = fish_data_sim)
  power_estimatesb[i] <- powersim(fit.cord_doub_theta, N = sample_sizes[i], coeff.alt, term = "Site.Type",
                                  nsim = 1000, ncrit = 4999, newdata = fish_data_sim)
}

Finally we can store the results in a data frame with our original power curve estimates to compare to, and plot the results.

powercurve_experiment_dat <- data.frame(sample_sizes = rep(sample_sizes, 3), 
                                        Method = factor(rep(c("Original", "No abund", "Double theta"), each = 13), 
                                                        levels = c("Double theta", "Original", "No abund")),
                                              Power = c(unlist(c(power_estimates, power_estimatesa, power_estimatesb))))


#plot power curve
ggplot(powercurve_experiment_dat, aes(x = sample_sizes, y = Power, colour = Method)) +
  geom_line(size = 1.1) + theme_bw() + geom_hline(yintercept = 0.8, linetype = "dashed", color = "red", size = 1) +
  xlab("N") + ylab("Power") + scale_x_continuous(breaks = sample_sizes) + theme(legend.text.align = 0) +
  scale_y_continuous(breaks = c(0.2, 0.4, 0.6, 0.8, 1.0)) +
  scale_colour_manual(
    values = c("#D55E00", "#56B4E9", "#009E73"),
    labels = c("Halve overdispersion", "Original", expression(paste("No " , italic(A.strigatus))))
  )

Here we can observe that halving the amount of overdispersion in our simulations increased power by as much as \(25\%\) and when we changed our most abundant species Atypichthys.strigatus from an ‘increaser’ to a ‘no effect’ taxa, power would in some cases decrease by as much as \(20\%\).

This demonstrates how sensitive our prior assumptions can be on resulting power, and in particular for abundant species.

For further details behind the methodology please go to Maslen et al. 2023, or for examples with more complicated effect size specifications go to the examples within ecopower.

If you use this package for sample size analysis in an associated paper or otherwise, please cite Maslen et al. 2023.