Package 'heplots'

Description

The heplots package provides functions for visualizing hypothesis tests in multivariate linear models (MANOVA, multivariate multiple regression, MANCOVA, and repeated measures designs). HE plots represent sums-of-squares-and-products matrices for linear hypotheses and for error using ellipses (in two dimensions), ellipsoids (in three dimensions), or by line segments in one dimension.

Details

The basic theory behind HE plots is described by Friendly (2007). See Fox, Friendly and Monette (2007) for a brief introduction; Friendly & Sigal (2016) for a tutorial on these methods; and Friendly, Monette and Fox (2013) for a general discussion of the role of elliptical geometry in statistical understanding.

Other topics now addressed here include robust MLMs, tests for equality of covariance matrices in MLMs, and chi square Q-Q plots for MLMs.

The package also provides a collection of data sets illustrating a variety of multivariate linear models of the types listed above, together with graphical displays.

Several tutorial vignettes are also included. See vignette(package="heplots").

The graphical functions contained here all display multivariate model effects in variable (data) space, for one or more response variables (or contrasts among response variables in repeated measures designs).

list(list("heplot"))

constructs two-dimensional HE plots for model terms and linear hypotheses for pairs of response variables in multivariate linear models.

list(list("heplot3d"))

constructs analogous 3D plots for triples of response variables.

list(list("pairs.mlm"))

constructs a “matrix” of pairwise HE plots.

list(list("heplot1d"))

constructs 1-dimensional analogs of HE plots for model terms and linear hypotheses for single response variables.

For repeated measure designs, between-subject effects and within-subject effects must be plotted separately, because the error terms (E matrices) differ. For terms involving within-subject effects, these functions carry out a linear transformation of the matrix Y of responses to a matrix Y M, where M is the model matrix for a term in the intra-subject design and produce plots of the H and E matrices in this transformed space. The vignette repeated describes these graphical methods for repeated measures designs.

The related car package calculates Type II and Type III tests of multivariate linear hypotheses using the Anova and linearHypothesis functions.

The candisc-package package provides functions for visualizing effects for MLM model terms in a low-dimensional canonical space that shows the largest hypothesis relative to error variation. The candisc package now also includes related methods for canonical correlation analysis.

The heplots package also contains a large number of multivariate data sets with examples of analyses and graphic displays. Use data(package="heplots") to see the current list.

Author(s)

Michael Friendly, John Fox, and Georges Monette

Maintainer: Michael Friendly, [email protected], http://datavis.ca

References

Friendly, M. (2006). Data Ellipses, HE Plots and Reduced-Rank Displays for Multivariate Linear Models: SAS Software and Examples. Journal of Statistical Software, 17(6), 1-42. % https://www.jstatsoft.org/v17/i06/, doi:10.18637/jss.v017.i06

Friendly, M. (2007). HE plots for Multivariate General Linear Models. Journal of Computational and Graphical Statistics, 16(2) 421-444. http://datavis.ca/papers/jcgs-heplots.pdf, doi:10.1198/106186007X208407

Fox, J., Friendly, M. & Monette, G. (2007). Visual hypothesis tests in multivariate linear models: The heplots package for R. DSC 2007: Directions in Statistical Computing. https://socialsciences.mcmaster.ca/jfox/heplots-dsc-paper.pdf

Friendly, M. (2010). HE Plots for Repeated Measures Designs. Journal of Statistical Software, 37(4), 1-40. doi:10.18637/jss.v037.i04.

Fox, J., Friendly, M. & Weisberg, S. (2013). Hypothesis Tests for Multivariate Linear Models Using the car Package. The R Journal, 5(1), https://journal.r-project.org/articles/RJ-2013-004/RJ-2013-004.pdf.

Friendly, M., Monette, G. & Fox, J. (2013). Elliptical Insights: Understanding Statistical Methods Through Elliptical Geometry. Statistical Science, 2013, 28 (1), 1-39, http://datavis.ca/papers/ellipses.pdf.

Friendly, M. & Sigal, M. (2014). Recent Advances in Visualizing Multivariate Linear Models. Revista Colombiana de Estadistica, 37, 261-283 doi:10.15446/rce.v37n2spe.47934.

Friendly, M. & Sigal, M. (2016). Graphical Methods for Multivariate Linear Models in Psychological Research: An R Tutorial. Submitted for publication.

See Also

\code{\link[car]{Anova}}, \code{\link[car]{linearHypothesis}} for Anova.mlm computations and tests

\code{\link[candisc]{candisc-package}} for reduced-rank views in canonical space

\code{\link[stats]{manova}} for a different approach to testing effects in MANOVA designs

Description

This data was taken from the National Longitudinal Study of Adolescent Health. It is a cross-sectional sample of participants from grades 7–12, described and analyzed by Warne (2014).

Format

A data frame with 4344 observations on the following 3 variables.

grade

an ordered factor with levels 7 < 8 < 9 < 10 < 11 < 12

depression

a numeric vector

anxiety

a numeric vector

Details

depression is the response to the question "In the last month, how often did you feel depressed or blue?"

anxiety is the response to the question "In the last month, how often did you have trouble relaxing?"

The responses for depression and anxiety were recorded on a 5-point Likert scale, with categories ⁠0="Never", 1="Rarely", 2="Occasionally", 3="Often", 4="Every day"⁠

Source

Warne, R. T. (2014). A primer on Multivariate Analysis of Variance (MANOVA) for Behavioral Scientists. Practical Assessment, Research & Evaluation, 19 (1).

Examples

data(AddHealth)

if(require(dplyr) & require(ggplot2)) {
# find means & std.errors by grade
means <- AddHealth |>
group_by(grade) |>
  summarise(
    n = n(),
    dep_se = sd(depression, na.rm = TRUE) / sqrt(n),
    anx_se = sd(anxiety, na.rm = TRUE) / sqrt(n),
    depression = mean(depression),
    anxiety = mean(anxiety) ) |> 
  relocate(depression, anxiety, .after = grade) |>
  print()
  
# plot means with std.error bars
ggplot(data = means, aes(x = anxiety, y = depression, 
color = grade)) +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = anxiety - anx_se, 
                     xmax = anxiety + anx_se)) +
  geom_errorbar(aes(ymin = depression - dep_se, 
                    ymax = depression + dep_se)) +
  geom_line(aes(group = 1), linewidth = 1.5) +
  geom_label(aes(label = grade), 
             nudge_x = -0.015, nudge_y = 0.02) +
  scale_color_discrete(guide = "none") +
  theme_bw(base_size = 15)
}

# fit mlm
AH.mod <- lm(cbind(anxiety, depression) ~ grade, data=AddHealth)

car::Anova(AH.mod)
summary(car::Anova(AH.mod))

heplot(AH.mod, hypotheses="grade.L", 
       fill=c(TRUE, FALSE),
       level = 0.4)

Description

Data are a subset from an observational, longitudinal, study on adopted children. Is child's intelligence related to intelligence of the biological mother and the intelligence of the adoptive mother?

Format

A data frame with 62 observations on the following 6 variables.

AMED

adoptive mother's years of education (proxy for her IQ)

BMIQ

biological mother's score on IQ test

Age2IQ

IQ of child at age 2

Age4IQ

IQ of child at age 4

Age8IQ

IQ of child at age 8

Age13IQ

IQ of child at age 13

Details

The child's intelligence was measured at age 2, 4, 8, and 13 for this sample. How does intelligence change over time, and how are these changes related to intelligence of the birth and adoptive mother?

Source

Ramsey, F.L. and Schafer, D.W. (2002). The Statistical Sleuth: A Course in Methods of Data Analysis (2nd ed), Duxbury.

This data set is identical to ex1605 in the Sleuth2 package.

References

Friendly, M. (2010). HE Plots for Repeated Measures Designs. Journal of Statistical Software, 37(4), 1-40. doi:10.18637/jss.v037.i04.

Skodak, M. and Skeels, H.M. (1949). A Final Follow-up Study of One Hundred Adopted Children, Journal of Genetic Psychology 75: 85–125.

See Also

ex1605

Examples

# Treat as multivariate regression problem
Adopted.mod <- lm(cbind(Age2IQ, Age4IQ, Age8IQ, Age13IQ) ~ AMED + BMIQ, 
                  data=Adopted)
Adopted.mod


require(car)
# test overall multivariate regression
print(linearHypothesis(Adopted.mod, c("AMED","BMIQ")), SSP=FALSE)

# show separate linear regressions
op <- par(mfcol=c(2,4), mar=c(4,4,1,1)+.1)
for (i in 3:6) {
	dataEllipse(as.matrix(Adopted[,c(1,i)]),
	            col="black", levels=0.68, ylim=c(70,140))
	abline(lm(Adopted[,i] ~ Adopted[,1]), col="red", lwd=2)

	dataEllipse(as.matrix(Adopted[,c(2,i)]),
	            col="black", levels=0.68, ylim=c(70,140))
	abline(lm(Adopted[,i] ~ Adopted[,2]), col="red", lwd=2)
	abline(a=0,b=1, lty=1, col="blue")
}
par(op)

# between-S (MMReg) plots
heplot(Adopted.mod, hypotheses=list("Reg"=c("AMED", "BMIQ")),
	main="IQ scores of adopted children: MMReg")

pairs(Adopted.mod, hypotheses=list("Reg"=c("AMED", "BMIQ")))

if(requireNamespace("rgl")){
heplot3d(Adopted.mod, hypotheses=list("Reg"=c("AMED", "BMIQ")),
	col = c("red", "blue", "black", "gray"), wire=FALSE)
}

# Treat IQ at different ages as a repeated measure factor
# within-S models & plots
Age <- data.frame(Age=ordered(c(2,4,8,13)))
car::Anova(Adopted.mod, idata=Age, idesign=~Age, test="Roy")

# within-S plots
heplot(Adopted.mod, idata=Age, idesign=~Age, iterm="Age",
	cex=1.25, cex.lab=1.4, fill=c(FALSE, TRUE),
	hypotheses=list("Reg"=c("AMED", "BMIQ"))
	)

Description

Draws a 3D arrow in an rgl scene with barbs at the arrow head

Usage

arrow3d(
  p0 = c(0, 0, 0),
  p1 = c(1, 1, 1),
  barblen,
  s = 0.05,
  theta = pi/6,
  n = 3,
  ...
)

Arguments

Value

Returns (invisibly): integer ID of the line added to the scene %%

Author(s)

Barry Rowlingson, posted to R-help, 1/10/2010

See Also

lines3d, segments3d,

Other 3D plotting: bbox3d(), cross3d(), ellipse3d.axes(), heplot3d()

Examples

arrow3d(c(0,0,0), c(2,2,2), barblen=.2, lwd=3, col="black")
arrow3d(c(0,0,0), c(-2,2,2), barblen=.2, lwd=3, col="red")

Description

This function extends bartlett.test to a multivariate response setting. It performs the Bartlett test of homogeneity of variances for each of a set of response variables, and prints a compact summary.

Bartlett's test is the univariate version of Box's M test for equality of covariance matrices. This function provides a univariate follow-up test to Box's M test to give one simple assessment of which response variables contribute to significant differences in variances among groups.

Usage

bartlettTests(y, ...)

## Default S3 method:
bartlettTests(y, group, ...)

## S3 method for class 'formula'
bartlettTests(y, data, ...)

## S3 method for class 'lm'
bartlettTests(y, ...)

Arguments

Value

An object of classes "anova" and "data.frame", with one observation for each response variable in y.

Author(s)

Michael Friendly

References

Bartlett, M. S. (1937). Properties of sufficiency and statistical tests. Proceedings of the Royal Society of London Series A, 160, 268-282.

See Also

boxM for Box's M test for all responses together.

Other homogeneity tests: leveneTests()

Examples

bartlettTests(iris[,1:4], iris$Species)

data(Skulls, package="heplots")
bartlettTests(Skulls[,-1], Skulls$epoch)

# formula method
bartlettTests(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls)

Description

Ellipsoids are created by rgl functions as meshes of points, segments, ... from coordinates in various forms. This function calculates the bounding box, defined as the range of the x, y, and z coordinates.

Usage

bbox3d(x, ...)

Arguments

Value

A 2 x 3 matrix, giving the minimum and maximum values in the rows and x, y, z coordinates
        in the columns.

See Also

Other 3D plotting: arrow3d(), cross3d(), ellipse3d.axes(), heplot3d()

Description

Pabalan, Davey and Packe (2000) studied the effects of captivity and maltreatment on reproductive capabilities of queen and worker bees in a complex factorial design.

Format

A data frame with 246 observations on the following 6 variables.

caste

a factor with levels Queen Worker

treat

a factor with levels "" CAP MAL

time

an ordered factor: time of treatment

Iz

an index of ovarian development

Iy

an index of ovarian reabsorption

trtime

a factor with levels 0 CAP05 CAP07 CAP10 CAP12 CAP15 MAL05 MAL07 MAL10 MAL12 MAL15

Details

Bees were placed in a small tube and either held captive (CAP) or shaken periodically (MAL) for one of 5, 7.5, 10, 12.5 or 15 minutes, after which they were sacrificed and two measures: ovarian development (Iz) and ovarian reabsorption (Iy), were taken. A single control group was measured with no such treatment, i.e., at time 0; there are n=10 per group.

The design is thus nearly a three-way factorial, with factors caste (Queen, Worker), treat (CAP, MAL) and time, except that there are only 11 combinations of Treatment and Time; we call these trtime below.

Models for the three-way factorial design, using the formula cbind(Iz,Iy) ~ caste*treat*time ignore the control condition at time==0, where treat==NA.

To handle the additional control group at time==0, while separating the effects of Treatment and Time, 10 contrasts can be defined for the trtime factor in the model cbind(Iz,Iy) ~ caste*trtime See demo(bees.contrasts) for details.

In the heplot examples below, the default size="evidence" displays are too crowded to interpret, because some effects are so highly significant. The alternative effect-size scaling, size="effect", makes the relations clearer.

Source

Pabalan, N., Davey, K. G. & Packe, L. (2000). Escalation of Aggressive Interactions During Staged Encounters in Halictus ligatus Say (Hymenoptera: Halictidae), with a Comparison of Circle Tube Behaviors with Other Halictine Species Journal of Insect Behavior, 13, 627-650.

References

Friendly, M. (2006). Data Ellipses, HE Plots and Reduced-Rank Displays for Multivariate Linear Models: SAS Software and Examples Journal of Statistical Software, 17, 1-42.

Examples

data(Bees)
require(car)

# 3-way factorial, ignoring 0 group
bees.mod <- lm(cbind(Iz,Iy) ~ caste*treat*time, data=Bees)
car::Anova(bees.mod)

op<-palette(c(palette()[1:4],"brown","magenta", "olivedrab","darkgray"))
heplot(bees.mod, 
    xlab="Iz: Ovarian development", 
    ylab="Iz: Ovarian reabsorption",
		main="Bees: ~caste*treat*time")

heplot(bees.mod, size="effect",
    xlab="Iz: Ovarian development", 
    ylab="Iz: Ovarian reabsorption",
    main="Bees: ~caste*treat*time", 
    )

# two-way design, using trtime
bees.mod1 <- lm(cbind(Iz,Iy) ~ caste*trtime, data=Bees)
Anova(bees.mod1)

# HE plots for this model, with both significance and effect size scaling

heplot(bees.mod1, 
    xlab="Iz: Ovarian development", 
    ylab="Iz: Ovarian reabsorption",
		main="Bees: ~caste*trtime")
heplot(bees.mod1, 
    xlab="Iz: Ovarian development", 
    ylab="Iz: Ovarian reabsorption",
    main="Bees: ~caste*trtime",
    size="effect")
palette(op)

# effect plots for separate responses
if(require(effects)) {
	bees.lm1 <-lm(Iy ~ treat*caste*time, data=Bees)
	bees.lm2 <-lm(Iz ~ treat*caste*time, data=Bees)
	
	bees.eff1 <- allEffects(bees.lm1)
	plot(bees.eff1,multiline=TRUE,ask=FALSE)
	
	bees.eff2 <- allEffects(bees.lm2)
	plot(bees.eff2,multiline=TRUE,ask=FALSE)
}

Description

boxM() performs the Box's (1949) M-test for homogeneity of covariance matrices obtained from multivariate normal data according to one or more classification factors. The test compares the product of the log determinants of the separate covariance matrices to the log determinant of the pooled covariance matrix, analogous to a likelihood ratio test. The test statistic uses a chi-square approximation.

Usage

boxM(Y, ...)

## S3 method for class 'formula'
boxM(Y, data, ...)

## S3 method for class 'lm'
boxM(Y, ...)

## Default S3 method:
boxM(Y, group, ...)

## S3 method for class 'boxM'
print(x, ...)

## S3 method for class 'boxM'
summary(
  object,
  digits = getOption("digits") - 2,
  cov = FALSE,
  quiet = FALSE,
  ...
)

Arguments

Details

As an object of class "boxM", a few methods are available: print.boxM(), summary.boxM() and plot.boxM().

There is no general provision as yet for handling missing data. Missing data are simply removed, with a warning.

As well, the computation assumes that the covariance matrix for each group is non-singular, so that $\log det(S_i)$ can be calculated for each group. At the minimum, this requires that $n > p$ for each group.

Box's M test for a multivariate linear model highly sensitive to departures from multivariate normality, just as the analogous univariate test. It is also affected adversely by unbalanced designs. Some people recommend to ignore the result unless it is very highly significant, e.g., p < .0001 or worse.

In general, heterogeneity of covariance matrices can be more easily seen and understood by plotting the covariance ellipses using covEllipses.

The summary method prints a variety of additional statistics based on the eigenvalues of the covariance matrices. These are returned invisibly, as a list containing the following components:

logDet

the vector of log determinants

eigs

eigenvalues of the covariance matrices

eigstats

statistics computed on the eigenvalues for each covariance matrix:

product

the product of eigenvalues, $\prod{\lambda_i}$

sum

the sum of eigenvalues, $\sum{\lambda_i}$

precision

the average precision of eigenvalues, $1/\sum(1/\lambda_i)$

max

the maximum eigenvalue, $\lambda_1$

Value

A list with class c("boxM", "htest") containing the following components:

Author(s)

The default method was taken from the biotools package, Anderson Rodrigo da Silva [email protected]

Generalized by Michael Friendly and John Fox

References

Box, G. E. P. (1949). A general distribution theory for a class of likelihood criteria. Biometrika, 36, 317-346.

Morrison, D.F. (1976) Multivariate Statistical Methods.

See Also

leveneTest carries out homogeneity of variance tests for univariate models with better statistical properties.

plot.boxM, a simple dot plot of the log determinants compared with that of the pooled covariance matrix, and also of other quantities computed from their eigenvalues

covEllipses plots covariance ellipses in variable space for several groups.

Examples

data(iris)

# default method, using `Y`, `group` 
res <- boxM(iris[, 1:4], iris[, "Species"])
res

# summary method gives details
summary(res)

# visualize (this is what is done in the plot method)
dets <- res$logDet
ng <- length(res$logDet)-1
dotchart(dets, xlab = "log determinant")
points(dets , 1:4, cex=c(rep(1.5, ng), 2.5), pch=c(rep(16, ng), 15),
       col= c(rep("blue", ng), "red"))

# plot method gives confidence intervals for logDet
plot(res, gplabel="Species")

# formula method
boxM( cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species,
      data=iris)

### Skulls data
data(Skulls)

# lm method
skulls.mod <- lm(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls)
skulls.boxM <- boxM(skulls.mod) |>
  print()
summary(skulls.boxM)

Description

Displays confidence ellipses for all parameters in an multivariate linear model, for a given pair of variables. As such, it is a generalization of confidenceEllipse.

Usage

coefplot(object, ...)

## S3 method for class 'mlm'
coefplot(
  object,
  variables = 1:2,
  parm = NULL,
  df = NULL,
  level = 0.95,
  intercept = FALSE,
  Scheffe = FALSE,
  bars = TRUE,
  fill = FALSE,
  fill.alpha = 0.2,
  labels = !add,
  label.pos = NULL,
  xlab,
  ylab,
  xlim = NULL,
  ylim = NULL,
  axes = TRUE,
  main = "",
  add = FALSE,
  lwd = 1,
  lty = 1,
  pch = 19,
  col = palette(),
  cex = 2,
  cex.label = 1.5,
  cex.lab = par("cex.lab"),
  lty.zero = 3,
  col.zero = 1,
  pch.zero = "+",
  verbose = FALSE,
  ...
)

Arguments

Value

Returns invisibly a list of the coordinates of the ellipses drawn

Author(s)

Michael Friendly

See Also

confidenceEllipse

Other multivariate linear models: glance.mlm()

Examples

rohwer.mlm <- lm(cbind(SAT,PPVT,Raven)~n+s+ns, data=Rohwer)

coefplot(rohwer.mlm, lwd=2, 
         main="Bivariate coefficient plot for SAT and PPVT", fill=TRUE)
coefplot(rohwer.mlm, add=TRUE, Scheffe=TRUE, fill=TRUE)

coefplot(rohwer.mlm, var=c(1,3))

mod1 <- lm(cbind(SAT,PPVT,Raven)~n+s+ns+na+ss, data=Rohwer)
coefplot(mod1, lwd=2, fill=TRUE, parm=(1:5),
	main="Bivariate 68% coefficient plot for SAT and PPVT", level=0.68)

Description

colDevs calculates the column deviations of data values from a central value (mean, median, etc.), possibly stratified by a grouping variable.

Usage

colDevs(x, group, center = mean, group.var = FALSE, ...)

Arguments

Details

Conceptually, the function is similar to a column-wise sweep, by group, allowing an arbitrary center function.

Non-numeric columns of x are removed, with a warning.

Value

By default, it returns a numeric matrix containing the deviations from the centering function. If levels==TRUE, it returns a data.frame containing the group factor prepended to the matrix of deviations.

Author(s)

Michael Friendly

See Also

colMeans for column means,

sweep

Examples

data(iris)

Species <- iris$Species
irisdev <- colDevs(iris[,1:4], Species, mean)

irisdev <- colDevs(iris[,1:4], Species, median)
# trimmed mean, using an anonymous function
irisdev <- colDevs(iris[,1:4], Species, function(x) mean(x, trim=0.25))

# include the group factor in output
irisdev <- colDevs(iris[,1:4], Species, group.var = "Species")
head(irisdev)

# no grouping variable: deviations from column grand means
# include all variables (but suppress warning for this doc)
irisdev <- suppressWarnings( colDevs(iris) )

# two-way design
colDevs(Plastic[,1:3], Plastic[,"rate"])
colDevs(Plastic[,1:3], Plastic[,"additive"])
# cell deviations
#' colDevs(Plastic[,1:3], interaction(Plastic[,c("rate", "additive")]))

Description

The function draws covariance ellipses for one or more groups and optionally for the pooled total sample. It uses either the classical product-moment covariance estimate, or a robust alternative, as provided by cov.rob. Provisions are provided to do this for more than two variables, in a scatterplot matrix format.

These plot methods provide one way to visualize possible heterogeneity of within-group covariance matrices in a one-way MANOVA design. When covariance matrices are nearly equal, their covariance ellipses should all have the same shape. When centered at a common mean, they should also all overlap.

They can also be used to visualize the difference between classical and robust covariance matrices by overlaying the two in a single plot (via add=TRUE).

Usage

covEllipses(x, ...)

## S3 method for class 'data.frame'
covEllipses(
  x,
  group,
  pooled = TRUE,
  method = c("classical", "mve", "mcd"),
  ...
)

## S3 method for class 'matrix'
covEllipses(
  x,
  group,
  pooled = TRUE,
  method = c("classical", "mve", "mcd"),
  ...
)

## S3 method for class 'formula'
covEllipses(x, data, ...)

## S3 method for class 'boxM'
covEllipses(x, ...)

## Default S3 method:
covEllipses(
  x,
  means,
  df,
  labels = NULL,
  variables = 1:2,
  level = 0.68,
  segments = 60,
  center = FALSE,
  center.pch = "+",
  center.cex = 2,
  col = getOption("heplot.colors", c("red", "blue", "black", "darkgreen", "darkcyan",
    "brown", "magenta", "darkgray")),
  lty = 1,
  lwd = 2,
  fill = FALSE,
  fill.alpha = 0.3,
  label.pos = 0,
  xlab,
  ylab,
  vlabels,
  var.cex = 2,
  main = "",
  xlim,
  ylim,
  axes = TRUE,
  offset.axes,
  add = FALSE,
  ...
)

Arguments

Details

The arguments labels, col, lty, lwd, fill, fill.alpha and label.pos are used to draw the ellipses for the groups and also for the pooled, within-group covariance, which is the last in a list when these are computed by the functions. These arguments are each taken in the order specified, and recycled as necessary.

Value

Nothing is returned. The function is used for its side-effect of producing a plot. %Returns invisibly an object of class "covEllipse", %which is a list of the coordinates for the ellipses drawn.

Author(s)

Michael Friendly

See Also

heplot, boxM,

cov.rob

Other covariance ellipses: ellipse.axes(), ellipse.box(), ellipse3d.axes()

Examples

data(iris)

# compare classical and robust covariance estimates
covEllipses(iris[,1:4], iris$Species)
# use `method="mve"`
covEllipses(iris[,1:4], iris$Species, 
    fill=TRUE, method="mve", add=TRUE, labels="")

# method for a boxM object	
iris.boxM <- boxM(iris[, 1:4], iris[, "Species"])
iris.boxM
# show the associated covariance ellipses
covEllipses(iris.boxM, fill=c(rep(FALSE,3), TRUE) )
# use centering
covEllipses(iris.boxM, fill=c(rep(FALSE,3), TRUE), center=TRUE, label.pos=1:4 )

# method for a list of covariance matrices
cov <- c(iris.boxM$cov, pooled=list(iris.boxM$pooled))
df <- c(table(iris$Species)-1, nrow(iris)-3)
covEllipses(cov, iris.boxM$means, df, label.pos=3, fill=c(rep(FALSE,3), TRUE))
 
covEllipses(cov, iris.boxM$means, df, label.pos=3, fill=c(rep(FALSE,3), TRUE), center=TRUE)

# scatterplot matrix version, specifying `variables`
covEllipses(iris[,1:4], iris$Species, 
	fill=c(rep(FALSE,3), TRUE), variables=1:4, 
	fill.alpha=.1)

Description

A chi square quantile-quantile plots show the relationship between data-based values which should be distributed as $\chi^2$ and corresponding quantiles from the $\chi^2$ distribution. In multivariate analyses, this is often used both to assess multivariate normality and check for or identify outliers.

For a data frame of numeric variables or a matrix supplied as the argument x, it uses the Mahalanobis squared distances ($D^2$) of observations $\mathbf{x}$ from the centroid $\bar{\mathbf{x}}$ taking the sample covariance matrix $\mathbf{S}$ into account,

$D^2 = (\mathbf{x} - \bar{\mathbf{x}})^\prime \; \mathbf{S}^{-1} \; (\mathbf{x} - \bar{\mathbf{x}}) \; .$

The method for "mlm" objects fit using lm for a multivariate response applies this to the residuals from the model.

Usage

cqplot(x, ...)

## S3 method for class 'mlm'
cqplot(x, ...)

## Default S3 method:
cqplot(
  x,
  method = c("classical", "mcd", "mve"),
  detrend = FALSE,
  pch = 19,
  col = palette()[1],
  cex = par("cex"),
  ref.col = "red",
  ref.lwd = 2,
  conf = 0.95,
  env.col = "gray",
  env.lwd = 2,
  env.lty = 1,
  env.fill = TRUE,
  fill.alpha = 0.2,
  fill.color = trans.colors(ref.col, fill.alpha),
  labels = if (!is.null(rownames(x))) rownames(x) else 1:nrow(x),
  id.n,
  id.method = "r",
  id.cex = 1,
  id.col = palette()[1],
  xlab,
  ylab,
  main,
  what = deparse(substitute(x)),
  ylim,
  ...
)

Arguments

Details

cqplot is a more general version of similar functions in other packages that produce chi square QQ plots. It allows for classical Mahalanobis squared distances as well as robust estimates based on the MVE and MCD; it provides an approximate confidence (concentration) envelope around the line of unit slope, a detrended version, where the reference line is horizontal, the ability to identify or label unusual points, and other graphical features.

Cases with any missing values are excluded from the calculation and graph with a warning.

Confidence envelope

In the typical use of QQ plots, it essential to have something in the nature of a confidence band around the points to be able to appreciate whether, and to what degree the observed data points differ from the reference distribution. For cqplot, this helps to assess whether the data are reasonably distributed as multivariate normal and also to flag potential outliers.

The calculation of the confidence envelope here follows that used in the SAS program, http://www.datavis.ca/sasmac/cqplot.html which comes from Chambers et al. (1983), Section 6.8.

The essential formula computes the standard errors as:

$\text{se} ( D^2_{(i)} ) = \frac{\hat{b}} {d ( q_i)} \times \sqrt{ p_i (1-p_i) / n }$

where $D^2_{(i)}$ is the i-th ordered value of $D^2$, $\hat{b}$ is an estimate of the slope of the reference line obtained from the ratio of the interquartile range of the $D^2$ values to that of the $\chi^2_p$ distribution and $d(q_i)$ is the density of the chi square distribution at the quantile $q_i$.

The pointwise confidence envelope of coverage conf = $1-\alpha$ is then calculated as $D^2_{(i)} \pm z_{1-\alpha/2} \text{se} ( D^2_{(i)} )$

Note that this confidence envelope applies only to the $D^2$ computed using the classical estimates of location ($\bar{\mathbf{x}}$) and scatter ($\mathbf{S}$). The qqPlot function provides for simulated envelopes, but only for a univariate measure. Oldford (2016) provides a general theory and methods for QQ plots.

Value

Returns invisibly a data.frame containing squared Mahalanobis distances (DSQ), their quantiles and p-values corresponding to the rows of x or the residuals of the model for the identified points, else NULL if no points are identified.

Author(s)

Michael Friendly

References

J. Chambers, W. S. Cleveland, B. Kleiner, P. A. Tukey (1983). Graphical methods for data analysis, Wadsworth.

R. W. Oldford (2016), "Self calibrating quantile-quantile plots", The American Statistician, 70, 74-90.

See Also

Mahalanobis for calculation of Mahalanobis squared distance;

qqplot; qqPlot can give a similar result for Mahalanobis squared distances of data or residuals; qqtest has many features for all types of QQ plots.

Other diagnostic plots: distancePlot(), eigstatCI(), plot.boxM()

Examples

cqplot(iris[, 1:4])

iris.mod <- lm(as.matrix(iris[,1:4]) ~ Species, data=iris)
out <- cqplot(iris.mod, id.n=3)
# show return value
out

# compare with car::qqPlot
car::qqPlot(Mahalanobis(iris[, 1:4]), dist="chisq", df=4)


# Adopted data
Adopted.mod <- lm(cbind(Age2IQ, Age4IQ, Age8IQ, Age13IQ) ~ AMED + BMIQ, 
                  data=Adopted)
cqplot(Adopted.mod, id.n=3)
cqplot(Adopted.mod, id.n=3, method="mve")


# Sake data
Sake.mod <- lm(cbind(taste, smell) ~ ., data=Sake)
cqplot(Sake.mod)
cqplot(Sake.mod, method="mve", id.n=2)

# SocialCog data -- one extreme outlier
data(SocialCog)
SC.mlm <-  lm(cbind(MgeEmotions,ToM, ExtBias, PersBias) ~ Dx,
               data=SocialCog)
cqplot(SC.mlm, id.n=1)

# data frame example: stackloss data
data(stackloss)
cqplot(stackloss[, 1:3], id.n=4)                # very strange
cqplot(stackloss[, 1:3], id.n=4, detrend=TRUE)
cqplot(stackloss[, 1:3], id.n=4, method="mve")
cqplot(stackloss[, 1:3], id.n=4, method="mcd")

Description

Draws a 3D cross or axis vectors in an rgl scene.

Usage

cross3d(centre = rep(0, 3), scale = rep(1, 3), ...)

Arguments

Value

Used for its side-effect, but returns (invisibly) a 6 by 3 matrix containing the end-points of three axes, in pairs.

Author(s)

Michael Friendly

See Also

segments3d

Other 3D plotting: arrow3d(), bbox3d(), ellipse3d.axes(), heplot3d()

Description

Find degrees of freedom for model terms

Usage

df.terms(model, term, ...)

## Default S3 method:
df.terms(model, term, ...)

Arguments

Description

Reaven and Miller (1979) examined the relationship among blood chemistry measures of glucose tolerance and insulin in 145 nonobese adults. They used the PRIM9 system at the Stanford Linear Accelerator Center to visualize the data in 3D, and discovered a peculiar pattern that looked like a large blob with two wings in different directions.

Format

A data frame with 145 observations on the following 6 variables.

relwt

relative weight, expressed as the ratio of actual weight to expected weight, given the person's height, a numeric vector

glufast

fasting plasma glucose level, a numeric vector

glutest

test plasma glucose level, a measure of glucose intolerance, a numeric vector

instest

plasma insulin during test, a measure of insulin response to oral glucose, a numeric vector

sspg

steady state plasma glucose, a measure of insulin resistance, a numeric vector

group

diagnostic group, a factor with levels Normal Chemical_Diabetic Overt_Diabetic

Details

After further analysis, the subjects were classified as subclinical (chemical) diabetics, overt diabetics and normals. This study was influential in defining the stages of development of Type 2 diabetes. Overt diabetes is the most advanced stage, characterized by elevated fasting blood glucose concentration and classical symptoms. Preceding overt diabetes is the latent or chemical diabetic stage, with no symptoms of diabetes but demonstrable abnormality of oral or intravenous glucose tolerance.

glutest was defined as the "area under the plasma glucose curve for the three hour oral glucose tolerance test." Reaven & Miller refer to this variable as "Glucose area".

instest was defined as the "area under the plasma insulin curve", and is referred to in the paper as "Insulin area".

This study was influential in defining the stages of development of Type 2 diabetes. Overt diabetes is the most advanced stage, characterized by elevated fasting blood glucose concentration and classical symptoms. Preceding overt diabetes is the latent or chemical diabetic stage, with no symptoms of diabetes but demonstrable abnormality of oral or intravenous glucose tolerance.

Source

Andrews, D. F. & Herzberg, A. M. (1985). Data: A Collection of Problems from Many Fields for the Student and Research Worker, Springer-Verlag, Ch. 36.

Friendly, M. (1991). SAS System for Statistical Graphics, Cary, NC: SAS Institute.

References

Reaven, G. M. and Miller, R. G. (1979). An attempt to define the nature of chemical diabetes using a multidimensional analysis. Diabetologia, 16, 17-24.

Examples

data(Diabetes)
col <- c("blue", "red", "darkgreen")[Diabetes$group]
pch <- c(16,15,17)[Diabetes$group]

# a perplexing plot, similar to Fig 2, but with a loess smooth 
plot(instest ~ glutest, data=Diabetes, pch=16,
	cex.lab=1.25,
	xlab="Glucose area (glutest)",
	ylab="Insulin area (instest)")
lines( loess.smooth(Diabetes$glutest, Diabetes$instest), col="blue", lwd=2) 

# scatterplot matrix, colored by group
plot(Diabetes[,1:5], col=col, pch=pch)

# covariance ellipses
covEllipses(Diabetes[,2:5], Diabetes$group, fill=TRUE, pooled=FALSE, 
	col=col)

covEllipses(Diabetes[,2:5], Diabetes$group, fill=TRUE, pooled=FALSE, 
	col=col, variables=1:4)

# Box's M test
diab.boxm <- boxM(Diabetes[,2:5], Diabetes$group)
diab.boxm
plot(diab.boxm)

# heplots
diab.mlm <- lm(cbind(glufast, glutest, instest, sspg) ~ group, data=Diabetes)

heplot(diab.mlm)
pairs(diab.mlm, fill=TRUE, fill.alpha=0.1)

Description

This plot, suggested by Rousseeuw & van Zomeren (1991), Rousseeu et al. (2004) typically plots Mahalanobis distances ($D$) of the Y response variables against the distances of the X variables in a multivariate linear model (MLM). When applied to a multivariate linear model itself, it plots the distances of the residuals for the Y variables against the predictor terms in the model.matrix X.

This diagnostic plot combines the information on regression outliers and leverage points, and often more useful than either distance separately.

Usage

distancePlot(X, Y, ...)

## Default S3 method:
distancePlot(
  X,
  Y,
  method = c("classical", "mcd", "mve"),
  level = 0.975,
  ids = rownames(X),
  pch = c(1, 16),
  col = c("black", "red"),
  label.pos = 2,
  xlab,
  ylab,
  verbose = FALSE,
  ...
)

## S3 method for class 'formula'
distancePlot(X, Y, data, ...)

## S3 method for class 'mlm'
distancePlot(X, ...)

Arguments

Details

Observations with "large" distances on X or Y are labeled with their ids. The cutoffs are calculated as $\sqrt{\chi^2_{k, \text{level}}}$.

Value

   Returns invisibly a data frame containing the distances, `distX`, `distY`

References

Rousseeuw P. J. & van Zomeren B. C. (1991). “Robust Distances: Simulation and Cutoff Values.” In W Stahel, S Weisberg (eds.), Directions in Robust Statistics and Diagnostics, Part II. Springer-Verlag, New York.

Rousseeuw, P. J., Van Driessen, K., Van Aelst, S., & Agullo, J. (2004). Robust multivariate regression. Technometrics, 46(3), 293–305. doi:10.1198/004017004000000329.

See Also

Mahalanobis

Other diagnostic plots: cqplot(), eigstatCI(), plot.boxM()

Examples

if(require("robustbase")) {
  # Examples from Rousseeuw etal (2004)
  data(pulpfiber, package="robustbase")
  # Figure 1
  distancePlot(pulpfiber[, 1:4], pulpfiber[, 5:8])   
  # Figure 3
  pulp.mod <- lm(cbind(Y1, Y2, Y3, Y4) ~ X1 + X2 + X3 + X4, data = pulpfiber)
  distancePlot(pulp.mod, method = "mcd")
}

# NLSY data
data(NLSY, package = "heplots")
NLSY.mlm <- lm(cbind(math, read) ~ income + educ + antisoc + hyperact,
               data = NLSY)

distancePlot(NLSY.mlm)

# gives the same result
distancePlot(NLSY[, 3:6], residuals(NLSY.mlm), level = 0.975)

distancePlot(NLSY.mlm, method ="mve")

# distancePlot(cbind(math, read) ~ income + educ + antisoc + hyperact,
#                data = NLSY)

# schooldata dataset
data(schooldata)
school.mod <- lm(cbind(reading, mathematics, selfesteem) ~ ., data=schooldata)
distancePlot(school.mod, cex = 1.5, cex.lab = 1.2)

data(Hernior)
Hern.mod <- lm(cbind(leave, nurse, los) ~
               age + sex +  pstat +  build + cardiac + resp, data=Hernior)
distancePlot(Hern.mod)

Description

A tiny hypothetical dataset to illustrate one-way MANOVA.

A dogfood manufacturer wanted to study preference for different dogfood formulas, two of their own ("Old", "New") and two from other manufacturers ("Major", "Alps"). In a between-dog design, 4 dogs were presented with a bowl of one formula and the time to start eating and amount eaten were recorded.

Usage

data("dogfood")

Format

A data frame with 16 observations on the following 3 variables.

formula

factor, a factor with levels Old, New, Major, Alps

start

numeric, time to start eating

amount

numeric, amount eaten

Details

In addition to testing the overall effects of formula, three useful (and orthogonal) contrasts can specified for this 3-df factor:

• Ours vs. Theirs, with weights c(1, 1, -1, -1)

• Major vs. Alps, with weights c(0, 0, 1, -1)

• Old vs. New, with weights c(1, -1, 0, 0)

Because these are orthogonal contrasts, they fully decompose the main effect of formula, in that their sum of squares add to the overall sum of squares.

Source

Used in my Psych 6140 lecture notes, http://friendly.apps01.yorku.ca/psy6140/

Examples

data(dogfood)
library(car)
library(candisc)

# make some boxplots
op <- par(mfrow = c(1,2))
boxplot(start ~ formula, data = dogfood)
points(start ~ formula, data = dogfood, pch=16, cex = 1.2)

boxplot(amount ~ formula, data = dogfood)
points(amount ~ formula, data = dogfood, pch=16, cex = 1.2)
par(op)

# setup contrasts to test interesting comparisons
C <- matrix(
       c( 1,  1, -1, -1,         #Ours vs. Theirs
          0,  0,  1, -1,           #Major vs. Alps
          1, -1,  0,  0),             #New vs. Old
       nrow=4, ncol=3)
# assign these to the formula factor
contrasts(dogfood$formula) <- C
# re-fit the model
dogfood.mod <- lm(cbind(start, amount) ~ formula, data=dogfood)

dogfood.mod <- lm(cbind(start, amount) ~ formula, data=dogfood)
Anova(dogfood.mod)

# data ellipses
covEllipses(cbind(start, amount) ~ formula, data=dogfood,
  fill = TRUE, fill.alpha = 0.1)

# test these contrasts with multivariate tests 
linearHypothesis(dogfood.mod, "formula1", title="Ours vs. Theirs")
linearHypothesis(dogfood.mod, "formula2", title="Old vs. New")
linearHypothesis(dogfood.mod, "formula3", title="Alps vs. Major")

heplot(dogfood.mod, fill = TRUE, fill.alpha = 0.1)

# display contrasts in the heplot 
hyp <- list("Ours/Theirs" = "formula1",
            "Old/New" = "formula2")
heplot(dogfood.mod, hypotheses = hyp,
       fill = TRUE, fill.alpha = 0.1)

dogfood.can <- candisc(dogfood.mod, data=dogfood)
heplot(dogfood.can, 
       fill = TRUE, fill.alpha = 0.1, 
       lwd = 2, var.lwd = 2, var.cex = 2)

Description

Computes bootstrap confidence intervals for statistics based on eigenvalues of grouped covariance matrices. This is intended for use with plot.boxM() to provide confidence intervals for measures beyond the default "logDet".

This is a convenience wrapper that extracts the necessary information from a boxM object, but still requires the original data and grouping variable since these are not stored in the boxM object.

Usage

eigstatCI(
  Y,
  group,
  which = c("product", "sum", "precision", "max"),
  R = 1000,
  conf = 0.95,
  type = c("perc", "bca", "norm", "basic"),
  parallel = FALSE,
  ncpus = 2,
  seed = NULL,
  ...
)

eigstatCI_boxM(boxm, Y, group, ...)

Arguments

Details

For each group (and the pooled data), this function performs nonparametric bootstrap resampling to estimate the sampling distribution of the specified eigenvalue-based statistic. Confidence intervals are computed using the percentile method or bias-corrected and accelerated (BCa) method.

Unlike logdetCI() which uses analytic approximations based on asymptotic theory, this function makes no distributional assumptions and can handle small to moderate sample sizes, though computational cost increases with the number of bootstrap replicates.

Value

A data frame with one row for each group plus the pooled data. Columns include:

group

Group name (factor level)

statistic

Observed value of the statistic

lower

Lower confidence limit

upper

Upper confidence limit

bias

Bootstrap estimate of bias (if available)

se

Bootstrap standard error (if available)

A data frame with bootstrap confidence intervals

Author(s)

Michael Friendly

References

Efron, B., & Tibshirani, R. J. (1994). An Introduction to the Bootstrap. CRC Press.

See Also

boxM, plot.boxM, logdetCI

Other diagnostic plots: cqplot(), distancePlot(), plot.boxM()

Examples

## Not run: 
library(boot)
data(iris)

# Bootstrap CI for product of eigenvalues
CI_prod <- eigstatCI(iris[,1:4], iris$Species, which="product", R=500)
CI_prod

# Bootstrap CI for sum of eigenvalues (= trace)
CI_sum <- eigstatCI(iris[,1:4], iris$Species, which="sum", R=500)
CI_sum

# Use with parallel processing for speed
CI_max <- eigstatCI(iris[,1:4], iris$Species, which="max",
                    R=1000, parallel=TRUE, ncpus=4)

## End(Not run)

## Not run: 
library(boot)
data(iris)

# Fit boxM
res <- boxM(iris[,1:4], iris$Species)

# Get bootstrap CIs (must provide original data again)
CI <- eigstatCI_boxM(res, Y = iris[,1:4], group = iris$Species,
                     which = "sum", R = 500)

## End(Not run)

Description

A function to draw the principal axes of a 2D ellipse from a correlation, covariance or sums of squares and cross products matrix in an existing plot.

Usage

ellipse.axes(
  x,
  centre = c(0, 0),
  center = centre,
  scale,
  which = 1:2,
  level = 0.95,
  radius = sqrt(qchisq(level, 2)),
  extend = 0,
  labels = TRUE,
  label.ends = c(2, 4),
  label.pos = c(2, 4, 1, 3),
  type = c("lines", "arrows"),
  ...
)

Arguments

Value

Invisibly returns a 4 x 2 matrix containing the end points of the axes in pairs (min, max) by rows.

Author(s)

Michael Friendly

See Also

lines, text

Other covariance ellipses: covEllipses(), ellipse.box(), ellipse3d.axes()

Examples

data(iris)
cov <- cov(iris[,1:2])
mu <- colMeans(iris[,1:2])

radius <- sqrt(qchisq(0.68, 2))
plot(iris[,1:2], asp=1)
car::ellipse(mu, cov, radius = radius)
res <- ellipse.axes(cov, center=mu, level = 0.68,
                    labels = TRUE)
res

# try some options 
plot(iris[,1:2], asp=1)
car::ellipse(mu, cov, radius = radius)
abline(h=mu[2], v=mu[1], col = "grey")
ellipse.axes(cov, centre=mu, level = 0.68,
             labels = "Dim", label.ends = 1:4,
             lwd = 2, lty = 2, col = "red",
             cex = 1.5)
             
# draw arrows rather than lines            
plot(iris[,1:2], asp=1)
car::ellipse(mu, cov, radius = radius)
ellipse.axes(cov, center=mu, level = 0.68,
             type = "arrows", extend = 0.2)

Description

Draw Conjugate Axes and Parallelogram Surrounding a Covariance Ellipse

Usage

ellipse.box(
  x,
  center = c(0, 0),
  which = 1:2,
  level = 0.95,
  radius = sqrt(qchisq(level, 2)),
  factor = c("cholesky", "pca"),
  draw = c("box", "diameters", "both"),
  ...
)

Arguments

Value

Invisibly returns a 2 column matrix containing the end points of lines.

See Also

Other covariance ellipses: covEllipses(), ellipse.axes(), ellipse3d.axes()

Examples

data(iris)
cov <- cov(iris[,3:4])
mu <- colMeans(iris[,3:4])

radius <- sqrt(qchisq(0.68, 2))
plot(iris[,3:4], asp=1)
car::ellipse(mu, cov, radius = radius)
ellipse.axes(cov, center=mu, level = 0.68,
            labels = TRUE)
ellipse.box(cov, center=mu, level = 0.68, 
            factor = "pca", 
            col = "red", lwd = 2 )

res <- ellipse.box(cov, center=mu, level = 0.68, factor = "chol", col = "green", lwd = 2 )
res

Description

A function to draw the major axes of a 3D ellipsoid from a correlation, covariance or sums of squares and cross products matrix.

Usage

ellipse3d.axes(
  x,
  centre = c(0, 0, 0),
  center = centre,
  scale = c(1, 1, 1),
  level = 0.95,
  t = sqrt(qchisq(level, 3)),
  which = 1:3,
  labels = TRUE,
  label.ends = c(2, 4, 6),
  ...
)

Arguments

Value

Returns a 6 x 3 matrix containing the end points of the three axis lines in pairs by rows.

Author(s)

Michael Friendly

See Also

segments3d, text3d, ellipse3d

Other covariance ellipses: covEllipses(), ellipse.axes(), ellipse.box()

Other 3D plotting: arrow3d(), bbox3d(), cross3d(), heplot3d()

Examples

data(iris)
iris3 <- iris[,1:3]
cov <- cov(iris3)
mu <- colMeans(iris3)
col <-c("blue", "green", "red")[iris$Species]

library(rgl)
plot3d(iris3, type="s", size=0.4, col=col, cex=2, box=FALSE, aspect="iso")
plot3d( ellipse3d(cov, centre=mu, level=0.68), col="gray", alpha=0.2,  add = TRUE)

axes <- ellipse3d.axes(cov, centre=mu, level=0.68, color="gray", lwd=2)

Description

This is an experimental function designed to separate internal code in link{heplot3d}.

Usage

Ellipsoid(x, ...)

## S3 method for class 'data.frame'
Ellipsoid(x, which = 1:3, method = c("classical", "mve", "mcd"), ...)

## Default S3 method:
Ellipsoid(
  x,
  center = c(0, 0, 0),
  which = 1:3,
  radius = 1,
  df = Inf,
  label = "",
  cex.label = 1.5,
  col = "pink",
  lwd = 1,
  segments = 40,
  shade = TRUE,
  alpha = 0.1,
  wire = TRUE,
  verbose = FALSE,
  warn.rank = FALSE,
  ...
)

Arguments

Value

returns the bounding box of the ellipsoid invisibly; otherwise used for it's side effect of drawing the ellipsoid

Examples

# none yet

Description

Calculates partial eta-squared for linear models or multivariate analogs of eta-squared (or R^2), indicating the partial association for each term in a multivariate linear model. For a multivariate model, this is a summary across all response variables.

Usage

etasq(x, ...)

## S3 method for class 'mlm'
etasq(x, ...)

## S3 method for class 'Anova.mlm'
etasq(x, anova = FALSE, ...)

## S3 method for class 'lm'
etasq(x, anova = FALSE, partial = TRUE, ...)

Arguments

Details

There is a different analog of $\eta^2$ for each of the four standard multivariate test statistics: Pillai's trace, Hotelling-Lawley trace, Wilks' Lambda and Roy's maximum root test.

For univariate linear models, classical $\eta^2$ = SSH / SST and partial $\eta^2$ = SSH / (SSH + SSE). These are identical in one-way designs.

Partial eta-squared describes the proportion of total variation attributable to a given factor, partialling out (excluding) other factors from the total nonerror variation. These are commonly used as measures of effect size or measures of (non-linear) strength of association in ANOVA models.

All multivariate tests are based on the $s=min(p, df_h)$ latent roots of $H E^{-1}$. The analogous multivariate partial $\eta^2$ measures are calculated as:

Pillai's trace (V)

$\eta^2 = V/s$

Hotelling-Lawley trace (T)

$\eta^2 = T/(T+s)$

Wilks' Lambda (L)

$\eta^2 = L^{1/s}$

Roy's maximum root (R)

$\eta^2 = R/(R+1)$

Value

When anova=FALSE, a one-column data frame containing the eta-squared values for each term in the model.

When anova=TRUE, a 5-column (lm) or 7-column (mlm) data frame containing the eta-squared values and the test statistics produced by print.Anova() for each term in the model.

Author(s)

Michael Friendly

References

Muller, K. E. and Peterson, B. L. (1984). Practical methods for computing power in testing the Multivariate General Linear Hypothesis Computational Statistics and Data Analysis, 2, 143-158.

Muller, K. E. and LaVange, L. M. and Ramey, S. L. and Ramey, C. T. (1992). Power Calculations for General Linear Multivariate Models Including Repeated Measures Applications. Journal of the American Statistical Association, 87, 1209-1226.

See Also

Anova

eta_squared for a function that calculates this effect size measure for each response variable separately.

Examples

library(car)
data(Soils, package="carData")
soils.mod <- lm(cbind(pH,N,Dens,P,Ca,Mg,K,Na,Conduc) ~ Block + Contour*Depth, data=Soils)
#Anova(soils.mod)
etasq(Anova(soils.mod))
etasq(soils.mod) # same
etasq(Anova(soils.mod), anova=TRUE)

etasq(soils.mod, test="Wilks")
etasq(soils.mod, test="Hotelling")

Description

Data collected as part of a preliminary study examining the relation between football helmet design and neck injuries. There are 30 subjects in each of three groups: High school football players, college players and non-football players.

Format

A data frame with 90 observations on the following 7 variables.

group

a factor with levels ⁠High school⁠ College Non-football

width

a numeric vector: head width at widest dimension

circum

a numeric vector: head circumference

front.back

a numeric vector: front to back distance at eye level

eye.top

a numeric vector: eye to top of head

ear.top

a numeric vector:ear to top of head

jaw

a numeric vector: jaw width

Source

Rencher, A. C. (1995), Methods of Multivariate Analysis, New York: Wiley, Table 8.3.

Examples

data(FootHead)
str(FootHead)
require(car)

# use Helmert contrasts for group
contrasts(FootHead$group) <- contr.helmert
contrasts(FootHead$group)

foot.mod <- lm(cbind(width, circum,front.back,eye.top,ear.top,jaw) ~ group, 
               data=FootHead)
Manova(foot.mod)

# show the HE plot for the first two variables
heplot(foot.mod, main="HE plot for width and circumference", fill=TRUE,
	col=c("red", "blue"))

# show it with tests of Helmert contrasts
heplot(foot.mod, hypotheses=list("group.1"="group1","group.2"="group2"),
	col=c("red", "blue", "green3", "green3" ),
	main="HE plot with orthogonal Helmert contrasts")

# show all pairwise HE plots
pairs(foot.mod)

# ... with tests of Helmert contrasts
pairs(foot.mod, hypotheses=list("group.1"="group1","group.2"="group2"),
	col=c("red", "blue", "green3", "green3"), hyp.labels=FALSE)

# see that the hypothesis for groups really is 2D
if(requireNamespace("rgl")){
heplot3d(foot.mod, variables=c(1,2,6),
	hypotheses=list("group.1"="group1","group.2"="group2"),
	col=c("red", "blue", "green3", "green3" ), wire=FALSE)
}

Description

This function takes an "mlm" object, fit by lm with a multivariate response. The goal is to return something analogous to glance.lm for a univariate response linear model.

Usage

## S3 method for class 'mlm'
glance(x, ...)

Arguments

Details

In the multivariate case, it returns a tibble with one row for each response variable, containing goodness of fit measures, F-tests and p-values.

Value

A tibble with one row for each response variable and the columns:

r.squared

R squared statistic, or the percent of variation explained by the model.

sigma

Estimated standard error of the residuals

fstatitic

Overall F statistic for the model

numdf

Numerator degrees of freedom for the overall test

dendf

Denominator degrees of freedom for the overall test

p.value

P-value corresponding to the F statistic

nobs

Number of observations used

See Also

Other multivariate linear models: coefplot()

Examples

iris.mod <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data=iris)
glance(iris.mod)

Description

gsorth uses sequential, orthogonal projections, as in the Gram-Schmidt method, to transform a matrix or numeric columns of a data frame into an uncorrelated set, possibly retaining the same column means and standard deviations as the original.

Usage

gsorth(y, order, recenter = TRUE, rescale = TRUE, adjnames = TRUE)

Arguments

Details

In statistical applications, interpretation depends on the order of the variables orthogonalized. In multivariate linear models, orthogonalizing the response, Y variables provides the equivalent of step-down tests, where Y1 is tested alone, and then Y2.1, Y3.12, etc. can be tested to determine their additional contributions over the previous response variables.

Similarly, orthogonalizing the model X variables provides the equivalent of Type I tests, such as provided by anova.

The method is equivalent to setting each of columns 2:p to the residuals from a linear regression of that column on all prior columns, i.e.,

z[,j] <- resid( lm( z[,j] ~ as.matrix(z[,1:(j-1)]), data=z) )

However, for accuracy and speed the transformation is carried out using the QR decomposition.

Value

Returns a matrix or data frame with uncorrelated columns. Row and column names are copied to the result.

Author(s)

Michael Friendly

See Also

qr,

Examples

GSiris <- gsorth(iris[,1:4])
GSiris <- gsorth(iris, order=1:4)   # same, using order
str(GSiris)
zapsmall(cor(GSiris))
colMeans(GSiris)
# sd(GSiris) -- sd(<matrix>) now deprecated
apply(GSiris, 2, sd)

# orthogonalize Y side
GSiris <- data.frame(gsorth(iris[,1:4]), Species=iris$Species)
iris.mod1 <- lm(as.matrix(GSiris[,1:4]) ~ Species, data=GSiris)
car::Anova(iris.mod1)

# orthogonalize X side
rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer)
car::Anova(rohwer.mod)

# type I tests for Rohwer data
Rohwer.orth <- cbind(Rohwer[,1:5], gsorth(Rohwer[, c("n", "s", "ns", "na", "ss")], adjnames=FALSE))

rohwer.mod1 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss, data=Rohwer.orth)
car::Anova(rohwer.mod1)
# compare with anova()
anova(rohwer.mod1)

# compare heplots for original Xs and orthogonalized, Type I
heplot(rohwer.mod)
heplot(rohwer.mod1)

Description

A study was conducted investigating the effectiveness of different kinds of psychological treatment on the sensitivity of headache sufferers to noise, described in Hand and Taylor (1987), Study E.

Format

A data frame with 98 observations on the following 6 variables.

type

Type of headache, a factor with levels Migrane Tension

treatment

Treatment group, a factor with levels T1 T2 T3 Control. See Details

u1

Noise level rated as Uncomfortable, initial measure

du1

Noise level rated as Definitely Uncomfortable, initial measure

u2

Noise level rated as Uncomfortable, final measure

du2

Noise level rated as Definitely Uncomfortable, final measure

Details

In a pre-post design, 98 patients were first assessed for the volume of noise which they found uncomfortable (U) and definitely uncomfortable (DU). They were then given relaxation training, where they listened to the noise at the DU level and given instruction breathing techniques and the use of visual imagery to distract them from discomfort. One of four treatments was then applied, and all patients were reassessed for the noise volume they considered uncomfortable (U) and definitely uncomfortable (DU).

The treatments are described as follows:

T1

Listened again to the tone at their initial DU level, for the same amount of time they were able to tolerate it before.

T2

Same as T1, with one additional minute exposure

T3

Same as T2, but were explicitly instructed to use the relaxation techniques

Control

These subject experienced no further exposure to the noise tone until the final sensitivity measures were taken

Hand and Taylor described several substantive hypotheses related to the differences among treatments. In the Headache data frame, these have been included as contrasts(Headache$treatment)

Source

D. J. Hand and C. C. Taylor (1987). Multivariate analysis of variance and repeated measures: a practical approach for behavioural scientists London: Chapman and Hall. ISBN: 0412258005. Table E.1.

Examples

library(car)
data(Headache)
str(Headache)

# basic MLM, specifying between-S effects
headache.mod <- lm(cbind(u1, du1, u2, du2) ~ type * treatment, data=Headache)

##############################
## between-S tests
##############################
Anova(headache.mod, test="Roy")

# test each contrast separately
print(linearHypothesis(headache.mod, hypothesis="treatment1", test="Roy"), SSP=FALSE)
print(linearHypothesis(headache.mod, hypothesis="treatment2", test="Roy"), SSP=FALSE)
print(linearHypothesis(headache.mod, hypothesis="treatment3", test="Roy"), SSP=FALSE)


heplot(headache.mod, variables=c(1,3),
	hypotheses=paste("treatment", 1:3, sep=""),
	hyp.labels=c("extra.exp", "no.inst", "explicit.inst"),
	xlab="u1: Initial sensitivity", ylab="u2: Final sensitivity",
	main="Headache data: Unpleasant levels")
abline(0, 1, lty=5, col="gray")

heplot(headache.mod, variables=c(2,4),
	hypotheses=paste("treatment", 1:3, sep=""),
	xlab="du1: Initial sensitivity", ylab="du2: Final sensitivity",
	main="Headache data: Definitely Unpleasant levels")
abline(0, 1, lty=5, col="gray")

pairs(headache.mod)

##############################
# between-S and within-S tests
##############################
idata = expand.grid(level=factor(c("U", "DU")), phase=factor(1:2))
Anova(headache.mod, idata=idata, idesign=~level*phase)

Description

This function plots ellipses representing the hypothesis and error sums-of-squares-and-products matrices for terms and linear hypotheses in a multivariate linear model. These include MANOVA models (all explanatory variables are factors), multivariate regression (all quantitative predictors), MANCOVA models, homogeneity of regression, as well as repeated measures designs treated from a multivariate perspective.

Usage

heplot(mod, ...)

## S3 method for class 'mlm'
heplot(
  mod,
  terms,
  hypotheses,
  term.labels = TRUE,
  hyp.labels = TRUE,
  err.label = "Error",
  label.pos = NULL,
  label.cex = par("cex"),
  variables = 1:2,
  error.ellipse = !add,
  factor.means = !add,
  grand.mean = !add,
  remove.intercept = TRUE,
  type = c("II", "III", "2", "3"),
  idata = NULL,
  idesign = NULL,
  icontrasts = c("contr.sum", "contr.poly"),
  imatrix = NULL,
  iterm = NULL,
  markH0 = !is.null(iterm),
  manova,
  size = c("evidence", "effect.size", "significance"),
  level = 0.68,
  alpha = 0.05,
  segments = 60,
  center.pch = "+",
  center.cex = 2,
  col = getOption("heplot.colors", c("red", "blue", "black", "darkgreen", "darkcyan",
    "magenta", "brown", "darkgray")),
  lty = 2:1,
  lwd = 1:2,
  fill = FALSE,
  fill.alpha = 0.3,
  xlab,
  ylab,
  main = "",
  xlim,
  ylim,
  axes = TRUE,
  offset.axes,
  add = FALSE,
  verbose = FALSE,
  warn.rank = FALSE,
  ...
)

Arguments

Details

The heplot function plots a representation of the covariance ellipses for hypothesized model terms and linear hypotheses (H) and the corresponding error (E) matrices for two response variables in a multivariate linear model (mlm).

The plot helps to visualize the nature and dimensionality response variation on the two variables jointly in relation to error variation that is summarized in the various multivariate test statistics (Wilks' Lambda, Pillai trace, Hotelling-Lawley trace, Roy maximum root). Roy's maximum root test has a particularly simple visual interpretation, exploited in the size="evidence" version of the plot. See the description of argument alpha.

For a 1 df hypothesis term (a quantitative regressor, a single contrast or parameter test), the H matrix has rank 1 (one non-zero latent root of $H E^{-1}$) and the H "ellipse" collapses to a degenerate line.

Typically, you fit a mlm with mymlm <- lm(cbind(y1, y2, y3, ...) ~ modelterms), and plot some or all of the modelterms with heplot(mymlm, ...). Arbitrary linear hypotheses related to the terms in the model (e.g., contrasts of an effect) can be included in the plot using the hypotheses argument. See linearHypothesis for details.

For repeated measure designs, where the response variables correspond to one or more variates observed under a within-subject design, between-subject effects and within-subject effects must be plotted separately, because the error terms (E matrices) differ. When you specify an intra-subject term (iterm), the analysis and HE plots amount to analysis of the matrix Y of responses post-multiplied by a matrix M determined by the intra-subject design for that term. See Friendly (2010) or the vignette("repeated") in this package for an extended discussion and examples.

The related candisc package provides functions for visualizing a multivariate linear model in a low-dimensional view via a generalized canonical discriminant analyses. heplot.candisc and heplot3d.candisc provide a low-rank 2D (or 3D) view of the effects for a given term in the space of maximum discrimination.

When an element of fill is TRUE, the ellipse outline is drawn using the corresponding color in col, and the interior is filled with a transparent version of this color specified in fill.alpha. To produce filled (non-degenerate) ellipses without the bounding outline, use a value of lty=0 in the corresponding position.

Value

The function invisibly returns an object of class "heplot", with coordinates for the various hypothesis ellipses and the error ellipse, and the limits of the horizontal and vertical axes. These may be useful for adding additional annotations to the plot, using standard plotting functions. (No methods for manipulating these objects are currently available.)

The components are:

H

a list containing the coordinates of each ellipse for the hypothesis terms

E

a matrix containing the coordinates for the error ellipse

center

x,y coordinates of the centroid

xlim

x-axis limits

ylim

y-axis limits

radius

the radius for the unit circles used to generate the ellipses

References

Friendly, M. (2006). Data Ellipses, HE Plots and Reduced-Rank Displays for Multivariate Linear Models: SAS Software and Examples Journal of Statistical Software, 17(6), 1–42. % https://www.jstatsoft.org/v17/i06/, DOI: 10.18637/jss.v017.i06

Friendly, M. (2007). HE plots for Multivariate General Linear Models. Journal of Computational and Graphical Statistics, 16(2) 421–444. http://datavis.ca/papers/jcgs-heplots.pdf

Friendly, Michael (2010). HE Plots for Repeated Measures Designs. Journal of Statistical Software, 37(4), 1-40. DOI: 10.18637/jss.v037.i04.

Fox, J., Friendly, M. & Weisberg, S. (2013). Hypothesis Tests for Multivariate Linear Models Using the car Package. The R Journal, 5(1), https://journal.r-project.org/articles/RJ-2013-004/RJ-2013-004.pdf.

Friendly, M. & Sigal, M. (2014) Recent Advances in Visualizing Multivariate Linear Models. Revista Colombiana de Estadistica, 37, 261-283. DOI: 10.15446/rce.v37n2spe.47934.

See Also

Anova, linearHypothesis for details on testing MLMs.

heplot1d, heplot3d, pairs.mlm, mark.H0 for other HE plot functions. coefplot.mlm for plotting confidence ellipses for parameters in MLMs.

trans.colors for calculation of transparent colors. label.ellipse for labeling positions in plotting H and E ellipses.

candisc, heplot.candisc for reduced-rank views of mlms in canonical space.

Other HE plot functions: heplot1d(), heplot3d(), pairs.mlm()

Examples

## iris data
contrasts(iris$Species) <- matrix(c(0,-1,1, 2, -1, -1), 3,2)
contrasts(iris$Species)

iris.mod <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~
Species, data=iris)

hyp <- list("V:V"="Species1","S:VV"="Species2")
heplot(iris.mod, hypotheses=hyp)
# compare with effect-size scaling
heplot(iris.mod, hypotheses=hyp, size="effect", add=TRUE)

# try filled ellipses; include contrasts
heplot(iris.mod, hypotheses=hyp, fill=TRUE, 
       fill.alpha=0.2, col=c("red", "blue"))
heplot(iris.mod, hypotheses=hyp, fill=TRUE, 
       col=c("red", "blue"), lty=c(0,0,1,1))

# vary label position and fill.alpha
heplot(iris.mod, hypotheses=hyp, fill=TRUE, fill.alpha=c(0.3,0.1), col=c("red", "blue"), 
       lty=c(0,0,1,1), label.pos=0:3)

# what is returned?
hep <-heplot(iris.mod, variables=c(1,3),  hypotheses=hyp)
str(hep)

# all pairs
pairs(iris.mod, hypotheses=hyp, hyp.labels=FALSE)


## Pottery data, from car package
data(Pottery, package = "carData")
pottery.mod <- lm(cbind(Al, Fe, Mg, Ca, Na) ~ Site, data=Pottery)
heplot(pottery.mod)
heplot(pottery.mod, terms=FALSE, add=TRUE, col="blue", 
  hypotheses=list(c("SiteCaldicot = 0", "SiteIsleThorns=0")),
  hyp.labels="Sites Caldicot and Isle Thorns")

## Rohwer data, multivariate multiple regression/ANCOVA
#-- ANCOVA, assuming equal slopes
rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ SES + n + s + ns + na + ss, data=Rohwer)
car::Anova(rohwer.mod)
col <- c("red", "black", "blue", "cyan", "magenta", "brown", "gray")
heplot(rohwer.mod, col=col)

# Add ellipse to test all 5 regressors
heplot(rohwer.mod, hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")), 
       col=col, fill=TRUE)
# View all pairs
pairs(rohwer.mod, hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")))
# or 3D plot

if(requireNamespace("rgl")){
col <- c("pink", "black", "blue", "cyan", "magenta", "brown", "gray")
heplot3d(rohwer.mod, hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")), col=col)
}

Description

This function plots a 1-dimensional representation of the hypothesis (H) and error (E) sums-of-squares-and-products matrices for terms and linear hypotheses in a multivariate linear model.

Usage

heplot1d(mod, ...)

## S3 method for class 'mlm'
heplot1d(
  mod,
  terms,
  hypotheses,
  term.labels = TRUE,
  hyp.labels = TRUE,
  variables = 1,
  error.ellipse = !add,
  factor.means = !add,
  grand.mean = !add,
  remove.intercept = TRUE,
  type = c("II", "III", "2", "3"),
  idata = NULL,
  idesign = NULL,
  icontrasts = c("contr.sum", "contr.poly"),
  imatrix = NULL,
  iterm = NULL,
  manova,
  size = c("evidence", "effect.size", "significance"),
  level = 0.68,
  alpha = 0.05,
  center.pch = "|",
  col = getOption("heplot.colors", c("red", "blue", "black", "darkgreen", "darkcyan",
    "magenta", "brown", "darkgray")),
  lty = 2:1,
  lwd = 1:2,
  xlab,
  main = "",
  xlim,
  axes = TRUE,
  offset.axes = 0.1,
  add = FALSE,
  verbose = FALSE,
  ...
)