| Title: | A Collection of Robust Statistical Methods |
|---|---|
| Description: | A collection of robust statistical methods based on Wilcox' WRS functions. It implements robust t-tests (independent and dependent samples), robust ANOVA (including between-within subject designs), quantile ANOVA, robust correlation, robust mediation, and nonparametric ANCOVA models based on robust location measures. |
| Authors: | Patrick Mair [cre, aut], Rand Wilcox [aut], Indrajeet Patil [ctb] |
| Maintainer: | Patrick Mair <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1-6 |
| Built: | 2024-10-16 05:16:16 UTC |
| Source: | https://github.com/r-forge/psychor |
This function computes robust ANCOVA for 2 independent groups and one covariate. It compares trimmed means. No parametric assumption (e.g. homogeneity) is made about the form of the regression lines. A running interval smoother is used. A bootstrap version which computes confidence intervals using a percentile t-bootstrap is provided as well.
ancova(formula, data, tr = 0.2, fr1 = 1, fr2 = 1, pts = NA, ...) ancboot(formula, data, tr = 0.2, nboot = 599, fr1 = 1, fr2 = 1, pts = NA, ...)ancova(formula, data, tr = 0.2, fr1 = 1, fr2 = 1, pts = NA, ...) ancboot(formula, data, tr = 0.2, nboot = 599, fr1 = 1, fr2 = 1, pts = NA, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
fr1 |
values of the span for the first group (1 means unspecified) |
fr2 |
values of the span for the second group (1 means unspecified) |
pts |
can be used to specify the design points where the regression lines are to be compared; if |
nboot |
number of bootstrap samples |
... |
currently ignored. |
Returns an object of class ancova containing:
evalpts |
covariate values (including points close to these values) where the test statistic is evaluated |
n1 |
number of subjects at evaluation point (first group) |
n2 |
number of subjects at evaluation point (first group) |
trDiff |
trimmed mean differences |
se |
standard errors for trimmed mean differences |
ci.low |
lower confidence limit for trimmed mean differences |
ci.hi |
upper confidence limit for trimmed mean differences |
test |
values of the test statistic |
crit.vals |
critical values |
p.vals |
p-values |
fitted.values |
fitted values from interval smoothing |
call |
function call |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
head(invisibility) ancova(mischief2 ~ cloak + mischief1, data = invisibility) ## specifying covariate evaluation points ancova(mischief2 ~ cloak + mischief1, data = invisibility, pts = c(3, 4, 8, 1)) ## bootstrap version ancboot(mischief2 ~ cloak + mischief1, data = invisibility)head(invisibility) ancova(mischief2 ~ cloak + mischief1, data = invisibility) ## specifying covariate evaluation points ancova(mischief2 ~ cloak + mischief1, data = invisibility, pts = c(3, 4, 8, 1)) ## bootstrap version ancboot(mischief2 ~ cloak + mischief1, data = invisibility)
binband compares two independent variables in terms of their probability function. discANOVA Tests the global hypothesis that for two or more independent groups, the corresponding discrete distributions are identical. That is, test the hypothesis that independent groups have identical multinomial distributions. discmcp provides multiple comparisons for J independent groups having discrete distributions. discstep implements the step-down multiple comparison procedure for comparing J independent discrete random variables.
binband(x, y, KMS = FALSE, alpha = 0.05, ADJ.P = FALSE, ...) discANOVA(formula, data, nboot = 500, ...) discmcp(formula, data, alpha = 0.05, nboot = 500, ...) discstep(formula, data, nboot = 500, alpha = 0.05, ...)binband(x, y, KMS = FALSE, alpha = 0.05, ADJ.P = FALSE, ...) discANOVA(formula, data, nboot = 500, ...) discmcp(formula, data, alpha = 0.05, nboot = 500, ...) discstep(formula, data, nboot = 500, alpha = 0.05, ...)
x |
an numeric vector of data values for group 1. |
y |
an numeric vector of data values for group 2. |
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
nboot |
number of bootstrap samples. |
alpha |
alpha level. |
KMS |
whether the Kulinskaya-Morgenthaler-Staudte method for comparing binomials should be used. |
ADJ.P |
whether the critical p-value should be adjusted to control FWE when the sample size is small (<50) |
... |
currently ignored. |
discANOVA returns an object of class "med1way" containing:
test |
value of the test statistic |
crit.val |
critical value |
p.value |
p-value |
call |
function call |
The remaining functions return an object of class "robtab" containing:
partable |
parameter table |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
Kulinskaya, E., Morgenthaler, S. and Staudte, R. (2010). Variance stabilizing the difference of two binomial proportions. American Statistician, 64, p. 350-356.
## Consider a study aimed at comparing two methods for reducing shoulder pain after surgery. ## For the first method, the shoulder pain measures are: x1 <- c(2, 4, 4, 2, 2, 2, 4, 3, 2, 4, 2, 3, 2, 4, 3, 2, 2, 3, 5, 5, 2, 2) ## and for the second method they are: x2 <- c(5, 1, 4, 4, 2, 3, 3, 1, 1, 1, 1, 2, 2, 1, 1, 5, 3, 5) fit1 <- binband(x1, x2) fit1 fit2 <- binband(x1, x2, KMS = TRUE, alpha = 0.01) fit2 ## More than two groups: discANOVA(libido ~ dose, viagra, nboot = 200) ## Multiple comparisons: discmcp(libido ~ dose, viagra) discstep(libido ~ dose, viagra)## Consider a study aimed at comparing two methods for reducing shoulder pain after surgery. ## For the first method, the shoulder pain measures are: x1 <- c(2, 4, 4, 2, 2, 2, 4, 3, 2, 4, 2, 3, 2, 4, 3, 2, 2, 3, 5, 5, 2, 2) ## and for the second method they are: x2 <- c(5, 1, 4, 4, 2, 3, 3, 1, 1, 1, 1, 2, 2, 1, 1, 5, 3, 5) fit1 <- binband(x1, x2) fit1 fit2 <- binband(x1, x2, KMS = TRUE, alpha = 0.01) fit2 ## More than two groups: discANOVA(libido ~ dose, viagra, nboot = 200) ## Multiple comparisons: discmcp(libido ~ dose, viagra) discstep(libido ~ dose, viagra)
In the TV show "I'm a celebrity, get me out of here" the celebrities had to eat things like stick insects, fish eyes, etc. This dataset records the time taken to retch when eating these things.
bushbush
A data frame with 5 variables and 8 observations:
participantparticipant ID
stick_insecttime taken to retch when eating a stick insect
kangaroo_testicletime taken to retch when eating a kangaroo testicle
fish_eyetime taken to retch when eating a fish eye
witchetty_grubtime taken to retch when eating a witchetty grub
Dataset from Field et al. book (p. 557).
Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. Sage.
bush summary(bush)bush summary(bush)
The bwtrim function computes a two-way between-within subjects ANOVA on the trimmed means. It is designed for one between-subjects variable and one within-subjects variable. The functions sppba, sppbb, and sppbi compute the main fixed effect, the main within-subjects effect, and the interaction effect only, respectively, using bootstrap. For these 3 functions the user can choose
an M-estimator for group comparisons.
bwtrim(formula, id, data, tr = 0.2, ...) tsplit(formula, id, data, tr = 0.2, ...) sppba(formula, id, data, est = "mom", avg = TRUE, nboot = 500, MDIS = FALSE, ...) sppbb(formula, id, data, est = "mom", nboot = 500, ...) sppbi(formula, id, data, est = "mom", nboot = 500, ...)bwtrim(formula, id, data, tr = 0.2, ...) tsplit(formula, id, data, tr = 0.2, ...) sppba(formula, id, data, est = "mom", avg = TRUE, nboot = 500, MDIS = FALSE, ...) sppbb(formula, id, data, est = "mom", nboot = 500, ...) sppbi(formula, id, data, est = "mom", nboot = 500, ...)
formula |
an object of class formula. |
id |
subject ID. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
est |
Estimate to be used for the group comparisons: either |
avg |
If |
nboot |
number of bootstrap samples. |
MDIS |
if |
... |
currently ignored. |
The tsplit function is doing exactly the same thing as bwtrim. It is kept in the package in order to be consistent with older versions of the Wilcox (2012) book.
For sppba, sppbb, and sppbi the analysis is carried out on the basis of all pairs of difference scores. The null hypothesis is that all such differences have a robust location value of zero. In the formula interface it is required to specify full model.
bwtrim returns an object of class "bwtrim" containing:
Qa |
first main effect |
A.p.value |
p-value first main effect |
A.df |
df F-distribution first main effect |
Qb |
second main effect |
B.p.value |
p-value second main effect |
B.df |
df F-distribution second main effect |
Qab |
interaction effect |
AB.p.value |
p-value interaction effect |
AB.df |
df F-distribution interaction |
call |
function call |
varnames |
variable names |
sppba, sppbb, and sppbi returns an object of class "spp" containing:
test |
value of the test statistic |
p.value |
p-value |
contrasts |
contrasts matrix |
Wilcox, R. (2017). Introduction to Robust Estimation and Hypothesis Testing (4th ed.). Elsevier.
## data need to be on long format pictureLong <- reshape(picture, direction = "long", varying = list(3:4), idvar = "case", timevar = c("pictype"), times = c("couple", "alone")) pictureLong$pictype <- as.factor(pictureLong$pictype) colnames(pictureLong)[4] <- "friend_requests" ## 2-way within-between subjects ANOVA bwtrim(friend_requests ~ relationship_status*pictype, id = case, data = pictureLong) ## between groups effect only (MOM estimator) sppba(friend_requests ~ relationship_status*pictype, case, data = pictureLong) ## within groups effect only (MOM estimator) sppbb(friend_requests ~ relationship_status*pictype, case, data = pictureLong) ## interaction effect only (MOM estimator) sppbi(friend_requests ~ relationship_status*pictype, case, data = pictureLong)## data need to be on long format pictureLong <- reshape(picture, direction = "long", varying = list(3:4), idvar = "case", timevar = c("pictype"), times = c("couple", "alone")) pictureLong$pictype <- as.factor(pictureLong$pictype) colnames(pictureLong)[4] <- "friend_requests" ## 2-way within-between subjects ANOVA bwtrim(friend_requests ~ relationship_status*pictype, id = case, data = pictureLong) ## between groups effect only (MOM estimator) sppba(friend_requests ~ relationship_status*pictype, case, data = pictureLong) ## within groups effect only (MOM estimator) sppbb(friend_requests ~ relationship_status*pictype, case, data = pictureLong) ## interaction effect only (MOM estimator) sppbi(friend_requests ~ relationship_status*pictype, case, data = pictureLong)
Originally from pepperjoe.com, this dataset contains the name, length, and heat of chiles. Heat is measured on a scale from 0-11. (0-2 ... for sissys, 3-4 ... sort of hot, 5-6 ... fairly hot, 7-8 ... real hot, 9.5-9 ... torrid, 9.5-11 ... nuclear).
chilechile
A data frame with 3 variables and 85 observations:
namename of the chile
lengthlength in cm
heatheat of the chile
Wright, D. B., & London, K. (2009). Modern Regression Techniques Using R. Sage.
summary(chile)summary(chile)
Weight loss is studied for three different types of diets.
bushbush
A data frame with 7 variables and 78 observations:
gendergender
ageage
heightbody height
diet.typethree types of diet
initial.weightinitial weight before diet
final.weightfinal weight after diet
weight.lossweight loss
Couturier, D. L., Nicholls, R., and Fernandes, M. (2018). ANOVA with R: analysis of the diet dataset. Retrieved online.
str(diet)str(diet)
These data are based on an educational TV show for children called "The Electric Company". In each of four grades, the classes were randomized into treated (TV show) and control groups (no TV show). At the beginning and at the end of the school year, students in all the classes were given a reading test. The average test scores per class were recorded.
electricelectric
A data frame with 5 variables and 192 observations:
CityFresno and Youngstown
Gradegrade
Pretestreading scores at the beginning of the semester
Posttestreading scores at the end of the semester
Grouptreatment vs. control
Gelman, A., & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press: New York, NY.
summary(electric)summary(electric)
This study looked at the effects of two forms of written corrective feedback on lexico-grammatical accuracy in the academic writing of English as a foreign language university students. It had a 3x4 within-by-between design with three groups (two treatment and one control) measured over four occasions (pretest/treatment, treatment, posttest, delayed posttest).
bushbush
A data frame with 4 variables and 120 observations:
idparticipant ID
groupcontrol, direct, indirect
essayfour measurement occasions
errorRatioerror ratio
McGrath, D. (2016). The Effects of Comprehensive Direct and Indirect Written Corrective Feedback on Accuracy in English as a Foreign Language Students' Writing (Unpublished master's thesis). Macquarie University, Sydney, Australia.
head(essays) summary(essays)head(essays) summary(essays)
Contains various team stats from five European soccer leagues (2008/09 season).
eurosoccereurosoccer
A data frame with 11 variables and 96 teams:
LeagueCountry
TeamTeam
GamesNumber of games
WonGames won
TiedGames tied
LostGames lost
GoalsScoredGoals scored
GoalsConcededGoals conceded
GoalDifferenceGoal difference
PointsFinal amount of points
GoalsGameGoal scored per game
head(eurosoccer)head(eurosoccer)
This dataset is about the effects of alcohol on mate selection in night-clubs. The hypothesis is that after alcohol had been consumed, subjective perceptions of physical attractiveness would become more inaccurate (beer-goggles effect). There are 48 participants: 24 males, 24 females. The reseacher took 3 groups of 8 participants to a night club. One group got no alcohol, one group 2 pints, and one group 4 pints. At the end of the evening the researcher took a photograph of the person the participant was chatting up. The attractiveness of the person on the photo was then evaluated by independent judges.
gogglesgoggles
A data frame with 3 variables and 48 observations:
gender24 male, 24 female students
alcoholamount of alcohol consumed
attractivenessattractiveness rating (0-100)
Dataset from Field et al. book (p. 501).
Field, A., Miles, J, & Field, Z. (2012). Discovering Statistics Using R. Sage.
goggles summary(goggles)goggles summary(goggles)
In a study on the effect of consuming alcohol, hangover symptoms were measured for two independent groups, with each subject consuming alcohol and being measured on three different occasions. One group consisted of sons of alcoholics and the other was a control group.
hangoverhangover
A data frame with 4 variables and 120 observations:
symptomsnumber of hangover symptoms
groupson of alcoholic vs. control
timemeasurement occasion
idsubject ID
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
summary(hangover)summary(hangover)
We are interested in the effect that wearing a cloak of invisibility has on people's tendency to mischief. 80 participants were placed in an enclosed community. Hidden cameras recorded mischievous acts. It was recorded how many mischievous acts were conducted in the first 3 weeks (mischief1). After 3 weeks 34 participants were told that the cameras were switched off so that no one would be able to see what they're getting up to. The remaining 46 subjects were given a cloak of invisibility. These people were told not to tell anyone else about their cloak and they could wear it whenever they liked. The number of mischievous acts were recorded over the next 3 weeks (mischief2).
invisibilityinvisibility
A data frame with 3 variables and 80 observations:
cloakfactor with 34 subjects in the no cloak condition, 46 in the cloak condition
mischief1number of mischievous acts during the first 3 weeks
mischief2number of mischievous acts during the second 3 weeks
Fictional dataset from Field et al. book (p. 485).
Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. Sage.
invisibility summary(invisibility)invisibility summary(invisibility)
TIn this dataset (n = 92) the relationship between how girls were raised by there own mother and their later feelings of maternal self-efficacy, i.e. our belief in our ability to succeed in specific situations. The other variable is self-esteem (can act as mediator). All variables are scored on a continuous scale from 1 to 4.
LeerkesLeerkes
A data frame with 3 variables and 92 observations:
MatCarematernal care
Efficacymaternal self-efficacy
Esteemself-esteen
Leerkes, E.M. & Crockenberg, S.C. (2002). The development of maternal self-efficacy and its impact on maternal behavior. Infancy, 3, 227-247.
Howell, D.C. (2012). Statistical Methods for Psychology (8th edition). Wadsworth, Belmont, CA.
summary(Leerkes)summary(Leerkes)
This function computes a one-way ANOVA for the medians. Homoscedasticity assumption not required. There shouldn't be too many ties.
med1way(formula, data, iter = 1000, ...)med1way(formula, data, iter = 1000, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
iter |
number of iterations to determine critical value. |
... |
currently ignored. |
Evaluating the test statistic using the df proposed in the literature can result in the actual level being less than the nominal level, (i.e., around 0.02-0.025 when testing at the 0.05 level and the sample size is small). A better strategy is to simulate the critical value and computing the p-value accordingly, as implemented in this function.
Returns an object of class med1way containing:
test |
F-value of the test statistic |
crit.val |
critical value |
p.value |
p-value |
call |
function call |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
set.seed(123) med1way(libido ~ dose, data = viagra, iter = 3000)set.seed(123) med1way(libido ~ dose, data = viagra, iter = 3000)
This function computes a two-way ANOVA medians with interactions effects.
med2way(formula, data, ...)med2way(formula, data, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
... |
currently ignored. |
There should not be too many ties in the data. The test statistics for the main effects in med2way are F-distributed, the (heteroscedastic) test for the interaction is chi-square distributed. Post hoc tests can be performed using
mcp2a.
Returns an object of class t2way containing:
Qa |
first main effect |
A.p.value |
p-value first main effect |
Qb |
second main effect |
B.p.value |
p-value second main effect |
Qab |
interaction effect |
AB.p.value |
p-value interaction effect |
call |
function call |
varnames |
variable names |
dim |
design dimensions |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
med2way(attractiveness ~ gender*alcohol, data = goggles) mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median")med2way(attractiveness ~ gender*alcohol, data = goggles) mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median")
Participants are randomly assigned to one of two groups. The first group watches a violent film, and the other watches a nonviolent film. Afterwards, the aggressive affect is measured, and it is desired to compare three groups, taking gender and degress into account as well.
moviemovie
A data frame with 4 variables and 68 observations:
degreeno degree vs. degree
gender36 males, 32 females
typeviolend vs. nonviolent
aggressiveaggressive affect
Artificial dataset from Wilcox book (p. 316).
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
movie summary(movie)movie summary(movie)
Tests whether a robust location measure (median, Huber Psi) differs from a null value and reports a 95% confidence interval based on percentile bootstrap.
onesampb(x, est = "onestep", nboot = 2000, nv = 0, alpha = 0.05, ...)onesampb(x, est = "onestep", nboot = 2000, nv = 0, alpha = 0.05, ...)
x |
a numeric vector. |
est |
robust estimator to be used ( |
nboot |
number of bootstrap samples. |
nv |
value under H0. |
alpha |
alpha level. |
... |
currently ignored. |
Returns an object of class "onesampb" containing:
ci |
95% confidence interval |
estimate |
robust location sample estimate |
p.value |
p-value |
n |
number of effective observations |
call |
function call |
Wilcox, R. (2017). Introduction to Robust Estimation and Hypothesis Testing (4th ed.). Elsevier.
set.seed(123) x <- rnorm(30) onesampb(x, nboot = 100) ## H0: Psi = 0 set.seed(123) x <- rlnorm(30) onesampb(x, est = "median", nv = 1) ## H0: median = 1set.seed(123) x <- rnorm(30) onesampb(x, nboot = 100) ## H0: Psi = 0 set.seed(123) x <- rlnorm(30) onesampb(x, est = "median", nv = 1) ## H0: median = 1
The pbcor function computes the percentage bend correlation coefficient, wincor the Winsorized correlation,
pball the percentage bend correlation matrix, winall the Winsorized correlation matrix.
pbcor(x, y = NULL, beta = 0.2, ci = FALSE, nboot = 500, alpha = 0.05, ...) pball(x, beta = 0.2, ...) wincor(x, y = NULL, tr = 0.2, ci = FALSE, nboot = 500, alpha = 0.05, ...) winall(x, tr = 0.2, ...)pbcor(x, y = NULL, beta = 0.2, ci = FALSE, nboot = 500, alpha = 0.05, ...) pball(x, beta = 0.2, ...) wincor(x, y = NULL, tr = 0.2, ci = FALSE, nboot = 500, alpha = 0.05, ...) winall(x, tr = 0.2, ...)
x |
a numeric vector, a matrix or a data frame. |
y |
a second numeric vector (for correlation functions). |
beta |
bending constant. |
tr |
amount of Winsorization. |
ci |
whether boostrap CI should be computed or not. |
nboot |
number of bootstrap samples for CI computation. |
alpha |
alpha level for CI computation. |
... |
currently ignored. |
It tested is whether the correlation coefficient equals 0 (null hypothesis) or not. Missing values are deleted pairwise. The tests are sensitive to heteroscedasticity. The test statistic H in pball tests the hypothesis that all correlations are equal to zero.
pbcor and wincor return an object of class "pbcor" containing:
cor |
robust correlation coefficient |
test |
value of the test statistic |
p.value |
p-value |
n |
number of effective observations |
cor_ci |
bootstrap confidence interval |
call |
function call |
pball and winall return an object of class "pball" containing:
pbcorm |
robust correlation matrix |
p.values |
p-values |
H |
H-statistic |
H.p.value |
p-value H-statistic |
cov |
variance-covariance matrix |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
x1 <- subset(hangover, subset = (group == "control" & time == 1))$symptoms x2 <- subset(hangover, subset = (group == "control" & time == 2))$symptoms pbcor(x1, x2) pbcor(x1, x2, beta = 0.1, ci = TRUE) wincor(x1, x2) wincor(x1, x2, tr = 0.1, ci = TRUE) require(reshape) hanglong <- subset(hangover, subset = group == "control") hangwide <- cast(hanglong, id ~ time, value = "symptoms")[,-1] pball(hangwide) winall(hangwide)x1 <- subset(hangover, subset = (group == "control" & time == 1))$symptoms x2 <- subset(hangover, subset = (group == "control" & time == 2))$symptoms pbcor(x1, x2) pbcor(x1, x2, beta = 0.1, ci = TRUE) wincor(x1, x2) wincor(x1, x2, tr = 0.1, ci = TRUE) require(reshape) hanglong <- subset(hangover, subset = group == "control") hangwide <- cast(hanglong, id ~ time, value = "symptoms")[,-1] pball(hangwide) winall(hangwide)
This dataset examines how the profile pictures on social network platforms affect the number of friend requests when females are in a relationship. The relationship status is a between-subject variable (part of the participants did set their status to relationship). For the first 3 weeks the subjects had a picture of their own in their profiles. For the following 3 weeks they posted a picture with a man.
picturepicture
A data frame with 4 variables and 40 observations:
casesubject id
relationship_statusRelationship status on social network platform
coupleamount of friend requests when profile picture as couple
aloneamount of friend requests when profile picture as single
Dataset from Field et al. book (p. 644).
Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. Sage.
picture summary(picture)picture summary(picture)
The Pygmalion effect is the phenomenon where higher expectations lead to an increase in performance. For instance, when teachers expect students to do well and show intellectual growth, they do; when teachers do not have such expectations, performance and growth are not so encouraged and may in fact be discouraged in a variety of ways. This dataset contains reasoning IQ scores of children. For the experimental group, positive expectancies had been suggested to teachers after the pretest. For the experimental group, no expectancies had been suggested after the pretest. For both groups we have reasoning IQ posttest scores. The dataset is taken from Elashoff and Snow (1970).
PygmalionPygmalion
A data frame with 3 variables and 114 observations:
Pretestpretest score
Posttestposttest score
Grouptreatment vs. control
Elashoff, J. D., & Snow, R. E. (1970). A case study in statistical inference: Reconsideration of the Rosenthal-Jacobson data on teacher expectancy. Technical Report No. 15, School of Education, Stanford University.
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
summary(Pygmalion)summary(Pygmalion)
One-way ANOVA based on quantiles. Only known method to work well when tied values are likely to occur.
Qanova(formula, data, q = 0.5, nboot = 600, ...)Qanova(formula, data, q = 0.5, nboot = 600, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
q |
quantile (or vector of quantiles) to be used. |
nboot |
number of bootstrap samples |
... |
currently ignored. |
Test global hypothesis that J independent groups have equal quantiles (default: median) using the Harrell-Davis estimator. Performs well when there are tied values.
Qanova an object of class "qanova" containing:
psihat |
value of the test statistics |
contrasts |
contrasts |
p.value |
p-value |
call |
function call |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
## median comparison set.seed(123) fitmed <- Qanova(libido ~ dose, viagra, nboot = 200) fitmed ## 1st, 3rd quartile comparison set.seed(123) fitquart <- Qanova(libido ~ dose, viagra, q = c(0.25, 0.75), nboot = 200) fitquart## median comparison set.seed(123) fitmed <- Qanova(libido ~ dose, viagra, nboot = 200) fitmed ## 1st, 3rd quartile comparison set.seed(123) fitquart <- Qanova(libido ~ dose, viagra, q = c(0.25, 0.75), nboot = 200) fitquart
The rmanova function computes a one-way repeated measures ANOVA for the trimmed means. Homoscedasticity assumption not required. Corresponding post hoc tests can be performed using rmmcp.
rmanova(y, groups, blocks, tr = 0.2, ...) rmmcp(y, groups, blocks, tr = 0.2, alpha = 0.05, ...)rmanova(y, groups, blocks, tr = 0.2, ...) rmmcp(y, groups, blocks, tr = 0.2, alpha = 0.05, ...)
y |
a numeric vector of data values (response). |
groups |
a vector giving the group of the corresponding elements of y. |
blocks |
a vector giving the block of the corresponding elements of y. |
tr |
trim level for the mean. |
alpha |
alpha level for post hoc comparisons. |
... |
currently ignored. |
rmanova an object of class "t1way" containing:
test |
value of the test statistic |
df1 |
degrees of freedom |
df2 |
degrees of freedom |
p.value |
p-value |
call |
function call |
rmmcp returns an object of class "mcp1" containing:
comp |
inference for all pairwise comparisons |
fnames |
names of the factor levels |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
head(WineTasting) rmanova(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster) ## post hoc rmmcp(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster) head(bush) require(reshape) bushLong <- melt(bush, id.var = "participant", variable_name = "food") rmanova(bushLong$value, bushLong$food, bushLong$participant) ## post hoc rmmcp(bushLong$value, bushLong$food, bushLong$participant)head(WineTasting) rmanova(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster) ## post hoc rmmcp(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster) head(bush) require(reshape) bushLong <- melt(bush, id.var = "participant", variable_name = "food") rmanova(bushLong$value, bushLong$food, bushLong$participant) ## post hoc rmmcp(bushLong$value, bushLong$food, bushLong$participant)
The rmanovab function computes a bootstrap version of the one-way repeated measures ANOVA for the trimmed means. Homoscedasticity assumption not required. Corresponding post hoc tests can be performed using pairdepb.
rmanovab(y, groups, blocks, tr = 0.2, nboot = 599, ...) pairdepb(y, groups, blocks, tr = 0.2, nboot = 599, ...)rmanovab(y, groups, blocks, tr = 0.2, nboot = 599, ...) pairdepb(y, groups, blocks, tr = 0.2, nboot = 599, ...)
y |
a numeric vector of data values (response). |
groups |
a vector giving the group of the corresponding elements of y. |
blocks |
a vector giving the block of the corresponding elements of y. |
tr |
trim level for the mean. |
nboot |
number of bootstrap samples. |
... |
currently ignored. |
rmanovab an object of class "rmanovab" containing:
test |
value of the test statistic |
crit |
critical value |
call |
function call |
pairdepb returns an object of class "mcp2" containing:
comp |
inference for all pairwise comparisons |
fnames |
names of the factor levels |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
set.seed(123) rmanovab(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster, nboot = 300) ## post hoc set.seed(123) pairdepb(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster, nboot = 300)set.seed(123) rmanovab(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster, nboot = 300) ## post hoc set.seed(123) pairdepb(WineTasting$Taste, WineTasting$Wine, WineTasting$Taster, nboot = 300)
The runmean implements a running interval smoother on the trimmed mean, rungen uses general M-estimators, runmbo performs
interval smoothing on M-estimators with bagging.
runmean(x, y, fr = 1, tr = 0.2, ...) rungen(x, y, fr = 1, est = "mom", ...) runmbo(x, y, fr = 1, est = "mom", nboot = 40, ...)runmean(x, y, fr = 1, tr = 0.2, ...) rungen(x, y, fr = 1, est = "mom", ...) runmbo(x, y, fr = 1, est = "mom", nboot = 40, ...)
x |
a numeric vector of data values (predictor) |
y |
a numeric vector of data values (response) |
fr |
smoothing factor (see details) |
tr |
trim level for the mean |
est |
type of M-estimator ( |
nboot |
number of bootstrap samples |
... |
currently ignored. |
The larger the smoothing factor, the stronger the smoothing. Often the choice fr = 1 gives good results; the general strategy is
to find the smallest constant so that the plot looks reasonably smooth.
Returns the fitted values.
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
## trimmed mean smoother fitmean <- runmean(Pygmalion$Pretest, Pygmalion$Posttest) ## MOM smoother fitmest <- rungen(Pygmalion$Pretest, Pygmalion$Posttest) ## median smoother fitmed <- rungen(Pygmalion$Pretest, Pygmalion$Posttest, est = "median") ## bagged onestep smoother fitbag <- runmbo(Pygmalion$Pretest, Pygmalion$Posttest, est = "onestep") ## plot smoothers plot(Pygmalion$Pretest, Pygmalion$Posttest, col = "gray", xlab = "Pretest", ylab = "Posttest", main = "Pygmalion Smoothing") orderx <- order(Pygmalion$Pretest) lines(Pygmalion$Pretest[orderx], fitmean[orderx], lwd = 2) lines(Pygmalion$Pretest[orderx], fitmest[orderx], lwd = 2, col = 2) lines(Pygmalion$Pretest[orderx], fitmed[orderx], lwd = 2, col = 3) lines(Pygmalion$Pretest[orderx], fitbag[orderx], lwd = 2, col = 4) legend("topleft", legend = c("Trimmed Mean", "MOM", "Median", "Bagged Onestep"), col = 1:4, lty = 1)## trimmed mean smoother fitmean <- runmean(Pygmalion$Pretest, Pygmalion$Posttest) ## MOM smoother fitmest <- rungen(Pygmalion$Pretest, Pygmalion$Posttest) ## median smoother fitmed <- rungen(Pygmalion$Pretest, Pygmalion$Posttest, est = "median") ## bagged onestep smoother fitbag <- runmbo(Pygmalion$Pretest, Pygmalion$Posttest, est = "onestep") ## plot smoothers plot(Pygmalion$Pretest, Pygmalion$Posttest, col = "gray", xlab = "Pretest", ylab = "Posttest", main = "Pygmalion Smoothing") orderx <- order(Pygmalion$Pretest) lines(Pygmalion$Pretest[orderx], fitmean[orderx], lwd = 2) lines(Pygmalion$Pretest[orderx], fitmest[orderx], lwd = 2, col = 2) lines(Pygmalion$Pretest[orderx], fitmed[orderx], lwd = 2, col = 3) lines(Pygmalion$Pretest[orderx], fitbag[orderx], lwd = 2, col = 4) legend("topleft", legend = c("Trimmed Mean", "MOM", "Median", "Bagged Onestep"), col = 1:4, lty = 1)
24 arachnophobes were used in all. 12 were asked to play with a big hairy tarantula spider with big fangs and an evil look. Their subsequent anxiety was measured. The remaining 12 were shown only picture of the same hairy tarantula. Again, the anxiety was measured.
spiderspider
A data frame with 2 variables and 24 observations:
Grouppicture vs. real spider
Anxietyanxiety measure
Dataset from Field et al. book (p. 362).
Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. Sage.
spiderspider
At a swimming team practice, all participants were asked to swim their best event as far as possible, but in each case the time that was reported was falsified to indicate poorer than expected performance (i.e., each swimmer was disappointed). 30 min later, they did the same performance. The authors predicted that on the second trial more pessimistic swimmers would do worse than on their first trial, whereas optimistic swimmers would do better. The response is ratio = Time1/Time2 (> 1 means that a swimmer did better in trial 2).
swimmingswimming
A data frame with 4 variables and 58 observations:
OptimOptimists and pessimists
SexGender of the swimmer
EventSwimming event: freestyle, breaststroke, backstroke
RatioRatio between the swimming times
Seligman, M. E. P., Nolen-Hoeksema, S., Thornton, N., & Thornton, C. M. (1990). Explanatory style as a mechanism of disappointing athletic performance. Psychological Science, 1, 143-146.
summary(swimming)summary(swimming)
The t1way function computes a one-way ANOVA on trimmed means. Homoscedasticity assumption not required. It uses a generalization of Welch's method. Corresponding post hoc tests can be performed using lincon.
t1way(formula, data, tr = 0.2, alpha = 0.05, nboot = 100, ...) lincon(formula, data, tr = 0.2, alpha = 0.05, method = "hochberg", ...)t1way(formula, data, tr = 0.2, alpha = 0.05, nboot = 100, ...) lincon(formula, data, tr = 0.2, alpha = 0.05, method = "hochberg", ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
alpha |
alpha level for CI computation. |
nboot |
number of bootstrap samples for effect size CI computation. |
method |
method to correct the p-value (see |
... |
currently ignored. |
In the post hoc computations, confidence intervals and p-values are adjusted to control FWE. The default for the p-values is to use Hochberg's 1988 sharper Bonferroni procedure.
t1way returns an object of class "t1way" containing:
test |
value of the test statistic (F-statistic) |
df1 |
degrees of freedom |
df2 |
degrees of freedom |
p.value |
p-value |
effsize |
explanatory measure of effect size |
effsize_ci |
boostrap effect size CI |
call |
function call |
lincon returns an object of class "mcp1" containing:
comp |
inference for all pairwise comparisons |
fnames |
names of the factor levels |
linconv2 returns an object of class "linconv2" containing:
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
set.seed(123) t1way(libido ~ dose, data = viagra) ## post hoc tests lincon(libido ~ dose, data = viagra)set.seed(123) t1way(libido ~ dose, data = viagra) ## post hoc tests lincon(libido ~ dose, data = viagra)
Test the hypothesis of equal trimmed means using a percentile t bootstrap method. Corresponding post hoc tests are provided in mcppb20.
t1waybt(formula, data, tr = 0.2, nboot = 599, ...) mcppb20(formula, data, tr = 0.2, nboot = 599, ...)t1waybt(formula, data, tr = 0.2, nboot = 599, ...) mcppb20(formula, data, tr = 0.2, nboot = 599, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
nboot |
number of bootstrap samples. |
... |
currently ignored. |
Returns an object of class t1waybt containing:
test |
value of the test statistic |
p.value |
p-value |
Var.Explained |
explained amount of variance |
Effect.Size |
effect size |
nboot.eff |
effective number of bootstrap samples |
call |
function call |
mcppb20 returns an object of class "mcp1" containing:
comp |
inference for all pairwise comparisons |
fnames |
names of the factor levels |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
t1waybt(libido ~ dose, data = viagra) ## post hoc mcppb20(libido ~ dose, data = viagra)t1waybt(libido ~ dose, data = viagra) ## post hoc mcppb20(libido ~ dose, data = viagra)
The t2way function computes a two-way ANOVA for trimmed means with interactions effects. Corresponding post hoc tests are in mcp2atm. pbad2way performs a two-way ANOVA using M-estimators for location with mcp2a for post hoc tests.
t2way(formula, data, tr = 0.2, ...) pbad2way(formula, data, est = "mom", nboot = 599, pro.dis = FALSE, ...) mcp2atm(formula, data, tr = 0.2, ...) mcp2a(formula, data, est = "mom", nboot = 599, ...)t2way(formula, data, tr = 0.2, ...) pbad2way(formula, data, est = "mom", nboot = 599, pro.dis = FALSE, ...) mcp2atm(formula, data, tr = 0.2, ...) mcp2a(formula, data, est = "mom", nboot = 599, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
est |
Estimate to be used for the group comparisons: either |
nboot |
number of bootstrap samples. |
pro.dis |
If |
... |
currently ignored. |
t2way does not report any degrees of freedom since an adjusted critical value is used.
pbad2way returns p-values only; if it happens that the variance-covariance matrix in the Mahalanobis distance computation is singular, it is suggested to use the projection distances by setting pro.dis = TRUE.
The functions t2way and pbad2way return an object of class t2way containing:
Qa |
first main effect |
A.p.value |
p-value first main effect |
Qb |
second main effect |
B.p.value |
p-value second main effect |
Qab |
interaction effect |
AB.p.value |
p-value interaction effect |
call |
function call |
varnames |
variable names |
dim |
design dimensions |
The functions mcp2atm and mcp2a return an object of class mcp containing:
effects |
list with post hoc comparisons for all effects |
contrasts |
design matrix |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
## 2-way ANOVA on trimmed means t2way(attractiveness ~ gender*alcohol, data = goggles) ## post hoc tests mcp2atm(attractiveness ~ gender*alcohol, data = goggles) ## 2-way ANOVA on MOM estimator pbad2way(attractiveness ~ gender*alcohol, data = goggles) ## post hoc tests mcp2a(attractiveness ~ gender*alcohol, data = goggles) ## 2-way ANOVA on medians pbad2way(attractiveness ~ gender*alcohol, data = goggles, est = "median") ## post hoc tests mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median") ## extract design matrix model.matrix(mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median"))## 2-way ANOVA on trimmed means t2way(attractiveness ~ gender*alcohol, data = goggles) ## post hoc tests mcp2atm(attractiveness ~ gender*alcohol, data = goggles) ## 2-way ANOVA on MOM estimator pbad2way(attractiveness ~ gender*alcohol, data = goggles) ## post hoc tests mcp2a(attractiveness ~ gender*alcohol, data = goggles) ## 2-way ANOVA on medians pbad2way(attractiveness ~ gender*alcohol, data = goggles, est = "median") ## post hoc tests mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median") ## extract design matrix model.matrix(mcp2a(attractiveness ~ gender*alcohol, data = goggles, est = "median"))
This function computes a three-way ANOVA for trimmed means with all interactions effects.
t3way(formula, data, tr = 0.2, ...)t3way(formula, data, tr = 0.2, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
... |
currently ignored. |
Returns an object of class t3way containing:
Qa |
first main effect |
A.p.value |
p-value first main effect |
Qb |
second main effect |
B.p.value |
p-value second main effect |
Qc |
third main effect |
C.p.value |
p-value third main effect |
Qab |
first two-way interaction effect |
AB.p.value |
p-value first two-way interaction effect |
Qac |
second two-way interaction effect |
AC.p.value |
p-value second two-way interaction effect |
Qbc |
third two-way interaction effect |
BC.p.value |
p-value third two-way interaction effect |
Qabc |
three-way interaction effect |
ABC.p.value |
p-value three-way interaction effect |
call |
function call |
varnames |
variable names |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
t3way(aggressive ~ degree*gender*type, data = movie)t3way(aggressive ~ degree*gender*type, data = movie)
Compute a 1-alpha confidence interval for the trimmed mean using a bootstrap percentile t method.
trimcibt(x, nv = 0, tr = 0.2, alpha = 0.05, nboot = 200, ...)trimcibt(x, nv = 0, tr = 0.2, alpha = 0.05, nboot = 200, ...)
x |
a numeric vector. |
nv |
value under H0. |
tr |
trim level for the mean. |
alpha |
alpha level. |
nboot |
number of bootstrap samples. |
... |
currently ignored. |
Returns an object of class "trimcibt" containing:
ci |
95% confidence interval |
estimate |
trimmed mean |
p.value |
p-value |
test.stat |
t-statistic |
tr |
trimming level |
n |
number of effective observations |
Wilcox, R. (2017). Introduction to Robust Estimation and Hypothesis Testing (4th ed.). Elsevier.
set.seed(123) x <- rnorm(30) trimcibt(x, nboot = 100) ## H0: Psi = 0set.seed(123) x <- rnorm(30) trimcibt(x, nboot = 100) ## H0: Psi = 0
The following functions for estimating robust location measures and their standard errors are provided: winmean
for the Winsorized mean, winse for its se, trimse for the trimmend mean se, msmedse for the median se,
mest for the M-estimator with se in mestse. The functions onestep and mom compute the one-step and
modified one-step (MOM) M-estimator. The Winsorized variance is implemented in winvar.
winmean(x, tr = 0.2, na.rm = FALSE, ...) winvar(x, tr = 0.2, na.rm = FALSE, STAND = NULL, ...) winse(x, tr = 0.2, ...) trimse(x, tr = 0.2, na.rm = FALSE, ...) msmedse(x, sewarn = TRUE, ...) mest(x, bend = 1.28, na.rm = FALSE, ...) mestse(x, bend = 1.28, ...) onestep(x, bend = 1.28, na.rm = FALSE, MED = TRUE, ...) mom(x, bend = 2.24, na.rm = TRUE, ...)winmean(x, tr = 0.2, na.rm = FALSE, ...) winvar(x, tr = 0.2, na.rm = FALSE, STAND = NULL, ...) winse(x, tr = 0.2, ...) trimse(x, tr = 0.2, na.rm = FALSE, ...) msmedse(x, sewarn = TRUE, ...) mest(x, bend = 1.28, na.rm = FALSE, ...) mestse(x, bend = 1.28, ...) onestep(x, bend = 1.28, na.rm = FALSE, MED = TRUE, ...) mom(x, bend = 2.24, na.rm = TRUE, ...)
x |
a numeric vector containing the values whose measure is to be computed. |
tr |
trim lor Winsorizing level. |
na.rm |
a logical value indicating whether NA values should be stripped before the computation proceeds. |
sewarn |
a logical value indicating whether warnings for ties should be printed. |
bend |
bending constant for M-estimator. |
MED |
if |
STAND |
no functionality, kept for WRS compatibility purposes. |
... |
currently ignored. |
The standard error for the median is computed according to McKean and Shrader (1984).
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
McKean, J. W., & Schrader, R. M. (1984). A comparison of methods for studentizing the sample median. Communications in Statistics - Simulation and Computation, 13, 751-773.
Dana, E. (1990). Salience of the self and salience of standards: Attempts to match self to standard. Unpublished PhD thesis, Department of Psychology, University of Southern California.
## Self-awareness data (Dana, 1990): Time persons could keep a portion of an ## apparatus in contact with a specified range. self <- c(77, 87, 88, 114, 151, 210, 219, 246, 253, 262, 296, 299, 306, 376, 428, 515, 666, 1310, 2611) mean(self, 0.1) ## .10 trimmed mean trimse(self, 0.1) ## se trimmed mean winmean(self, 0.1) ## Winsorized mean (.10 Winsorizing amount) winse(self, 0.1) ## se Winsorized mean winvar(self, 0.1) ## Winsorized variance median(self) ## median msmedse(self) ## se median mest(self) ## Huber M-estimator mestse(self) onestep(self) ## one-step M-estimator mom(self) ## modified one-step M-estimator## Self-awareness data (Dana, 1990): Time persons could keep a portion of an ## apparatus in contact with a specified range. self <- c(77, 87, 88, 114, 151, 210, 219, 246, 253, 262, 296, 299, 306, 376, 428, 515, 666, 1310, 2611) mean(self, 0.1) ## .10 trimmed mean trimse(self, 0.1) ## se trimmed mean winmean(self, 0.1) ## Winsorized mean (.10 Winsorizing amount) winse(self, 0.1) ## se Winsorized mean winvar(self, 0.1) ## Winsorized variance median(self) ## median msmedse(self) ## se median mest(self) ## Huber M-estimator mestse(self) onestep(self) ## one-step M-estimator mom(self) ## modified one-step M-estimator
The twopcor function tests whether the difference between two Pearson correlations is 0. The twocor function performs the same test on a robust correlation coefficient (percentage bend correlation or Winsorized correlation).
twopcor(x1, y1, x2, y2, nboot = 599, ...) twocor(x1, y1, x2, y2, corfun = "pbcor", nboot = 599, tr = 0.2, beta = 0.2, ...)twopcor(x1, y1, x2, y2, nboot = 599, ...) twocor(x1, y1, x2, y2, corfun = "pbcor", nboot = 599, tr = 0.2, beta = 0.2, ...)
x1 |
a numeric vector. |
y1 |
a numeric vector. |
x2 |
a numeric vector. |
y2 |
a numeric vector. |
nboot |
number of bootstrap samples. |
corfun |
Either |
tr |
amount of Winsorization. |
beta |
bending constant. |
... |
currently ignored. |
It is tested whether the first correlation coefficient (based on x1 and y1) equals to the second correlation coefficient (based on x2 and y2). Both approaches return percentile bootstrap CIs.
twopcor and twocor return an object of class "twocor" containing:
r1 |
robust correlation coefficient |
r2 |
value of the test statistic |
ci |
confidence interval |
p.value |
p-value |
call |
function call |
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
ct1 <- subset(hangover, subset = (group == "control" & time == 1))$symptoms ct2 <- subset(hangover, subset = (group == "control" & time == 2))$symptoms at1 <- subset(hangover, subset = (group == "alcoholic" & time == 1))$symptoms at2 <- subset(hangover, subset = (group == "alcoholic" & time == 2))$symptoms set.seed(111) twopcor(ct1, ct2, at1, at2) set.seed(123) twocor(ct1, ct2, at1, at2, corfun = "pbcor", beta = 0.15) set.seed(224) twocor(ct1, ct2, at1, at2, corfun = "wincor", tr = 0.15, nboot = 50)ct1 <- subset(hangover, subset = (group == "control" & time == 1))$symptoms ct2 <- subset(hangover, subset = (group == "control" & time == 2))$symptoms at1 <- subset(hangover, subset = (group == "alcoholic" & time == 1))$symptoms at2 <- subset(hangover, subset = (group == "alcoholic" & time == 2))$symptoms set.seed(111) twopcor(ct1, ct2, at1, at2) set.seed(123) twocor(ct1, ct2, at1, at2, corfun = "pbcor", beta = 0.15) set.seed(224) twocor(ct1, ct2, at1, at2, corfun = "wincor", tr = 0.15, nboot = 50)
Participants were assigned randomly to three viagra dosages (placebo, low dosage, high dosage). The dependent variable was an objective measure of libido.
viagraviagra
A data frame with 2 variables and 15 observations:
doseviagra dosage
libidoobjective measure of libido
Artificial dataset from Field et al. book (p. 401).
Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. Sage.
viagraviagra
In this hypothetical dataset we have three types of wine (A, B and C). We asked 22 friends to taste each of the three wines (in a blind fold fashion), and then to give a grade of 1 to 7. We asked them to rate the wines 5 times each, and then averaged their results to give a number for a persons preference for each wine.
WineTastingWineTasting
A data frame with 3 variables and 66 observations:
TasteTaste Rating
WineWine (A, B, C)
TasterTaster (index)
WineTasting summary(WineTasting)WineTasting summary(WineTasting)
Compute an AKP-type effect size for dependent sample ANOVA
wmcpAKP(x, tr = 0.2, nboot = 200, ...)wmcpAKP(x, tr = 0.2, nboot = 200, ...)
x |
data frame in wide format (no missing values allowed) |
tr |
trim level for the means. |
nboot |
number of bootstrap samples. |
... |
currently ignored. |
The computation is based on a modification of the Algina-Keselman-Penfield effect size for J dependent samples.
Algina, J., Keselman, H.J., & Penfield, R.D. (2005). An alternative to Cohen's standardized mean difference effect size: A robust parameter and confidence interval in the two independent groups case. Psychological Methods, 10, 317-328.
## Not run: require(reshape) WineLong <- cast(WineTasting, Taster ~ Wine, value = "Taste")[,-1] set.seed(123) effsize <- wmcpAKP(WineLong, nboot = 20) effsize ## End(Not run)## Not run: require(reshape) WineLong <- cast(WineTasting, Taster ~ Wine, value = "Taste")[,-1] set.seed(123) effsize <- wmcpAKP(WineLong, nboot = 20) effsize ## End(Not run)
The function yuen performs Yuen's test for trimmed means, yuenbt is a bootstrap version of it. akp.effect and yuen.effect.ci can be used for effect size computation. The pb2gen function performs a t-test based on various robust estimators, medpb2 compares two independent groups using medians, and qcomhd compares arbitrary quantiles.
yuen(formula, data, tr = 0.2, ...) yuenbt(formula, data, tr = 0.2, nboot = 599, side = TRUE, ...) akp.effect(formula, data, EQVAR = TRUE, tr = 0.2, nboot = 200, alpha = 0.05, ...) yuen.effect.ci(formula, data, tr = 0.2, nboot = 400, alpha = 0.05, ...) pb2gen(formula, data, est = "mom", nboot = 599, ...) medpb2(formula, data, nboot = 2000, ...) qcomhd(formula, data, q = c(0.1, 0.25, 0.5, 0.75, 0.9), nboot = 2000, alpha = 0.05, ADJ.CI = TRUE, ...)yuen(formula, data, tr = 0.2, ...) yuenbt(formula, data, tr = 0.2, nboot = 599, side = TRUE, ...) akp.effect(formula, data, EQVAR = TRUE, tr = 0.2, nboot = 200, alpha = 0.05, ...) yuen.effect.ci(formula, data, tr = 0.2, nboot = 400, alpha = 0.05, ...) pb2gen(formula, data, est = "mom", nboot = 599, ...) medpb2(formula, data, nboot = 2000, ...) qcomhd(formula, data, q = c(0.1, 0.25, 0.5, 0.75, 0.9), nboot = 2000, alpha = 0.05, ADJ.CI = TRUE, ...)
formula |
an object of class formula. |
data |
an optional data frame for the input data. |
tr |
trim level for the mean. |
nboot |
number of bootstrap samples. |
side |
|
est |
estimate to be used for the group comparisons: either |
q |
quantiles to be used for comparison. |
alpha |
alpha level. |
ADJ.CI |
whether CIs should be adjusted. |
EQVAR |
whether variances are assumed to be equal across groups. |
... |
currently ignored. |
If yuenbt is used, p-value computed only when side = TRUE. medpb2 is just a wrapper function for pb2gen with the median
as M-estimator. It is the only known method to work well in simulations when tied values are likely to occur.qcomhd returns p-values and critical p-values based on Hochberg's method.
Returns objects of classes "yuen" or "pb2" containing:
test |
value of the test statistic (t-statistic) |
p.value |
p-value |
conf.int |
confidence interval |
df |
degress of freedom |
diff |
trimmed mean difference |
effsize |
explanatory measure of effect size |
call |
function call |
qcomhd returns an object of class "robtab" containing:
partable |
parameter table |
Algina, J., Keselman, H.J., & Penfield, R.D. (2005). An alternative to Cohen's standardized mean difference effect size: A robust parameter and confidence interval in the two independent groups case. Psychological Methods, 10, 317-328.
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
Wilcox, R., & Tian, T. (2011). Measuring effect size: A robust heteroscedastic approach for two or more groups. Journal of Applied Statistics, 38, 1359-1368.
Yuen, K. K. (1974). The two sample trimmed t for unequal population variances. Biometrika, 61, 165-170.
set.seed(123) ## Yuen's test yuen(Anxiety ~ Group, data = spider) ## Bootstrap version of Yuen's test (symmetric CIs) yuenbt(Anxiety ~ Group, data = spider) ## Robust Cohen's delta akp.effect(Anxiety ~ Group, data = spider) ## Using an M-estimator pb2gen(Anxiety ~ Group, data = spider, est = "mom") pb2gen(Anxiety ~ Group, data = spider, est = "mean") pb2gen(Anxiety ~ Group, data = spider, est = "median") ## Using the median medpb2(Anxiety ~ Group, data = spider) ## Quantiles set.seed(123) qcomhd(Anxiety ~ Group, data = spider, q = c(0.8, 0.85, 0.9), nboot = 500)set.seed(123) ## Yuen's test yuen(Anxiety ~ Group, data = spider) ## Bootstrap version of Yuen's test (symmetric CIs) yuenbt(Anxiety ~ Group, data = spider) ## Robust Cohen's delta akp.effect(Anxiety ~ Group, data = spider) ## Using an M-estimator pb2gen(Anxiety ~ Group, data = spider, est = "mom") pb2gen(Anxiety ~ Group, data = spider, est = "mean") pb2gen(Anxiety ~ Group, data = spider, est = "median") ## Using the median medpb2(Anxiety ~ Group, data = spider) ## Quantiles set.seed(123) qcomhd(Anxiety ~ Group, data = spider, q = c(0.8, 0.85, 0.9), nboot = 500)
The function yuend performs Yuen's test on trimmed means for dependent samples. Dqcomhd compares the quantiles of the marginal distributions associated with two dependent groups via hd estimator. Tied values are allowed.
dep.effect computes various effect sizes and confidence intervals for two dependent samples (see Details).
yuend(x, y, tr = 0.2, ...) Dqcomhd(x, y, q = c(1:9)/10, nboot = 1000, na.rm = TRUE, ...) dep.effect(x, y, tr = 0.2, nboot = 1000, ...)yuend(x, y, tr = 0.2, ...) Dqcomhd(x, y, q = c(1:9)/10, nboot = 1000, na.rm = TRUE, ...) dep.effect(x, y, tr = 0.2, nboot = 1000, ...)
x |
an numeric vector of data values (e.g. for time 1). |
y |
an numeric vector of data values (e.g. for time 2). |
tr |
trim level for the means. |
q |
quantiles to be compared. |
nboot |
number of bootstrap samples. |
na.rm |
whether missing values should be removed. |
... |
currently ignored. |
The test statistic is a paired samples generalization of Yuen's independent samples t-test on trimmed means.
dep.effect computes the following effect sizes:
AKP: robust standardized difference similar to Cohen's d
QS: Quantile shift based on the median of the distribution of difference scores,
QStr: Quantile shift based on the trimmed mean of the distribution of X-Y
SIGN: P(X<Y), probability that for a random pair, the first is less than the second.
yuend returns an object of class "yuen" containing:
test |
value of the test statistic (t-statistic) |
p.value |
p-value |
conf.int |
confidence interval |
df |
degress of freedom |
diff |
trimmed mean difference |
call |
function call |
Dqcomhd returns an object of class "robtab" containing:
partable |
parameter table |
dep.effect returns a matrix with the null value of the effect size, the estimated effect size, small/medium/large conventions, and lower/upper CI bounds.
Wilcox, R. (2012). Introduction to Robust Estimation and Hypothesis Testing (3rd ed.). Elsevier.
## Cholesterol data from Wilcox (2012, p. 197) before <- c(190, 210, 300,240, 280, 170, 280, 250, 240, 220) after <- c(210, 210, 340, 190, 260, 180, 200, 220, 230, 200) yuend(before, after) set.seed(123) Dqcomhd(before, after, nboot = 200, q = c(0.25, 0.5, 0.75)) set.seed(123) dep.effect(before, after)## Cholesterol data from Wilcox (2012, p. 197) before <- c(190, 210, 300,240, 280, 170, 280, 250, 240, 220) after <- c(210, 210, 340, 190, 260, 180, 200, 220, 230, 200) yuend(before, after) set.seed(123) Dqcomhd(before, after, nboot = 200, q = c(0.25, 0.5, 0.75)) set.seed(123) dep.effect(before, after)
Performs a robust mediation test according to the method proposed by Zu & Yuan (2010).
ZYmediate(x, y, med, nboot = 2000, alpha = 0.05, kappa = 0.05, ...)ZYmediate(x, y, med, nboot = 2000, alpha = 0.05, kappa = 0.05, ...)
x |
vector with predictor values. |
y |
vector with response values. |
med |
vector with mediator values. |
nboot |
number of bootstrap samples. |
alpha |
alpha level. |
kappa |
the percent of cases to be controlled when robust method is used. |
... |
currently ignored. |
ZYmediate returns an object of class "robmed" containing:
a.est |
effect of predictor on mediator) |
b.est |
effect of mediator on response (in multiple regression model which includes the predictor as well) |
ab.est |
indirect effect (mediation effect) |
CI.ab |
confidence interval mediation effect |
p.value |
p-value mediation test |
call |
function call |
Zu J., Yuan, K.-H. (2010). Local influence and robust procedures for mediation analysis. Multivariate Behavioral Research, 45, 1-44.
## Leerkes data: ## Y: Efficacy ## X: MatCare ## M: Esteem ## fitting robust mediator regressions require(MASS) summary(rlm(Efficacy ~ MatCare, data = Leerkes)) summary(rlm(Esteem ~ MatCare, data = Leerkes)) summary(rlm(Efficacy ~ MatCare + Esteem, data = Leerkes)) ## robust testing of mediating effect (indirect effect) with(Leerkes, ZYmediate(MatCare, Efficacy, Esteem))## Leerkes data: ## Y: Efficacy ## X: MatCare ## M: Esteem ## fitting robust mediator regressions require(MASS) summary(rlm(Efficacy ~ MatCare, data = Leerkes)) summary(rlm(Esteem ~ MatCare, data = Leerkes)) summary(rlm(Efficacy ~ MatCare + Esteem, data = Leerkes)) ## robust testing of mediating effect (indirect effect) with(Leerkes, ZYmediate(MatCare, Efficacy, Esteem))