Tidy Approaches to Missing Data

class: center, middle, inverse, title-slide

# Tidy Approaches to Missing Data
## WOMBAT Conference
### Nicholas Tierney, Monash University
### Friday, 2nd June, 2017

---

class: center, middle
# Question

### Who has dealt with missing data

???

- Question: introduce missing data by polling the room.

Then briefly say something along the lines of:

> Missing data is important to identify because it often creates bias.

> It is prevalent is so many pieces of data analysis, for example, in survey design, individuals might not wish to provide answers to sensitive questions regarding drug and alcohol consumption.

---
class: inverse, center, middle
# Motivation

### Why would I want to deal with missing data?

---
class: center, middle
# Example

### Australian weather data 
### || Date | Wind | Temp | Rain ||
### Using the R package: [`bomrang`](https://github.com/ToowoombaTrio/bomrang)

---

## Weather data from Melbourne Olympic Park

<div id="htmlwidget-ec1c9ec6386bef060ab4" style="width:100%;height:80%;" class="leaflet html-widget"></div>
<script type="application/json" data-for="htmlwidget-ec1c9ec6386bef060ab4">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["//{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"maxNativeZoom":null,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"continuousWorld":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":null,"unloadInvisibleTiles":null,"updateWhenIdle":null,"detectRetina":false,"reuseTiles":false,"attribution":"© <a href=\"http://openstreetmap.org\">OpenStreetMap<\/a> contributors, <a href=\"http://creativecommons.org/licenses/by-sa/2.0/\">CC-BY-SA<\/a>"}]},{"method":"addMarkers","args":[-37.82463,144.9813,null,null,null,{"clickable":true,"draggable":false,"keyboard":true,"title":"","alt":"","zIndexOffset":0,"opacity":1,"riseOnHover":false,"riseOffset":250},null,null,null,null,null,null,null]},{"method":"addProviderTiles","args":["Hydda.RoadsAndLabels",null,null,{"errorTileUrl":"","noWrap":false,"zIndex":null,"unloadInvisibleTiles":null,"updateWhenIdle":null,"detectRetina":false,"reuseTiles":false}]}],"limits":{"lat":[-37.82463,-37.82463],"lng":[144.9813,144.9813]}},"evals":[],"jsHooks":[]}</script>

---
# Identifying Missingness Mechanisms

???

> Now, A brief example of this might be if you have some data from a weather station might have some malfunctions when there are extreme weather events, it might not record data, as the sensor could shut off.

Just show

And this is a pretty well defined situation here - the answer in fact seems obvious:

> The weather station turns off when there is extremely low and extremely high temperatures.

But, like all good explanations, this one is simple - but the process to get there, where we descri probably was not.

In order to get to a position where you could generate this graphic, and show this data, you probably had to spend more time than you would have liked developing exploratory data analyses and models

It might seem obvious to point out the precise mechanism for generating the missingness, but this is kind of difficult when there is a lot of missing data, or when there are many variables.

---
# Identifying Missingness Mechanisms

---
class: inverse, center, middle
# Wait, What?

---
background-image: url(https://i.redd.it/wku7d6zhk3ry.jpg)
background-size: contain 
background-position: 50% 50%
class: center, bottom

---

---

.center[
# The Questions
]

.large[
1. How to START looking at missing data? 
2. How to explore missingness mechanisms? 
3. Can we model missingness?
]

.center[
# The Answer(s)?
]

.center[
.large[
Tools that facilitate exploration of missing values
]
]
???

To look at these questions I am going to discuss some tools that  I have built with Di Cook and Miles McBain, to facilitate the highlighting and exploration of missing values.

I'll then move discuss the way that we can handle and manipulate missing data.

I will then move to discuss the problems with current methods for exploring and visualising missings, and what we can do about it.

---
class: inverse, middle, center
# Back to the beginning

```
## # A tibble: 155 x 5
##                 ymd_hm     time gust_kmh air_temp rain_trace
##                  <chr>    <chr>   <fctr>    <chr>      <chr>
##  1 2017-06-01 14:30:00 14:30:00       19     13.1          0
##  2 2017-06-01 14:00:00 14:00:00       19     13.1          0
##  3 2017-06-01 13:30:00 13:30:00       13     13.3          0
##  4 2017-06-01 13:00:00 13:00:00       15     13.4          0
##  5 2017-06-01 12:30:00 12:30:00       13     13.2          0
##  6 2017-06-01 12:00:00 12:00:00       15     12.5          0
##  7 2017-06-01 11:30:00 11:30:00       20     12.2          0
##  8 2017-06-01 11:00:00 11:00:00        0     11.8          0
##  9 2017-06-01 10:30:00 10:30:00        7     11.2          0
## 10 2017-06-01 10:00:00 10:00:00       13     10.9          0
## # ... with 145 more rows
```

---
# Data Exploration: First Steps

```r
summary(melb_data)
```

```
##     ymd_hm              time              gust_kmh    air_temp        
##  Length:155         Length:155         13     :26   Length:155        
##  Class :character   Class :character   15     :16   Class :character  
##  Mode  :character   Mode  :character   11     :15   Mode  :character  
##                                        17     :13                     
##                                        19     :12                     
##                                        20     :12                     
##                                        (Other):61                     
##   rain_trace       
##  Length:155        
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 
```

???

Before you start looking at missing data, you'll need to look at the data.

Some useful tools to start getting a look at your data are things like:

---
# Data Exploration: First Steps

```r
str(melb_data)
```

```
## Classes 'tbl_df', 'tbl' and 'data.frame':	155 obs. of  5 variables:
##  $ ymd_hm    : chr  "2017-06-01 14:30:00" "2017-06-01 14:00:00" "2017-06-01 13:30:00" "2017-06-01 13:00:00" ...
##  $ time      : chr  "14:30:00" "14:00:00" "13:30:00" "13:00:00" ...
##  $ gust_kmh  : Factor w/ 17 levels "0","7","9","11",..: 8 8 5 6 5 6 9 1 2 5 ...
##  $ air_temp  : chr  "13.1" "13.1" "13.3" "13.4" ...
##  $ rain_trace: chr  "0" "0" "0" "0" ...
```

---
# Data Exploration: First Steps

```r
glimpse(melb_data)
```

```
## Observations: 155
## Variables: 5
## $ ymd_hm     <chr> "2017-06-01 14:30:00", "2017-06-01 14:00:00", "2017...
## $ time       <chr> "14:30:00", "14:00:00", "13:30:00", "13:00:00", "12...
## $ gust_kmh   <fctr> 19, 19, 13, 15, 13, 15, 20, 0, 7, 13, 11, 9, 13, 1...
## $ air_temp   <chr> "13.1", "13.1", "13.3", "13.4", "13.2", "12.5", "12...
## $ rain_trace <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "...
```

---
# Data Exploration: First Steps

- `skimr` (developed at rOpenSci unconf)

```r
skimr::skim(melb_data)
```

```
## 
## Character Variables
## # A tibble: 4 x 9
##          var      type complete missing empty     n   min   max n_unique
## *      <chr>     <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>
## 1   air_temp character      143      12     0   155     1     4       62
## 2 rain_trace character      155       0     0   155     1     3       18
## 3       time character      155       0     0   155     8     8       59
## 4     ymd_hm character      155       0     0   155    19    19      155
## 
## Factor Variables
## # A tibble: 1 x 7
##        var   type complete missing     n n_unique
##      <chr>  <chr>    <dbl>   <dbl> <dbl>    <dbl>
## 1 gust_kmh factor      155       0   155       17
## # ... with 1 more variables: stat <chr>
```

???
Developed at the rOpenSci unconf

---

# Data Exploration: Pre-exploratory visualisation?

- Printed summary = Manual parsing of text

- Lots of Squinting

???

However, these summary methods still require you to parse a lot of text, and read through summaries. These are very useful, but they fall over when you find, for example, that all of your data might be encoded as characters and factors, instead of dates, integers, and numerics.

---
# Data Exploration: Pre-exploratory visualisation?

---
class: inverse, middle, center
# `visdat`

Visualise whole data frames at once

???
can use two functions from the `visdat` package:

- `vis_dat`, and
- `vis_miss`

---
.center[
## `vis_dat`
]

```r
vis_dat(melb_data)
```

???

These visualise the whole dataframe at once, and provide some useful information about whether the data is missing or not.

In particular, `vis_miss` provides a summary of whether the data is all present, and also provides small summaries in each column as to the amount of missings in each of the columns.

---
class: middle, center
# Cleaning...

---
class: center
## `vis_dat(melb_data)`

---
class: center
### `vis_miss(melb_data)`

---
# Exploring missingness relationships

.large[

We can identify key variables that are missing

For further exploration, we need to explore the relationship amongst:

]

- Air temperature,
- Gust (Km/hr), and 
- Rainfall.

???

This process of visualising the data as a heatmap allows you identify what is missing and where. For here though, are interested in exploring these missingness relationships.

---

## Exploring missingness relationships

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_point()
```

```
## Warning: Removed 12 rows containing missing values (geom_point).
```

???

Typically, when exploring this data, you would do something like this:

The problem with this is that ggplot does not handle missings be default, and removes the missing values. This makes it hard to explore the missing values.

---

## Exploring missingness relationships

How do you visualise something that is not there?

One approach, (MANET, ggobi) is to shift missing values below the minimum value, so that they can be plotted on the same axis.

.pull-left[

```
## # A tibble: 11 x 1
##    air_temp
##       <dbl>
##  1      9.3
##  2      9.7
##  3      9.0
##  4      8.7
##  5      8.9
##  6      9.1
##  7      9.1
##  8       NA
##  9       NA
## 10       NA
## 11     11.9
```
]

.pull-right[

```
## # A tibble: 11 x 2
##    air_temp air_temp_shift
##       <dbl>          <dbl>
##  1      9.3       9.300000
##  2      9.7       9.700000
##  3      9.0       9.000000
##  4      8.7       8.700000
##  5      8.9       8.900000
##  6      9.1       9.100000
##  7      9.1       9.100000
##  8       NA       8.328230
##  9       NA       8.431627
## 10       NA       8.388872
## 11     11.9      11.900000
```
]

---

# Exploring missingness relationships

However, doing this for each variable would typically involve some code like this:

```r
melb_data %>%  
mutate(miss_cat = miss_cat(., "rain_trace", "air_temp")) %>%   
ggplot(data = .,
       aes(x = shadow_shift(rain_trace),
           y = shadow_shift(air_temp),
           colour = miss_cat)) +
   geom_point()
```

Which is rather involved, and so not ideal

-- 
.large[
.center[
So where do we go?
]
]
---
background-image: url(https://68.media.tumblr.com/c9eaec31221bcb2c6b7c7b3d9e464954/tumblr_nbgwierwKV1rcb0d2o6_400.jpg)

background-size: contain 
background-position: 50% 50%
class: center, bottom, inverse

---
class: inverse middle, center

# Narnia

### An R package to simplify dealing with missing data

---
.center[
# What does it do?
]

--
.pull-left[

## Visualisations

- `geom_missing_point()`
- `gg_missing_var()`
- `bind_shadow()` patterns

]

.pull-right[

## Numerical Summaries

- `n_miss()`
- `n_complete()`
- `miss_case_pct()`
- `miss_case_table()`
- `miss_case_summary()`
- `miss_var_pct()`
- `miss_var_summary()`
- `miss_var_table()`

]

---
.center[
## `geom_missing_point()`
]

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_point()
```

```
## Warning: Removed 12 rows containing missing values (geom_point).
```

---
.center[
## `geom_missing_point()`
]

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_missing_point()
```

---
.center[
## `geom_missing_point()`
]

```r
ggplot(messy_melb_data, 
       aes(x = gust_kmh, 
           y = air_temp)) + 
  geom_missing_point()
```

???

Instead, we have create a ggplot2 geom, `geom_missing_point()`.

This `geom` allows for missing values to be displayed, and also works with the rest of ggplot2 - themes, and facets as well.

---
## Facets!

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_missing_point() + 
  facet_wrap(~lubridate::month(messy_melb_data$ymd_hm))
```

---
## Themes!

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_missing_point() + 
  facet_wrap(~lubridate::month(messy_melb_data$ymd_hm)) +
  theme_dark()
```

---
class: inverse, middle, center
# The shadowlands

### Representing Missing values

---
## The shadowlands

`as_shadow(melb_data)`

```
## # A tibble: 21 x 5
##    ymd_hm_NA time_NA gust_kmh_NA air_temp_NA rain_trace_NA
##       <fctr>  <fctr>      <fctr>      <fctr>        <fctr>
##  1       !NA     !NA         !NA         !NA           !NA
##  2       !NA     !NA         !NA         !NA           !NA
##  3       !NA     !NA         !NA         !NA           !NA
##  4       !NA     !NA         !NA         !NA           !NA
##  5       !NA     !NA         !NA         !NA           !NA
##  6       !NA     !NA         !NA         !NA           !NA
##  7       !NA     !NA         !NA         !NA           !NA
##  8       !NA     !NA         !NA          NA           !NA
##  9       !NA     !NA         !NA          NA           !NA
## 10       !NA     !NA         !NA          NA           !NA
## # ... with 11 more rows
```

---

## The shadowlands

```r
melb_data_shadow <- bind_shadow(melb_data)

glimpse(melb_data_shadow)
```

```
## Observations: 155
## Variables: 10
## $ ymd_hm        <dttm> 2017-06-01 14:30:00, 2017-06-01 14:00:00, 2017-...
## $ time          <time> 14:30:00, 14:00:00, 13:30:00, 13:00:00, 12:30:0...
## $ gust_kmh      <int> 19, 19, 13, 15, 13, 15, 20, 0, 7, 13, 11, 9, 13,...
## $ air_temp      <dbl> 13.1, 13.1, 13.3, 13.4, 13.2, 12.5, 12.2, 11.8, ...
## $ rain_trace    <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0...
## $ ymd_hm_NA     <fctr> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N...
## $ time_NA       <fctr> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N...
## $ gust_kmh_NA   <fctr> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N...
## $ air_temp_NA   <fctr> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N...
## $ rain_trace_NA <fctr> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N...
```

---

## The shadowlands

```r
ggplot(melb_data_shadow,
       aes(x = rain_trace,
           colour = air_temp_NA)) + 
  geom_density()
```

---

## The shadowlands

```r
ggplot(melb_data_shadow,
       aes(x = gust_kmh,
           colour = air_temp_NA)) + 
  geom_density()
```

???

- Provides a consistent way to set the existing location of missing values.

- This has additional great benefits when combined with imputation

- Also allows for additional visualisations

---

# What about imputation?

```r
library(simputation)
  simputation::impute_lm(air_temp ~ rain_trace + gust_kmh, 
                         dat = melb_data) %>%
  ggplot(aes(x = ymd_hm,
             y = air_temp)) + 
  geom_point()
```

---
class: inverse, middle, center

## Yeah, what about it?

---

## What is imputed, stays imputed

**it's invisible!!**

---

# Shadows make things clearer!

```r
  simputation::impute_lm(air_temp ~ rain_trace + gust_kmh,
                         dat = melb_data_shadow) %>%
  ggplot(aes(x = ymd_hm,
             y = air_temp,
             colour = air_temp_NA)) + 
  geom_point()
```

---
## Narnia Numerical Summaries

```r
n_miss(melb_data)
```

```
## [1] 12
```

```r
n_complete(melb_data)
```

```
## [1] 763
```

```r
miss_case_pct(melb_data)
```

```
## [1] 7.741935
```

```r
miss_var_pct(melb_data)
```

```
## [1] 20
```

---

## Summaries

**for each case**

```r
miss_case_summary(melb_data)
```

```
## # A tibble: 155 x 3
##     case n_missing percent
##    <int>     <int>   <dbl>
##  1     1         0       0
##  2     2         0       0
##  3     3         0       0
##  4     4         0       0
##  5     5         0       0
##  6     6         0       0
##  7     7         0       0
##  8     8         0       0
##  9     9         0       0
## 10    10         0       0
## # ... with 145 more rows
```

---

## Summaries

**for each variable**

```r
miss_var_summary(melb_data)
```

```
## # A tibble: 5 x 3
##     variable n_missing  percent
##        <chr>     <int>    <dbl>
## 1   air_temp        12 7.741935
## 2     ymd_hm         0 0.000000
## 3       time         0 0.000000
## 4   gust_kmh         0 0.000000
## 5 rain_trace         0 0.000000
```

---

## Tabulations

**for each case**

```r
miss_case_table(melb_data)
```

```
## # A tibble: 2 x 3
##   n_missing_in_case n_cases   percent
##               <int>   <int>     <dbl>
## 1                 0     143 92.258065
## 2                 1      12  7.741935
```

---

## Tabulations

**for each variable**

```r
miss_var_table(melb_data)
```

```
## # A tibble: 2 x 3
##   n_missing_in_var n_vars percent
##              <int>  <int>   <dbl>
## 1                0      4      80
## 2               12      1      20
```

---
class: inverse, middle, center
# Back to our questions

---
# 1. How to START looking at missing data?

```r
vis_dat(melb_data)
```

---
# 1. How to START looking at missing data?

```r
vis_miss(melb_data)
```

---
# 2. How to explore missingness mechanisms?

```r
ggplot(messy_melb_data, 
       aes(x = rain_trace, 
           y = air_temp)) + 
  geom_missing_point()
```

---
# 2. How to explore missingness mechanisms?

```r
ggplot(melb_data_shadow,
       aes(x = gust_kmh,
           colour = air_temp_NA)) + 
  geom_density()
```

---
# 3. Can we model missingness?

You can model missing data using methods from Tierney et el. (2015).

Predict the proportion of missingness - use `add_prop_miss()` from `narnia`.

```r
melb_data %>%
  add_prop_miss() %>%
  slice(99:104)
```

```
## # A tibble: 6 x 6
##                ymd_hm     time gust_kmh air_temp rain_trace prop_miss
##                <dttm>   <time>    <int>    <dbl>      <dbl>     <dbl>
## 1 2017-05-30 18:11:00 18:11:00       17      8.9        1.2       0.0
## 2 2017-05-30 18:00:00 18:00:00       15      9.1        0.0       0.0
## 3 2017-05-30 17:30:00 17:30:00       13      9.1        0.0       0.0
## 4 2017-05-30 17:00:00 17:00:00       32       NA        0.0       0.2
## 5 2017-05-30 16:30:00 16:30:00       33       NA        0.0       0.2
## 6 2017-05-30 16:00:00 16:00:00       32       NA        0.0       0.2
```

---
# 3. Can we model missingness?

```r
library(rpart)
library(rpart.plot)

messy_melb_data %>%
  add_prop_miss() %>%
  rpart(prop_miss ~ ., data = .) %>%
  prp(type = 4, extra = 101, prefix = "Prop. Miss = ")
```

```r
# library(visdat)
# vis_miss(messy_melb_data)
```

---

# So what?

.large[
- These tools help us identify missingness, while maintaining a tidy workflow

Why do we care about these mechanisms or these patterns?

- It helps us understand potential mechanisms, such as:
	- Equipment failures
		- Future work and maintenance
	- Possible solutions?
	  - Impute based upon nearby station records
	  
]

---

## Missing Data: Future development

.large[
- Continue making `narnia` play well with `dplyr` and `tidyr`
- Big Data tools (sparklyr, sparklingwater)
- Expand ggplot `geom_missing_*` family
- Develop methods for handling and visualising multiple imputation
]

---
# Thank you

.large[
- Di Cook

- Rob Hyndman

- Miles McBain
]

---
# References

.pull-left[
## R packages
- github.com/njtierney/visdat
- github.com/njtierney/narnia
- github.com/ToowoombaTrio/BOMRang

## My details

- twitter.com/nj_tierney
- github.com/njtierney
- njtierney.com
]

.pull-right[

## Software / Papers
- MANET: http://www.rosuda.org/MANET/
- ggobi: http://www.ggobi.org/
- Tierney NJ et al, BMJ Open 2015, doi: 10.1136/bmjopen-2014-007450
]