Provides p-values for omnibus tests based on permutations for factorial and repeated measures ANOVA. This function produces the F statistics, parametric p-values (based, on Gaussian and sphericity assumptions) and p-values based on the permutation methods that handle nuisance variables.

aovperm(
  formula,
  data = NULL,
  np = 5000,
  method = NULL,
  type = "permutation",
  ...
)

Arguments

formula

A formula object. The formula for repeated measures ANOVA should be written using the same notation as aov by adding +Error(id/within), where id is the factor that identify the subjects and within is the within factors.

data

A data frame or matrix.

np

The number of permutations. Default value is 5000.

method

A character string indicating the method used to handle nuisance variables. Default is NULL and will change if to "freedman_lane" for the fixed effects model and "Rd_kheradPajouh_renaud" for the random effect models. See Details for other methods.

type

A character string to specify the type of transformations: "permutation" and "signflip" are available. Is overridden if P is given. See help from Pmat.

...

Futher arguments, see details.

Value

A lmperm object containing most of the objects given in an lm object, an ANOVA table with parametric and permutation p-values, the test statistics and the permutation distributions.

Details

The following methods are available for the fixed effects model defined as \(y = D\eta + X\beta + \epsilon\). If we want to test \(\beta = 0\) and take into account the effects of the nuisance variables \(D\), we transform the data :

method argument\(y*\)\(D*\)\(X*\)
"draper_stoneman"\(y\)\(D\)\(PX\)
"freedman_lane"\((H_D+PR_D)y\)\(D\)\(X\)
"manly"\(Py\)\(D\)\(X\)
"terBraak"\((H_{X,D}+PR_{X,D})y\)\(D\)\(X\)
"kennedy"\(PR_D y\)\(R_D X\)
"huh_jhun"\(PV'R_Dy\)\(V'R_D X\)
"dekker"\(y\)\(D\)\(PR_D X\)

The following methods are available for the random effects model \(y = D\eta + X\beta + E\kappa + Z\gamma+ \epsilon\). If we want to test \(\beta = 0\) and take into account the effect of the nuisance variable \(D\) we can transform the data by permutation:

method argument\(y*\)\(D*\)\(X*\)\(E*\)\(Z*\)
"Rd_kheradPajouh_renaud"\(PR_D y\)\(R_D X\)\(R_D E\)\(R_D Z\)
"Rde_kheradPajouh_renaud"\(PR_{D,E}y\)\(R_{D,E} X\)\(R_{D,E}Z\)

Other arguments could be pass in ... :

P : a matrix, of class matrix or Pmat, containing the permutations (for the reproductibility of the results). The first column must be the identity permutation (not checked). P overwrites np argument.
rnd_rotation : a random matrix of size \(n \times n\) to compute the rotation used for the "huh_jhun" method. coding_sum : a logical set to TRUE defining the coding of the design matrix to contr.sum to test the main effects. If it is set to FALSE the design matrix is computed with the coding defined in the dataframe. The tests of simple effets are possible with a coding of the factors of the dataframe set to contr.treatment.

See also

lmperm plot.lmperm

Other main function: clusterlm(), lmperm()

Author

jaromil.frossard@unige.ch

Examples

## data
data("emergencycost")

## centrering the covariate to the mean
emergencycost$LOSc <- scale(emergencycost$LOS, scale = FALSE)

## ANCOVA
## Warning : np argument must be greater (recommendation: np>=5000)
mod_cost_0 <- aovperm(cost ~ LOSc*sex*insurance, data = emergencycost, np = 2000)
mod_cost_0
#> Anova Table
#> Resampling test using freedman_lane to handle nuisance variables and 2000 permutations.
#>                           SS  df        F parametric P(>F) resampled P(>F)
#> LOSc               2.162e+09   1 483.4422          0.00000          0.0005
#> sex                1.463e+07   1   3.2714          0.07229          0.0840
#> insurance          6.184e+05   1   0.1383          0.71048          0.6785
#> LOSc:sex           8.241e+06   1   1.8427          0.17646          0.1630
#> LOSc:insurance     2.911e+07   1   6.5084          0.01163          0.0265
#> sex:insurance      1.239e+05   1   0.0277          0.86801          0.8620
#> LOSc:sex:insurance 1.346e+07   1   3.0091          0.08463          0.0870
#> Residuals          7.514e+08 168                                          

## same analysis but with signflip
## Warning : np argument must be greater (recommendation: np>=5000)
mod_cost_0s <- aovperm(cost ~ LOSc*sex*insurance, data = emergencycost, type="signflip", np = 2000)
mod_cost_0s
#> Anova Table
#> Resampling test using freedman_lane to handle nuisance variables and 2000 signflips.
#>                           SS  df        F parametric P(>F) resampled P(>F)
#> LOSc               2.162e+09   1 483.4422          0.00000          0.0005
#> sex                1.463e+07   1   3.2714          0.07229          0.2015
#> insurance          6.184e+05   1   0.1383          0.71048          0.7390
#> LOSc:sex           8.241e+06   1   1.8427          0.17646          0.1740
#> LOSc:insurance     2.911e+07   1   6.5084          0.01163          0.0385
#> sex:insurance      1.239e+05   1   0.0277          0.86801          0.8850
#> LOSc:sex:insurance 1.346e+07   1   3.0091          0.08463          0.1360
#> Residuals          7.514e+08 168                                          

## Testing at 14 days
emergencycost$LOS14 <- emergencycost$LOS - 14

mod_cost_14 <- aovperm(cost ~ LOS14*sex*insurance, data = emergencycost, np = 2000)
mod_cost_14
#> Anova Table
#> Resampling test using freedman_lane to handle nuisance variables and 2000 permutations.
#>                            SS  df        F parametric P(>F) resampled P(>F)
#> LOS14               2.162e+09   1 483.4422          0.00000          0.0005
#> sex                 2.760e+07   1   6.1703          0.01397          0.0210
#> insurance           9.864e+05   1   0.2206          0.63923          0.6075
#> LOS14:sex           8.241e+06   1   1.8427          0.17646          0.1490
#> LOS14:insurance     2.911e+07   1   6.5084          0.01163          0.0215
#> sex:insurance       7.722e+05   1   0.1727          0.67829          0.6640
#> LOS14:sex:insurance 1.346e+07   1   3.0091          0.08463          0.0760
#> Residuals           7.514e+08 168                                          

## Effect of sex within the public insured
contrasts(emergencycost$insurance) <- contr.treatment
contrasts(emergencycost$sex) <- contr.sum
emergencycost$insurance <- relevel(emergencycost$insurance, ref = "public")

mod_cost_se <- aovperm(cost ~ LOSc*sex*insurance, data = emergencycost,
                        np = 2000, coding_sum = FALSE)
mod_cost_se
#> Anova Table
#> Resampling test using freedman_lane to handle nuisance variables and 2000 permutations.
#>                           SS  df         F parametric P(>F) resampled P(>F)
#> LOSc               9.512e+09   1 2126.7539        0.000e+00          0.0005
#> sex                6.092e+07   1   13.6210        3.019e-04          0.0005
#> insurance          6.184e+05   1    0.1383        7.105e-01          0.6920
#> LOSc:sex           1.510e+08   1   33.7708        3.039e-08          0.0005
#> LOSc:insurance     2.911e+07   1    6.5084        1.163e-02          0.0215
#> sex:insurance      1.239e+05   1    0.0277        8.680e-01          0.8600
#> LOSc:sex:insurance 1.346e+07   1    3.0091        8.463e-02          0.0785
#> Residuals          7.514e+08 168                                           


## Repeated measures ANCOVA
## data
data(jpah2016)

## centrering the covariate
jpah2016$bmic <- scale(jpah2016$bmi, scale = FALSE)

## Warning : np argument must be greater (recommendation: np>=5000)
mod_jpah2016 <- aovperm(iapa ~ bmic*condition*time+ Error(id/(time)),
                    data = jpah2016, method = "Rd_kheradPajouh_renaud")
#> Warning: The data are not balanced, the results may not be exact.
mod_jpah2016
#> 
#> Resampling test using Rd_kheradPajouh_renaud to handle nuisance variables and 5000 permutations.
#>                          SSn dfn    SSd dfd     MSEn  MSEd        F
#> bmic                   18.68   1 106884  13    18.68  8222 0.002272
#> condition           27878.20   2 106884  13 13939.10  8222 1.695381
#> bmic:condition      89238.48   2 106884  13 44619.24  8222 5.426938
#> time                  268.84   1 167305  13   268.84 12870 0.020889
#> bmic:time             366.49   1 167305  13   366.49 12870 0.028477
#> condition:time      21159.77   2 167305  13 10579.89 12870 0.822083
#> bmic:condition:time 29145.72   2 167305  13 14572.86 12870 1.132347
#>                     parametric P(>F) resampled P(>F)
#> bmic                         0.96271          0.9652
#> condition                    0.22169          0.2206
#> bmic:condition               0.01934          0.0214
#> time                         0.88730          0.8804
#> bmic:time                    0.86859          0.8644
#> condition:time               0.46112          0.4428
#> bmic:condition:time          0.35209          0.3526