R/aovperm.R
aovperm.RdProvides 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",
...
)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.
A data frame or matrix.
The number of permutations. Default value is 5000.
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.
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.
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.
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.
## 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