ETC5521: Exploratory Data Analysis

.info-box.w-50.bg-white[
These slides are viewed best by Chrome or Firefox and occasionally need to be refreshed if elements did not load properly. See <a href=lecture-05B.pdf>here for the PDF <i class="fas fa-file-pdf"></i></a>. 
]

<br>

---

# .monash-blue[ETC5521: Exploratory Data Analysis]

<br>

<h2 style="font-weight:900!important;">Working with a single variable, making transformations, detecting outliers, using robust statistics</h2>

.bottom_abs.width100[

Lecturer: *Di Cook*

<i class="fas fa-envelope"></i>  ETC5521.Clayton-x@monash.edu

<i class="fas fa-calendar-alt"></i> Week 5 - Session 2

<br>

]

---

# Categorical variables

This lecture is based on Chapter 4  of <br><br>Unwin (2015) Graphical Data Analysis with R

---

# There are two types of categorical variables

.monash-blue[**Nominal**] where there is no intrinsic ordering to the categories<br>
**E.g.** blue, grey, black, white.

<br>

.monash-blue[**Ordinal**] where there is a clear order to the categories.<Br>
**E.g.** Strongly disagree, disagree, neutral, agree, strongly agree.

---

# Categorical variables in R

```r
data <- c(2, 2, 1, 1, 3, 3, 3, 1)
factor(data)
```

```
## [1] 2 2 1 1 3 3 3 1
## Levels: 1 2 3
```
]
* You can easily change the labels of the variables:
.f4[

```r
factor(data, labels = c("I", "II", "III"))
```

```
## [1] II  II  I   I   III III III I  
## Levels: I II III
```
]
]
.item.f4[

]
]

* Order of the factors are determined by the input:

```r
*# numerical input are ordered in increasing order
factor(c(1, 3, 10))
```

```
## [1] 1  3  10
## Levels: 1 3 10
```

```r
*# character input are ordered alphabetically
factor(c("1", "3", "10"))
```

```
## [1] 1  3  10
## Levels: 1 10 3
```

```r
*# you can specify order of levels explicitly
factor(c("1", "3", "10"),
  levels = c("1", "3", "10")
)
```

```
## [1] 1  3  10
## Levels: 1 3 10
```

---

# Numerical factors in R

```r
x <- factor(c(10, 20, 30, 10, 20))
mean(x)
```

```
## Warning in mean.default(x): argument is not numeric or logical: returning NA
```

```
## [1] NA
```

<i class="fas fa-exclamation-triangle"></i> `as.numeric` function returns the internal integer values of the factor

```r
mean(as.numeric(x))
```

```
## [1] 1.8
```

You probably want to use:

```r
mean(as.numeric(levels(x)[x]))
```

```
## [1] 18
```

]
.w-50[

```r
mean(as.numeric(as.character(x)))
```

```
## [1] 18
```

]
]

---
# Numerical summaries: counts, proportions, percentages and odds

```
## # A tibble: 22 × 7
##    country   iso3   year count      p   pct  odds
##    <chr>     <chr> <dbl> <dbl>  <dbl> <dbl> <dbl>
##  1 Australia AUS    2000   982 0.0522  5.22 1    
##  2 Australia AUS    2001   953 0.0507  5.07 0.970
##  3 Australia AUS    2002  1008 0.0536  5.36 1.03 
##  4 Australia AUS    2003   926 0.0493  4.93 0.943
##  5 Australia AUS    2004  1036 0.0551  5.51 1.05 
##  6 Australia AUS    2005  1030 0.0548  5.48 1.05 
##  7 Australia AUS    2006  1127 0.0600  6.00 1.15 
##  8 Australia AUS    2007  1081 0.0575  5.75 1.10 
##  9 Australia AUS    2008  1182 0.0629  6.29 1.20 
## 10 Australia AUS    2009  1176 0.0626  6.26 1.20 
## 11 Australia AUS    2010  1146 0.0610  6.10 1.17 
## 12 Australia AUS    2011  1202 0.0640  6.40 1.22 
## 13 Australia AUS    2012  1259 0.0670  6.70 1.28 
## 14 Australia AUS    2013   512 0.0272  2.72 0.521
## 15 Australia AUS    2014   474 0.0252  2.52 0.483
## 16 Australia AUS    2015   438 0.0233  2.33 0.446
## 17 Australia AUS    2016   481 0.0256  2.56 0.490
## 18 Australia AUS    2017   524 0.0279  2.79 0.534
## 19 Australia AUS    2018   502 0.0267  2.67 0.511
## 20 Australia AUS    2019   554 0.0295  2.95 0.564
## 21 Australia AUS    2020   609 0.0324  3.24 0.620
## 22 Australia AUS    2021   593 0.0316  3.16 0.604
```
]]
]
.item.w-40[

For qualitative data, compute

- count/frequency, 
- proportion/percentage
- and sometimes, an **odds ratio**. Here we have used ratio relative to the count in year 2000.

<br><br>
**Note:** For exploration, no rounding of digits was done, but to report you would need to make the numbers pretty.
]]
---

# .orange[Revisiting Case study] .circle.bg-orange.white[1] 2019 Australian Federal Election

.flex[
.w-50[
<img src="images/week4B/aus-election-plot1-1.png" width="432" style="display: block; margin: auto;" />

]
.w-50[
<img src="images/week4B/aus-election-plot2-1.png" width="432" style="display: block; margin: auto;" />
]

]

]
.panel[.panel-name[data]
.scroll-sign[.s400.f4[

```r
df1 <- read_csv(here::here("data/HouseFirstPrefsByCandidateByVoteTypeDownload-24310.csv"),
  skip = 1,
  col_types = cols(
    .default = col_character(),
    OrdinaryVotes = col_double(),
    AbsentVotes = col_double(),
    ProvisionalVotes = col_double(),
    PrePollVotes = col_double(),
    PostalVotes = col_double(),
    TotalVotes = col_double(),
    Swing = col_double()
  )
)
tdf3 <- df1 %>%
  group_by(DivisionID) %>%
  summarise(
    DivisionNm = unique(DivisionNm),
    State = unique(StateAb),
    votes_GRN = TotalVotes[which(PartyAb == "GRN")],
    votes_total = sum(TotalVotes)
  ) %>%
  mutate(perc_GRN = votes_GRN / votes_total * 100)
```
]]]
.panel[.panel-name[R]
.f4[

```r
tdf3 %>%
  ggplot(aes(perc_GRN, State)) +
  ggbeeswarm::geom_quasirandom(groupOnX = FALSE, varwidth = TRUE) +
  labs(
    x = "Percentage of first preference votes per division",
    y = "State",
    title = "First preference votes for the Greens party"
  )
tdf3 %>%
  mutate(State = fct_reorder(State, perc_GRN)) %>%
  ggplot(aes(perc_GRN, State)) +
  ggbeeswarm::geom_quasirandom(groupOnX = FALSE, varwidth = TRUE) +
  labs(
    x = "Percentage of first preference votes per division",
    y = "State",
    title = "First preference votes for the Greens party"
  )
```
]]
]

--
<br><br>
<center>Sorting levels is (almost) always better when plotting</center>

---

# Order nominal variables meaningfully

<i class="fas fa-code"></i> **Coding tip**: use below functions to easily change the order of factor levels

```r
stats::reorder(factor, value, mean)
forcats::fct_reorder(factor, value, median)
forcats::fct_reorder2(factor, value1, value2, func)
```

---

# .orange[Case study] .circle.bg-orange.white[6] Aspirin use after heart attack

.grid[
.item[
<img src="images/week4B/meta-plot1-1.png" width="432" style="display: block; margin: auto;" />

]
.item[

* Meta-analysis is a statistical analysis that combines the results of multiple scientific studies.
* This data studies the use of aspirin for death prevention after myocardial infarction, or in plain terms, a heart attack.
* The ISIS-2 study has more patients than all other studies combined.
* You could consider lumping the categories with low frequencies together.
]

]

]
.panel[.panel-name[data]
.h300.f4.scroll-sign[

```r
data("Fleiss93", package = "meta")
df6 <- Fleiss93 %>%
  mutate(total = n.e + n.c)
skimr::skim(df6)
```

```
## ── Data Summary ────────────────────────
##                            Values
## Name                       df6   
## Number of rows             7     
## Number of columns          7     
## _______________________          
## Column type frequency:           
##   character                1     
##   numeric                  6     
## ________________________         
## Group variables            None  
## 
## ── Variable type: character ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate min max empty n_unique whitespace
## 1 study                 0             1   3   6     0        7          0
## 
## ── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate  mean      sd   p0    p25  p50   p75  p100 hist 
## 1 year                  0             1 1979.    4.39 1974 1978.  1979 1980   1988 ▇▇▇▁▃
## 2 event.e               0             1  304   563.     32   46.5   85  174   1570 ▇▁▁▁▁
## 3 n.e                   0             1 2027. 2959.    317  686.   810 1550.  8587 ▇▂▁▁▂
## 4 event.c               0             1  327.  618.     38   58     67  172.  1720 ▇▁▁▁▁
## 5 n.c                   0             1 1974. 2993.    309  515    771 1554.  8600 ▇▂▁▁▂
## 6 total                 0             1 4000. 5950.    626 1228.  1529 3103  17187 ▇▂▁▁▂
```
]]
.panel[.panel-name[R]
.f4[

```r
df6 %>%
  mutate(study = fct_reorder(study, desc(total))) %>%
  ggplot(aes(study, total)) +
  geom_col() +
  labs(x = "", y = "Frequency") +
  guides(x = guide_axis(n.dodge = 2))

df6 %>%
  mutate(
    study = ifelse(total < 2000, "Other", study),
    study = fct_reorder(study, desc(total))
  ) %>%
  ggplot(aes(study, total)) +
  geom_col() +
  labs(x = "", y = "Frequency")
```
]]
]

.f5.footnote[
Fleiss JL (1993): The statistical basis of meta-analysis. *Statistical Methods in Medical Research* **2** 121–145<br>
Balduzzi S, Rücker G, Schwarzer G (2019), How to perform a meta-analysis with R: a practical tutorial, Evidence-Based Mental Health.
]

---

# Consider combining factor levels with low frequencies

<i class="fas fa-code"></i> **Coding tip**: the following family of functions help to easily lump factor levels together:

```r
forcats::fct_lump()
forcats::fct_lump_lowfreq()
forcats::fct_lump_min()
forcats::fct_lump_n()
forcats::fct_lump_prop()
# if conditioned on another variable
ifelse(cond, "Other", factor)
dplyr::case_when(
  cond1 ~ "level1",
  cond2 ~ "level2",
  TRUE ~ "Other"
)
```

---

# .orange[Case study] .circle.bg-orange.white[7] Anorexia

.grid[
.item[
<img src="images/week4B/anorexia-plot1-1.png" width="432" style="display: block; margin: auto;" />

<table class=" lightable-classic" style='font-family: "Arial Narrow", "Source Sans Pro", sans-serif; width: auto !important; margin-left: auto; margin-right: auto;'>
 <thead>
  <tr>
   <th style="text-align:left;"> Treatment </th>
   <th style="text-align:right;"> Frequency </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> CBT </td>
   <td style="text-align:right;"> 29 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Cont </td>
   <td style="text-align:right;"> 26 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> FT </td>
   <td style="text-align:right;"> 17 </td>
  </tr>
</tbody>
</table>

]
.item[

**Table or Plot?**

]

]
.panel[.panel-name[data]
.h200.f4.scroll-sign[

```r
data(anorexia, package = "MASS")

df9tab <- table(anorexia$Treat) %>%
  as.data.frame() %>%
  rename(Treatment = Var1, Frequency = Freq)

skimr::skim(anorexia)
```

```
## ── Data Summary ────────────────────────
##                            Values  
## Name                       anorexia
## Number of rows             72      
## Number of columns          3       
## _______________________            
## Column type frequency:             
##   factor                   1       
##   numeric                  2       
## ________________________           
## Group variables            None    
## 
## ── Variable type: factor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate ordered n_unique top_counts              
## 1 Treat                 0             1 FALSE          3 CBT: 29, Con: 26, FT: 17
## 
## ── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate mean   sd   p0  p25  p50  p75  p100 hist 
## 1 Prewt                 0             1 82.4 5.18 70   79.6 82.3 86    94.9 ▂▅▇▆▁
## 2 Postwt                0             1 85.2 8.04 71.3 79.3 84.1 91.6 104.  ▆▇▅▆▂
```
]]
.panel[.panel-name[R]
.f4[

```r
ggplot(anorexia, aes(Treat)) +
  geom_bar() +
  labs(x = "", y = "Frequency")
```

```r
ggplot(anorexia, aes(Treat)) +
  stat_count(geom = "point", size = 4) +
  stat_count(geom = "line", group = 1) +
  labs(y = "Frequency", x = "")
```

]]
]

.f6.footnote[

Hand, D. J., Daly, F., McConway, K., Lunn, D. and Ostrowski, E. eds (1993) A Handbook of Small Data Sets. Chapman & Hall, Data set 285 (p. 229) <br>
Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0

]

* Table for accuracy, plot for visual communication

**Why not a point or line?**

<div class="f4">
<ul>
<li>This can be appropriate depending on what you want to communicate </li>
<li>A barplot occupies more area compared to a point and the area does a better job of communicating size</li>
<li>A line is suggestive of a trend </li>
</ul>
</div>

---

# .orange[Case study] .circle.bg-orange.white[8] Titanic

.flex[
.w-50[
.flex[
.w-50[
<img src="images/week4B/titanic-plot1-1.png" width="288" style="display: block; margin: auto;" />

]
.w-50[
<img src="images/week4B/titanic-plot3-1.png" width="216" style="display: block; margin: auto;" />

]
]

]
.w-50[

**What does the graphs for each categorical variable tell us?**

* There were more crews than 1st to 3rd class passengers
* There were far more males on ship; possibly because majority of crew members were male. You can further explore this by constructing two-way tables or graphs that consider both variables.
* Most passengers were adults. 
* More than two-thirds of passengers died.

]

]
.panel[.panel-name[data]
.h350.f4.scroll-sign[

```r
df9 <- as_tibble(Titanic)
skimr::skim(df9)
```

```
## ── Data Summary ────────────────────────
##                            Values
## Name                       df9   
## Number of rows             32    
## Number of columns          5     
## _______________________          
## Column type frequency:           
##   character                4     
##   numeric                  1     
## ________________________         
## Group variables            None  
## 
## ── Variable type: character ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate min max empty n_unique whitespace
## 1 Class                 0             1   3   4     0        4          0
## 2 Sex                   0             1   4   6     0        2          0
## 3 Age                   0             1   5   5     0        2          0
## 4 Survived              0             1   2   3     0        2          0
## 
## ── Variable type: numeric ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate mean   sd p0  p25  p50 p75 p100 hist 
## 1 n                     0             1 68.8 136.  0 0.75 13.5  77  670 ▇▁▁▁▁
```
]]
.panel[.panel-name[R]
.f4.scroll-sign[.s500[

```r
df9 %>%
  group_by(Class) %>%
  summarise(total = sum(n)) %>%
  ggplot(aes(Class, total)) +
  geom_col(fill = "#ee64a4") +
  labs(x = "", y = "Frequency")

df9 %>%
  group_by(Sex) %>%
  summarise(total = sum(n)) %>%
  ggplot(aes(Sex, total)) +
  geom_col(fill = "#746FB2") +
  labs(x = "", y = "Frequency")

df9 %>%
  group_by(Age) %>%
  summarise(total = sum(n)) %>%
  ggplot(aes(Age, total)) +
  geom_col(fill = "#C8008F") +
  labs(x = "", y = "Frequency")

df9 %>%
  group_by(Survived) %>%
  summarise(total = sum(n)) %>%
  ggplot(aes(Survived, total)) +
  geom_col(fill = "#795549") +
  labs(x = "Survived", y = "Frequency")
```
]]
]]

.f5.footnote[
British Board of Trade (1990), Report on the Loss of the ‘Titanic’ (S.S.). British Board of Trade Inquiry Report (reprint). Gloucester, UK: Allan Sutton Publishing
]

---

# Coloring bars

* Colour here doesn't add information as the x-axis already tells us about the categories, but colouring bars can make it more visually appealing, and provide subtle clues. 
* If you have too many categories colour won't work well to differentiate the categories.

---

# .orange[Case study] .circle.bg-orange.white[9] Opinion poll in Ireland Aug 2013

.grid[
.item[
<img src="images/week4B/poll-plot1-1.png" width="504" style="display: block; margin: auto;" />

]
.item[

* Pie chart is popular in mainstream media but are not generally recommended as people are generally poor at comparing angles.
* 3D pie charts should definitely be avoided!
* Here you can see that there are many people that are "Undecided" for which political party to support and failing to account for this paints a different picture.

]

]
.panel[.panel-name[data]
.f4[

```r
df9 <- tibble(
  party = c(
    "Fine Gael", "Labour", "Fianna Fail",
    "Sinn Fein", "Indeps", "Green", "Undecided"
  ),
  nos = c(181, 51, 171, 119, 91, 4, 368)
)
df9v2 <- df9 %>% filter(party != "Undecided")
df9
```

```
## # A tibble: 7 × 2
##   party         nos
##   <chr>       <dbl>
## 1 Fine Gael     181
## 2 Labour         51
## 3 Fianna Fail   171
## 4 Sinn Fein     119
## 5 Indeps         91
## 6 Green           4
## 7 Undecided     368
```
]]
.panel[.panel-name[R]
.f4[

```r
g9 <- df9 %>%
  ggplot(aes("", nos, fill = party)) +
  geom_col(color = "black") +
  labs(y = "", x = "") +
  coord_polar("y") +
  theme(
    axis.line = element_blank(),
    axis.line.y = element_blank(),
    axis.text = element_blank(),
    panel.grid.major = element_blank()
  ) +
  scale_fill_discrete_qualitative(name = "Party")
g9
g9 %+% df9v2 +
  # below is needed to keep the same color scheme as before
  scale_fill_manual(values = qualitative_hcl(7)[1:6])
```
]]]

---

# Piechart is a stacked barplot just with a transformed coordinate system

```r
df <- data.frame(var = c("A", "B", "C"), perc = c(40, 40, 20))
g <- ggplot(df, aes("", perc, fill = var)) +
  geom_col()
g
```

<img src="images/week4B/barplot-1.png" width="216" style="display: block; margin: auto;" />
--

```r
g + coord_polar("y")
```

---

# Roseplot is a barplot just with a transformed coordinate system

```r
dummy <- data.frame(
  var = LETTERS[1:20],
  n = round(rexp(20, 1 / 100))
)
g <- ggplot(dummy, aes(var, n)) +
  geom_col(fill = "pink", color = "black")
g
```

<img src="images/week4B/nonstacked-barplot-1.png" width="720" style="display: block; margin: auto;" />
--

```r
g + coord_polar("x") + theme_void()
```

---
# Visual inference

.item.w-50[

Typical plot description:

```r
ggplot(data, aes(x=var1)) +
  geom_col()

ggplot(data, aes(x=var1)) +
  geom_bar()
```

<br><br>
*Is the distribution consistent with a sample from a binomial distribution with a given p?*

]
.item.w-5[
.white[space]
]
.item.w-45[

Potential simulation method from binomial

```r
# Only one option
null_dist("var1", "binom", 
  list(size=n, p=phat))
```
]
]

---
# Lineup of tuberculosis count between sexes

.panelset[
.panel[.panel-name[📊]
<img src="images/week4B/tb-lineup-1.png" width="80%" style="display: block; margin: auto;" />

.f5[**Note**: This is nothing more than you can learn from a conventional hypothesis test of `$$H_0: p=0.5$$`. .monash-orange2[**Stay tuned for more interesting visual inference lineups in coming weeks!**]]
]
.panel[.panel-name[R]
.scroll-sign[.s500.f5[

```r
library(nullabor)
set.seed(252)
tb_oz_2012 <- tb %>%
  filter(iso3 == "AUS",
         year == 2012) %>%
  select(iso3, year, new_sp_m04:new_ep_m65) %>%
  pivot_longer(cols=new_sp_m04:new_ep_m65, 
               names_to = "var", 
               values_to = "count") %>%
  separate(var, into=c("new", "type", "sexage")) %>%
  select(-new) %>%
  filter(!(sexage %in% c("sexunk014", "sexunk04", 
                         "sexunk15plus", "sexunk514",
                         "f15plus", "m15plus"))) %>%
  group_by(sexage) %>%
  summarise(count = sum(count, na.rm=TRUE)) %>%
  mutate(sex=str_sub(sexage, 1, 1),
         age=str_sub(sexage, 2, str_length(sexage))) %>%
  group_by(sex) %>%
  summarise(count=sum(count)) %>%
  ungroup() %>%
  mutate(sex01 = ifelse(sex=="m", 0, 1)) %>%
  select(-sex)

ggplot(lineup(null_dist("count", "binom", 
                        list(size=sum(tb_oz_2012$count),
                               p=0.5)), 
              tb_oz_2012, n=10),
       aes(x=sex01, y=count)) +
  geom_col() +
  facet_wrap(~.sample, ncol=5) +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid.major = element_blank())
```
]]
]
]

---

# Take away messages

<ul class="fa-ul">
{{content}}
</ul>

]
]

<li><span class="fa-li"><i class="fas fa-paper-plane"></i></span>Again, be prepared to do multiple plots</li>
{{content}}
--

<li><span class="fa-li"><i class="fas fa-paper-plane"></i></span>Changing bins or binwidth/bandwidth in histogram, violin or density plots can paint a different picture</li>
{{content}}
--

<li><span class="fa-li"><i class="fas fa-paper-plane"></i></span>Consider different representations of categorical variables (reordering meaningfully, lumping low frequencies together, plot or table, pie or barplot, missing categories)</li>

---
# Resources and Acknowledgement

- Slides originally written by Emi Tanaka and constructed with [`xaringan`](https://github.com/yihui/xaringan), [remark.js](https://remarkjs.com), [`knitr`](http://yihui.name/knitr), and [R Markdown](https://rmarkdown.rstudio.com).

---

background-size: cover
class: title-slide
background-image: url("images/bg-01.png")

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a>.

.bottom_abs.width100[

Lecturer: *Di Cook*

<i class="fas fa-envelope"></i>  ETC5521.Clayton-x@monash.edu

<i class="fas fa-calendar-alt"></i> Week 5 - Session 2

<br>

]