powersim: Provide power estimates for multivariate abundance models

Description Usage Arguments Details Value Functions See Also Examples

View source: R/generic_functions.R

Description

powersim returns a power estimate for a cord object for a given sample size N and effect size of interest.

Usage

 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
powersim(
  object,
  coeffs,
  term,
  N = nrow(object$obj$data),
  coeffs0 = effect_null(object$obj, term),
  nsim = 999,
  test = "score",
  alpha = 0.05,
  newdata = NULL,
  n_replicate = NULL,
  ncores = detectCores() - 1,
  show.time = TRUE
)

## S3 method for class 'cord'
powersim(
  object,
  coeffs,
  term,
  N = nrow(object$obj$data),
  coeffs0 = effect_null(object$obj, term),
  nsim = 999,
  test = "score",
  alpha = 0.05,
  newdata = NULL,
  n_replicate = NULL,
  ncores = detectCores() - 1,
  show.time = TRUE
)

Arguments

object

objects of class cord, typically the result of a call to cord.

coeffs

Coefficient matrix for a manyglm object that characterises the size of effects to be simulated. See effect_alt for help in producing this matrix.

term

Name of predictor of interest in quotes.

N

Number of samples for power estimate. Defaults to the number of observations in the original sample.

coeffs0

Coefficient matrix under the null hypothesis. Defaults to being specified by effect_null.

nsim

Number of simulations for power estimate to be based upon. Defaults to 999.

test

Test statistic for power estimate to based upon. Defaults to "score", however "wald" is also allowed.

alpha

Type I error rate for power estimate, defaults to 0.05.

newdata

Data frame of the same size as the original data frame from the cord object (object$obj$data), that specifies a different design of interest.

n_replicate

Number of unique replicates of the original data frame. Defaults to NULL, overwrites N if specified.

ncores

Number of cores for parallel computing. Defaults to the total number of cores available on the machine minus 1.

show.time

Logical. Displays time elapsed. Defaults to TRUE.

Details

powersim takes a cord object, sample size N and coefficient matrix coeffs which specifies an effect size of interest and returns a power estimate.

The power estimate is obtained by first parsing the cord object into extend, nsim times with an effect size specified by coeffs. Next, the cord object is parsed into extend an additional nsim times with a null effect, which is defined by default by effect_null. This effectively simulates nsim manyglm models under both the null and alternative hypothesis.

For each simulated manyglm object a test statistic test is obtained. A critical test statistic is then obtained as the upper 1 - alpha quantile of simulated test statistics under the null hypothesis. Power is then estimated as the proportion of times the test statistics simulated under the alternative hypothesis exceed the critical test statistic under the null.

To improve computation time, simulations are computed in parallel using the "socket" approach, which by default uses all available cores minus 1 for clustering. Using 1 less than the number of available cores for your machine (detectCores()-1) is recommended to better avoid error relating to clustering or nodes.

Value

Power estimate result, and;

power

power.

Functions

See Also

effect_alt, effect_null, extend

Examples

 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
library(ecoCopula)
library(mvabund)
data(spider)
spiddat = mvabund(spider$abund)
X = data.frame(spider$x)

# Specify increasers and decreasers
increasers = c("Alopacce", "Arctlute", "Arctperi", "Pardnigr", "Pardpull")
decreasers = c("Alopcune", "Alopfabr", "Zoraspin")

# Find power for continuous predictor at effect_size=1.5
fit.glm = manyglm(spiddat~bare.sand, family="negative.binomial", data=X)
effect_mat = effect_alt(fit.glm, effect_size=1.5,
       increasers, decreasers, term="bare.sand")
fit.cord = cord(fit.glm)
powersim(fit.cord, coeffs=effect_mat, term="bare.sand", nsim=99, ncores=2)

# Find power for categorical predictor with 4 levels at effect_size=1.5
X$Treatment = rep(c("A","B","C","D"),each=7)
fit_factors.glm = manyglm(spiddat~Treatment, family="negative.binomial", data=X)
effect_mat = effect_alt(fit_factors.glm, effect_size=1.5,
       increasers, decreasers, term="Treatment")
fit_factors.cord = cord(fit_factors.glm)
powersim(fit_factors.cord, coeffs=effect_mat, term="Treatment", nsim=99, ncores=2)

# Change effect size parameterisation
effect_mat = effect_alt(fit_factors.glm, effect_size=1.5,
                         increasers, decreasers, term="Treatment",
                         K=c(3,1,2))
powersim(fit_factors.cord, coeffs=effect_mat, term="Treatment", nsim=99, ncores=2)

ecopower documentation built on Sept. 6, 2021, 9:09 a.m.