Title: | Welch-James Statistic for Robust Hypothesis Testing under Heterocedasticity and Non-Normality |
---|---|
Description: | Implementation of Johansen's general formulation of Welch-James's statistic with Approximate Degrees of Freedom, which makes it suitable for testing any linear hypothesis concerning cell means in univariate and multivariate mixed model designs when the data pose non-normality and non-homogeneous variance. Some improvements, namely trimmed means and Winsorized variances, and bootstrapping for calculating an empirical critical value, have been added to the classical formulation. The code departs from a previous SAS implementation by L.M. Lix and H.J. Keselman, available at <http://supp.apa.org/psycarticles/supplemental/met_13_2_110/SAS_Program.pdf> and published in Keselman, H.J., Wilcox, R.R., and Lix, L.M. (2003) <DOI:10.1111/1469-8986.00060>. |
Authors: | Pablo J. Villacorta <[email protected]> |
Maintainer: | Pablo J. Villacorta <[email protected]> |
License: | LGPL (>= 3) |
Version: | 0.3.2 |
Built: | 2024-12-02 06:44:36 UTC |
Source: | CRAN |
The data (Keselman et al., 2003) represent the reaction times in milliseconds of children with attention-deficit hyperactivity (ADHD) and normal children when they are presented four kinds of inputs: a target alone or an arrow stimuli incongruent, congruent and neutral to the target. According to the authors, the dataset was artificially generated from the summary measures given in the original study by Jonkman et al. (1999), in groups of 20 and 10 children to create an unbalanced design.
adhdData
adhdData
A data frame with 30 rows and 5 variables:
whether the child has ADHD or is healty (normal)
.
Reaction time (milliseconds) to a target alone.
Reaction time (milliseconds) to a congruent stimulus.
Reaction time (milliseconds) to a neutral stimulus.
Reaction time (milliseconds) to an incongruent stimulus.
H. J. Keselman, R. R. Wilcox, and L. M. Lix. A generally robust approach to hypothesis testing in independent and correlated groups designs. Psychophyisiology, 40:586-596, 2003. (Data displayed in page 593).
L. Jonkman, C. Kemner, M. Verbaten, H. van Engeland, J. Kenemans, G. Camfferman, J. Buitelaar, and H. Koelega. Perceptual and response interference in children with attention-deficit hyperactivity disorder, and the effects of methylphenidate. 36(4):419-429, 1999.
# Assuming a multivariate response omnibus_LSM_multi <- welchADF.test(adhdData, response = c("TargetAlone", "Congruent", "Neutral", "Incongruent"), between.s = "Group", within.s = "multivariate", contrast = "omnibus") # The same using the S3 method for class formula omnibus_LSM_multi_form <- welchADF.test(cbind(TargetAlone, Congruent, Neutral, Incongruent) ~ Group, data = adhdData) # Pairwise comparisons of the implicit within-subjects effect present in the multivariate response pairwise_LSM_multi <- update(omnibus_LSM_multi, contrast = "all.pairwise", effect = "multivariate") summary(omnibus_LSM_multi) summary(pairwise_LSM_multi)
# Assuming a multivariate response omnibus_LSM_multi <- welchADF.test(adhdData, response = c("TargetAlone", "Congruent", "Neutral", "Incongruent"), between.s = "Group", within.s = "multivariate", contrast = "omnibus") # The same using the S3 method for class formula omnibus_LSM_multi_form <- welchADF.test(cbind(TargetAlone, Congruent, Neutral, Incongruent) ~ Group, data = adhdData) # Pairwise comparisons of the implicit within-subjects effect present in the multivariate response pairwise_LSM_multi <- update(omnibus_LSM_multi, contrast = "all.pairwise", effect = "multivariate") summary(omnibus_LSM_multi) summary(pairwise_LSM_multi)
Exactly the same data explained in "adhdData" but reshaped as follows.
adhdData2
adhdData2
A data frame with 120 rows and 4 variables:
whether the child has ADHD or is healty (normal)
.
the stimulus to which the reaction time in this row corresponds.
an integer ID that identifies the subject to which the reaction time corresponds.
reaction time (milliseconds) to of the aforementioned subject to the aforementioned stimulus.
# Omnibus test of a mixed between x within subjects model, # using trimmed means and Winsorized variances omnibus_trimmed <- welchADF.test(adhdData2, response = "Milliseconds",between.s = "Group", within.s = "Stimulus", subject = "Subject", contrast = "omnibus", trimming = TRUE) # The same using S3 method for class formula. The data can be in separate variables in # the scope of the call, not necessarily in a data.frame millis <- adhdData2$Milliseconds gr <- adhdData2$Group st <- adhdData2$Stimulus sbj <- adhdData2$Subject omnibus_trimmed_formula <- welchADF.test(millis ~ gr*st + (st|sbj), contrast = "omnibus", trimming = TRUE) summary(omnibus_trimmed_formula) # Pairwise contrasts of the effects pairwise_LSM <- welchADF.test(adhdData2, response = "Milliseconds", between.s = "Group", within.s = "Stimulus", subject = "Subject", contrast = "all.pairwise", effect = "Stimulus") pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE) # Bootstrapping to obtain an empirical critical value ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap =TRUE, seed = 123456) summary(pairwise_trimmed_boot) ## End(Not run)
# Omnibus test of a mixed between x within subjects model, # using trimmed means and Winsorized variances omnibus_trimmed <- welchADF.test(adhdData2, response = "Milliseconds",between.s = "Group", within.s = "Stimulus", subject = "Subject", contrast = "omnibus", trimming = TRUE) # The same using S3 method for class formula. The data can be in separate variables in # the scope of the call, not necessarily in a data.frame millis <- adhdData2$Milliseconds gr <- adhdData2$Group st <- adhdData2$Stimulus sbj <- adhdData2$Subject omnibus_trimmed_formula <- welchADF.test(millis ~ gr*st + (st|sbj), contrast = "omnibus", trimming = TRUE) summary(omnibus_trimmed_formula) # Pairwise contrasts of the effects pairwise_LSM <- welchADF.test(adhdData2, response = "Milliseconds", between.s = "Group", within.s = "Stimulus", subject = "Subject", contrast = "all.pairwise", effect = "Stimulus") pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE) # Bootstrapping to obtain an empirical critical value ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap =TRUE, seed = 123456) summary(pairwise_trimmed_boot) ## End(Not run)
For the effects indicated by the user, extracts confidence intervals
for the effect size from an object of class welchADFt
returned by a call to welchADF.test
.
## S3 method for class 'welchADFt' confint(object, parm, level = 0.95, ...)
## S3 method for class 'welchADFt' confint(object, parm, level = 0.95, ...)
object |
a object of (S3 class) welchADF.t |
parm |
a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
confidence level (completely ignored, since the confidence threshold
had been in the call to |
... |
additional argument(s) for methods. |
If effect.size
was set to TRUE
in the call to welchADF.test
,
a matrix (or vector) with columns giving lower and upper confidence limits for each parameter.
These will be labelled as (1-level)/2 and 1 - (1-level)/2 in
If effect.size
was set to FALSE
, NULL
is returned.
Wild strain house mice were, at birth, cross fostered onto house mouse (Mus), deer mouse (Peromyscus) or rat (Rattus) nursing mothers. Ten days after weaning, each subject was tested in an apparatus that allowed it to enter four different tunnels: one scented with clean pine shavings, and the other three tunnels with shavings bearing the scent of Mus, Peromyscus, or Rattus respectively. Three variables were measured for each tunnel: the number of visits to the tunnel during a twenty minute test, the time spent by each subject in each of the four tunnels and the latency to first visit of each tunnel.
miceData
miceData
A data frame with 144 rows and 6 variables:
an ID that identifies the mouse to which the measurements of the row correspond.
type of nursing mother this mouse had at birth.
type of tunnel to which the measurements of this row correspond.
number of visits to the aforementioned tunnel.
total time spent by the mouse at this tunnel.
seconds passed before the mouse first visited this tunnel.
http://core.ecu.edu/psyc/wuenschk/SPSS/TUNNEL4b.sav
K. L.Wuensch. Fostering house mice onto rats and deer mice: Effects on response to species odors. Animal Learning & Behavior, 20(3):253 - 258, 1992
omnibus_LSM <- welchADF.test(miceData, response = c("visits", "time", "latency"), between.s = "nurs", within.s = "tunnel", subject = "Subject", contrast = "omnibus") # Formula interface, using cbind() to specify a multivariate response. omnibus_LSM_formula <- welchADF.test( cbind(visits, time, latency) ~ nurs*tunnel + (tunnel|Subject), data = miceData) pairwise_LSM_nurs <- welchADF.test(miceData, response = c("visits", "time", "latency"), between.s = "nurs", within.s = "tunnel", subject = "Subject", effect = "nurs", contrast = "all.pairwise") pairwise_LSM_tunnel <- update(pairwise_LSM_nurs, effect = "tunnel") ## Not run: pairwise_nurs_trimmed_boot <- update(pairwise_LSM_nurs, trimming = TRUE, bootstrap = TRUE) pairwise_tunnel_trimmed_boot <- update(pairwise_nurs_trimmed_boot, effect = "tunnel") summary(pairwise_nurs_trimmed_boot) ## End(Not run) summary(omnibus_LSM) summary(pairwise_LSM_nurs) summary(pairwise_LSM_tunnel)
omnibus_LSM <- welchADF.test(miceData, response = c("visits", "time", "latency"), between.s = "nurs", within.s = "tunnel", subject = "Subject", contrast = "omnibus") # Formula interface, using cbind() to specify a multivariate response. omnibus_LSM_formula <- welchADF.test( cbind(visits, time, latency) ~ nurs*tunnel + (tunnel|Subject), data = miceData) pairwise_LSM_nurs <- welchADF.test(miceData, response = c("visits", "time", "latency"), between.s = "nurs", within.s = "tunnel", subject = "Subject", effect = "nurs", contrast = "all.pairwise") pairwise_LSM_tunnel <- update(pairwise_LSM_nurs, effect = "tunnel") ## Not run: pairwise_nurs_trimmed_boot <- update(pairwise_LSM_nurs, trimming = TRUE, bootstrap = TRUE) pairwise_tunnel_trimmed_boot <- update(pairwise_nurs_trimmed_boot, effect = "tunnel") summary(pairwise_nurs_trimmed_boot) ## End(Not run) summary(omnibus_LSM) summary(pairwise_LSM_nurs) summary(pairwise_LSM_tunnel)
An artificial dataset created by Lix et al. which recreates data reported by a real study on perception and concentration, on which 42 students were given several puzzles to be solved. The students are divided into three balanced groups as they had previously been asked to imagine solving puzzles in the distant future, near future, or not to imagine anything at all (control group).
perceptionData
perceptionData
A data frame with 42 rows and 2 variables:
group of the student (what the student was asked to imagine)
number of puzzles the student was able to solve, out of 12
Forster, J., Liberman, N., & Friedman, R.S. (2004). Temporal construal effects on abstract and concrete thinking: consequences for insight and creative cognition. Journal of Personality and Social Psychology, 87, 2, 177-189.
omnibus_LSM <- welchADF.test(perceptionData, response = "y", between.s = "Group") omnibus_trimmed <- update(omnibus_LSM, trimming = TRUE) pairwise_LSM <- update(omnibus_LSM, effect = "Group", contrast = "all.pairwise") pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE, effect.size = TRUE) summary(omnibus_LSM) ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap = TRUE, seed = 12345, numsim_b = 600) summary(pairwise_trimmed_boot, digits = 6) # digits defaults to max(4, getOption("digits")) ## End(Not run)
omnibus_LSM <- welchADF.test(perceptionData, response = "y", between.s = "Group") omnibus_trimmed <- update(omnibus_LSM, trimming = TRUE) pairwise_LSM <- update(omnibus_LSM, effect = "Group", contrast = "all.pairwise") pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE, effect.size = TRUE) summary(omnibus_LSM) ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap = TRUE, seed = 12345, numsim_b = 600) summary(pairwise_trimmed_boot, digits = 6) # digits defaults to max(4, getOption("digits")) ## End(Not run)
summary
and print
methods for class "welchADF"
## S3 method for class 'welchADFt' summary(object, verbose = FALSE, digits = max(4, getOption("digits")), ...) ## S3 method for class 'welchADFt' print(x, digits = getOption("digits"), ...)
## S3 method for class 'welchADFt' summary(object, verbose = FALSE, digits = max(4, getOption("digits")), ...) ## S3 method for class 'welchADFt' print(x, digits = getOption("digits"), ...)
object |
an object of class "welchADFt" returned by a call to
|
verbose |
whether the summary table should be preceded by an extended description of the welch ADF test performed, including the kind of means (trimmed or not), kind of test (omnibus or pairwise) and whether bootstrap was used or not. Defaults to FALSE. |
digits |
the number of significant digits to use when printing. |
... |
further arguments passed to or from other methods. |
x |
an object of class "welchADFt". |
a string which summarizes the info, with class "summary.welchADFt"
ready to be printed by specific method print.summary.welchADFt
.
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.
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, ...)
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, ...)
formula |
A data frame or a formula or an lm object returned by |
... |
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 |
response |
A string or vector of strings with the name(s) of the column(s) of |
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 |
effect |
If |
correction |
The type of p-value correction when applying multiple pariwise comparisons (i.e. when |
trimming |
Boolean to control the use of robust estimators.
|
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. |
numsim_b |
If |
effect.size |
Boolean; whether effect size estimates should be computed. |
numsim_es |
Positive integer defining the number of bootstrap samples to generate a CI for the effect size estimate (defaults to 999). Ignored if |
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. |
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 - |
seed |
Initial seed (positive integer) for the (R default) random number generator. Defaults to NULL (seed taken from the current time). |
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
formula
: Specialized method that accepts a formula
lm
: Specialized method that accepts a linear model object of class lm
aov
: Specialized method that accepts an Analysis of Variance Model of class aov
lmer
: Specialized method that accepts a Linear Mixed-Effects Model of class lmer
default
: Default method that accepts a data.frame (see argument formula
) and where
factors are passed as strings corresponding to column names
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.
print.welchADFt
, summary.welchADFt
p.adjust.methods
perceptionData
adhdData
adhdData2
womenStereotypeData
miceData
# 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)
# 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)
An artificial dataset created by Lix et al. from summary data presented by Wicherts et al. (2005) (see the vignette). These authors examined the effects of stereotype threat on women's mathematics ability. Originally there were four different tests administered to study participants (arithmetic, number series, word problems, and sums tests). The dataset contains only scores for the arithmetic test because these scores exhibited a greater magnitude of variance heterogeneity than scores for the other tests. It is an unbalanced design with cell sizes ranging from 45 to 50 participants, and a total sample size of 283.
womenStereotypeData
womenStereotypeData
A data frame with 283 rows and 3 variables:
test condition (control, nullified, stereotype threat)
the individual's sex (male, female)
score on the arithmetic test, out of 40
J.Wicherts, C. Dolan, and D. Hessen. Stereotype threat and group differences in test performance: a question of measurement invariance. Journal of Personality and Social Psychology, 89(5):696-716, 2005.
omnibus_LSM <- welchADF.test(womenStereotypeData, response = "y", between.s = c("condition", "sex"), contrast = "omnibus", effect = "condition") omnibus_trimmed <- update(omnibus_LSM, trimming = TRUE, effect = NULL) # unset value of "effect" pairwise_LSM <- update(omnibus_LSM, contrast = "all.pairwise", effect = c("condition", "sex")) pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE) summary(omnibus_LSM) summary(omnibus_trimmed) summary(pairwise_trimmed) ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap = TRUE, seed = 12345) summary(pairwise_trimmed_boot) ## End(Not run)
omnibus_LSM <- welchADF.test(womenStereotypeData, response = "y", between.s = c("condition", "sex"), contrast = "omnibus", effect = "condition") omnibus_trimmed <- update(omnibus_LSM, trimming = TRUE, effect = NULL) # unset value of "effect" pairwise_LSM <- update(omnibus_LSM, contrast = "all.pairwise", effect = c("condition", "sex")) pairwise_trimmed <- update(pairwise_LSM, trimming = TRUE) summary(omnibus_LSM) summary(omnibus_trimmed) summary(pairwise_trimmed) ## Not run: pairwise_trimmed_boot <- update(pairwise_trimmed, bootstrap = TRUE, seed = 12345) summary(pairwise_trimmed_boot) ## End(Not run)