```{r, include = FALSE}
# Load libraries 
source("libraries.R")

# Set up chunk for all slides
source("chunk_options_and_themes.R")
```

## Outline

- Wrangling using verbs
- Plotting data using a grammar
- Mapping for data analysis
- Statistical inference with data visualisation
- Perceptual principles for good graphic design
- Representing uncertainty

## dplyr verbs

There are five primary dplyr **verbs**, representing distinct data analysis tasks:

- `filter`: Extract specified rows of a data frame
- `arrange`: Reorder the rows of a data frame
- `select`: Select specified columns of a data frame
- `mutate`: Add new or transform a column of a data frame
- `summarise`: Create collapsed summaries of a data frame
- (`group_by`: Introduce structure to a data frame)
 
 
## Filter

:::: {.columns}
::: {.column width=30%}

::: {style="font-size: 70%;"}
Data from an experiment on the taste of french fries, cooked using three different oils over a 10 week period. There were 12 people tasting two fries (replicates) from each batch.
:::

:::
::: {.column width=70%}

Get a subset of the observations (horizontal slice)

```{r}
#| label: select_a_subset_of_the_observations
#| code-line-highlight: "3"
load("data/french_fries.rda")
french_fries |>
    filter(subject == 3, time == 1) 
ff_long <- french_fries |> 
  pivot_longer(potato:painty, 
               names_to = "type", 
               values_to = "rating")
```
:::
::::


## Arrange

Order the observations (hierarchically)

```{r order the observations}
#| code-line-numbers: true
#| code-line-highlight: "2"
french_fries |>
    arrange(desc(rancid)) |> 
    head()
```

## Select

Select a subset of the variables (vertical selection)

```{r select a subset of the variables}
#| code-line-numbers: "2"
french_fries |>
    select(time, treatment, subject, rep, potato) |> 
    head()
```

## Summarise

Summarise observations into a (set of) one-number statistic(s)

```{r summarize observations into one-number statistic}
#| code-line-numbers: "2,5"
french_fries |>
    summarise( 
      mean_rancid = mean(rancid, na.rm=TRUE), 
      sd_rancid = sd(rancid, na.rm = TRUE)
      ) 
```

## Summarise and group_by

```{r summarise and group_by}
french_fries |>
    group_by(time, treatment) |>
    summarise(mean_rancid = mean(rancid), sd_rancid = sd(rancid))
```

## Let's use these tools

To answer these french fry experiment questions:

- Is the design complete?
- Are replicates like each other? 
- How do the ratings on the different scales differ?
- Are raters giving different scores on average?
- Do ratings change over the weeks?

## Completeness 

- If the data is complete it should be 12 x 10 x 3 x 2, that is, 6 records for each person in each week. 

- *To check*: tabulate  number of records for each subject, time and treatment. See `n()`.

:::: {.columns}
::: {.column}

```{r}
#| label:  checking design completeness
french_fries |> 
  group_by(subject) |> 
  summarize(n = n()) 
```
:::
::: {.column}

*Three subjects skipped a week.* 

::: {.fragment style="font-size: 70%;"}

<br>
**Other nice short cuts**

instead of `group_by(subject) |> summarize(n = n())` we can use:

- `group_by(subject) |> tally()`
- `count(subject)`

:::
:::
::::

## Counts for subject by time

```{r counts for subject by time}
options(width=120)
french_fries |>
  count(subject, time) |>
  pivot_wider(names_from="time", values_from="n")
```

## How do scores change over time?

::: {.columns}

::: {.column width=30%}

```{r do-scores-change-over-time, fig.show='hide'}
ggplot(data=ff_long, aes(x=time, y=rating, colour=treatment)) +
  geom_point() +
  facet_grid(subject~type) 
```

:::
::: {.column width=70%}

```{r ref.label="do-scores-change-over-time", echo=FALSE, fig.width=12, fig.height=9, out.width="100%"}
```

:::
::::

## Plot differently

:::: {.columns}
::: {.column width=30%}

Get summary of ratings over replicates and connect the dots

::: 
::: {.column width=70%}

```{r ff line plots, fig.width=12, fig.height=9, out.width="80%"}
#| code-fold: true
ff_av <- ff_long |> 
  group_by(subject, time, type, treatment) |>
  summarise(rating=mean(rating))

ggplot(data=ff_long, aes(x=time, y=rating, colour=treatment)) + 
  facet_grid(subject~type) +
  geom_line(data=ff_av, aes(group=treatment))
```

:::
::::

##  {.information-slide .center style="text-align: center;"}

With tidy data, wrangling is a logical and organised process. 

## Plotting data using a grammar {.transition-slide .center style="text-align: center;"}

## Grammar of graphics

:::: {.columns}
::: {.column}
What is a data plot?  

- **data**
- **aesthetics: mapping of variables to graphical elements**
- **geom**: type of plot structure to use
- **transformations**: log scale, ...
- **layers**: multiple geoms, multiple data sets, annotation
- **facets**: show subsets in different plots
- **themes**: modifying style
:::
::: {.column width=10%}
:::
::: {.column width=40%}

::: {.fragment}

<img src="images/ggplot2.png" width="20%" />

```
ggplot(data = <DATA>) + 
  <GEOM_FUNCTION>(
     mapping = aes(<MAPPINGS>),
     stat = <STAT>, 
     position = <POSITION>
  ) +
  <COORDINATE_FUNCTION> +
  <FACET_FUNCTION> +
  <THEME>
```
:::
:::
::::

## Why use a grammar?

- Change the script a little, make all sorts of types of plots.

- Instead of plots being like unique animals in a zoo, the script is like having the genetic footprint to see how different plots are related.

- There is a tight connection with statistical thinking. This makes it possible to do statistical inference and make valid decisions from data plots.


## Example: Tuberculosis data

(Current) TB case notifications data from [WHO](http://www.who.int/tb/country/data/download/en/). Also available via R package [`getTBinR`](https://github.com/seabbs/getTBinR). Here we have Indonesia incidence 1995-2012.

```{r}
#| label: read-and-wrangle-TB-data
#| echo: false
tb <- readr::read_csv("data/TB_notifications_2025-07-22.csv") |>
  dplyr::select(country, iso3, year, new_sp_m04:new_sp_fu) |>
  pivot_longer(new_sp_m04:new_sp_fu, names_to="stuff", values_to="count") |>
  separate(stuff, c("stuff1", "stuff2", "sexage")) |>
  dplyr::select(-stuff1, -stuff2) |>
  mutate(sex=substr(sexage, 1, 1), 
         age=substr(sexage, 2, length(sexage))) |>
  dplyr::select(-sexage)

# Filter data to get the US
tb_idn <- tb |> 
  filter(iso3 == "IDN") |>
  filter(!(age %in% c("04", "014", "514", "u"))) |>
  filter(year > 2000, year < 2013) |>
  mutate(
    age_group = factor(age, labels = c("15-24", "25-34", "35-44", "45-54", "55-64", "65+"))
  )
```

```{r}
#| label: make-a-barchart-of-US-TB-incidence
#| code-fold: true
#| out-width: 80%
#| fig-width: 10
#| fig-height: 3
ggplot(tb_idn, aes(x = year, 
                  y = count, 
                  fill = sex)) +
  geom_bar(stat = "identity") +
  facet_grid(~age) 
```


- What do you learn about tuberculosis incidences in Indonesia from this plot?
- Give three changes to the plot that would improve it. 

## Compare males and females

```{r}
#| label: compare-proportions-of-males-females
#| out-width: 80%
#| fig-width: 10
#| fig-height: 3
#| code-line-numbers: "3,4"
#| code-fold: true
ggplot(tb_idn, aes(x=year, y=count, fill=sex)) +  
  geom_bar(stat="identity", position="fill") + 
  ylab("proportion") + 
  facet_grid(~age_group) +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10"))
```

::: {style="font-size: 80%;"}
- What do we learn about the data that is different from the previous plot?
- What is easier and what is harder or impossible to learn from this arrangement?
:::


::: notes
- Focus is now on proportions of male and female each year, within age group
- Proportions are similar across year 
- Roughly equal proportions at young and old age groups, more male incidence in middle years

:::

## Separate plots

:::: {.columns}
::: {.column width="70%"}

```{r}
#| label: compare-counts-of-males-females
#| out-width: 100%
#| code-fold: true
#| fig-height: 5
#| fig-width: 10
#| code-line-numbers: "5"
# Make separate plots for males and females, focus on counts by category
ggplot(tb_idn, aes(x=year, y=count, fill=sex)) +
  geom_bar(stat="identity") +
  facet_grid(sex~age_group) + 
  scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"))
```
:::
::: {.column width=30%}
Focus is now separately on the temporal trend for each sex and age group. What do we see?

::: {.fragment}
Trend was tapering in younger age groups, but it is steadily increasing for over 65.
:::

:::
::::


## Make a pie

:::: {.columns}
::: {.column width="70%"}

```{r}
#| label: rose-plot-of-males-females
#| code-fold: true
#| fig-height: 5
#| fig-width: 10
#| code-line-numbers: "6"
# How to make a pie instead of a barchart - not straight forward
ggplot(tb_idn, aes(x=year, y=count, fill=sex)) +
  geom_bar(stat="identity") + 
  facet_grid(sex~age_group) + 
  coord_polar() + 
  scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"))
```

:::
::: {.column width=30%}

::: {.fragment}
This isn't a pie, it's a [rose plot](https://datavizcatalogue.com/methods/nightingale_rose_chart.html)!
:::
:::
::::


## Stacked bar

:::: {.columns}
::: {.column width="70%"}

```{r}
#| label: stacked-barchart-of-males-females
#| code-fold: true
#| fig-height: 5
#| fig-width: 8
#| out-width: 100%
#| code-line-numbers: "3"
# Step 1 to make the pie
ggplot(tb_idn, aes(x = 1, y = count, fill = factor(year))) +
  geom_bar(stat="identity", position="fill") + 
  facet_grid(sex~age_group) +
  scale_fill_viridis_d("", option="inferno") 
```
:::
::: {.column width=30%}

Step 1: turn it into a stacked bar.

:::
::::

## Pie chart

:::: {.columns}
::: {.column width="70%"}

```{r}
#| label: pie-chart-of-males-females
#| out-width: 100%
#| code-fold: true
#| fig-height: 6
#| fig-width: 10
#| code-line-numbers: "3,4,7"
# Now we have a pie, note the mapping of variables
# and the modification to the coord_polar
ggplot(tb_idn, aes(x = 1, y = count, fill = factor(year))) + 
  geom_bar(stat="identity", position="fill") + 
  facet_grid(sex~age_group) +
  scale_fill_viridis_d("", option="inferno") +
  coord_polar(theta = "y") 
```

:::
::: {.column width=30%}

- What are the pros, and cons, of using the pie chart for this data?
- Would it be better if the pies used age for the segments, and facetted by year (and sex)?

```{r}
#| eval: false
#| out-width: 60%
#| code-fold: true
#| fig-height: 5
#| code-line-numbers: "1,3"
# age for segments, facet by year and sex
ggplot(tb_idn, aes(x = 1, y = count, fill = factor(age_group))) + 
  geom_bar(stat="identity", position="fill") + 
  facet_grid(sex~year) + 
  scale_fill_viridis_d("", option="inferno") +
  coord_polar(theta = "y") +
  theme(legend.position="bottom")
```


:::
::::

## Line plot vs barchart

:::: {.columns}
::: {.column width="70%"}

```{r use a line plot instead of bar, fig.width=10, fig.height=3, out.width="100%"}
#| code-fold: true
ggplot(tb_idn, aes(x=year, y=count, colour=sex)) +
  geom_line() + 
  geom_point() +
  facet_grid(~age_group) +
  ylim(c(0,NA)) + 
  scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"),
                     limits=c(2000, 2012))
```
:::

::: {.column width=30%}

- We can read counts for both sexes
- Males and females can be directly compared
- Temporal trend is more visible
:::
::::

## Line plot vs barchart

:::: {.columns}
::: {.column width="70%"}

```{r use a line plot of proportions, fig.width=10, fig.height=3, out.width="100%"}
#| code-fold: true
tb_idn |> group_by(year, age_group) |> 
  summarise(p = count[sex=="m"]/sum(count)) |>
  ggplot(aes(x=year, y=p)) +
    geom_hline(yintercept = 0.50, colour="grey50", linewidth=2) +
    geom_line() + 
    geom_point() +
    facet_grid(~age_group) +
    ylab("Proportion of Males") + 
    scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"),
                     limits=c(2000, 2012))
```
:::

::: {.column width=30%}

- Attention is forced to proportion of males
- Direct comparison of counts within year and age
- Equal proportion guideline provides a baseline for comparison

:::
::::

## Avoid "connect the dots" plots

:::: {.columns}
::: {.column width="70%"}

```{r use a smoother instead of bar, fig.width=10, fig.height=3, out.width="100%"}
#| code-fold: true
ggplot(tb_idn, aes(x=year, y=count, colour=sex)) +
  geom_point() +
  geom_smooth(se=F) + 
  facet_grid(~age_group) +
  ylim(c(0,NA)) + 
  scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"),
                     limits=c(2000, 2012))
```
:::

::: {.column width=30%}

- Focus on the trend
- Actual counts might vary a little
- Temporal trend is even clearer
:::
::::

## Avoid "connect the dots" plots

:::: {.columns}
::: {.column width="70%"}

```{r use a smoother of proportions, fig.width=10, fig.height=3, out.width="100%"}
#| code-fold: true
tb_idn |> group_by(year, age_group) |> 
  summarise(p = count[sex=="m"]/sum(count)) |>
  ggplot(aes(x=year, y=p)) +
    geom_hline(yintercept = 0.50, colour="grey50", linewidth=2) +
    geom_point() +
    geom_smooth(se=F) + 
    facet_grid(~age_group) +
    ylab("Proportion of Males") + 
    scale_x_continuous("year", breaks=seq(2000, 2012, 5), 
                     labels=c("00", "05", "10"),
                     limits=c(2000, 2012))
```
:::

::: {.column width=30%}

- Attention is trend of proportion of males

:::
::::

## {.information-slide .center style="text-align: center;"}

With a small change in code, many different plots of the same data were generated.  It allowed us to see different patterns and relationships in the data.

## `gg` extensions

::: {.smaller}
- There are more than 150 extensions to `ggplot2`
- Many adhere to the grammar, to define new types of displays, like 
    - `ggdist`: including representation of uncertainty and error
    - `gganimate`: specifies animations as layers
- Others are supporting packages, such as 
    - `patchwork` for laying out multiple displays
    - `ggthemes` for styling plots
    
:::

[https://exts.ggplot2.tidyverse.org/gallery/](https://exts.ggplot2.tidyverse.org/gallery/)

## Exploring distributions and uncertainty

`ggdist`

:::: {.columns}
::: {.column width="10%"}

<img src="https://mjskay.github.io/ggdist/logo.svg" width="80%">

:::

::: {.column width="45%"}

::: { .smaller}

```{r}
#| label: tb-comparisons1
#| fig-width: 5
#| fig-height: 4
#| code-fold: true
tb_inc_100k <- read_csv("data/TB_burden_countries_2025-07-22.csv") |>
  filter(iso3 %in% c("USA", "AUS"))
ggplot(tb_inc_100k, aes(y = iso3, 
                        x = e_inc_100k)) +
  stat_gradientinterval(fill = "darkorange") +
  ylab("") +
  xlab("Inc per 100k") +
  theme_ggdist()
```
:::
:::

::: {.column width="45%"}

::: {.smaller}
 
```{r}
#| label: tb-comparisons2
#| fig-width: 5
#| fig-height: 4
#| code-fold: true
ggplot(tb_inc_100k, aes(y = iso3, 
                        x = e_inc_100k)) +
  stat_halfeye(side = "right") +
  geom_dots(side="left", 
                    fill = "darkorange", color = "darkorange") +
  ylab("") +
  xlab("Inc per 100k") +
  theme_ggdist()
```

:::
:::
::::

## Distribution, variation, uncertainty [(1/3)]{.smallest}

:::: {.columns}

::: {.column}

::: {.callout-note appearance="simple"}

::: {style="font-size: 150%;"}

The most valuable way to show uncertainty is to show all the data.

:::

:::

Plot shows first preference % for greens in the 2019 Australian Federal Election, for the 150 electorates. 

Plot of choice is the **jittered dotplot**, where points are spread vertically according to density.


:::
::: {.column}

::: {style="font-size: 70%;"}
```{r}
#| label: election1
#| code-fold: true
election <- read_csv("data/election2019.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()
  )
)
e_grn <- election |>
  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)

e_grn |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(x=perc_GRN, y=State)) +
    geom_quasirandom(groupOnX = FALSE, varwidth = TRUE) +
    labs(
      x = "First preference votes %",
      y = ""
    ) +
  xlim(c(0,50))

```
:::

:::
::::

## Distribution, variation, uncertainty  [(2/3)]{.smallest}

:::: {.columns}

::: {.column}

```{r}
#| label: election1
#| echo: false
```

:::
::: {.column}

What do we learn?

::: {.incremental}
- Different number of observations in each state
- One outlier in Vic
- As a group, ACT has higher %'s
- Vic has a small cluster of points with higher %'s
- %'s are mostly very low
:::

::: {.fragment}
This plot ONLY shows uncertainty! 
:::

:::
::::

## Distribution, variation, uncertainty  [(3/3)]{.smallest}

:::: {.columns}

::: {.column}

What would be other common ways to display this data?

::: {.incremental}
- Side-by-side boxplots
- Side-by-side violin
- On a map of electorates
:::

::: {.fragment}
For each plot think about 

- what is uncertainty, and what is estimate
- what the plot shows or hides

:::

:::
::: {.column}


::: {.panel-tabset}

## Dotplot

```{r}
#| label: election1
#| echo: false
#| out-width: 80%
```

## Boxplot

::: {style="font-size: 70%;"}
```{r}
#| label: election2
#| out-width: 80%
#| code-fold: true
e_grn |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(x=perc_GRN, y=State)) +
    geom_boxplot(varwidth = TRUE) +
    labs(
      x = "First preference votes %",
      y = ""
    ) +
  xlim(c(0,50))
```
:::

## Violin

::: {style="font-size: 70%;"}
```{r}
#| label: election3
#| out-width: 80%
#| code-fold: true
e_grn |>
  mutate(State = fct_reorder(State, perc_GRN)) |>
  ggplot(aes(x=perc_GRN, y=State)) +
    geom_violin(draw_quantiles = c(0.25, 0.5, 0.75),
      fill="#006dae", alpha=0.5) +
    labs(
      x = "First preference votes %",
      y = ""
    ) +
  xlim(c(0,50))
```
:::

## Map

::: {style="font-size: 70%;"}
```{r}
#| label: election4
#| fig-width: 8
#| fig-height: 6
#| out-width: 100%
#| code-fold: true
oz_states <- ozmaps::ozmap_states %>% filter(NAME != "Other Territories")
oz_votes <- rmapshaper::ms_simplify(ozmaps::abs_ced)
oz_votes_grn <- full_join(oz_votes, e_grn, by=c("NAME"="DivisionNm"))

ggplot(oz_votes_grn, aes(fill=perc_GRN)) +
  geom_sf(colour="white") +
  scale_fill_viridis_c(direction=-1, trans = "log", 
    guide = "colourbar", 
    labels = scales::label_number(accuracy = 0.1)) +
  theme_map() +
  theme(legend.position = "right", 
    legend.title = element_blank())
```
:::


:::

:::
::::

## Models and uncertainty 

::: {.panel-tabset}

## model

:::: {.columns}
::: {.column width=20%}

Plotting the fitted model alone

:::
::: {.column}

::: {style="font-size: 70%;"}

```{r}
#| label: model
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
#| code-fold: true
data("wages")
wages_fct <- wages |>
  select(id, ln_wages, xp, high_grade) |>
  mutate(high_grade = factor(high_grade))
wages_fit <- lmer(ln_wages~xp + high_grade + (xp|id), data=wages_fct)
wages_fe <- summary(wages_fit)$coefficients
wages_fe_d <- tibble(xp = rep(seq(0, 13, 1), 7),
     high_grade = rep(c(6, 7, 8, 9, 10, 11, 12), rep(14, 7))) |>
  mutate(ln_wages = case_when(
    high_grade == 6 ~ wages_fe[1,1] + wages_fe[2,1]*xp,
    high_grade == 7 ~ wages_fe[1,1] + wages_fe[3,1] + wages_fe[2,1]*xp,
    high_grade == 8 ~ wages_fe[1,1] + wages_fe[4,1]  + wages_fe[2,1]*xp,
    high_grade == 9 ~ wages_fe[1,1] + wages_fe[5,1]  + wages_fe[2,1]*xp,
    high_grade == 10 ~ wages_fe[1,1] + wages_fe[6,1]  + wages_fe[2,1]*xp,
    high_grade == 11 ~ wages_fe[1,1] + wages_fe[7,1]  + wages_fe[2,1]*xp,
    high_grade == 12 ~ wages_fe[1,1] + wages_fe[8,1]  + wages_fe[2,1]*xp)
  ) |>
  mutate(high_grade = factor(high_grade))
ggplot(wages_fe_d) + 
  geom_line(aes(x=xp, 
                y=ln_wages, 
                colour=high_grade, 
                group=high_grade)) +
  scale_colour_discrete_divergingx(palette = "Zissou 1") +
  labs(x="Experience (years)", y="Wages (ln)", colour="Grade") 
```
:::

:::
::::

## and data

:::: {.columns}
::: {.column width=40%}

Adding the data makes the model look less impressive

:::
::: {.column}

::: {style="font-size: 70%;"}

```{r}
#| label: modelanddata
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
#| code-fold: true
ggplot() + 
  geom_line(data=wages_fct, aes(x=xp, y=ln_wages, group=id), alpha=0.1) +
  geom_line(data=wages_fe_d, aes(x=xp, 
                y=ln_wages, 
                colour=high_grade, 
                group=high_grade)) +
  scale_colour_discrete_divergingx(palette = "Zissou 1") +
  labs(x="Experience (years)", y="Wages (ln)", colour="Grade") 
```

:::
:::
::::

## and SE

:::: {.columns}
::: {.column width=40%}

Standard errors as produced by the model fit. 
:::


::: {.column }

::: {style="font-size: 70%;"}

```{r}
#| label: SE
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
#| code-fold: true
wages_fe_d <- wages_fe_d |>
  mutate(ln_wages_l = case_when(
    high_grade == 6 ~ wages_fe[1,1] - wages_fe[1,2] +
                      (wages_fe[2,1]-wages_fe[2,2])*xp ,
    high_grade == 7 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[3,1] - wages_fe[3,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp,
    high_grade == 8 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[4,1] - wages_fe[4,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp,
    high_grade == 9 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[5,1] - wages_fe[5,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp,
    high_grade == 10 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[6,1] - wages_fe[6,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp,
    high_grade == 11 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[7,1] - wages_fe[7,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp,
    high_grade == 12 ~ wages_fe[1,1] - wages_fe[1,2] + 
                      wages_fe[8,1] - wages_fe[8,2] + 
                      (wages_fe[2,1]-wages_fe[2,2])*xp)
  ) |>
  mutate(ln_wages_u = case_when(
    high_grade == 6 ~ wages_fe[1,1] + wages_fe[1,2] +
                      (wages_fe[2,1]+wages_fe[2,2])*xp ,
    high_grade == 7 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[3,1] + wages_fe[3,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp,
    high_grade == 8 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[4,1] + wages_fe[4,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp,
    high_grade == 9 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[5,1] + wages_fe[5,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp,
    high_grade == 10 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[6,1] + wages_fe[6,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp,
    high_grade == 11 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[7,1] + wages_fe[7,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp,
    high_grade == 12 ~ wages_fe[1,1] + wages_fe[1,2] + 
                      wages_fe[8,1] + wages_fe[8,2] + 
                      (wages_fe[2,1]+wages_fe[2,2])*xp)
  ) 

ggplot() + 
  geom_ribbon(data=wages_fe_d, 
            aes(x=xp, 
                ymin=ln_wages_l,
                ymax=ln_wages_u,
                fill=high_grade), colour=NA, alpha=0.1) +
  geom_line(data=wages_fe_d, 
            aes(x=xp, 
                y=ln_wages, 
                colour=high_grade, 
                group=high_grade)) +
  scale_fill_discrete_divergingx(palette = "Zissou 1") +
  scale_colour_discrete_divergingx(palette = "Zissou 1") +
  labs(x="Experience (years)", y="Wages (ln)", colour="Grade", fill="Grade") 

```

:::

:::
::::

## individual fits

:::: {.columns}
::: {.column width=40%}

Examine individual fits. Too many to show all - sample multiple times to get through them all. 

:::

::: {.column }

::: {style="font-size: 70%;"}

```{r}
#| label: indiv-fits
#| fig-width: 8
#| fig-height: 5
#| out-width: 100%
#| code-fold: true
wages_full <- wages_fct |>
  add_predictions(wages_fit, 
                  var = "pred") |>
  add_residuals(wages_fit, 
                var = "res")
set.seed(1222)
wages_full |> add_n_obs() |> filter(n_obs > 4) |>
  sample_n_keys(size = 12) |>
  ggplot() + 
  geom_line(aes(x = xp, y = pred, group = id, 
             colour = factor(id))) + 
  geom_point(aes(x = xp, y = ln_wages, 
                 colour = factor(id))) + 
  facet_wrap(~id, ncol=4)  +
  scale_x_continuous("Experience (years)", 
    breaks=seq(0, 12, 2)) +
  ylab("Wages (ln)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```
:::

:::
::::


## which is honest?

::: {style="font-size: 70%;"}

- Plotting the model alone suggests higher education increases wages (by XXX), particularly if the student has 12 years (full high school).
- Adding the individual level observations shows the large individual-to-individual variability, which swamps the education differences. It also suggests the individual observations might not have linear increase.
- Adding representation of the standard error shows the "overlap" between education levels, that the estimate for wage increase with 12 years of education is not substantially better than 10 years, but it is better than 6 years. 
- Displaying the observed and fitted values separately by individual, examines individual level uncertainty.

Note that, this particular model fit makes the assumption that wages increase linearly with years of experience. 
:::

:::


## Common measures and representations 

:::: {.columns}
::: {.column width=15%}

::: {style="font-size: 70%;"}

Melbourne pedestrian counts at Southern Cross Station, Sunday Aug 31, 2025.

:::

:::

::: {.column width=85%}

::: {.panel-tabset}

```{r}
#| eval: false
#| echo: false
ped <- rwalkr::melb_walk(from = as.Date("2025-08-01"), 
                         to = as.Date("2025-08-31"))
ped <- ped |>
  mutate(wday = wday(Date, label=TRUE, 
      abbr=TRUE, week_start=1)) |>
  mutate(Count = as.numeric(Count))
save(ped, file="data/ped_Aug2025.rda")
```

```{r}
#| label: pedestrians
#| echo: false
load("data/ped_Aug2025.rda")
ped_sc <- ped |>
  filter(Sensor == "Southern Cross Station") |>
  filter(Date == ymd("2025-08-31")) |>
  group_by(Time) |>
  summarise(Count = sum(Count), .groups = "drop") |>
  mutate(se = sqrt(Count))
b1 <- ggplot(ped_sc, aes(x=Time, y=Count)) +
  geom_col(fill = "#20794D") +
  xlab("Hour")
b2 <- ggplot(ped_sc, aes(x=Time, y=Count)) +
  geom_col(fill = "#b9ca4a") +
  geom_errorbar(aes(ymin = Count - se, ymax = Count + se),
    width=0.5, colour="#20794D") +
  xlab("Hour")
b3 <- ggplot(ped_sc, aes(x=Time,
    ydist=distributional::dist_normal(Count, se))) +
  stat_pointinterval(colour = "#20794D") +
  xlab("Hour") + ylab("Count")
b4 <- ggplot(ped_sc, aes(x=Time,
    ydist=distributional::dist_normal(Count, se))) +
  stat_gradientinterval(colour = NA, fill="#20794D", 
    .width=1) +
  geom_line(aes(x=Time, y=Count), colour="#20794D") +
  xlab("Hour") + ylab("Count")
b5 <- ggplot(ped_sc, aes(x=Time, y=Count)) +
  geom_ribbon(aes(ymin = Count - qnorm(0.975)*se, 
                  ymax = Count + qnorm(0.975)*se),
    fill = "#b9ca4a") +
  geom_line(colour="#20794D") +
  xlab("Hour")
ped_sc_ci <- ped_sc |>
  mutate(l50 = Count - qnorm(0.75)*se,
         u50 = Count + qnorm(0.75)*se,
         l80 = Count - qnorm(0.9)*se,
         u80 = Count + qnorm(0.9)*se,
         l99 = Count - qnorm(0.995)*se,
         u99 = Count + qnorm(0.995)*se
  ) |>
  pivot_longer(cols=l50:u99, names_to = "intprob", 
    values_to="value") |>
  mutate(bound = str_sub(intprob, 1, 1),
       prob = str_sub(intprob, 2, 3)) |>
  select(Time, Count, se, prob, bound, value) |>
  pivot_wider(names_from = bound, values_from = value)
b6 <- ggplot(ped_sc_ci, aes(x=Time, y=Count)) +
  geom_lineribbon(aes(ymin = l, ymax = u, fill = prob)) +
  labs(x="Hour", fill="Confidence") +
  scale_fill_discrete_sequential(palette = "Greens", 
    rev=FALSE, n=5)  
b7 <- ggplot(ped_sc, aes(x=Time, y=Count)) +
  geom_smooth(colour = "#20794D", fill = "#b9ca4a") +
  geom_point(colour = "#20794D") +
  xlab("Hour")
```

## [Bars]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b1
```

## [bar+CI]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b2
```

## [CI]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b3
```

## [Gradient]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b4
```

## [Ribbon]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b5
```

## [LineRibbon]{.smaller} 

```{r}
#| echo: false
#| fig-width: 9
#| fig-height: 5
#| out-width: 70%
b6
```

## [Loess]{.smaller} 

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 5
#| out-width: 70%
b7
```

## [Code]{.smaller} 

```{r}
#| label: pedestrians
#| echo: true
#| eval: false
```


:::

:::
::::

## Forecasting

:::: {.columns}
::: {.column}

This is a forecast of business trips to Melbourne for 2018-2022, based on data from 2000-2017. How is the uncertainty represented?

::: {style="font-size: 70%;"}

```{r}
#| label: forecast1
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
#| code-fold: true
tourism_melb <- tourism %>%
  filter(Region == "Melbourne", Purpose == "Business")
fit <- tourism_melb %>%
  model(
    ets = ETS(Trips ~ trend("A"))
  )
fc <- fit |>
  forecast(h = "5 years")
fc |>
  autoplot(tourism_melb) +
    theme(aspect.ratio = 0.6)
```

:::

:::
::: {.column}

::: {.fragment}

Representation is by simulating and displaying a number of possible forecasts.

::: {style="font-size: 70%;"}

```{r}
#| label: forecast2
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
#| code-fold: true
fc_b <- fit |>
  forecast(h = "5 years", bootstrap = TRUE)
fc_b_samples <- fit |>
  generate(h = 20, times = 50, bootstrap = TRUE)
fc_b_samples |>
  ggplot() +
    geom_line(aes(x = Quarter, y = .sim, group = .rep),
      colour = "#027EB6", alpha=0.1) +
    geom_line(data=fc_b, aes(x = Quarter, y = .mean), 
      colour = "#027EB6") +
    autolayer(tourism_melb, Trips) +
    ylab("Trips") +
    theme(aspect.ratio = 0.6)
``` 

:::

:::
:::
::::


## Mapping and data {.transition-slide .center style="text-align: center;"}

## Map thinning

The big change from working with maps in a GIS and maps for data analysis is the SIZE of the map data.

We are going to demonstrate how you need to change your approach, using the code in `map.R`.

The shapefiles for this example are downloaded from [https://imcarto.webflow.io/gdb](https://drive.google.com/file/d/1pYYMFtB6rZXiuZtysahr78Foy3n3KRQ3/view). It is the "Province" data for Indonesian Administrative boundaries.

## Displaying statistics with a map

:::: {.columns}
::: {.column style="font-size: 80%;"}

A choropleth map is used to show a [measured variable]{.monash-blue2} associated with a political or geographic [region]{.monash-blue2}. Polygons for the region are filled with colour. 

The purpose is to examine the spatial distribution of a variable.


```{r}
#| label: setup-choro
#| echo: false
library(sugarbag)

invthm <- theme_map() + 
  theme(
    panel.background = element_rect(fill = "black", colour = NA), 
    plot.background = element_rect(fill = "black", colour = NA),
    legend.background = element_rect(fill = "transparent", colour = NA),
    legend.key = element_rect(fill = "transparent", colour = NA),
    text = element_text(colour = "white"),
    axis.text = element_blank()
  )

# function to allocate colours to regions
aus_colours <- function(sir_p50){
  value <- case_when(
    sir_p50 <  0.74 ~ "#33809d",
    sir_p50 >= 0.74 & sir_p50 < 0.98 ~ "#aec6c7",
    sir_p50 >= 0.98 & sir_p50 < 1.05 ~ "#fff4bc",
    sir_p50 >= 1.05 & sir_p50 < 1.45 ~ "#ff9a64",
    sir_p50 >= 1.45 ~ "#ff3500",
    TRUE ~ "#FFFFFF")
  return(value)
}
```

```{r }
#| label: thyroiddata
#| code-fold: true
#| eval: false
sa2 <- strayr::read_absmap("sa22011") |> 
  filter(!st_is_empty(geometry)) |> 
  filter(!state_name_2011 == "Other Territories") |> 
  filter(!sa2_name_2011 == "Lord Howe Island")
sa2 <- sa2 |> rmapshaper::ms_simplify(keep = 0.5, keep_shapes = TRUE) # Simplify the map!!!
SIR <- read_csv("data/SIR Downloadable Data.csv") |> 
  filter(SA2_name %in% sa2$sa2_name_2011) |> 
  dplyr::select(Cancer_name, SA2_name, Sex_name, p50) |> 
  filter(Cancer_name == "Thyroid", Sex_name == "Females")
ERP <- read_csv("data/ERP.csv") |>
  filter(REGIONTYPE == "SA2", Time == 2011, Region %in% SIR$SA2_name) |> 
  dplyr::select(Region, Value)
# Alternative maps
# Join with sa2 sf object
sa2thyroid_ERP <- SIR |> 
  left_join(sa2, ., by = c("sa2_name_2011" = "SA2_name")) |>
  left_join(., ERP |> 
              dplyr::select(Region, 
              Population = Value), by = c("sa2_name_2011"= "Region")) |> 
  filter(!st_is_empty(geometry))
sa2thyroid_ERP <- sa2thyroid_ERP |> 
  #filter(!is.na(Population)) |> 
  filter(!sa2_name_2011 == "Lord Howe Island") |> 
  mutate(SIR = map_chr(p50, aus_colours)) |> 
  st_as_sf() 
save(sa2, file="data/sa2.rda")
save(sa2thyroid_ERP, file="data/sa2thyroid_ERP.rda")
```

:::
::: {.column style="font-size: 80%;"}

```{r}
#| label: choro
#| code-fold: true
#| fig-width: 10
#| fig-height: 8
#| out-width: 100%
# Plot the choropleth
load("data/sa2thyroid_ERP.rda")
aus_ggchoro <- ggplot(sa2thyroid_ERP) + 
  geom_sf(aes(fill = SIR), size = 0.1) + 
  scale_fill_identity() + invthm
aus_ggchoro
```

:::
::::

## The problem with choropleth maps

The problem is that [high density population areas may be very small geographically]{.monash-blue2}. They can [disappear]{.monash-orange2} in a choropleth map, which means that we get a biased sense of the spatial distribution of a variable.


## Cartograms

:::: {.columns}
::: {.column width=40%}

A [cartogram](https://www.r-graph-gallery.com/cartogram.html) transforms the [geographic shape to match the value of a statistic or the population]{.monash-blue2}. Its a useful exploratory technique for examining the spatial distribution of a measured variable. 

<br>
<br>
BUT they [don't work for Australia]{.monash-orange2}.
:::
::: {.column width=60% .smallest}

```{r}
#| label: cartogram
#| code-fold: true
#| fig-width: 9
#| fig-height: 10
#| out-width: 70%
# transform to NAD83 / UTM zone 16N
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet=TRUE)

nc <- nc |>
  mutate(lBIR79 = log(BIR79))
nc_utm <- st_transform(nc, 26916)

orig <- ggplot(nc) + 
  geom_sf(aes(fill = lBIR79)) +
  ggtitle("choropleth") +
  theme_map() +
  theme(legend.position = "none")

nc_utm_carto <- cartogram_cont(nc_utm, weight = "BIR74", itermax = 5)

carto <- ggplot(nc_utm_carto) + 
  geom_sf(aes(fill = lBIR79)) +
  ggtitle("cartogram") +
  theme_map() +
  theme(legend.position = "none")

nc_utm_dorl <- cartogram_dorling(nc_utm, weight = "BIR74")

dorl <- ggplot(nc_utm_dorl) + 
  geom_sf(aes(fill = lBIR79)) +
  ggtitle("dorling") +
  theme_map() +
  theme(legend.position = "none")

orig + carto + dorl + plot_layout(ncol=1)
```

:::
::::

## Hexagon tile

:::: {.columns}
::: {.column}
A hexagon tile map represents every spatial polygon with an equal sized hexagon. In dense areas these will be tesselated, but separated hexagons are placed at centroids of the remote spatial regions.

<br>

Now the higher thyroid incidence in Perth suburbs, some Melbourne suburbs, and Sydney are more visible.
:::
::: {.column}

```{r}
#| label: hexmap
#| code-fold: true
#| fig-width: 10
#| fig-height: 8
#| out-width: 100%
if (!file.exists("data/aus_hexmap.rda")) {
  
## Create centroids set
centroids <- sa2 |> 
  create_centroids(., "sa2_name_2011")
## Create hexagon grid
grid <- create_grid(centroids = centroids,
                    hex_size = 0.2,
                    buffer_dist = 5)
## Allocate polygon centroids to hexagon grid points
aus_hexmap <- allocate(
  centroids = centroids,
  hex_grid = grid,
  sf_id = "sa2_name_2011",
  ## same column used in create_centroids
  hex_size = 0.2,
  ## same size used in create_grid
  hex_filter = 10,
  focal_points = capital_cities,
  width = 35,
  verbose = FALSE
)
save(aus_hexmap, 
     file = "data/aus_hexmap.rda")
}

load("data/aus_hexmap.rda")
## Prepare to plot
fort_hex <- fortify_hexagon(data = aus_hexmap,
                            sf_id = "sa2_name_2011",
                            hex_size = 0.2) |> 
            left_join(sa2thyroid_ERP |> select(sa2_name_2011, SIR, p50))
## Make a plot
aus_hexmap_plot <- ggplot() +
  geom_sf(data=sa2thyroid_ERP, fill=NA, colour="grey60", size=0.1) +
  geom_polygon(data = fort_hex, aes(x = long, y = lat, group = hex_id, fill = SIR)) +
  scale_fill_identity() +
  invthm 
aus_hexmap_plot  
```

:::
::::

## Plots and statistical inference {.transition-slide .center style="text-align: center;"}

## Tidy data and random variables 


::: columns
::: column

- Tidy data mirrors elementary statistics
- Tabular form puts variables in columns and observations in rows
- Not all tabular data is in this form
- In this form, we can think about  $X_1 \sim N(0,1), ~~X_2 \sim \text{Exp}(1) ...$
:::

::: column

$$\begin{align}X &= \left[ \begin{array}{rrrr}
           X_1 & X_2 & ... & X_p 
           \end{array} \right] \\
  &= \left[ \begin{array}{rrrr}
           X_{11} & X_{12} & ... & X_{1p} \\
           X_{21} & X_{22} & ... & X_{2p} \\
           \vdots & \vdots & \ddots& \vdots \\
           X_{n1} & X_{n2} & ... & X_{np}
           \end{array} \right]\end{align}$$
:::
:::

## Grammar of graphics and statistics 

::: {.incremental}

- A statistic is a function on the values of items in a sample, e.g. for $n$ iid random variates $\bar{X}_1=\displaystyle\sum_{i=1}^n X_{i1}$, $s_1^2=\displaystyle\frac{1}{n-1}\displaystyle\sum_{i=1}^n(X_{i1}-\bar{X}_1)^2$
- We study the behaviour of the statistic over all possible samples of size $n$. 
- The grammar of graphics is the mapping of (random) variables to graphical elements, making plots of data into statistics

:::

## Making inference from data

:::: {.columns}
::: {.column}

Traditional statistical inference

You need a:

- statistic computed from the data
- null and alternative hypothesis
- reference distribution on which to measure the statistic
    - if it is extreme on this scale, reject the null

:::
::: {.column}

Inference with *data plots*

You need a:

- plot description provided by the grammar [(a statistic)]{.blue .fragment}
    - This implies one or more null hypotheses
- null generating mechanism, e.g. permutation, simulation from a distribution or model
- visual evaluation: is one plot in the array different?

:::
::::

## Example using TB [(1/2)]{.smallest}

:::: {.columns}
::: {.column width=40%}
We made various statements about patterns. Let's take the trend plot comparing temporal trend of counts males and females in different age groups. We said:

- Counts for men are higher across all age groups

We will make comparison plots by simulating sex from a binomial distribution, for each age group and year.

:::
::: {.column width=60%}

```{r}
#| label: tb-lineup1
#| code-fold: true
#| fig-width: 8
#| fig-height: 6
#| out-width: 100%
# Need a binary variable to use simulation
tb_idn_long <- tb_idn |>
  filter(year > 2000, year < 2013) |>
  mutate(sex01 = ifelse(sex=="m", 0, 1)) |>
  select(year, age_group, sex01, count) |>
  uncount(count)

# Generate the null samples
set.seed(733)
pos = sample(1:4)
tb_in_s1 <- tb_idn_long |>
  group_by(year, age_group) |>
  mutate(sex01 = rbinom(n(), 1, 0.5)) |>
  ungroup() |>
  count(year, age_group, sex01) |>
  rename(count = n) |>
  mutate(.sample = pos[1])
tb_in_s2 <- tb_idn_long |>
  group_by(year, age_group) |>
  mutate(sex01 = rbinom(n(), 1, 0.5)) |>
  ungroup() |>
  count(year, age_group, sex01) |>
  rename(count = n) |>
  mutate(.sample = pos[2])
tb_in_s3 <- tb_idn_long |>
  group_by(year, age_group) |>
  mutate(sex01 = rbinom(n(), 1, 0.5)) |>
  ungroup() |>
  count(year, age_group, sex01) |>
  rename(count = n) |>
  mutate(.sample = pos[3])
tb_idn_d <- tb_idn_long |>
  count(year, age_group, sex01) |>
  rename(count = n) |>
  mutate(.sample = pos[4])

# Make the lineup  
tb_l <- bind_rows(tb_idn_d, tb_in_s1, tb_in_s2, tb_in_s3) 
  
ggplot(tb_l, aes(x=year, y=count, colour=factor(sex01))) +
  geom_point() +
  geom_smooth(se=F) + 
  facet_grid(.sample~age_group) +
  ylim(c(0,NA)) + 
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "none",
        strip.text.x = element_blank())
```

:::
::::

## Example using TB [(2/2)]{.smallest}

:::: {.columns}
::: {.column width=40%}
Now lets check the temporal trend of proportion of males and females

We will make comparison plots by randomising the proportion across years for each age group.

:::
::: {.column width=60%}

```{r}
#| label: tb-lineup2
#| code-fold: true
#| fig-width: 8
#| fig-height: 6
#| out-width: 100%
set.seed(854)
pos = sample(1:4)

tb_idn_prop <- tb_idn |> group_by(year, age_group) |> 
  summarise(p = count[sex=="m"]/sum(count)) 

tb_idn_prop_s1 <- tb_idn_prop |>
  group_by(age_group) |>
  mutate(p = sample(p)) |>
  mutate(.sample = pos[1])
tb_idn_prop_s2 <- tb_idn_prop |>
  group_by(age_group) |>
  mutate(p = sample(p)) |>
  mutate(.sample = pos[2])
tb_idn_prop_s3 <- tb_idn_prop |>
  group_by(age_group) |>
  mutate(p = sample(p)) |>
  mutate(.sample = pos[3])
tb_idn_prop <- tb_idn_prop |>
  mutate(.sample = pos[4])

tb_idn_prop_l <- bind_rows(tb_idn_prop, tb_idn_prop_s1, 
  tb_idn_prop_s2, tb_idn_prop_s3)

ggplot(tb_idn_prop_l, aes(x=year, y=p)) +
  geom_hline(yintercept = 0.50, colour="grey50", linewidth=2) +
  geom_point() +
  geom_smooth(se=F) + 
  facet_grid(.sample~age_group) +
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "none",
        strip.text.x = element_blank())
```

:::
::::

## Inference process

:::: {.columns}
::: {.column}

**Null-generating mechanisms**

- Permutation: randomizing the order of one of the variables breaks association, but keeps marginal distributions the same
- Simulation: from a given distribution, or model. Assumption is that the data comes from that model.

:::
::: {.column}

**Evaluation**

- Compute $p$-value, essentially following a binomial formula with multiple people reading the plots
- Power $=$ signal strength, based people picking the data plot with one plot design more than another

:::
::::



## Using the nullabor package

![](images/nullabor_hex.png){.absolute bottom=30 left=30 width="150px"}

:::: {.columns}
::: {.column width=40%}

```{r}
#| label: lineup 1
#| echo: true
# Make a lineup of mtcars data
# 20 plots, one data, 19 nulls
# Which one is different?
set.seed(913)
library(ggplot2)
l <- ggplot(
  lineup(
    null_permute('mpg'), 
    mtcars), 
  aes(mpg, wt)
) +
  geom_point() +
  facet_wrap(~ .sample) +
  theme(axis.text = element_blank(),
        axis.title = element_blank())
```

```{r}
# Compute the p-value for 10 observers
pvisual(5, 10, 20)
```

:::
::: {.column width=60%}

```{r}
#| echo: false
#| fig-height: 6
#| fig-width: 6
#| out-width: 80%
l
```
:::
::::

## Very important

- $H_0$ (underlying the null sample generation) is determined based on the plot type

- $H_0$ is **not** based on the structure seen in the data set

- Null data creation method does not match characteristics of original sample other than that in $H_0$



## Perceptual principles {.transition-slide .center style="text-align: center;"}

## Game: Which plot wears it better? {.smaller}

:::: columns
::: column

Coming up: **2 different plots** of 2012 TB incidence in Indonesia, based on variables `sex` and `age`.

- In arrangement A, separate plots are made for age, and sex is mapped to the x axis. 
- Conversely, in arrangement B, separate plots are made for sex, and age is mapped to the x axis. 


:::
::: column


> *At which age(s) are the counts for males and females relatively the same?*

Which plot makes this question easier to answer?

:::
::::

::: {layout-ncol=2}

```{r focus on one year gender side-by-side bars of males females, fig.height=3, echo=FALSE}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=sex, y=count, fill=sex)) +
    geom_bar(stat="identity", position="dodge") + 
    facet_wrap(~age, ncol=6) +
    ggtitle("Arrangement A")
```

```{r focus on one year age side-by-side bars of age group, fig.height=3, fig.width=8, echo=FALSE}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=age, y=count, fill=age)) +
    geom_bar(stat="identity", position="dodge") + 
    facet_wrap(~sex, ncol=6) +
    scale_fill_brewer("", palette="Dark2") +
    ggtitle("Arrangement B")
```
:::


## Three Variables {.smaller}

:::: columns
::: column
Next, we have **two different plots** of based on three variables: `age`, `sex`, `year`.

- In plot type A, a line plot of counts is drawn separately by age and sex, and year is mapped to the x axis. 
- Conversely, in plot type B, counts for sex, and age are stacked into a bar chart, separately by age and sex, and year is mapped to the x axis

:::

::: column

> *Is the trend for females increasing over time for all agre groups?* 

Which plot makes this easier? 

:::
::::


::: {layout-ncol=2}

```{r use a line plot instead of bar 2, echo = F}
#| fig-width: 10
#| fig-height: 3
ggplot(tb_idn, aes(x=year, y=count, colour=sex)) +
  geom_point() +
  geom_smooth(se=F, colour="black") + 
  facet_wrap(~age, ncol=6) +
  ggtitle("Type A") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10"))
```

```{r colour and axes fixes, echo = F}
#| fig-width: 10
#| fig-height: 3
ggplot(tb_idn, aes(x=year, y=count, fill=sex)) +
  geom_bar(stat="identity") + 
  facet_wrap(~age, ncol=6) +
  ggtitle("Type B") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10"))
```
:::


## Perceptual principles

- Hierarchy of mappings
- *Pre-attentive*: some elements are noticed before you even realise it.
- *Color palettes*: qualitative, sequential, diverging. Make sure the mapping is appropriate.
- *Proximity*: Place elements for primary comparison close together. 
- *Change blindness*: When focus is interrupted differences may not be noticed.
- *Aspect ratio*: the relative width to height makes some patterns easier or harder to read


## Hierarchy of mappings

:::: {.columns}
::: {.column width="60%"}

1. Position - common scale (BEST)
2. Position - nonaligned scale
3. Length, direction, angle
4. Area
5. Volume, curvature
6. Shading, color (WORST)

[(Cleveland, 1984; Heer and Bostock, 2009)]{.smaller}

:::

::: {.column width="40%"}

::: {.fragment}

1. scatterplot, barchart
2. side-by-side boxplot, stacked barchart
3. piechart, rose plot, gauge plot, donut, wind direction map, starplot
4. treemap, bubble chart, mosaicplot
5. chernoff face
6. choropleth map

:::

:::
::::

## Color palettes 

:::: {.columns}

::: {.column width="40%"}

- sequential 
- diverging 
- qualitative

[Color Brewer](http://colorbrewer2.org/#type=sequential&scheme=BuGn&n=3) annotates palettes with attributes.

:::
::: {.column width="60%"}

```{r show different types of color palettes, fig.height=7, fig.width=12, out.width="100%"}
display.brewer.all()
```

:::
::::


## Which type of palette is best?

:::: {.columns}
::: {.column width="33%"}

```{r palette1, fig.height=6, fig.width=6, echo = F}
ggplot(tb_idn, aes(x=year, y=count, colour=age)) +
  geom_line() + geom_point() +
  facet_wrap(~sex, ncol=6) +
  scale_colour_brewer("", palette="Dark2") +
  ggtitle("qualitative") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10")) +
  theme(legend.position = "bottom")
```

:::
::: {.column width=33%}

```{r palette2, fig.height=6, fig.width=6, echo = F}
ggplot(tb_idn, aes(x=year, y=count, colour=age)) +
  geom_line() + geom_point() +
  facet_wrap(~sex, ncol=6) +
  scale_color_manual(values = brewer.pal(8, "YlGnBu")[3:8]) +
  ggtitle("sequential") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10")) +
  theme(legend.position = "bottom")
```

:::
::: {.column width=33%}

```{r palette3, fig.height=6, fig.width=6, echo = F}
ggplot(tb_idn, aes(x=year, y=count, colour=age)) +
  geom_line() + geom_point() +
  facet_wrap(~sex, ncol=6) +
  scale_color_brewer("", palette = "PiYG") +
  ggtitle("diverging") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10")) +
  theme(legend.position = "bottom")
```

:::
::::

## Color blind-proofing

```{r using the dichromat package to check color blind appearance, echo=TRUE, eval = F}
clrs <- hue_pal()(9)
d + theme(legend.position = "none")

clrs <- dichromat(hue_pal()(9))
d + 
  scale_colour_manual("", values=clrs) + 
  theme(legend.position = "none")
```

- Online checking tool [coblis](https://www.color-blindness.com/coblis-color-blindness-simulator/): upload an image and it will re-map the colors for different colour perception issues. 
- The package `colorblind` has color blind friendly palettes    [(Susan: but the colours are awful `r emo::ji("sob")`).]{.smaller}
- The package `colorspace` also has good tools for proofing for colorblind proof palettes

[[Color deficiency friendly color schemes](https://davidmathlogic.com/colorblind/#%23005AB5-%23DC3220)]{.bottom}


## Color blind simulation

:::: {.columns}
::: {.column}

Original colours

```{r show the default colour scheme, echo=FALSE, fig.width=4, fig.height=4, out.width="80%"}
clrs <- hue_pal()(9)
dsamp <- diamonds |>
  sample_n(1000)
d <- ggplot(
  dsamp, aes(carat, price)) +
  geom_point(aes(
    colour = clarity))
p1 <- d +  
  scale_color_discrete() + 
  theme(legend.position = "none")  
p1
```
:::
::: {.column}

Color blind view

```{r show the dichromat adjusted colors, echo=FALSE, fig.width=4, fig.height=4, out.width="80%"}
clrs <- dichromat(hue_pal()(9))
p2 <- d + 
  scale_color_manual("", values=clrs) + 
  theme(legend.position = "none")

p2
```

:::
::::


## Pre-attentive

Can you find the odd one out?

```{r is shape preattentive, echo=FALSE, fig.width=4, fig.height=4, out.width="80%"}
set.seed(1106)
df <- data.frame(x=runif(100), y=runif(100), cl=sample(c(rep("A", 1), rep("B", 99))))
ggplot(data=df, aes(x, y, shape=cl)) + theme_bw() + 
  geom_point(size=3) +
  theme(legend.position="None", aspect.ratio=1, axis.text = element_blank(), axis.ticks=element_blank(), axis.title = element_blank()) 
```

## Pre-attentive

Is it easier now?

```{r is color preattentive, echo=FALSE, fig.width=4, fig.height=4, out.width="80%"}
ggplot(data=df, aes(x, y, colour=cl)) + 
  geom_point(size=3) +
  theme_bw() + 
  scale_colour_brewer(palette="Set1") +
  theme(legend.position="None", aspect.ratio=1, axis.text = element_blank(), axis.ticks=element_blank(), axis.title = element_blank()) 
```


## Proximity

Place elements that you want to compare close to each other. 
If there are multiple comparisons to make, you need to decide which one is most important.

::: {layout-ncol=2}

```{r a line plot on sex, fig.height=3, fig.width = 8, echo = F}
ggplot(tb_idn, aes(x=year, y=count, colour=sex)) +
  geom_line() + geom_point() +
  facet_wrap(~age, ncol=6) +
  ggtitle("Arrangement A") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10"))
```

```{r a line plot on age, fig.height=3, fig.width=8, echo = F}
ggplot(tb_idn, aes(x=year, y=count, colour=age)) +
  geom_line() + geom_point() +
  facet_wrap(~sex, ncol=6) +
  scale_colour_brewer("", palette="Dark2") +
  ggtitle("Arrangement B") +
  scale_x_continuous("year", breaks = seq(2000, 2012, 5), 
                     labels = c("00", "05", "10"))
```
:::


## Mapping and proximity

:::: {.columns}
::: {.column width="20%"}
Same proximity is used, but different geoms. 

*Which is better to compare the relative ratios of males to females by age?*
:::
::: {.column width="80%"}
```{r side-by-side bars of males females, fig.width=8, fig.height=3, echo = F, out.width="70%"}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=sex, y=count, fill=sex)) +
    geom_bar(stat="identity", position="dodge") + 
    facet_wrap(~age, ncol=6) +
    scale_y_continuous("count ('000)", breaks=seq(0, 28000, 5000),
      labels=c("0", "5", "10", "15", "20", "25")) +
    ggtitle("Position - common scale ")
```

```{r piecharts of males females, fig.width=10, fig.height=2.5, echo = F, out.width="100%"}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=1, y=count, fill=sex)) +
    geom_bar(stat="identity", position="fill") + 
    facet_wrap(~age, ncol=6) +
    ggtitle("Angle") + xlab("") + ylab("") +
    coord_polar(theta = "y")
```
:::
::::

## Mapping and proximity

:::: {.columns}
::: {.column width="20%"}
Same proximity is used, but different geoms. 

*Which is better to compare the relative ratios of ages by sex?*
:::

::: {.column width="40%"}
```{r side-by-side bars of age, fig.width=8, fig.height=4, echo = F, out.width="100%"}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=age, y=count, fill=age)) +
    geom_bar(stat="identity", position="dodge") + 
    facet_wrap(~sex, ncol=6) +
    scale_fill_brewer("", palette="Dark2") +
    scale_y_continuous("count ('000)", breaks=seq(0, 28000, 5000),
      labels=c("0", "5", "10", "15", "20", "25")) +
    ggtitle("Position - common scale ")
```

```{r piecharts of age, fig.width=8, fig.height=4, echo = F, out.width="100%"}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=1, y=count, fill=age)) +
    geom_bar(stat="identity", position="fill") + 
    facet_wrap(~sex, ncol=6) +
    scale_fill_brewer("", palette="Dark2") +
    ggtitle("Angle") + xlab("") + ylab("") +
    coord_polar(theta = "y")
```

:::

::: {.column width="40%"}

```{r stacked bars of age, fig.width=8, fig.height=6, echo = F, out.width="100%"}
tb_idn |> filter(year == 2012) |>
  ggplot(aes(x=1, y=count, fill=age)) +
    geom_bar(stat="identity", position="fill") + 
    facet_wrap(~sex, ncol=6) +
    scale_fill_brewer("", palette="Dark2") +
    ggtitle("Position - nonaligned") + xlab("") + ylab("")
```

:::
::::

## Change blindness


```{r facetting plots can result in change blindness, echo=TRUE, out.width="50%", fig.width=8, fig.height=5}
#| code-fold: true
ggplot(dsamp, aes(x=carat, y=price, colour = clarity)) +
  geom_point() +
  geom_smooth(se=FALSE) +
  scale_color_brewer(palette="Set1") +
  scale_y_continuous("price ('000)", breaks=seq(0, 18000, 5000),
    labels=c("0", "5", "10", "15")) +
  facet_wrap(~clarity, ncol=4)
```

*Which has the steeper slope, VS1 or VS2?*

## Change blindness

:::: columns
::: column
Making comparisons across plots requires the eye to jump from one focal point to another. 

It may result in missing important differences. 
:::

::: column

```{r averlaying makes comparisons easier, echo=TRUE, out.width="100%", fig.width=7, fig.height=4}
#| code-fold: true
ggplot(dsamp, aes(x=carat, y=price, 
                  colour = clarity)) +
  geom_point() +
  geom_smooth(se=FALSE) +
  scale_y_continuous("price ('000)", breaks=seq(0, 18000, 5000),
    labels=c("0", "5", "10", "15")) +
  scale_color_brewer(palette="Set1") 
```
:::
::::


## Core principles [1/2]{.smallest} {.information-slide .center style="text-align: center;"}

- Make a plot of your **data**! 
    - The hierarchy matters if the structure is weak or differences b/w groups are small. 
- Knowing how to use proximity is a valuable and rare skill
- Use of colour: don't over use 
    - Too many colours
    - Mapping cts variable to colour to add another dimension

## Core principles [2/2]{.smallest} {.information-slide .center style="text-align: center;"}

- Show the data! 
    - Statistics are good if there's too much data
    - Always plot the data for yourself to see the variability
- One plot is never enough
    - Plot the data in different ways
    - Understand the relationships between variables


## Resources

- [R for Data Science (2e)](https://r4ds.hadley.nz)
- [ggplot2: Elegant Graphics for Data Analysis, Hadley Wickham](https://ggplot2-book.org)
- [ggplot2 web site](https://ggplot2.tidyverse.org)
- [R Graphics Cookbook, Winston Chang](http://www.cookbook-r.com/Graphs/)
- [Data Visualization, Kieran Healy](https://socviz.co)
- [Data Visualization with R, Rob Kabacoff](https://rkabacoff.github.io/datavis/index.html)
- [Fundamentals of Data Visualization, Claus O. Wilke](https://serialmentor.com/dataviz/)
- [Naomi Robbins, Creating More Effective Graphs](http://www.nbr-graphs.com)
- [Cleveland, McGill (1984) Graphical perception: Theory, experimentation](https://www.tandfonline.com/doi/abs/10.1080/01621459.1984.10478080)
- [Heer, Bostock (2010) Crowdsourcing graphical perception](http://vis.stanford.edu/files/2010-MTurk-CHI.pdf)
- [Antony Unwin, Graphical Data Analysis with R](https://www.crcpress.com/Graphical-Data-Analysis-with-R/Unwin/9781498715232)
- Wickham, H., Cook, D., Hofmann, H. and Buja, A. (2010) Graphical Inference for Infovis,  [http://doi.ieeecomputersociety.org/10.1109/TVCG.2010.161](http://doi.ieeecomputersociety.org/10.1109/TVCG.2010.161). 
- Hofmann, H., Follett, L., Majumder, M. and Cook, D. (2012) Graphical Tests for Power Comparison of Competing Designs, [http://doi.ieeecomputersociety.org/10.1109/TVCG.2012.230](http://doi.ieeecomputersociety.org/10.1109/TVCG.2012.230).


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