title: "4. Mosaic plots"

author: "Michael Friendly"

date: "`r Sys.Date()`"

package: vcdExtra

output:

rmarkdown::html_vignette:

fig_caption: yes

bibliography: ["vcd.bib", "vcdExtra.bib"]

csl: apa.csl

vignette: >

```{r setup, include = FALSE}

Mosaic plots provide an ideal method both for visualizing contingency tables and for

visualizing the fit--- or more importantly--- **lack of fit** of a loglinear model.

For a two-way table, `mosaic()`, by default, fits a model of independence, $[A][B]$

or `~A + B` as an R formula. The `vcdExtra` package extends this to models fit

using `glm(..., family=poisson)`, which can include specialized models for

ordered factors, or square tables that are intermediate between the saturated model,

$[A B]$ = `A * B`, and the independence model $[A][B]$.

For $n$-way tables, `vcd::mosaic()` can fit any loglinear model, and can also be

used to plot a model fit with `MASS:loglm()`. The `vcdExtra` package extends this

to models fit using `stats::glm()` and, by extension, to non-linear models fit

using the [gnm package](https://cran.r-project.org/package=gnm).

See @vcd:Friendly:1994, @vcd:Friendly:1999 for the statistical ideas behind these

uses of mosaic displays in connection with loglinear models. Our book @FriendlyMeyer:2016:DDAR

gives a detailed discussion of mosaic plots and many more examples.

The essential ideas are to:

* recursively sub-divide a unit square into rectangular "tiles" for the

cells of the table, such that the area of each tile is proportional to the cell frequency. Tiles are split in a sequential order:

+ First according to the **marginal** proportions of a first variable, V1

+ Next according to the **conditional** proportions of a 2nd variable, V2 | V1

+ Next according to the **conditional** proportions of a 3rd variable, V3 | {V1, V2}

+ ...

* For a given loglinear model, the tiles can then be shaded in various ways to reflect

the residuals (lack of fit) for a given model.

* The pattern of residuals can then

be used to suggest a better model or understand *where* a given model fits or

does not fit.

`mosaic()` provides a wide range of options for the directions of splitting,

the specification of shading, labeling, spacing, legend and many other details.

It is actually implemented as a special case of a more general

class of displays for $n$-way tables called `strucplot`, including

sieve diagrams, association plots, double-decker plots as well as mosaic

plots.

For details, see `help(strucplot)` and the "See also" links therein,

and also @vcd:Meyer+Zeileis+Hornik:2006b, which is available as

an R vignette via `vignette("strucplot", package="vcd")`.

***Example***:

A mosaic plot for the Arthritis treatment data fits the model of independence,

`~ Treatment + Improved` and displays the association in the pattern of

residual shading. The goal is to visualize the difference in the proportions

of `Improved` for the two levels of `Treatment` : "Placebo" and "Treated".

The plot below is produced with the following call to `mosaic()`.

With the first split by `Treatment` and the shading used, it is easy to see

that more people given the placebo experienced no improvement, while more

people given the active treatment reported marked improvement.

`gp = shading_max` specifies that color in the plot signals a significant residual at a 90% or 99% significance level,

with the more intense shade for 99%.

Note that the residuals for the independence model were not large

(as shown in the legend),

yet the association between `Treatment` and `Improved`

is highly significant.

```{r, art1}

In contrast, one of the other shading schemes, from @vcd:Friendly:1994

(use: `gp = shading_Friendly`),

uses fixed cutoffs of $\pm 2, \pm 4$,

to shade cells which are *individually* significant

at approximately $\alpha = 0.05$ and $\alpha = 0.001$ levels, respectively.

The plot below uses `gp = shading_Friendly`.

Mosaic plots using tables or frequency data frames as input typically take the levels of the

table variables in the order presented in the dataset. For character variables, this is often

alphabetical order. That might be helpful for looking up a value, but is unhelpful for seeing

and understanding the pattern of association.

It is usually much better to order the levels of the row and column variables to help reveal

the nature of their association. This is an example of **effect ordering for data display**

[@FriendlyKwan:02:effect].

***Example***:

Data from @Glass:54 gave this 5 x 5 table on the occupations of 3500 British fathers and their sons, where the occupational categories are listed in alphabetic order.

```{r glass}

The mosaic display shows very strong association, but aside from the

diagonal cells, the pattern is unclear. Note the use of

`set_varnames` to give more descriptive labels for the variables

and abbreviate the occupational category labels.

and `interpolate` to set the shading levels for the mosaic.

```{r glass-mosaic1}

The occupational categories differ in **status**, and can be reordered correctly as follows,

from `Professional` down to `Unskilled`.

```{r glass-order}

The revised mosaic plot can be produced by indexing the rows and columns of

the table using `ord`.

```{r glass-mosaic2}

From this, and for the examples in the next section, it is useful to re-define

`father` and `son` as **ordered** factors in the original `Glass` frequency data.frame.

```{r glass-ord}

For mobility tables such as this, where the rows and columns refer to the same

occupational categories it comes as no surprise that there is a strong association

in the diagonal cells: most often, sons remain in the same occupational categories

as their fathers.

However, the re-ordered mosaic display also reveals something subtler:

when a son differs in occupation from the father, it is more likely that

he will appear in a category one-step removed than more steps removed.

The residuals seem to decrease with the number of steps from the diagonal.

For such tables, specialized loglinear models provide interesting cases

intermediate between the independence model, [A] [B],

and the saturated model, [A B]. These can be fit using `glm()`, with the data

in frequency form,

where `assoc` is a special term to handle a restricted form of association, different from

`A:B` which specifies the saturated model in this notation.

* **Quasi-independence**: Asserts independence, but ignores the diagonal cells by

fitting them exactly. The loglinear model is:

$\log m_{ij} = \mu + \lambda^A_i + \lambda^B_j + \delta_i I(i = j)$,

where $I()$ is the indicator function.

* **Symmetry**: This model asserts that the joint distribution of the row and column variables

is symmetric, that is $\pi_{ij} = \pi_{ji}$:

A son is equally likely to move from their father's occupational category $i$ to another

category, $j$, as the reverse, moving from $j$ to $i$.

Symmetry is quite strong, because it also implies

**marginal homogeneity**, that the marginal probabilities of the row and column variables

are equal,

$\pi{i+} = \sum_j \pi_{ij} = \sum_j \pi_{ji} = \pi_{+i}$ for all $i$.

* **Quasi-symmetry**: This model uses the standard main-effect terms in the loglinear

model, but asserts that the association parameters are symmetric,

$\log m_{ij} = \mu + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij}$,

where $\lambda^{AB}_{ij} = \lambda^{AB}_{ji}$.

The [gnm package](https://cran.r-project.org/package=gnm) provides a variety of these functions:

`gnm::Diag()`, `gnm::Symm()` and `gnm::Topo()` for an interaction factor as specified by an array of levels, which may be arbitrarily structured.

For example, the following generates a term for a diagonal factor in a $4 \times 4$

table. The diagonal values reflect parameters fitted for each diagonal cell. Off-diagonal

values, "." are ignored.

```{r diag}

`Symm()` constructs parameters for symmetric cells. The particular values don't matter.

All that does matter is that the same value, e.g., `1:2` appears in both the (1,2) and

(2,1) cells.

```{r symm}

***Example***:

To illustrate, we fit the four models below, starting with the independence model

`Freq ~ father + son` and then adding terms to reflect the restricted

forms of association, e.g., `Diag(father, son)` for diagonal terms and

`Symm(father, son)` for symmetry.

```{r glass-models}

We can visualize these using the `vcdExtra::mosaic.glm()` method, which extends

mosaic displays to handle fitted `glm` objects. *Technical note*: for

models fitted using `glm()`, standardized residuals, `residuals_type="rstandard"`

have better statistical properties than the default Pearson residuals in mosaic

plots and analysis.

```{r glass-quasi}

Mosaic plots for the other models would give further visual assessment of these models,

however we can also test differences among them. For nested models, `anova()` gives tests

of how much better a more complex model is compared to the previous one.

```{r glass-anova}

Alternatively, `vcdExtra::LRstats()` gives model summaries for a collection of

models, not necessarily nested, with AIC and BIC statistics reflecting

model parsimony.

```{r glass-lrstats}

By all criteria, the model of quasi symmetry fits best. The residual deviance $G^2

is not significant. The mosaic is largely unshaded, indicating a good fit, but there

are a few shaded cells that indicate the remaining positive and negative residuals.

For comparative mosaic displays, it is sometimes useful to show the $G^2$ statistic

in the main title, using `vcdExtra::modFit()` for this purpose.

```{r glass-qsymm}

When natural orders for row and column levels are not given a priori, we can find

orderings that make more sense using correspondence analysis.

The general ideas are that:

* Correspondence analysis assigns scores to the row and column variables to best account for the association in 1, 2, ... dimensions

* The first CA dimension accounts for largest proportion of the Pearson $\chi^2$

* Therefore, permuting the levels of the row and column variables by the CA Dim1 scores gives a more coherent mosaic plot, more clearly showing the nature of the association.

* The [seriation package](https://cran.r-project.org/package=seriation) now has a method to order variables in frequency tables using CA.

***Example***:

As an example, consider the `HouseTasks` dataset, a 13 x 4 table of frequencies of household tasks performed by couples, either by the `Husband`, `Wife`, `Alternating` or `Jointly`.

You can see from the table that some tasks (Repairs) are done largely by the husband; some (laundry, main meal)

are largely done by the wife, while others are done jointly or alternating between husband and

wife. But the `Task` and `Who` levels are both in alphabetical order.

```{r housetasks}

The naive mosaic plot for this dataset is shown below, splitting first by `Task` and then by `Who`. Due to the length of the factor labels, some features of `labeling` were used to make the display more readable.

```{r housetasks-mos1}

Correspondence analysis, using the [ca package](https://cran.r-project.org/package=ca),

shows that nearly 89% of the $\chi^2$ can be accounted for in two dimensions.

```{r housetasks-ca}

The CA plot has a fairly simple interpretation: Dim1 is largely the distinction between

tasks primarily done by the wife vs. the husband. Dim2 distinguishes tasks that are done

singly vs. those that are done jointly.

```{r housetasks-ca-plot}

So, we can use the `CA` method of `seriation::seriate()` to find the order of permutations of

`Task` and `Who` along the CA dimensions.

```{r housetasks-seriation}

Now, use `seriation::permute()` to use `order` for the permutations of `Task` and `Who`,

and plot the resulting mosaic:

```{r housetasks-mos2}

It is now easy to see the cluster of tasks (laundry and cooking) done largely by the wife

at the top, and those (repairs, driving) done largely by the husband at the bottom.