---
title: "Multiway Frequency Tables"
---

## Packages

```{r}
library(tidyverse)
```

## Multi-way frequency analysis

-   A study of gender and eyewear-wearing finds the following
    frequencies:
      
```         
gender contacts glasses none
female 121      32      129
male   42       37      85
```

-   Is there association between eyewear and gender?

-   Normally answer this with chisquare test (based on observed and
    expected frequencies from null hypothesis of no association).

-   Two categorical variables and a frequency.

-   We assess in way that generalizes to more categorical variables.

## The data file

```         
gender contacts glasses none
female 121      32      129
male   42       37      85
```

-   This is *not tidy!*

-   Two variables are gender and *eyewear*, and those numbers all
    frequencies.

```{r bMultiway-1, message=F}
my_url <- "http://ritsokiguess.site/datafiles/eyewear.txt"
(eyewear <- read_delim(my_url, " "))
```

## Tidying the data

```{r bMultiway-2, size="footnotesize"}
eyewear %>%
  pivot_longer(contacts:none, names_to="eyewear", 
               values_to="frequency") -> eyes
eyes
```


## Modelling

-   Predict frequency from other factors and combos.
-   `glm` with `poisson` family.

```{r bMultiway-4 }
eyes.1 <- glm(frequency ~ gender * eyewear,
  data = eyes,
  family = "poisson"
)
```

-   Called **log-linear model**.

## What can we get rid of?

\small

```{r bMultiway-5 }
drop1(eyes.1, test = "Chisq")
```

nothing!

## Conclusions

-   `drop1` says what we can remove at this step. Significant = must
    stay.

-   Cannot remove anything.

-   Frequency depends on `gender-eye wear` *combination*, cannot be
    simplified further.

-   Gender and eyewear are *associated*.

-   For modelling, stop here.

## Making a graph

```{r}
#| fig-height: 5
ggplot(eyes, aes(x = gender, y = frequency, fill = eyewear)) + 
  geom_col(position = "fill")
```


## Conclusions

- Females are more likely to wear contacts than males are.
- Females are *less* likely to wear glasses than males are.
- The previous two comments are the reasons for the significant association.


## Code comments 1/2

- The code again:

```{r}
#| eval: false
ggplot(eyes, aes(x = gender, y = frequency, fill = eyewear)) + 
  geom_col(position = "fill")
```

- Variation on two-variable bar chart that we saw in C32.
- Comparing (most easily) *proportions*, so `fill` clearer than `dodge`.
- Each row of dataframe represents many people (the number in `frequency`), so use `geom_col` rather than `geom_bar`.
- `geom_col` takes a `y` that should be the frequency.

## Code comments 2/2

```{r}
#| eval: false
ggplot(eyes, aes(x = gender, y = frequency, fill = eyewear)) + 
  geom_col(position = "fill")
```


- Often in this work, one variable in association is explanatory rather than response. Have that as `x` (here `gender`); `eyewear` is response and goes in `fill`.
- Interpretation: out of each category of explanatory ("out of females"), what proportion in each response category and where do they differ?

## No association

-   Suppose table had been as shown below:


```{r bMultiway-8, message=F}
my_url <- "http://ritsokiguess.site/datafiles/eyewear2.txt"
eyewear2 <- read_table(my_url)
eyewear2
eyewear2 %>% 
   pivot_longer(contacts:none, names_to = "eyewear", 
                values_to = "frequency") -> eyes2
```


## Comments

-   Females and males wear contacts and glasses *in same proportions*
    -   though more females and more contact-wearers.
-   No *association* between gender and eyewear.

## Analysis for revised data

```{r bMultiway-9}
eyes.2 <- glm(frequency ~ gender * eyewear,
  data = eyes2,
  family = "poisson"
)
drop1(eyes.2, test = "Chisq")
```

No longer any association. Take out interaction.

## No interaction


```{r bMultiway-10 }
eyes.3 <- update(eyes.2, . ~ . - gender:eyewear)
drop1(eyes.3, test = "Chisq")
```

-   More females (gender effect) over all eyewear

-   fewer glasses-wearers (eyewear effect) over both genders

-   no association (no interaction).

## Graph shows no association

```{r}
#| fig-height: 5
ggplot(eyes2, aes(x = gender, y = frequency, fill = eyewear)) + 
  geom_col(position = "fill")
```


## Chest pain, being overweight and being a smoker

-   In a hospital emergency department, 176 subjects who attended for
    acute chest pain took part in a study.

-   Each subject had a normal or abnormal electrocardiogram reading
    (ECG), were overweight (as judged by BMI) or not, and were a smoker
    or not.

-   How are these three variables related, or not?

## The data

In modelling-friendly format:

```         

ecg bmi smoke count
abnormal overweight yes 47
abnormal overweight no 10
abnormal normalweight yes 8 
abnormal normalweight no 6
normal overweight yes 25 
normal overweight no 15 
normal normalweight yes 35
normal normalweight no 30
```

## First step

\small

```{r bMultiway-11, message=F}
my_url <- "http://ritsokiguess.site/datafiles/ecg.txt"
chest <- read_delim(my_url, " ")
chest.1 <- glm(count ~ ecg * bmi * smoke,
  data = chest,
  family = "poisson"
)
drop1(chest.1, test = "Chisq")
```

\normalsize

That 3-way interaction comes out.

## Removing the 3-way interaction

\small

```{r bMultiway-12}
chest.2 <- update(chest.1, . ~ . - ecg:bmi:smoke)
drop1(chest.2, test = "Chisq")
```

\normalsize

At $\alpha=0.05$, `bmi:smoke` comes out.

## Removing `bmi:smoke`


```{r bMultiway-13}
chest.3 <- update(chest.2, . ~ . - bmi:smoke)
drop1(chest.3, test = "Chisq")
```


- `ecg:smoke` has become significant. So we have to stop.

- `ecg` is associated with both `bmi` and `smoke`, but separately (it doesn't depend on the combination of `bmi` and `smoke`).

## Understanding the final model

- For each of the significant associations, make a bar chart (here, two-variable because two-way interactions)
- Here, `ecg` is response (patients came into the study being smokers or overweight) so use as `fill` in both graphs.
- `y` is the frequency column.

## `ecg:bmi`

```{r}
#| fig-height: 5
ggplot(chest, aes(x = bmi, y = count, fill = ecg)) + 
  geom_col(position = "fill")
```

## Comment

- Most of the normal weight people had a normal ECG as well, but for the overweight people, a small majority had an abnormal ECG.

## `ecg:smoke`

```{r}
#| fig-height: 5
ggplot(chest, aes(x = smoke, y = count, fill = ecg)) + 
  geom_col(position = "fill")
```

## Comments

-   Most nonsmokers have a normal ECG, but smokers are about 50--50
    normal and abnormal ECG.

-   Don't look at `smoke:bmi` since not significant.

## Simpson's paradox: the airlines example

```         
               Alaska Airlines       America West
Airport       On time    Delayed   On time    Delayed
Los Angeles      497        62        694       117
Phoenix          221        12       4840       415
San Diego        212        20        383        65
San Francisco    503       102        320       129
Seattle         1841       305        201        61

Total           3274       501       6438       787
```

Use `status` as variable name for "on time/delayed".

-   Alaska: 13.3% flights delayed ($501/(3274+501)$).

-   America West: 10.9% ($787/(6438+787)$).

-   America West more punctual, right?

## Arranging the data

-   Can only have single thing in columns, so we have to construct
    column names like this: \small

```         

airport    aa_ontime aa_delayed aw_ontime aw_delayed
LosAngeles   497          62       694        117
Phoenix      221          12      4840        415
SanDiego     212          20       383         65
SanFrancisco 503         102       320        129
Seattle     1841         305       201         61
```

\normalsize

-   Read in:

```{r bMultiway-16, message=F}
my_url <- "http://ritsokiguess.site/datafiles/airlines.txt"
airlines <- read_table(my_url)
```

## Data, as read in

```{r}
airlines
```


## Tidying

-   Some tidying gets us the right layout, with frequencies all in one
    column and the airline and delayed/on time status separated out.
    This uses one of the fancy versions of `pivot_longer`:

```{r bMultiway-17, message=F, size="small"}
airlines %>%
   pivot_longer(-airport, 
                names_to = c("airline", "status"), 
                names_sep = "_", 
                values_to = "freq" ) -> punctual
```

## The data frame `punctual`

```{r bMultiway-18, echo=F}
punctual
```

## Proportions delayed by airline

```{r}
#| fig-height: 5
ggplot(punctual, aes(x = airline, y = freq, fill = status)) + 
  geom_col(position = "fill")
```

## Shrinking the $y$-axis

```{r}
#| fig-height: 5
ggplot(punctual, aes(x = airline, y = freq, fill = status)) + 
  geom_col(position = "fill") +
  coord_cartesian(ylim = c(0.75, 1))
```


## Comment

- Most flights are on time, but Alaska Airlines is late a little more often.


## Proportion delayed by airport, for each airline

We now have *three* categorical variables, so use one of the explanatories (for me, airport) as facets:

## The graph(s)

```{r}
#| fig-height: 5
ggplot(punctual, aes(x = airline, y = freq, fill = status)) +
  geom_col(position = "fill") + facet_wrap(~ airport)
```

## Zoom in on the $y$-scale

```{r}
#| fig-height: 5
ggplot(punctual, aes(x = airline, y = freq, fill = status)) +
  geom_col(position = "fill") + facet_wrap(~ airport) +
  coord_cartesian(ylim = c(0.6, 1))
```


## Simpson's Paradox

-   America West more punctual overall,

-   but worse at *every single* airport!

-   How is that possible?

-   Log-linear analysis sheds some light.

## Model 1 and output

```{r bMultiway-22}
punctual.1 <- glm(freq ~ airport * airline * status,
  data = punctual, family = "poisson"
)
drop1(punctual.1, test = "Chisq")
```

## Remove 3-way interaction

\footnotesize

```{r bMultiway-23}
punctual.2 <- update(punctual.1, ~ . - airport:airline:status)
drop1(punctual.2, test = "Chisq")
```

\normalsize

Stop here, and draw graphs to understand significant results.

##   `airline:status`:

```{r}
#| fig-height: 5

ggplot(punctual, aes(x = airline, y = freq, fill = status)) + 
  geom_col(position = "fill")
```

## Comments

- We did this one before.

- Slightly more of Alaska Airlines' flights delayed overall.

##   `airport:status`:

```{r}
#| fig-height: 5

ggplot(punctual, aes(x = airport, y = freq, fill = status)) + 
  geom_col(position = "fill")
```

## Comments

-   Flights into San Francisco (and maybe Seattle) are often late, and
    flights into Phoenix are usually on time.

-   Considerable variation among airports.

##   `airport:airline`:

```{r}
#| fig-height: 5

ggplot(punctual, aes(x = airport, y = freq, fill = airline)) + 
  geom_col(position = "fill")
```

## Comments

-   What fraction of each airline's flights are to each airport.

-   Most of Alaska Airlines' flights to Seattle and San Francisco.

-   Most of America West's flights to Phoenix.

## The resolution

-   Most of America West's flights to Phoenix, where it is easy to be on
    time.

-   Most of Alaska Airlines' flights to San Francisco and Seattle, where
    it is difficult to be on time.

-   Overall comparison looks bad for Alaska because of this.

-   But, *comparing like with like*, if you compare each airline's
    performance *to the same airport*, Alaska does better.

-   Aggregating over the very different airports was a (big) mistake:
    that was the cause of the Simpson's paradox.

-   Alaska Airlines is *more* punctual when you do the proper
    comparison.

## Ovarian cancer: a four-way table

-   Retrospective study of ovarian cancer done in 1973.

-   Information about 299 women operated on for ovarian cancer 10 years
    previously.

-   Recorded:

    -   stage of cancer (early or advanced)

    -   type of operation (radical or limited)

    -   X-ray treatment received (yes or no)

    -   10-year survival (yes or no)

-   Survival looks like response (suggests logistic regression).

-   Log-linear model finds any associations at all.

## The data

after tidying:

\scriptsize

```         

stage operation xray survival freq
early radical no no 10
early radical no yes 41
early radical yes no 17
early radical yes yes 64
early limited no no 1
early limited no yes 13
early limited yes no 3
early limited yes yes 9
advanced radical no no 38
advanced radical no yes 6
advanced radical yes no 64
advanced radical yes yes 11
advanced limited no no 3
advanced limited no yes 1
advanced limited yes no 13
advanced limited yes yes 5

```

\normalsize

## Reading in data

\small

```{r bMultiway-27, message=F}
my_url <- "http://ritsokiguess.site/datafiles/cancer.txt"
cancer <- read_delim(my_url, " ")
cancer %>% slice(1:6)
```

\normalsize

## Model 1

hopefully looking familiar by now:

```{r bMultiway-28}
cancer.1 <- glm(freq ~ stage * operation * xray * survival,
  data = cancer, family = "poisson")
```

## Output 1

See what we can remove:

\scriptsize

```{r bMultiway-29}
drop1(cancer.1, test = "Chisq")
```

\normalsize

Non-significant interaction can come out.

## Model 2

\scriptsize

```{r bMultiway-30}
cancer.2 <- update(cancer.1, . ~ . - stage:operation:xray:survival)
drop1(cancer.2, test = "Chisq")
```

\normalsize

Least significant term is `stage:xray:survival`: remove.

## Take out `stage:xray:survival`

\scriptsize

```{r bMultiway-31}
cancer.3 <- update(cancer.2, . ~ . - stage:xray:survival)
drop1(cancer.3, test = "Chisq")
```

\normalsize

`operation:xray:survival` comes out next.

## Remove `operation:xray:survival`

\scriptsize

```{r bMultiway-32}
cancer.4 <- update(cancer.3, . ~ . - operation:xray:survival)
drop1(cancer.4, test = "Chisq")
```

\normalsize

## Comments

-   `stage:operation:xray` has now become significant, so won't remove
    that.

-   Shows value of removing terms one at a time.

-   There are no higher-order interactions containing both `xray` and
    `survival`, so now we get to test (and remove) `xray:survival`.

## Remove `xray:survival`

\scriptsize

```{r bMultiway-33}
cancer.5 <- update(cancer.4, . ~ . - xray:survival)
drop1(cancer.5, test = "Chisq")
```

\normalsize

## Remove `stage:operation:survival`

\scriptsize

```{r bMultiway-34}
cancer.6 <- update(cancer.5, . ~ . - stage:operation:survival)
drop1(cancer.6, test = "Chisq")
```

\normalsize

## Last step?

Remove `operation:survival`.

\footnotesize

```{r bMultiway-35, size="footnotesize"}
cancer.7 <- update(cancer.6, . ~ . - operation:survival)
drop1(cancer.7, test = "Chisq")
```

\normalsize

Finally done!

## Conclusions

-   What matters is things associated with `survival` (`survival` is
    "response").

-   Only significant such term is `stage:survival`.


## The graph

```{r}
#| fig-height: 5

ggplot(cancer, aes(x = stage, y = freq, fill = survival)) + 
  geom_col(position = "fill")
```

## Comments

-   Most people in early stage of cancer survived, and most people in
    advanced stage did not survive.

-   This true *regardless* of type of operation or whether or not X-ray
    treatment was received. These things have no impact on survival.

## What about that other interaction?

```{r}
#| fig-height: 5

ggplot(cancer, aes(x = stage, y = freq, fill = xray)) + 
  geom_col(position = "fill") + facet_wrap(~ operation)
```

## Comments

-   The association is between `stage` and `xray` *only for those who
    had the limited operation*.

-   For those who had the radical operation, there was no association
    between `stage` and `xray`.

-   This is of less interest than associations with `survival`.

## General procedure

-   Start with "complete model" including all possible interactions.

-   `drop1` gives highest-order interaction(s) remaining, remove least
    non-significant.

-   Repeat as necessary until everything significant.

-   Look at graphs of significant interactions.

-   Main effects not usually very interesting.

-   Interactions with "response" usually of most interest: show
    association with response.
