welchADF.test: General Approximate Degrees of Freedom Solution for Inference...

Description Usage Arguments Value Methods (by class) References See Also Examples

View source: R/welch-wrapper.R

Description

Computes the Welch-James statistic with Approximate Degrees of Freedom, based on the SAS code by H.J. Keselman, R.R. Wilcox and L.M. Lix.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
welchADF.test(formula, ...)

## S3 method for class 'formula'
welchADF.test(formula, data, subset, ...)

## S3 method for class 'lm'
welchADF.test(formula, ...)

## S3 method for class 'aov'
welchADF.test(formula, ...)

## S3 method for class 'lmer'
welchADF.test(formula, ...)

## Default S3 method:
welchADF.test(formula, response, between.s = NULL,
  within.s = NULL, subject = NULL, contrast = c("omnibus",
  "all.pairwise"), effect = NULL, correction = c("hochberg", "holm"),
  trimming = FALSE, per = 0.2, bootstrap = FALSE, numsim_b = 999,
  effect.size = FALSE, numsim_es = 999, scaling = TRUE,
  standardize.effsz = TRUE, alpha = 0.05, seed = 0, ...)

Arguments

formula

A data frame or a formula or an lm object returned by lm. If formula is a data frame, it must have the format described for argument data of welchADF.test for class formula. If formula is a formula, it must be a two-sided linear formula object describing the between-subjects effects of the model, with the response variable(s) (surrounded by cbind() if it is a multivariate response, as cbind(r1, r2, r3)) on the left of the ~ operator and the terms on the right. Within-subject factors must be specified as a parenthesis term on the rhs with the form (w1*w2*...|subject) with 'subject' being the name of the subject ID column, which must be present in the data (cannot be omitted like in the data.frame version) e.g. cbind(y1,y2) ~ A*B*C + (B*C|subjet). In this example, A would be a between-subject factor because it is the only variable not present again on the left side of the parenthesis term (where the within-subjects factors are specified). Multi-level or customized structures are not supported, e.g. those having strata or more than one random factor. The formula does not affect the results of the tests as no 'model' is actually fit to the data. It only affects the terms displayed as a result of an omnibus contrast when no effect argument is given (which means run an omnibus test for every effect or interaction (up to second order)). In such case, the factors or interactions not specified in the formula will not be applied an omnibus test. In case of a pairwise contrast, the structure in the formula is ignored beyond interpreting which variables are between-subject factors and which are within-subjects factors.

...

Further arguments to be passed to specialized methods. Should be named arguments.

data

A data.frame with the data, formatted as follows: one row per observation of the response variable(s), and as many columns as needed to indicate the factor combination to which the observation corresponds. If necessary, an extra column with the subject ID for designs having within-subjects factors (can be omitted if there is only one within-subjects factor, see the vignette).

subset

A specification of the rows to be used as in model.frame: defaults to all rows. This can be any valid indexing vector (see [.data.frame) for the rows of data or if that is not supplied, a data frame made up of the variables used in formula.

response

A string or vector of strings with the name(s) of the column(s) of data corresponding to the response variable(s). If a vector of strings is provided, the responses are taken as a set of repeated measurements or dependent variables.

between.s

Vector of strings with the columns that correspond to between-subjects factors, if any (defaults to NULL).

within.s

Vector of strings with the columns that correspond to within-subjects factors, if any (defaults to NULL).

subject

Name of the column (if any) containing the subject ID (defaults to NULL).

contrast

One of "omnibus", "all.pairwise" indicating the type of contrast to be performed with the data. "omnibus" stands for omnibus contrasts on either all the variables found in the data, including the within-subjects factor if any, or on a single effect specified in the effect argument. "all.pairwise" stands for pairwise comparisons of the response(s) variable(s) over all combinations of values of one effect, indicated in the effect parameter. Defaults to "omnibus".

effect

If contrast = "omnibus", then effect must be a string or vector of strings with the name of the effect(s) for which the "omnibus" contrasts will be applied. If contrast = "all.pairwise", then effect must be a tagged list in which the names are the factors, and the elements must be vectors of two strings each, indicating If contrast = "omnibus" but effect is not specified, then an omnibus contrast will be done separately for each of the effects and also all possible interactions. If contrast = "omnibus" and effect is a vector of strings, then only one contrast will be done, involving the simultaneous interaction of all effects indicated in the effect vector (i.e. no other interaction or main effect contrasts will be performed).

correction

The type of p-value correction when applying multiple pariwise comparisons (i.e. when contrast = "all.pairwise"). Defaults to "hochberg".

trimming

Boolean to control the use of robust estimators. FALSE indicates the ADF test statistic will be computed with the usual Least-Squares Estimators. TRUE indicates that trimmed means and Winsorized variances/covariances will be used to compute the ADF test statistic.

per

(Only if trimming is TRUE) The proportion of trimming in EACH tail of the data distribution. Must be between 0 and .49 (i.e., 49% trimming in each tail). Recommended per = 0.2 symmetric trimming (i.e., 20% of the observations in each tail are trimmed).

bootstrap

Boolean; whether bootstrap should be used to compute a critical value for the test statistic produced by WJGLM. FALSE means a theoretical critical value for the test statistic TWJ/c (or TWJt/c) will be computed. TRUE indicates the non-parametric percentile bootstrap method will be used to compute a critical value. When combined with contrast = "all.pairwise", then Family-Wise Error Rate (FWER) control is achieved using bootstrapping to compute an adjusted empirical critical value (refer to Keselman, Wilcox & Lix 2003, page 589 for details).

numsim_b

If bootstrap = TRUE, this is a positive integer that defines the number of bootstrap samples to generate a critical value (defaults to 999). Carpenter and Bithell (2000) recommend between 1000 and 2000 samples for setting confidence intervals around a parameter estimate, but Efron and Tibshirani (1986) note that for other applications of the bootstrap, such as computing the standard error, as few as 100 bootstrap samples will be sufficient. Choosing a value of B such that (B + 1)α is an integer value is recommended, because it avoids the need for interpolation.

effect.size

Boolean; whether effect size estimates should be computed. TRUE means effect size estimates together with its confidence intervals will be computed (defaults to FALSE).

numsim_es

Positive integer defining the number of bootstrap samples to generate a CI for the effect size estimate (defaults to 999). Ignored if effect.size = FALSE.

scaling

Boolean; whether a scaling factor for the effect size estimator should be used (0.642 for 20% symmetric trimming) when robust estimators are adopted. FALSE means no scaling factor will be used. TRUE indicates that a scaling factor will be adopted. Defaults to TRUE

standardize.effsz

Boolean: whether the effect size should be standardized by the average of variances or not. Defaults to TRUE.

alpha

Significance threshold for the confidence interval calculation, where 1 - alpha is the probability coverage. Defaults to 0.05.

seed

Initial seed (positive integer) for the (R default) random number generator. Defaults to NULL (seed taken from the current time).

Value

An object of class "welchADFt" which is a list of lists (one sub-list per effect, even if there is only one). There are methods print.welchADFt and summary.welchADFt

Methods (by class)

References

Villacorta, P.J. (2017). The welchADF Package for Robust Hypothesis Testing in Unbalanced Multivariate Mixed Models with Heteroscedastic and Non-normal Data. The R Journal, 9:2, 309 - 328.

Website with the original SAS code and examples: http://supp.apa.org/psycarticles/supplemental/met_13_2_110/met_13_2_110_supp.html

Keselman, H. J., Wilcox, R. R., & Lix, L. M. (2003). A generally robust approach to hypothesis testing in independent and correlated groups designs. Psychophysiology, 40, 586-596.

Lix, L. M., & Keselman, H. J. (1995). Approximate degrees of freedom tests: A unified perspective on testing for mean equality. Psychological Bulletin, 117, 547-560.

Carpenter, J., & Bithell, J. (2000). Bootstrap confidence intervals: When, which, what? A practical guide for medical statisticians. Statistics in Medicine, 19, 1141-1164.

Efron, B., & Tibshirani, R. (1986). Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statistical Science, 1, 54-75.

See Also

print.welchADFt, summary.welchADFt

p.adjust.methods perceptionData adhdData adhdData2 womenStereotypeData miceData

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
31
32
33
34
35
# Omnibus contrast only of those effects included, namely condition and sex (no interaction)
omnibus_LSM_formula <- welchADF.test(y ~ condition + sex, data = womenStereotypeData)
# Works well with update.default() method
omnibus_interact_formula <- update(omnibus_LSM_formula, . ~ condition*sex)
summary(omnibus_LSM_formula)
summary(omnibus_interact_formula)

# Fit a linear model using the built-in function stats::lm
lm.women <- lm(y ~ condition + sex, womenStereotypeData)

# Fit an Analysis of Variance model using the built-in function stats::aov
aov.women <- aov(lm.women)

# Now use the this object to apply a welchADF test to the same formula and data
omnibus_no_interact <- welchADF.test(lm.women, contrast = "omnibus")
omnibus_no_interactB <- welchADF.test(aov.women) # omnibus as well

# Integrates well with the update.default() method
omnibus_interact <- update(omnibus_no_interact, . ~ condition*sex)
summary(omnibus_no_interact)
summary(omnibus_interact)
# Two-way factorial design using the default interface. See the vignette.
omnibus_LSM <- welchADF.test(womenStereotypeData, response = "y", between.s = 
  c("condition", "sex"), contrast = "omnibus")
# Method update() also works with the welchADF.test.default interface
omnibus_trimmed <- update(omnibus_LSM, effect = "condition", trimming = TRUE)
pairwise_LSM <- update(omnibus_LSM, contrast = "all.pairwise", effect = c("condition", "sex"))
pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE) # just trimming
summary(omnibus_LSM)
pairwise_LSM
## Not run: 
  pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap = TRUE) # trimming and bootstrapping
  summary(pairwise_trimmed_boot)

## End(Not run)

Example output

Call:
   welchADF.test(formula = y ~ condition + sex, data = womenStereotypeData)

          WJ statistic Numerator DF Denominator DF Pr(>WJ)  
condition        2.151            2          154.7 0.11986  
sex              2.933            1          216.4 0.08824 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
   welchADF.test(formula = y ~ condition + sex + condition:sex, 
    data = womenStereotypeData)

              WJ statistic Numerator DF Denominator DF Pr(>WJ)  
condition            2.151            2          154.7 0.11986  
sex                  2.933            1          216.4 0.08824 .
condition:sex        2.521            2          154.7 0.08368 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
   welchADF.test(formula = y ~ condition + sex, data = womenStereotypeData, 
    contrast = "omnibus")

          WJ statistic Numerator DF Denominator DF Pr(>WJ)  
condition        2.151            2          154.7 0.11986  
sex              2.933            1          216.4 0.08824 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
   welchADF.test(formula = y ~ condition + sex + condition:sex, 
    data = womenStereotypeData, contrast = "omnibus")

              WJ statistic Numerator DF Denominator DF Pr(>WJ)  
condition            2.151            2          154.7 0.11986  
sex                  2.933            1          216.4 0.08824 .
condition:sex        2.521            2          154.7 0.08368 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
   welchADF.test(formula = womenStereotypeData, response = "y", 
    between.s = c("condition", "sex"), contrast = "omnibus")

              WJ statistic Numerator DF Denominator DF Pr(>WJ)  
condition            2.151            2          154.7 0.11986  
sex                  2.933            1          216.4 0.08824 .
condition:sex        2.521            2          154.7 0.08368 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
   welchADF.test(formula = womenStereotypeData, response = "y", 
    between.s = c("condition", "sex"), contrast = "all.pairwise", 
    effect = c("condition", "sex"))

Welch-James Approximate DF Test (Least squares means & variances)
Multiple tetrad interaction contrasts with respect to condition x sex interaction
using Hochberg method for adjusting the p-values
         control:stereotype x female:male nullified:stereotype x female:male
adj.pval                        0.8111564                           0.154972
         control:nullified x female:male
adj.pval                        0.154972 
Call:
   welchADF.test(formula = womenStereotypeData, response = "y", 
    between.s = c("condition", "sex"), contrast = "all.pairwise", 
    effect = c("condition", "sex"), trimming = TRUE, bootstrap = TRUE)

                                   WJ statistic Numerator DF Denominator DF
control:stereotype x female:male         0.4662            1          97.49
nullified:stereotype x female:male       1.9846            1          79.63
control:nullified x female:male          5.7662            1          88.55
                                   significant?
control:stereotype x female:male             no
nullified:stereotype x female:male           no
control:nullified x female:male              no

Bootstrap critical value: 5.827

welchADF documentation built on Sept. 8, 2019, 9:02 a.m.