| Title: | Power Estimates and Equivalence Testing for Multivariate Data |
|---|---|
| Description: | Estimates power by simulation for multivariate abundance data to be used for sample size estimates. Multivariate equivalence testing by simulation from a Gaussian copula model. The package also provides functions for parameterising multivariate effect sizes and simulating multivariate abundance data jointly. The discrete Gaussian copula approach is described in Popovic et al. (2018) <doi:10.1016/j.jmva.2017.12.002>. |
| Authors: | Ben Maslen [aut], Michelle Lim [aut, cre] |
| Maintainer: | Michelle Lim <[email protected]> |
| License: | LGPL (>= 2.1) |
| Version: | 0.2.0 |
| Built: | 2024-12-16 07:01:06 UTC |
| Source: | CRAN |
Dataset of fish abundances recorded at crayweed reference and restored sites.
data(crayweed)data(crayweed)
An object of class "list" containing:
A matrix with 27 observations of abundance of 34 fish species.
A data frame with treatment and time variables.
The matrix abund has the following species abundances:
Abudefduf.sp.
Acanthopagrus.australis
Acanthurus.nigrofuscus
Achoerodus.viridis
Aplodactylus.lophodon
Atypichthys.strigatus
Cheilodactylus.fuscus
Chromis.hypsilepis
Crinodus.lophodon
Girella.elevata
Girella.tricuspidata
Hypoplectrodes.maccullochi
Needle.fish.unidentified
Notolabrus.gymnogenis
Odax.cyanomelas
Olisthops.cyanomelas
Ophthalmolepis.lineolatus
Parma.microlepis
Parma.unifasciata
Pempheris.compressa
Pempheris.multiradiata
Pictilabrus.laticlavius
Prionurus.microlepidotus
Pseudocaranx.dentex
Pseudojuloides.elongatus
Pseudolabrus.gymnogenis
Sardinops.neopilchardus
Scorpis.lineolatus
Seriola.lalandi
Sphyraena.obtusata
Tetractenos.hamiltoni
Trachinops.taeniatus
Trachurus.novaezelandiae
Upeneichthyes.lineatus
The data frame X has the following variables:
treatment - reference / restored
time - sample period with seven time points
Data attributed to the crayweed restoration project (http://www.operationcrayweed.com/).
data(crayweed) head(crayweed$abund) head(crayweed$X)data(crayweed) head(crayweed$abund) head(crayweed$X)
effect_alt returns a coefficient matrix to be parsed to extend, powersim
and equivtest to specify an effect size of interest.
## S3 method for class 'manyglm' effect_alt(object, effect_size, increasers, decreasers, term, K = NULL) effect_alt(object, effect_size, increasers, decreasers, term, K = NULL)## S3 method for class 'manyglm' effect_alt(object, effect_size, increasers, decreasers, term, K = NULL) effect_alt(object, effect_size, increasers, decreasers, term, K = NULL)
object |
objects of class |
effect_size |
An effect size of interest, see details for interpretation. |
increasers |
A vector list of responses which increase relative to the control group/intercept. |
decreasers |
A vector list of responses which decrease relative to the control group/intercept. |
term |
Name of predictor of interest in quotes. |
K |
A vector of length |
effect_alt helps users to create interpretable multivariate effect sizes to be parsed into extend,
powersim and equivtest, so that researchers can investigate the relationship between effect size, power and
sample size in a complicated multivariate abundance setting.
effect_alt creates an effect of size log(effect_size) for a predictor of interest (term), for responses
who have been specified to increase (increasers) and -log(effect_size) for responses who have been specified
to decrease (decreasers). Responses that have not been specified in the increasers or decreasers vectors
are specified to have no effect with a coefficient of 0. The effect has been logged to make the effect size
interpretable within the coefficient matrix.
For poisson regression family=poisson() and negative binomial regression family="negative.binomial" the effect
size is interpreted for a categorical variable as the multiplicative change in mean abundance in the treatment
group relative to the control group, whilst for a continuous variable it is interpreted as the multiplicative
change in abundance for a 1 unit increase in the predictor of interest.
For logit regression family=binomial("logit") the effect size is interpreted as an odds ratio. For a categorical
variable this is the change in odds of obtaining outcome 1 when being in the treatment group relative to the
control group. Whilst for continuous variables, this is interpreted as the change in odds of obtaining outcome 1
with a 1 unit increase in the predictor of interest.
For cloglog regression family=binomial("cloglog") the effect size is interpreted similarly to poisson and
negative binomial regression. For a categorical variable it is interpreted as the multiplicative change in
the mean of the underlying count in the treatment group relative to the control. Whilst for a continuous
variable it is interpreted as the multiplicative change in the mean of the underlying count for a 1 unit
increase in the predictor of interest.
For categorical variables, the intercept is also changed to be the group mean intercept by taking the intercept of a model without the categorical predictor of interest. This is done to avoid messy comparisons of null control groups.
For categorical variables with more than two levels, effect size is changed to effect_size^K[i] where K defaults
to be c(1,2,...,nlevels - 1), where nlevels are the number of levels of the categorical variable and is
specified along the order of the levels. To change this, specify a vector K with length of nlevels - 1.
To change the control group, this must be done prior to specifying the manyglm object
using relevel (which can also change the order of the levels).
Note that if the predictor of interest is a categorical variable it must be classed either as a factor or character otherwise results may be misleading.
A coefficient matrix with the specified effect size.
effect_alt(): Specify multivariate effect sizes
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") # Obtain an effect matrix of effect_size=3 spid.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) effect_mat = effect_alt(spid.glm, effect_size=3, increasers, decreasers, term="soil.dry") # Obtain an effect matrix of effect_size=1.5 X$Treatment = rep(c("A","B","C","D"),each=7) spid.glm = manyglm(spiddat~Treatment, family="negative.binomial", data=X) effect_mat = effect_alt(spid.glm, effect_size=1.5, increasers, decreasers, term="Treatment") # Change effect size parameterisation effect_mat = effect_alt(spid.glm, effect_size=1.5, increasers, decreasers, term="Treatment", K=c(3,1,2))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") # Obtain an effect matrix of effect_size=3 spid.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) effect_mat = effect_alt(spid.glm, effect_size=3, increasers, decreasers, term="soil.dry") # Obtain an effect matrix of effect_size=1.5 X$Treatment = rep(c("A","B","C","D"),each=7) spid.glm = manyglm(spiddat~Treatment, family="negative.binomial", data=X) effect_mat = effect_alt(spid.glm, effect_size=1.5, increasers, decreasers, term="Treatment") # Change effect size parameterisation effect_mat = effect_alt(spid.glm, effect_size=1.5, increasers, decreasers, term="Treatment", K=c(3,1,2))
effect_null returns a coefficient matrix to be parsed to powersim by default
to specify a null effect.
## S3 method for class 'manyglm' effect_null(object, term) effect_null(object, term)## S3 method for class 'manyglm' effect_null(object, term) effect_null(object, term)
object |
objects of class |
term |
Name of predictor of interest in quotes. |
effect_null produces a coefficient matrix with a null effect that is specified by setting the parameter
estimates of a predictor of interest term to 0. This function is used by default in powersim.
Note that intercept values are parameterised as in effect_alt.
A coefficient matrix with the null effect.
effect_null(): Specify null effects for multivariate abundance data
library(mvabund) data(spider) spiddat = mvabund(spider$abund) X = data.frame(spider$x) # Find null effect size for continuous predictor spid.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) coeffs0 = effect_null(spid.glm, term="soil.dry")library(mvabund) data(spider) spiddat = mvabund(spider$abund) X = data.frame(spider$x) # Find null effect size for continuous predictor spid.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) coeffs0 = effect_null(spid.glm, term="soil.dry")
equivtest takes in a copula model fitted to data and a matrix of effect sizes to execute a
a multivariate equivalence test.
## S3 method for class 'cord' equivtest( object, coeffs, term = NULL, object0 = NULL, stats = NULL, test = "LR", nsim = 999, ncores = detectCores() - 1, show.time = TRUE ) equivtest( object, coeffs, term = NULL, object0 = NULL, stats = NULL, test = "LR", nsim = 999, ncores = detectCores() - 1, show.time = TRUE )## S3 method for class 'cord' equivtest( object, coeffs, term = NULL, object0 = NULL, stats = NULL, test = "LR", nsim = 999, ncores = detectCores() - 1, show.time = TRUE ) equivtest( object, coeffs, term = NULL, object0 = NULL, stats = NULL, test = "LR", nsim = 999, ncores = detectCores() - 1, show.time = TRUE )
object |
objects of class |
coeffs |
Coefficient matrix for a |
term |
Name of predictor of interest in quotes. Defaults to |
object0 |
object of class |
stats |
Statistics simulated under the null hypothesis. Optional, defaults to |
test |
Test statistic for computing p-value. Defaults to |
nsim |
Number of simulations for p-value estimate to be based upon. Defaults to |
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 |
equivtest takes a cord object and a coefficient matrix coeffs which specifies an effect size of
interest to perform an equivalence test.
First, marginal parameters of the data are obtained from a manyglm object. Next, a copula model is fitted
using cord to estimate the factor analytic covariance structure of the data. The cord function uses two
factors by default. The p-value is then obtained by parsing the cord object into extend,
nsim times with an effect size specified by coeffs.
The test statistics are simulated under the hypothesis that the effect size equals a certain threshold. The p-value is computed as the proportion of times the simulated test statistics are less than the observed statistic. Equivalence is declared if the estimated effect is less than the threshold.
equivtest can handle any user-defined null hypothesis, so only the fitted null model (object0) or the predictor of
interest (term) needs to be specified. If both object0 and term are NULL, equivtest will
automatically set the predictor of interest as the last term in the fitted object model or drop the only term in the model
to obtain the intercept model.
Simulations are computed in parallel using the "socket" approach, which uses all available cores minus 1 for clustering
to improve computation efficiency. Using 1 less than the number of available cores for your
machine (detectCores()-1) is recommended to leave one core available for other computer processes.
Equivalence test results, and;
p |
p-value; |
stat_obs |
observed statistic; |
stats |
simulated statistics. |
equivtest(): Multivariate equivalence testing
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") # Equivalence test for continuous predictor at effect_size=1.5 fit.glm = manyglm(spiddat~bare.sand, family="negative.binomial", data=X) threshold = effect_alt(fit.glm, effect_size=1.5, increasers, decreasers, term="bare.sand") fit.cord = cord(fit.glm) equivtest(fit.cord, coeffs=threshold, term="bare.sand", nsim=99, ncores=2) # Equivalence test 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) threshold = effect_alt(fit_factors.glm, effect_size=1.5, increasers, decreasers, term="Treatment") fit_factors.cord = cord(fit_factors.glm) equivtest(fit_factors.cord, coeffs=threshold, term="Treatment", nsim=99, ncores=2) # Specify object0 object0.glm = manyglm(spiddat~1, family="negative.binomial") object0.cord = cord(object0.glm) equivtest(fit_factors.cord, coeffs=threshold, object0=object0.cord, nsim=99, ncores=2)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") # Equivalence test for continuous predictor at effect_size=1.5 fit.glm = manyglm(spiddat~bare.sand, family="negative.binomial", data=X) threshold = effect_alt(fit.glm, effect_size=1.5, increasers, decreasers, term="bare.sand") fit.cord = cord(fit.glm) equivtest(fit.cord, coeffs=threshold, term="bare.sand", nsim=99, ncores=2) # Equivalence test 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) threshold = effect_alt(fit_factors.glm, effect_size=1.5, increasers, decreasers, term="Treatment") fit_factors.cord = cord(fit_factors.glm) equivtest(fit_factors.cord, coeffs=threshold, term="Treatment", nsim=99, ncores=2) # Specify object0 object0.glm = manyglm(spiddat~1, family="negative.binomial") object0.cord = cord(object0.glm) equivtest(fit_factors.cord, coeffs=threshold, object0=object0.cord, nsim=99, ncores=2)
extend returns a simulated response matrix or a manyglm object with N observations
and simulated response matrix that utilises the existing correlation structure of the data.
## S3 method for class 'cord' extend( object, N = nrow(object$obj$data), coeffs = coef(object$obj), newdata = NULL, n_replicate = NULL, do.fit = FALSE, seed = NULL ) extend( object, N = nrow(object$obj$data), coeffs = coef(object$obj), newdata = NULL, n_replicate = NULL, do.fit = FALSE, seed = NULL )## S3 method for class 'cord' extend( object, N = nrow(object$obj$data), coeffs = coef(object$obj), newdata = NULL, n_replicate = NULL, do.fit = FALSE, seed = NULL ) extend( object, N = nrow(object$obj$data), coeffs = coef(object$obj), newdata = NULL, n_replicate = NULL, do.fit = FALSE, seed = NULL )
object |
objects of class |
N |
Number of samples to be extended. Defaults to the number of observations in the original sample. |
coeffs |
Coefficient matrix for a |
newdata |
Data frame of same size as the original X covariates from the fitted |
n_replicate |
Number of unique replicates of the original data frame. Defaults to |
do.fit |
Logical. If |
seed |
Random number seed, defaults to a random seed number. |
extend takes a cord object and returns a new simulated response matrix or an "extended" manyglm object
with N observations and the new simulated response matrix. Response abundances are simulated through a Gaussian
copula model that utilises a coefficient matrix coeffs, the specified cord model and the joint
correlation structure exhibited between the response variables. To help with the specification of
coeffs, see effect_alt which simplifies this process.
Response variables are simulated through a copula model by first extracting Gaussian copular scores
as Dunn-Smyth residuals (Dunn & Smyth 1996), which are obtained from abundances with marginal distributions
which have been specified via the original manyglm model (fit.glm; see examples);
These scores then follow a multivariate Gaussian distribution with zero mean and covariance structure ,
To avoid estimating a large number pairwise correlations within , factor analysis is utilised
with two latent factor variables, which can be interpreted as an unobserved environmental covariate.
Thus, in order to simulate new multivariate abundances we simulate new copula scores and back transform them to
abundances as , where the coefficient matrix coeffs specifies the
effect size within the new marginal distributions .
The data frame is also extended in a manner that preserves the original design structure. This is done by first
repeating the design matrix until the number of samples exceeds N, then randomly removing rows from the last
repeated data frame until the number of samples equals N. Alternatively, a balanced design structure can be
obtained by specifying the number of replicates.
newdata can be utilised if a different data frame is wanted for simulation.
If users are interested in obtaining a manyglm model, do.fit=TRUE can be used to obtain a manyglm
object from the simulated responses.
Simulated data or manyglm object.
extend(): Simulate or extend multivariate abundance data
Dunn, P.K., & Smyth, G.K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics 5, 236-244.
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") # Simulate data fit.glm = manyglm(spiddat~1, family="negative.binomial") fit.cord = cord(fit.glm) simData = extend(fit.cord) # Simulate data with N=20 fit.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) fit.cord = cord(fit.glm) simData = extend(fit.cord, N=20) # Obtain a manyglm fit from simulated data with N=10 and 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) newFit.glm = extend(fit_factors.cord, N=10, coeffs=effect_mat, do.fit=TRUE) # Change sampling design X_new = X X_new$Treatment[6:7] = c("B","B") simData = extend(fit_factors.cord, N=NULL, coeffs=effect_mat, newdata=X_new, n_replicate=5)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") # Simulate data fit.glm = manyglm(spiddat~1, family="negative.binomial") fit.cord = cord(fit.glm) simData = extend(fit.cord) # Simulate data with N=20 fit.glm = manyglm(spiddat~soil.dry, family="negative.binomial", data=X) fit.cord = cord(fit.glm) simData = extend(fit.cord, N=20) # Obtain a manyglm fit from simulated data with N=10 and 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) newFit.glm = extend(fit_factors.cord, N=10, coeffs=effect_mat, do.fit=TRUE) # Change sampling design X_new = X X_new$Treatment[6:7] = c("B","B") simData = extend(fit_factors.cord, N=NULL, coeffs=effect_mat, newdata=X_new, n_replicate=5)
Dataset of fish abundances and associated environmental variables.
data(fish)data(fish)
An object of class "data.frame" containing:
A data frame with nine observations of abundance of 34 fish species, and five site-related variables.
The matrix abund has the following species abundances:
Abudefduf.sp
Acanthurus.nigrofuscus
Achoerodus.viridis
Aplodactylus.lophodon
Atypichthys.strigatus
Baitfish
Brachaluteres.jacksonianus
Chaeotodon.auriga
Cheilodactylus.fuscus
Chromis.hypsilepis
Girella.elevata
Girella.tricuspidata
Heterodontus.portusjacksoni
Kyphosus.sydneyanus
Latropiscis.purpurissatus
Meuschenia.spp
Microcanthus.strigatus
Naso.unicornis
Notolabrus.gymnogenis
Olisthops.cyanomelas
Ophthalmolepis.lineolatus
Pagrus.auratus
Parma.microlepis
Parma.unifasciata
Parupeneus.signatus
Pempheris.compressa
Pictilabrus.laticlavius
Prionurus.maculatus
Prionurus.microlepidotus
Scorpis.lineolatus
Seriola.lalandi
Seriola.sp
Trachinops.taeniatus
Unidentified.wrasse
Site.Type - control / reference / restored
Site.Name - location of site
Viz - visibility in metres
Temp - temperature in degrees Celsius
Depth - depth in metres
Data attributed to the crayweed restoration project (http://www.operationcrayweed.com/).
data(fish) head(fish)data(fish) head(fish)
powersim returns a power estimate for a cord object for a given sample size N
and effect size of interest.
## S3 method for class 'cord' powersim( object, coeffs, term, N = nrow(object$obj$data), coeffs0 = effect_null(object$obj, term), nsim = 1000, ncrit = 999, test = "score", alpha = 0.05, newdata = NULL, ncores = detectCores() - 1, show.time = TRUE, long_power = FALSE, n.samp = 10, nlv = 2 ) powersim( object, coeffs, term, N = nrow(object$obj$data), coeffs0 = effect_null(object$obj, term), nsim = 999, ncrit = nsim, test = "score", alpha = 0.05, newdata = NULL, ncores = detectCores() - 1, show.time = TRUE, long_power = FALSE, n.samp = 10, nlv = 2 )## S3 method for class 'cord' powersim( object, coeffs, term, N = nrow(object$obj$data), coeffs0 = effect_null(object$obj, term), nsim = 1000, ncrit = 999, test = "score", alpha = 0.05, newdata = NULL, ncores = detectCores() - 1, show.time = TRUE, long_power = FALSE, n.samp = 10, nlv = 2 ) powersim( object, coeffs, term, N = nrow(object$obj$data), coeffs0 = effect_null(object$obj, term), nsim = 999, ncrit = nsim, test = "score", alpha = 0.05, newdata = NULL, ncores = detectCores() - 1, show.time = TRUE, long_power = FALSE, n.samp = 10, nlv = 2 )
object |
objects of class |
coeffs |
Coefficient matrix for a |
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 |
nsim |
Number of simulated test statistics under the specified effect size ( |
ncrit |
Number of simulated test statistics under the null effect to estimate the critical value. Defaults to |
test |
Test statistic for power estimate to based upon. Defaults to |
alpha |
Type I error rate for power estimate, defaults to |
newdata |
Data frame of the same size as the original data frame from the |
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 |
long_power |
Logical. Whether to estimate power using separate critical test statistics for each |
n.samp |
integer, number of sets of residuals for importance sampling for the copula model with cord. Defaults to |
nlv |
number of latent variables (default = 2) for the copula model with cord, recommend setting this lower for smaller sample sizes |
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 ncrit times with a null effect, which is defined by default by
effect_null. This effectively simulates nsim + ncrit 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 errors relating to clustering or nodes.
Power estimate result, and;
power |
power. |
powersim(): Provide power estimates for multivariate abundance models
Ben Maslen <[email protected]>.
Maslen, B., Popovic, G., Lim, M., Marzinelli, E., & Warton, D. (2023). How many sites? Methods to assist design decisions when collecting multivariate data in ecology. Methods in Ecology and Evolution, 14(6), 1564-1573. Popovic, G. C., Hui, F. K., & Warton, D. I. (2018). A general algorithm for covariance modeling of discrete data. Journal of Multivariate Analysis, 165, 86-100.
effect_alt, effect_null, extend
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, ncrit=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, ncrit=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, ncrit=99, ncores=2)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, ncrit=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, ncrit=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, ncrit=99, ncores=2)