```{r, include = FALSE}
# Load libraries
source("libraries.R")
# Set up chunk for all slides
source("chunk_options_and_themes.R")
```
```{r}
#| label: setting_up
#| echo: false
load("data/french_fries.rda")
```
## Outline
- Framing problems
- Tidy models formulation
- Case study: hotel bookings
- Unsupervised learning
- Temporal data and forecasting
## Framing the problem
1. **Supervised learning**
a. **classification**: categorical $y_i$ is [available]{.monash-orange2} for all $x_i$
b. **regression**: continuous $y_i$ is [available]{.monash-orange2} for all $x_i$
2. **Unsupervised learning**: $y_i$ [unavailable]{.monash-orange2} for all $x_i$
```{r fig.width=9, fig.height=4, out.width="100%", fig.align='center', echo=FALSE}
#| label: overview-methods
load("data/flea.rda")
p1 <- ggplot(flea, aes(x=tars1, y=aede1, colour = species)) +
geom_point() +
scale_colour_brewer(palette = "Dark2") +
xlab("Var 1") + ylab("Var 2") +
ggtitle("Classification") +
theme(legend.position="None") +
theme(aspect.ratio=1)
p2 <- ggplot(mtcars, aes(x=hp, y=mpg)) +
geom_point() +
geom_smooth(se=F) +
ggtitle("Regression") +
theme(aspect.ratio=0.7)
p3 <- ggplot(flea, aes(x=tars1, y=aede1)) +
geom_point() + xlab("Var 1") + ylab("Var 2") +
ggtitle("Clustering") +
theme(aspect.ratio=1)
p1 + p2 + p3 + plot_layout(ncol=3)
```
## What type of problem is this? [(1/3)]{.smallest}
Food servers' tips in restaurants may be influenced by many factors, including the nature of the restaurant, size of the party, and table locations in the restaurant. Restaurant managers need to know which factors matter when they assign tables to food servers. For the sake of staff morale, they usually want to avoid either the substance or the appearance of unfair treatment of the servers, for whom tips (at least in restaurants in the United States) are a major component of pay.
In one restaurant, a food server recorded the following data on all
customers they served during an interval of two and a half months in
early 1990. The restaurant, located in a suburban shopping mall, was
part of a national chain and served a varied menu. In observance of
local law the restaurant offered seating in a non-smoking section to
patrons who requested it. Each record includes a day and time, and
taken together, they show the server's work schedule.
What is $y$? What is $x$?
## What type of problem is this? [(2/3)]{.smallest}
Every person monitored their email for a week and recorded information about each email message; for example, whether it was spam, and what day of the week and time of day the email arrived. We want to use this information to build a spam filter, a classifier that will catch spam with high probability but will never classify good email as spam.
What is $y$? What is $x$?
## What type of problem is this? [(3/3)]{.smallest}
A health insurance company collected the following information about households:
- Total number of doctor visits per year
- Total household size
- Total number of hospital visits per year
- Average age of household members
- Total number of gym memberships
- Use of physiotherapy and chiropractic services
- Total number of optometrist visits
The health insurance company wants to provide a small range of products, containing different bundles of services and for different levels of cover, to market to customers.
What is $y$? What is $x$?
## Model form
All (data-centric) models have a fitted values and residuals.
$$y = f(x_1, x_2, ..., x_p) + \varepsilon$$
where $y$ is the observed response, $x_1, ..., x_p$ are the observed values of $p$ predictors and $\varepsilon$ is the error. We conventionally use $n$ to specfify the sample size.
- We might not be able to exactly specify $f$ in some forms
- The residuals have some ideal properties, such as uniform over data space, symmetric around the fitted value, and some models impose a distribution.
## Accuracy vs interpretability
::: {.columns}
::: {.column}
[Predictive accuracy]{.monash-blue2}
The primary purpose is to be able to [predict]{.monash-orange2} $\widehat{Y}$ for new data. And we'd like to do that well! That is, [accurately]{.monash-orange2}.
{style="font-size: 50%;"}
:::
::: {.column}
[Interpretability]{.monash-blue2}
Almost equally important is that we want to [understand the relationship]{.monash-orange2} between ${\mathbf X}$ and $Y$. The simpler model that is (almost) as accurate is the one we choose, always.
> Person: *Why did you predict 42 for this value?*
> Computer: **Awkward silence**
[From Interpretable Machine Learning](https://christophm.github.io/interpretable-ml-book/){ style="font-size: 50%;"}
:::
:::
## Parametric vs non-parametric
::: {.columns}
::: {.column}
[Parametric]{.monash-purple2} methods
- [Assume]{.monash-purple2} that the model takes a specific form
- Fitting then is a matter of [estimating the parameters]{.monash-purple2} of the model
- Generally considered to be [less flexible]{.monash-purple2}
- If assumptions are [wrong]{.monash-purple2}, performance likely to be [poor]{.monash-purple2}
:::
::: {.column}
[Non-parametric]{.monash-pink2} methods
- [No]{.monash-pink2} specific assumptions
- Allow the [data to specify]{.monash-pink2} the model form, without being too rough or wiggly
- More [flexible]{.monash-pink2}
- Generally needs [more observations]{.monash-pink2}, and not [too many variables]{.monash-pink2}
- Easier to [over-fit]{.monash-pink2}
:::
:::
##
::: {.columns}
::: {.column width=30%}
{width=500 style="font-size: 50%;"}
:::
::: {.column width=33%}
```{r}
#| echo: false
#| eval: false
#| label: sine-curve-data
# Generate the sine-curve data
set.seed(1259)
x1 <- runif(340)
x2 <- runif(340)
y <- 3*x1+sin(x1*15)
y <- (y-min(y))/(max(y)-min(y))
d <- tibble(x1, x2, y)
d$cl <- ifelse(x2 > y, "A", "B")
d$cl[sample(1:340, 25)] <- sample(c("A", "B"), 25, replace=TRUE)
write_csv(d, file="data/sine-curve.csv")
# Test set
x1 <- runif(212)
x2 <- runif(212)
y <- 3*x1+sin(x1*15)
y <- (y-min(y))/(max(y)-min(y))
d <- tibble(x1, x2, y)
d$cl <- ifelse(x2 > y, "A", "B")
d$cl[sample(1:212, 18)] <- sample(c("A", "B"), 18, replace=TRUE)
write_csv(d, file="data/sine-curve-test.csv")
```
```{r}
#| label: sine-curve-plot
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 100%
w <- read_csv("data/sine-curve.csv") |>
mutate(cl = factor(cl))
ggplot(w, aes(x=x1, y=x2, colour = cl)) +
geom_point(size=2.5, alpha=0.8) +
geom_line(aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Observed data") +
theme(legend.position = "none")
```
[Black line [is true boundary]{.monash-blue2}.]{.smaller}
[[Grids]{.monash-pink2} (right) show [boundaries]{.monash-pink2} for two different models.]{.smaller}
:::
::: {.column width=30%}
```{r}
#| label: sin-curve-models
#| echo: false
#| fig-width: 4
#| fig-height: 8
library(geozoo)
w_grid <- cube.solid.random(p=2)$points
colnames(w_grid) <- c("x1", "x2")
w_grid <- w_grid |> as_tibble() |>
mutate(x1 = x1*(max(w$x1)-min(w$x1)+min(w$x1)),
x2 = x2*(max(w$x2)-min(w$x2)+min(w$x2)))
w_lda <- lda(cl~x1+x2, data=w)
w_lda_pred <- predict(w_lda, w_grid)$class
w_rf <- randomForest(cl~x1+x2, data=w, ntree=2)
w_rf_pred <- predict(w_rf, w_grid)
w_grid <- w_grid |>
mutate(plda = w_lda_pred,
prf = w_rf_pred)
p1 <- ggplot(w_grid, aes(x=x1, y=x2, colour = plda)) +
geom_point(alpha=0.5, size=2) +
geom_line(data=w, aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Parametric") +
theme(legend.position = "none",
axis.text = element_blank())
p2 <- ggplot(w_grid, aes(x=x1, y=x2, colour = prf)) +
geom_point(alpha=0.5, size=2) +
geom_line(data=w, aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Non-parametric") +
theme(legend.position = "none",
axis.text = element_blank())
p1 + p2 + plot_layout(ncol=1)
```
:::
:::
## Reducible vs irreducible error
::: {.columns}
::: {.column width=10%}
:::
::: {.column width=30%}
```{r}
#| label: model-errors1
#| echo: false
#| fig-width: 5
#| fig-height: 5
w <- w |>
mutate(plda = predict(w_lda, w)$class,
prf = predict(w_rf, w)) |>
mutate(lda_err = ifelse(cl != plda, "1", "0"),
rf_err = ifelse(cl != prf, "1", "0"))
ggplot(w, aes(x=x1, y=x2,
colour = cl,
shape=lda_err)) +
geom_point(size=5, alpha=0.8) +
geom_line(aes(x=x1, y=y),
inherit.aes = FALSE,
colour="black") +
scale_color_discrete_qualitative() +
scale_shape_manual(values=c(1, 19)) +
ggtitle("Parametric errors") +
theme(legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank())
```
:::
::: {.column width=10%}
:::
::: {.column width=30%}
```{r}
#| label: model-errors2
#| echo: false
#| fig-width: 5
#| fig-height: 5
ggplot(w, aes(x=x1, y=x2,
colour = cl,
shape=rf_err)) +
geom_point(size=5, alpha=0.8) +
geom_line(aes(x=x1, y=y),
inherit.aes = FALSE,
colour="black") +
scale_color_discrete_qualitative() +
scale_shape_manual(values=c(1, 19)) +
ggtitle("Non-parametric errors") +
theme(legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank())
```
:::
::: {.column width=10%}
:::
:::
If the [model form is incorrect]{.monash-umber2}, the error (solid circles) may arise from wrong shape, and is thus [reducible]{.monash-umber2}. [Irreducible]{.monash-olive2} means that we have got the right model and mistakes (solid circles) are [random noise]{.monash-olive2}.
## Flexible vs inflexible
```{r}
#| label: flexible-model
#| echo: false
#| fig-width: 9
#| fig-height: 3
p1 <- ggplot(w_grid, aes(x=x1, y=x2, colour = plda)) +
geom_point(size=2, alpha=0.5) +
geom_line(data=w, aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Less flexible") +
theme(legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank())
w_rf <- randomForest(cl~x1+x2, data=w, ntree=4)
w_rf_pred <- predict(w_rf, w_grid)
w_grid <- w_grid |>
mutate(prf = w_rf_pred)
p2 <- ggplot(w_grid, aes(x=x1, y=x2, colour = prf)) +
geom_point(size=2, alpha=0.5) +
geom_line(data=w, aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("<--------------->") +
theme(legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank())
w_rf <- randomForest(cl~x1+x2, data=w, ntree=100)
w_rf_pred <- predict(w_rf, w_grid)
w_grid <- w_grid |>
mutate(prf = w_rf_pred)
p3 <- ggplot(w_grid, aes(x=x1, y=x2, colour = prf)) +
geom_point(size=2, alpha=0.5) +
geom_line(data=w, aes(x=x1, y=y), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("More flexible") +
theme(legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank())
p1 + p2 + p3 + plot_layout(ncol=3)
```
[Parametric]{.monash-orange2} models tend to be [less flexible]{.monash-orange2} but [non-parametric]{.monash-blue2} models can be flexible or less flexible depending on [parameter settings]{.monash-blue2}.
## Bias vs variance
::: {.columns}
::: {.column}
[Bias]{.monash-orange2} is the error that is introduced by modeling a
complicated problem by a simpler problem.
- For example, linear regression assumes a linear relationship and perhaps the relationship is not exactly linear.
- In general, but not always, the [more flexible]{.monash-orange2} a method is, the [less bias]{.monash-orange2} it will have because it can fit a complex shape better.
:::
::: {.column}
::: {.fragment}
[Variance]{.monash-orange2}
refers to how much your estimate would change if you had different training data. Its measuring how much your model depends on the data you have, to the neglect of future data.
- In general, the [more flexible]{.monash-orange2} a method is, the [more variance]{.monash-orange2} it has.
- The [size]{.monash-orange2} of the training data can impact on the variance.
:::
:::
:::
## Bias
::: {.columns}
::: {.column}
```{r}
#| label: bias1
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
p1 + ggtitle("Large bias")
```
When you impose too many assumptions with a parametric model, or use an inadequate non-parametric model, such as not letting an algorithm converge fully.
:::
::: {.column}
```{r}
#| label: bias2
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
p3 + ggtitle("Small bias")
```
When the model closely captures the true shape, with a parametric model or a flexible model.
:::
:::
## Variance
::: {.columns}
::: {.column}
```{r}
#| label: variance1
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
p1 + ggtitle("Small variance")
```
This fit will be virtually identical even if we had a different training sample.
:::
::: {.column}
```{r}
#| label: variance2
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
p3 + ggtitle("Large variance")
```
Likely to get a very different model if a different training set is used.
:::
:::
## Training vs test splits
::: {.info}
When data are reused for multiple tasks, instead of carefully *spent* from the finite data budget, certain risks increase, such as the risk of accentuating bias or compounding effects from methodological errors. [*Julia Silge*](https://www.tmwr.org/splitting)
:::
- [Training set]{.monash-blue2}: Used to fit the model, might be also broken into a validation set for frequent assessment of fit.
- [Test set]{.monash-blue2}: Purely used to assess final models performance on future data.
## Training vs test splits
::: {.columns}
::: {.column}
```{r}
#| label: balanced-data
d_bal <- tibble(y=c(rep("A", 6), rep("B", 6)),
x=c(runif(12)))
d_bal$y
set.seed(130)
d_bal_split <- initial_split(d_bal, prop = 0.70)
training(d_bal_split)$y
testing(d_bal_split)$y
```
:::
::: {.column}
```{r}
#| label: unbalanced-data
d_unb <- tibble(y=c(rep("A", 2), rep("B", 10)),
x=c(runif(12)))
d_unb$y
set.seed(132)
d_unb_split <- initial_split(d_unb, prop = 0.70)
training(d_unb_split)$y
testing(d_unb_split)$y
```
::: {.fragment}
Always [stratify splitting]{.monash-orange2} by sub-groups, especially response variable classes.
```{r}
#| label: unbalanced-split
d_unb_strata <- initial_split(d_unb, prop = 0.70, strata=y)
training(d_unb_strata)$y
testing(d_unb_strata)$y
```
:::
:::
:::
## Measuring accuracy for categorical response
Compute $\widehat{y}$ from [training data]{.monash-orange2}, $\{(y_i, {\mathbf x}_i)\}_{i = 1}^n$. The error rate ([fraction of misclassifications]{.monash-orange2}) to get the [Training Error Rate]{.monash-green2}
$$\text{Error rate} = \frac{1}{n}\sum_{i=1}^n I(y_i \ne \widehat{y}({\mathbf x}_i))$$
A better estimate of future [accuracy]{.monash-orange2} is obtained using [test data]{.monash-orange2} to get the [Test Error Rate]{.monash-green2}.
::: {.info}
Training error will usually be [smaller]{.monash-orange2} than test error. When it is much smaller, it indicates that the model is too well fitted to the training data to be accurate on future data (over-fitted).
:::
## Confusion (misclassification) matrix
::: {.columns}
::: {.column}
| | | predicted |
| | | `1` | `0` |
| true | `1` | `a` | `b` |
| | `0` | `c` | `d` |
Consider `1`=positive (P), `0`=negative (N).
- True positive (TP): `a`
- True negative (TN): `d`
- False positive (FP): `c` (Type I error)
- False negative (FN): `b` (Type II error)
:::
::: {.column}
- Sensitivity, recall, hit rate, or true positive rate (TPR): TP/P = `a`/(`a`+`b`)
- Specificity, selectivity or true negative rate (TNR): TN/N = `d`/(`c`+`d`)
- Prevalence: P/(P+N) = (`a`+`b`)/(`a`+`b`+`c`+`d`)
- Accuracy: (TP+TN)/(P+N) = (`a`+`d`)/(`a`+`b`+`c`+`d`)
- Balanced accuracy: (TPR + TNR)/2 (or average class errors)
:::
:::
## Confusion (misclassification) matrix: computing
::: {.columns}
::: {.column}
Two classes
```{r}
#| echo: false
#| label: predictive-class-example
a2 <- tibble(y = c(rep("bilby", 12),
rep("quokka", 15)),
pred = c(rep("bilby", 9),
rep("quokka", 13),
rep("bilby", 5)),
bilby = c(0.9, 0.8, 0.9,
0.7, 0.6, 0.8,
0.9, 0.7, 0.6,# true
0.4, 0.3, 0.4,# error
0.2, 0.4, 0.1,# true
0.4, 0.1, 0.3,
0.2, 0.4, 0.3,# true
0.4,
0.6, 0.7,# error
0.6, 0.9, 0.7)) |>
mutate(quokka = 1-bilby)
a3 <- a2 |>
bind_rows(tibble(
y = rep("numbat", 8),
pred = c(rep("numbat", 6),
rep("quokka", 2))))
a2 <- a2 |>
mutate(y = factor(y),
pred = factor(pred))
a3 <- a3 |>
mutate(y = factor(y),
pred = factor(pred))
```
```{r}
#| eval: false
#| echo: false
# tidymodels has it transposed
#| label: confusion-matrix
cm <- conf_mat(a2, y, pred)
autoplot(cm)
# Make it show in right direction
conf_mat(a2, pred, y, dnn=c("Truth", "Pred"))
```
```{r}
# Write out the confusion matrix in standard form
#| label: oconfusion-matrix-tidy
cm <- a2 |> count(y, pred) |>
group_by(y) |>
mutate(cl_acc = n[pred==y]/sum(n))
cm |>
pivot_wider(names_from = pred,
values_from = n) |>
select(y, bilby, quokka, cl_acc)
```
```{r}
accuracy(a2, y, pred) |> pull(.estimate)
bal_accuracy(a2, y, pred) |> pull(.estimate)
sens(a2, y, pred) |> pull(.estimate)
specificity(a2, y, pred) |> pull(.estimate)
```
:::
::: {.column}
More than two classes
```{r}
# Write out the confusion matrix in standard form
cm3 <- a3 |> count(y, pred) |>
group_by(y) |>
mutate(cl_err = n[pred==y]/sum(n))
cm3 |>
pivot_wider(names_from = pred,
values_from = n, values_fill=0) |>
select(y, bilby, quokka, numbat, cl_err)
```
```{r}
accuracy(a3, y, pred) |> pull(.estimate)
bal_accuracy(a3, y, pred) |> pull(.estimate)
```
:::
:::
## Receiver Operator Curves (ROC)
::: {.columns}
::: {.column}
The balance of getting it right, without predicting everything as positive.
{style="font-size: 50%;"}
Need [predictive probabilities]{.monash-orange2}, probability of being each class.
:::
::: {.column}
```{r}
#| label: roc-curve
a2 |> slice_head(n=3)
roc_curve(a2, y, bilby) |>
autoplot()
```
:::
:::
## Preparing to fit model: IDA {.transition-slide .center style="text-align: center;"}
##
::: {.column width=60%}
*The first thing to do with data is to [look at them]{.monash-blue2} .... usually means [tabulating]{.monash-blue2} and [plotting]{.monash-blue2} the data in many different ways to [see what's going on]{.monash-blue2}. With the wide availability of computer packages and graphics nowadays there is no excuse for ducking the labour of this preliminary phase, and it may save some* [**red faces**]{.monash-red2} *later.*
[[Crowder, M. J. & Hand, D. J. (1990) "Analysis of Repeated Measures"](https://doi.org/10.1201/9781315137421)]{.smallest}
:::
## Pre-processing steps
- **Handle missing values**: most modeling methods require complete data.
- **Decide on an appropriate model type**: Explore relationships between response and predictors: nonlinear, clustered. Different shapes indicate which model will likely fit better
- **Check model assumptions**
- **Feature engineering**: Do you have the predictors needed, or could you create better predictors? Should some variables be transformed.
- **Check for anomalies**: Are there some observations that are unusal and might affect the fit? (We'll do this at the end, too.)
- **Split data**: Decide on how much of the data should be used for training a model, validation, and testing.
- **Reduce the number of predictors**: Too many predictors can result in over-fitting the training data.
## Check model assumptions
- **Statistical models** tend to be less flexible. They impose strict assumptions, such as *linear form*, and *Gaussian errors*. Check by plotting response vs predictors to check the relationship is what is assumed, and spread is uniform.
- **Non-parametric models, and algorithms**, don't explicitly impose assumptions, but the *process imposes constraints*, that means fit will fail if not satisfied. For example, tree-based methods mostly use a single variable for any split. If the relationship with the response is best explained by a combination of variables, this will result in a poor fit.
## Feature engineering
Creating new variables to get better fits is a special skill! Sometimes automated by the method. All are transformations of the original variables. (See `tidymodels` steps.)
- scaling, centering, sphering (`step_pca`)
- log or square root or box-cox transformation (`step_log`)
- ratio of values (`step_ratio`)
- polynomials or splines: $x_1^2, x_1^3$ (`step_ns`)
- dummy variables: categorical predictors expanded into multiple new binary variables (`step_dummy`)
- Convolutional Neural Networks: neural networks but with pre-processing of images to combine values of neighbouring pixels; flattening of images
## After fitting: diagnostics {.transition-slide .center style="text-align: center;"}
## Diagnosing the fit
::: {.columns}
::: {.column}
Compute and examine the [usual diagnostics]{.monash-blue2}, some methods have more
- fit statistics: accuracy, sensitivity, specificity
- errors/misclassifications
- variable importance
- plot residuals, examine the misclassifications
- [check test set is similar to training]{.monash-pink2}
[Go beyond ...]{.monash-orange2} Look at the data and the model together!
[[Wickham et al (2015) Removing the Blindfold](http://onlinelibrary.wiley.com/doi/10.1002/sam.11271/abstract)]{.smaller}
:::
::: {.column}
::: {.fragment}
[*Training - plusses; Test - dots*]{.smaller .center}
```{r}
#| label: training-test
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 80%
w_test <- read_csv("data/sine-curve-test.csv")
ggplot() +
geom_point(data=w_grid, aes(x=x1, y=x2,
colour = prf),
alpha=0.05, size=4) +
geom_point(data=w, aes(x=x1, y=x2, colour = cl),
shape=3, size=3) +
geom_point(data=w_test, aes(x=x1, y=x2, colour = cl),
shape=19, size=3) +
scale_color_discrete_qualitative() +
theme(legend.position = "none",
axis.text = element_blank())
```
:::
:::
:::
## Model-in-the-data-space
:::: {.columns}
::: {.column width=60%}
{width=450 style="font-size: 50%;"}
:::
::: {.column width=40%}
We plot the [model on the data]{.monash-orange2} to assess whether it [fits]{.monash-orange2} or is a [misfit]{.monash-orange2}!
::: {.fragment}
Doing this in high dimensions is considered difficult!
:::
::: {.fragment}
So it is common to only plot the [*data-in-the-model-space*]{.monash-blue2}.
:::
::: {.fragment}
NOTE, WE CAN PLOT THESE THINGS IN HIGH DIMENSIONS. But this is for another day.
:::
:::
::::
## Data-in-the-model-space
:::: {.columns}
::: {.column}
```{r}
#| label: data-in-model-space1
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
#| code-fold: true
w <- read_csv("data/sine-curve.csv") |>
mutate(cl = factor(cl))
w_rf <- randomForest(cl~x1+x2, data=w, ntree=200)
w_rf_pred <- predict(w_rf, w, type="prob") |>
as_tibble(.name_repair="unique")
w_rf_pred |>
bind_cols(w) |>
select(`A`, cl) |>
ggplot(aes(x=`A`, colour=cl, fill=cl)) +
geom_density(alpha=0.7) +
xlab("Prob class A") +
ggtitle("Random forest") +
theme(legend.position="none")
```
:::
::: {.column}
```{r}
#| label: data-in-model-space2
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
#| code-fold: true
w_lda <- lda(cl~x1+x2, data=w)
w_lda_pred <- predict(w_lda, w, method="predictive")$posterior |>
as_tibble(.name_repair="unique")
w_lda_pred |>
bind_cols(w) |>
select(`A`, cl) |>
ggplot(aes(x=`A`, colour=cl, fill=cl)) +
geom_density(alpha=0.7) +
xlab("Prob class A") +
ggtitle("Linear DA") +
theme(legend.position="none")
```
:::
::::
Predictive probabilities are aspects of the model. It is useful to plot. What do we learn here?
::: {.fragment}
[But it doesn't tell you why there is a difference.]{.monash-orange2}
:::
## Model-in-the-data-space
:::: {.columns}
::: {.column}
```{r}
#| label: model-in-the-data-space1
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
library(geozoo)
w_grid <- cube.solid.random(p=2)$points
colnames(w_grid) <- c("x1", "x2")
w_grid <- w_grid |> as_tibble() |>
mutate(x1 = x1*(max(w$x1)-min(w$x1)+min(w$x1)),
x2 = x2*(max(w$x2)-min(w$x2)+min(w$x2)))
w_lda <- lda(cl~x1+x2, data=w)
w_lda_pred <- predict(w_lda, w_grid)$class
w_rf <- randomForest(cl~x1+x2, data=w, ntree=200)
w_rf_pred <- predict(w_rf, w_grid)
w_grid <- w_grid |>
mutate(plda = w_lda_pred,
prf = w_rf_pred)
ggplot(w_grid, aes(x=x1, y=x2, colour = plda)) +
geom_point(alpha=0.5, size=2) +
geom_text(data=w, aes(x=x1, y=x2, label=cl), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Linear DA") +
theme(legend.position = "none",
axis.text = element_blank())
```
:::
::: {.column}
```{r}
#| label: model-in-the-data-space2
#| echo: false
#| fig-width: 4
#| fig-height: 4
#| out-width: 60%
ggplot(w_grid, aes(x=x1, y=x2, colour = prf)) +
geom_point(alpha=0.5, size=2) +
geom_text(data=w, aes(x=x1, y=x2, label=cl), colour="black") +
scale_color_discrete_qualitative() +
ggtitle("Random forest") +
theme(legend.position = "none",
axis.text = element_blank())
```
:::
::::
::: {style="font-size: 85%;"}
Model is displayed, as a grid of predicted points in the original variable space. Data is overlaid, using text labels. What do you learn?
::: {.fragment}
[One model has a linear boundary, and the other has the highly non-linear boundary, which matches the class cluster better. Also ...]{.monash-orange2}
:::
:::
## Plotting residuals
:::: {.columns}
::: {.column}
```{r}
#| label: surreal1
#| code-fold: true
#| fig-width: 10
#| fig-height: 2.5
#| out-width: 100%
# Load the Jack-o'-Lantern data
d <- jackolantern_surreal_data[sample(1:5395, 2000),]
ggduo(d, c("x1", "x2", "x3", "x4"), "y")
```
Data looks uninteresting: weak linear relationships between response and predictors.
:::
::: {.column}
::: {.fragment}
```{r}
#| label: surreal2
#| code-fold: true
#| fig-width: 5
#| fig-height: 5
#| out-width: 70%
# Fit a linear model to the surreal Jack-o'-Lantern data
d_lm <- lm(y ~ ., data = d)
# Plot the residuals to reveal the hidden image
d_all <- augment(d_lm, d)
ggplot(d_all, aes(x=.fitted, y=.resid)) +
geom_point() +
theme(aspect.ratio=1)
```
Residuals have a *spooky* pattern!
:::
:::
::::
## Reading residual plots using a lineup
:::: {.columns}
::: {.column width="30%"}
- Reading residual plots can be extremely hard!
- You have to decide whether there is NO PATTERN.
- This is a good place to use the lineup developed to do inference for data visualisation.
:::
::: {.column width="70%"}
```{r}
#| label: tips-lineup
#| fig-height: 10
#| fig-width: 10
#| out-width: 60%
#| code-fold: true
x <- lm(tip ~ total_bill, data = tips)
tips.reg <- data.frame(tips, .resid = residuals(x),
.fitted = fitted(x))
library(ggplot2)
ggplot(lineup(null_lm(tip ~ total_bill,
method = 'rotate'),
tips.reg)) +
geom_point(aes(x = total_bill,
y = .resid)) +
facet_wrap(~ .sample) +
theme(axis.text = element_blank(),
axis.title = element_blank())
```
:::
::::
## Re-sampling statistics
:::: {.columns}
::: {.column}
Techniques for assessing behaviour with new samples, and measuring the variability of models and predictions.
- Leave-one-out (algebraic measures typically) - outmoded
- **Cross-validation**: Multiple splits to repeatedly use one for training and the other for testing, to decide on best model parameters.
- **Bootstrap**: sample the sample with replacement. Mostly used to compute confidence intervals for estimates.
:::
::: {.column}
- **Permutation**: sample the values on one variable to break associations with other variables while keeping marginals fixed. Mostly used for conducting hypothesis tests.
- **Simulation**: sample from know distribution, or a fitted model. Used to check goodness-of-fit.
:::
::::
## Tidy models {.transition-slide .center style="text-align: center;"}
## Tidying model output
:::: {.columns}
::: {.column}
The package `broom` gets model results into a tidy format at different levels. (And `broom.mixed` does this for mixed models.)
- **model level**: `broom::glance` (values such as AIC, BIC, model fit, ...)
- **coefficients in the model**: `broom::tidy` (estimate, confidence interval, significance level, ...)
- **observation level**: `broom::augment` (fitted values, residuals, predictions, influence, ...)
:::
::: {.column}
```{r}
#| code-fold: true
glance(d_lm)
tidy(d_lm)
```
:::
::::
## `tidymodels`
`tidymodels` is a set of R packages that make it easier to build, tune, and evaluate statistical and machine learning models — all using the same consistent, tidyverse-style syntax. It provides a framework for the modeling workflow in R:
1. Preprocess data (with `recipes`)
2. Specify models (with `parsnip`)
3. Resample or cross-validate (with `rsample`)
4. Tune hyperparameters (with `tune`)
5. Evaluate and compare models (with `yardstick`, `workflows`, `broom`)
## Case study: hotel bookings [(1/10)]{.smallest}
Source: [A predictive modeling case study](https://www.tidymodels.org/start/case-study/)
:::: columns
::: column
**About the data**: The hotel bookings data from [Antonio, Almeida, and Nunes (2019)](https://doi.org/10.1016/j.dib.2018.11.126) to predict which hotel stays included children and/or babies, based on the other characteristics of the stays such as which hotel the guests stay at, how much they pay, etc.
**Model**: We will build a model to predict which actual hotel stays included children and/or babies, and which did not. Our outcome variable `children` is a factor variable with two levels.
:::
::: column
::: {style="font-size: 70%;"}
```{r}
#| code-fold: true
library(tidymodels)
load("data/hotels.rda")
glimpse(hotels)
```
```{r}
#| code-fold: true
hotels |>
count(children) |>
mutate(prop = n/sum(n))
```
:::
:::
::::
## Case study: hotel bookings [(2/10)]{.smallest}
:::: columns
::: column
**Data splitting and re-sampling**
Reserve 25% of the stays to the test set. The outcome variable children is pretty imbalanced so use a stratified random sample.
:::
::: column
```{r}
#| code-fold: true
set.seed(656)
splits <- initial_split(hotels, strata = children)
hotel_other <- training(splits)
hotel_test <- testing(splits)
# training set proportions by children
hotel_other |>
count(children) |>
mutate(prop = n/sum(n))
# test set proportions by children
hotel_test |>
count(children) |>
mutate(prop = n/sum(n))
```
:::
::::
## Case study: hotel bookings [(3/10)]{.smallest}
:::: columns
::: column
**Data splitting and re-sampling**
The "other" set is further divided into training and validation sets. These two are used to build the model.
The test set is only used when the final model is decided.
:::
::: column
```{r}
#| code-fold: true
set.seed(703)
hotel_val_set <- validation_split(hotel_other,
strata = children,
prop = 0.80)
hotel_val_set
```
:::
::::
## Case study: hotel bookings [(4/10)]{.smallest}
:::: columns
::: column
**IDA**
- Check for distribution of each variable
- Check association between predictors
:::
::: column
```{r}
#| code-fold: true
p1 <- ggplot(hotel_other, aes(x=children, y=lead_time)) +
geom_quasirandom()
p2 <- ggplot(hotel_other, aes(x=children, y=stays_in_weekend_nights)) +
geom_quasirandom()
p1 + p2 + plot_layout(ncol=2)
```
:::
::::
## Case study: hotel bookings [(5/10)]{.smallest}
:::: columns
::: column
**IDA**
- Check for distribution of each variable
- Check association between predictors
:::
::: column
```{r}
#| code-fold: true
#| fig-width: 10
#| fig-height: 8
p1 <- ggplot(hotel_other, aes(x=children, y=lead_time)) +
geom_boxplot()
p2 <- ggplot(hotel_other, aes(x=children, y=stays_in_weekend_nights)) +
geom_boxplot()
p3 <- ggplot(hotel_other, aes(x=children, y=average_daily_rate)) +
geom_boxplot()
p4 <- ggplot(hotel_other,
aes(x=reserved_room_type,
fill=children)) +
geom_bar(position="fill")
p1 + p2 + p3 + p4 + plot_layout(ncol=2)
```
:::
::::
## Case study: hotel bookings [(6/10)]{.smallest}
:::: columns
::: column
**Model setup**
Since the outcome variable is categorical, logistic regression is a reasonable choice. This fits an S-curve which essentially models proportion of classes (children v none) for different values of the predictors. Using `glmnet` allows for feature selection while training the model.
:::
::: column
::: {.panel-tabset}
## model
```{r}
#| code-fold: false
lr_mod <-
logistic_reg(penalty = tune(), mixture = 1) |>
set_engine("glmnet")
```
## recipe
```{r}
#| code-fold: true
holidays <- c("AllSouls", "AshWednesday", "ChristmasEve", "Easter",
"ChristmasDay", "GoodFriday", "NewYearsDay", "PalmSunday")
lr_recipe <-
recipe(children ~ ., data = hotel_other) |>
step_date(arrival_date) |>
step_holiday(arrival_date, holidays = holidays) |>
step_rm(arrival_date) |>
step_dummy(all_nominal_predictors()) |>
step_zv(all_predictors()) |>
step_normalize(all_predictors())
```
## workflow
```{r}
#| code-fold: true
lr_workflow <-
workflow() |>
add_model(lr_mod) |>
add_recipe(lr_recipe)
```
## tune
```{r}
#| code-fold: true
lr_reg_grid <- tibble(penalty = 10^seq(-4, -1, length.out = 30))
lr_reg_grid |> top_n(-5) # lowest penalty values
lr_reg_grid |> top_n(5) # highest penalty values
```
:::
:::
::::
## Case study: hotel bookings [(7/10)]{.smallest}
:::: columns
::: column
**Model training**
Better performance for smaller penalty, means that all predictors have some importance for the model fit.
:::
::: column
::: {.panel-tabset}
## tuning
```{r}
#| code-fold: true
lr_res <-
lr_workflow |>
tune_grid(hotel_val_set,
grid = lr_reg_grid,
control = control_grid(save_pred = TRUE),
metrics = metric_set(roc_auc))
```
## summary
```{r}
#| code-fold: true
lr_plot <-
lr_res |>
collect_metrics() |>
ggplot(aes(x = penalty, y = mean)) +
geom_point() +
geom_line() +
ylim(c(0.75, 0.9)) +
ylab("Area under the ROC Curve") +
scale_x_log10(labels = scales::label_number())
lr_plot
```
## best
```{r}
#| code-fold: true
lr_best <-
lr_res |>
collect_metrics() |>
arrange(penalty) |>
slice(12)
lr_best
```
## ROC
```{r}
#| code-fold: true
lr_auc <-
lr_res |>
collect_predictions(parameters = lr_best) |>
roc_curve(children, .pred_children) |>
mutate(model = "Logistic Regression")
autoplot(lr_auc)
```
:::
:::
::::
## Case study: hotel bookings [(8/10)]{.smallest}
:::: columns
::: column
**Alternative model**
Random forest models are generalist fits. It combines the predictions from many trees, each one fitted with a random selection of variables (features) and a random sample of observations.
Because each tree splits on single variables, it is also robust to missing values by using the complete data efficiently. It also handles the predictors as factors easily.
:::
::: column
::: {style="font-size: 70%;"}
::: {.panel-tabset}
## model
```{r}
#| code-fold: true
cores <- parallel::detectCores()
cores
```
```{r}
#| code-fold: false
rf_mod <-
rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |>
set_engine("ranger", num.threads = cores) |>
set_mode("classification")
```
## recipe
```{r}
#| code-fold: false
rf_recipe <-
recipe(children ~ ., data = hotel_other) |>
step_date(arrival_date) |>
step_holiday(arrival_date) |>
step_rm(arrival_date)
```
## workflow
```{r}
#| code-fold: true
rf_workflow <-
workflow() |>
add_model(rf_mod) |>
add_recipe(rf_recipe)
```
```{r}
#| code-fold: true
rf_mod
extract_parameter_set_dials(rf_mod)
```
## tune
```{r}
#| code-fold: true
#| eval: false
# This takes a while so the results are saved
set.seed(1012)
rf_res <-
rf_workflow |>
tune_grid(hotel_val_set,
grid = 25,
control = control_grid(save_pred = TRUE),
metrics = metric_set(roc_auc))
save(rf_res, file="data/rf_res.rda")
```
```{r}
#| code-fold: true
load("data/rf_res.rda")
rf_res |>
show_best(metric = "roc_auc")
```
```{r}
#| code-fold: true
#| out-width: 50%
autoplot(rf_res)
```
## summary
```{r}
#| code-fold: true
rf_best <-
rf_res |>
select_best(metric = "roc_auc")
rf_best
```
```{r}
#| code-fold: true
rf_res |>
collect_predictions()
```
```{r}
#| code-fold: true
rf_auc <-
rf_res |>
collect_predictions(parameters = rf_best) |>
roc_curve(children, .pred_children) |>
mutate(model = "Random Forest")
```
:::
:::
:::
::::
## Case study: hotel bookings [(9/10)]{.smallest}
:::: columns
::: column
**Choose the best model**
For the same false positive rate (predicting the observation to the target class, incorrectly), the random forest model has a higher true positive rate (correctly predicting the target class) than the logistic regression model.
Most important variables are `average_daily_room_rate`, `reserved_room_type`, `assigned_room_type`, `lead_time`.
:::
::: column
::: {.panel-tabset}
## ROC
```{r}
#| code-fold: true
bind_rows(rf_auc, lr_auc) |>
ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) +
geom_path(lwd = 1.5, alpha = 0.8) +
geom_abline(lty = 3) +
coord_equal() +
scale_color_viridis_d(option = "plasma", end = .6)
```
## best
```{r}
#| code-fold: true
# the last model
last_rf_mod <-
rand_forest(mtry = 8, min_n = 7, trees = 1000) |>
set_engine("ranger", num.threads = cores, importance = "impurity") |>
set_mode("classification")
# the last workflow
last_rf_workflow <-
rf_workflow |>
update_model(last_rf_mod)
```
```{r}
#| code-fold: true
#| eval: false
# the last fit
set.seed(345)
last_rf_fit <-
last_rf_workflow |>
last_fit(splits)
last_rf_fit
save(last_rf_fit, file="data/last_rf_fit.rda")
```
```{r}
#| code-fold: true
load("data/last_rf_fit.rda")
last_rf_fit |>
collect_metrics()
```
## vip
```{r}
#| code-fold: true
last_rf_fit |>
extract_fit_parsnip() |>
vip(num_features = 20)
```
## final
```{r}
#| code-fold: true
last_rf_fit |>
collect_predictions() |>
roc_curve(children, .pred_children) |>
autoplot()
```
:::
:::
::::
## Case study: hotel bookings [(10/10)]{.smallest}
:::: columns
::: column
**Diagnostics and summaries**
- Plot the training data to understand the model fit. (Not done here)
- Plot the test data to understand where it likely errs on new data.
- Examine where the model misclassified the observation, to check if there is anything systematic.
:::
::: column
::: {.panel-tabset}
## Predictions
```{r}
#| code-fold: true
hotel_test_pred <- tibble(hotel_test) |>
mutate(
pchild_cl = last_rf_fit$.predictions[[1]]$.pred_class,
pchild_pc = last_rf_fit$.predictions[[1]]$.pred_children,
pchild_pn = last_rf_fit$.predictions[[1]]$.pred_none)
hotel_test_pred |>
select(pchild_pc, pchild_pn, pchild_cl) |>
ggplot(aes(x=pchild_pc,
fill=pchild_cl,
colour=pchild_cl)) +
geom_density()
```
## Errors
```{r}
#| code-fold: true
hotel_test_pred <- hotel_test_pred |>
mutate(error = if_else(children != pchild_cl,
"yes", "no"))
ggplot(hotel_test_pred, aes(x=children,
y=average_daily_rate,
colour = children)) +
geom_boxplot() +
scale_x_discrete("", labels = c("C", "N")) +
facet_grid(error~reserved_room_type) +
theme(legend.position = "none")
```
:::
:::
::::
## {.information-slide .center style="text-align: center;"}
Tidy models makes it easy to compare and contrast fits, and to follow all the steps in pre-processing and fitting. Many case studies and examples to follow on the web site.
More work is needed to check assumptions, and to summarise the fit well.
## Unsupervised learning {.transition-slide .center style="text-align: center;"}
## Dimension reduction
:::: {.columns}
::: {.column}
1. Decide on necessary variables to keep, and unusable variables to drop. (All variables have to be **numeric**.)
2. Handle missing values - all methods require **complete data**.
3. Possibly transform variables to **symmetric**. Outliers can affect dimension reduction. (Similarly, clusters and nonlinear relationships can affect dimension reduction.)
4. Choose method, typically **principal component analysis** for simplicity and interpretability.
:::
::: {.column}
::: {.fragment}
::: {style="font-size: 80%;"}
Principal component analysis (PCA) produces a low-dimensional representation of a
dataset. It finds a sequence of linear combinations of the
variables that have [maximal variance]{.monash-orange2}, and are [mutually uncorrelated]{.monash-orange2}. It is an unsupervised learning method.
Use it, when:
- You have too many predictors for a regression. Instead, we can use the first few principal components.
- Need to understand relationships between variables.
- To make plots summarising the variation in a large number of variables.
*Make sure you use standardised variables!*
:::
:::
:::
::::
## Principal component analysis (PCA) [(1/3)]{.smallest}
:::: {.columns}
::: {.column}
- **new variables**: principal components are linear combinations (eigenvectors) of the original variables
- **summary of variance**: total variance, variance of new variables (eigenvalues), proportion of total variance
- choosing number of new variables
- summary of new variables
- checking what is left over
:::
::: {.column}
`hotels` data:
::: {style="font-size: 60%;"}
```{r}
#| code-fold: true
glimpse(hotels)
```
Usable variables: `lead_time`, `stays_XXX`, `adults`, `is_repeated_guest`, `previous_XXX`, `booking_changes`,
`days_in_waiting_list`, `average_daily_rate`, `total_of_special_requests`.
```{r}
#| code-fold: true
hotels_dr <- hotels |>
select(lead_time:adults,
is_repeated_guest:previous_bookings_not_canceled,
booking_changes, days_in_waiting_list,
average_daily_rate, total_of_special_requests)
```
:::
:::
::::
## Principal component analysis (PCA) [(2/3)]{.smallest}
:::: {.columns}
::: {.column}
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
hotels_pca <- prcomp(hotels_dr, scale=TRUE)
hotels_pca$sdev^2
mulgar::ggscree(hotels_pca, q=11)
```
:::
::: {.column}
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 6
#| out-width: 80%
ggfortify:::autoplot.prcomp(hotels_pca,
loadings = TRUE,
loadings.label = TRUE)
```
:::
::::
## Principal component analysis (PCA) [(3/3)]{.smallest}
:::: {.columns}
::: {.column}
PCA doesn't help to reduce dimension for the `hotels` data.
However, it could be helpful to understand some of the other patterns in the data, e.g. outliers, clusters.
:::
::: {.column}
```{r}
#| code-fold: true
#| fig-width: 6
#| fig-height: 6
#| out-width: 80%
ggfortify:::autoplot.prcomp(hotels_pca, x=6, y=7,
loadings = TRUE,
loadings.label = TRUE)
```
:::
::::
## Cluster analysis
:::: {.columns}
::: {.column}
- The aim of cluster analysis is to [group cases (objects) according to their similarity]{.monash-blue2} on the variables. It is also often called unsupervised classification, meaning that classification is the ultimate goal, but the classes (groups) are not known ahead of time.
- Hence the first task in cluster analysis is to construct the class information. To determine closeness we start with [measuring the interpoint distances]{.monash-blue2}.
:::
::: {.column}
Common methods:
- Hierarchical clustering, summarised by a dendrogram
- $k$-means clustering
- Model-based clustering
and many variations.
Clustering is one of the hardest analytical methods because there is no known truth. Deciding on a solution requires external knowledge and understanding of geometry.
:::
::::
## Case study: Australian tourism [(1/3)]{.smallest}
:::: {.columns}
::: {.column}
- Survey of 563 Australian tourists, see [Dolnicar S, Grün B, Leisch F (2018)](https://link.springer.com/book/10.1007/978-981-10-8818-6)
- Six different types of risks: recreational, health, career, financial, social and safety
- Rated on a scale from 1 (never) to 5 (very often)
Here we will use $k$-means, and examine solutions for $k=2,3,4,5$ clusters. We will use the [lionfish](https://mmedl94.github.io/lionfish/) R package, that is mostly written in python for highly reactive interactive graphics.
:::
::: {.column}
Clustering tends to break up data according to it's shape.
::: {layout-ncol=3}
{fig-alt="Whole apple and knife on cutting board."}
{fig-alt="Apple in two halves and knife on cutting board."}
{fig-alt="Apple sliced into eighths and knife on cutting board."}
{fig-alt="Whole banana and knife on cutting board."}
{fig-alt="Banana in half and knife on cutting board."}
{fig-alt="Banana cut into eight coin-shaped pieces and knife on cutting board."}
:::
:::
::::
## Case study: Australian tourism [(2/3)]{.smallest}
:::: {.columns}
::: {.column width=33%}
{width=450 fig-alt="Screenshot of the lionfish interface showing the two cluster result. The projection shown is where there is a fairly clean line separating the two clusters."}
{width=450 fig-alt="Screenshot of the lionfish interface showing the four cluster result. The projection shown is where there is a fairly clear separation the four clusters, along the direction of the largest spread of points, but the fattest part of the pear shape has been divided into two. So the four clusters are spread along the main direction of spread, with two side-by-side in the fat part of the pear."}
:::
::: {.column width=33%}
{width=450 fig-alt="Screenshot of the lionfish interface showing the three cluster result. The projection shown is where there is a fairly clear separation the three clusters, along the direction of the largest spread of points."}
{width=450 fig-alt="Screenshot of the lionfish interface showing the five cluster result. The projection shown is where there is a fairly clear separation the four clusters, along the direction of the largest spread of points, but the fattest part of the pear shape has been divided into three. The clusters at the bottom and the top of the pear shape have been faded so we can focus on the three clusters in the fat part of the pear."}
:::
::: {.column width=33%}
::: {style="font-size: 80%;"}
- Up to 3 clusters, clustering is along the main axis of variation (i.e. PC 1), which includes most of the different types of risks.
- With 4-5 clusters it divides up medium values of some activities.
:::
:::
::::
## Case study: Australian tourism [(3/3)]{.smallest}
:::: {.columns}
::: {.column width=80%}
::: {.f50}
```{r}
#| code-fold: true
#| fig-width: 8
#| fig-height: 6
#| out-width: 60%
library(lionfish)
data("risk")
colnames(risk) <- c("Rec", "Hea", "Car", "Fin", "Saf", "Soc")
risk_d <- apply(risk, 2, function(x) (x-mean(x))/sd(x))
# Two clusters
nc <- 3
set.seed(1145)
r_km <- kmeans(risk_d, centers=nc,
iter.max = 500, nstart = 5)
r_km_d <- risk |>
as_tibble() |>
mutate(cl = factor(r_km$cluster))
r_km_d |>
pivot_longer(Rec:Soc, names_to = "var", values_to = "val") |>
ggplot(aes(x=val, fill=cl)) +
geom_bar() +
facet_wrap(~var, ncol=3, scales="free_y") +
scale_fill_manual(values = c("#377EB8", "#FF7F00", "#4DAF4A")) +
xlab("") + ylab("") +
theme_minimal() +
theme(legend.position = "none",
axis.text = element_blank())
```
:::
:::
::: {.column width=20%}
All the activities contribute to the segmentation into three clusters.
:::
::::
## {.information-slide .center style="text-align: center;"}
Unsupervised learning provides a way to organise observations or variables into smaller sets. It can help to make big noisy data more digestible - *if done well*.
## Temporal data {.transition-slide .center style="text-align: center;"}
## `tsibble`: R temporal object {background-image="https://tsibble.tidyverts.org/reference/figures/logo.png" background-position="5% 15%" background-size="10%"}
The `tsibble` package provides a data infrastructure for tidy temporal data with wrangling tools. Adapting the tidy data principles, `tsibble` is a data- and model-oriented object. In `tsibble`:
- [Index]{.monash-blue2} is a variable with inherent ordering from past to present.
- [Key]{.monash-blue2} is a set of variables that define observational units over time.
- Each observation should be [uniquely identified]{.monash-blue2} by index and key.
- Each observational unit should be measured at a [common interval]{.monash-blue2}, if regularly spaced.
## Regular vs irregular
:::: {.columns}
::: {.column}
::: {style="font-size: 80%l"}
The Melbourne pedestrian sensor data has a [regular]{.monash-blue2} period. Counts are provided for every hour, at numerous locations.
```{r}
#| label: ped-reg
options(width=55)
pedestrian
```
:::
:::
::: {.column}
::: {style="font-size: 80%l"}
::: {.fragment}
In contrast, the US flights data, below, is [irregular]{.monash-blue2}.
```{r}
#| label: nycflights
options(width=55)
flights_ts <- flights |>
mutate(dt = ymd_hm(paste(paste(year, month, day, sep="-"),
paste(hour, minute, sep=":")))) |>
as_tsibble(index = dt, key = c(origin, dest, carrier, tailnum), regular = FALSE)
flights_ts
```
:::
:::
:::
::::
::: {.fragment}
[Is pedestrian traffic really regular?]{.monash-orange2}
:::
## Getting started
::: {.info}
Wrangling prior to analysing temporal data includes:
- **aggregate** by temporal unit.
- **construct temporal units** to study seasonality, such as month, week, day of the week, quarter, ...
- checking and imputing **missings**.
:::
For the airlines data, you can aggregate by multiple quantities, eg number of arrivals, departures, average hourly arrival delay and departure delays.
## Aggregating
:::: {.columns}
::: {.column}
The US flights data already has some temporal components created, so aggregating by these is easy. Here is departure delay.
```{r}
#| code-fold: true
flights_mth <- flights_ts |>
as_tibble() |>
group_by(month, origin) |>
summarise(dep_delay = mean(dep_delay, na.rm=TRUE)) |>
as_tsibble(key=origin, index=month)
ggplot(flights_mth, aes(x=month, y=dep_delay, colour=origin)) +
geom_point() +
geom_smooth(se=F) +
scale_x_continuous("", breaks = seq(1, 12, 1),
labels=c("J","F","M","A","M","J",
"J","A","S","O","N","D")) +
scale_y_continuous("av dep delay (mins)", limits=c(0, 25)) +
theme(aspect.ratio = 0.5)
```
:::
::: {.column}
Aggregate by month, but examine arrival delays.
```{r}
#| code-fold: true
flights_mth_arr <- flights_ts |>
as_tibble() |>
group_by(month, origin) |>
summarise(arr_delay = mean(arr_delay, na.rm=TRUE)) |>
as_tsibble(key=origin, index=month)
ggplot(flights_mth_arr, aes(x=month, y=arr_delay, colour=origin)) +
geom_point() +
geom_smooth(se=F) +
scale_x_continuous("", breaks = seq(1, 12, 1),
labels=c("J","F","M","A","M","J",
"J","A","S","O","N","D")) +
scale_y_continuous("av arr delay (mins)", limits=c(0, 25)) +
theme(aspect.ratio = 0.5)
```
:::
::::
## Constructing temporal units
:::: {.columns}
::: {.column style="font-size: 80%;"}
Week day vs weekend would be expected to have different patterns of delay, but this is not provided.
```{r}
#| code-fold: true
#| fig-height: 8
#| fig-width: 5
#| out-width: 48%
flights_wk <- flights_ts |>
as_tibble() |>
mutate(wday = wday(dt, label=TRUE, week_start = 1)) |>
group_by(wday, origin) |>
summarise(dep_delay = mean(dep_delay, na.rm=TRUE)) |>
mutate(weekend = ifelse(wday %in% c("Sat", "Sun"), "yes", "no")) |>
as_tsibble(key=origin, index=wday)
ggplot(flights_wk, aes(x=wday, y=dep_delay, fill=weekend)) +
geom_col() +
facet_wrap(~origin, ncol=1, scales="free_y") +
xlab("") +
ylab("av dep delay (mins)") +
theme(aspect.ratio = 0.5, legend.position = "none")
```
:::
::: {.column style="font-size: 80%;"}
::: {.fragment}
Be careful of times!
```{r}
#| code-fold: true
flights_airtm <- flights |>
mutate(dep_min = dep_time %% 100,
dep_hr = dep_time %/% 100,
arr_min = arr_time %% 100,
arr_hr = arr_time %/% 100) |>
mutate(dep_dt = ymd_hm(paste(paste(year, month, day, sep="-"),
paste(dep_hr, dep_min, sep=":")))) |>
mutate(arr_dt = ymd_hm(paste(paste(year, month, day, sep="-"),
paste(arr_hr, arr_min, sep=":")))) |>
mutate(air_time2 = as.numeric(difftime(arr_dt, dep_dt)))
fp <- flights_airtm |>
sample_n(3000) |>
ggplot(aes(x=air_time, y=air_time2, label = paste(origin, dest))) +
geom_abline(intercept=0, slope=1) +
geom_point()
ggplotly(fp, width=500, height=500)
```
Why is this not what we expect?
:::
:::
::::
## Checking and filling missings [(1/4)]{.smallest}
:::: {.columns}
::: {.column}
```{r}
#| label: missings-simple
set.seed(328)
harvest <- tsibble(
year = c(2010, 2011, 2013, 2011,
2012, 2013),
fruit = rep(c("kiwi", "cherry"),
each = 3),
kilo = sample(1:10, size = 6),
key = fruit, index = year
)
harvest
```
:::
::: {.column}
```{r}
#| label: missing-gaps
has_gaps(harvest, .full = TRUE)
```
[Can you see the gaps in time?]{.monash-orange2}
Both levels of the key have missings.
:::
::::
## Checking and filling missings [(2/4)]{.smallest}
:::: {.columns}
::: {.column}
```{r}
#| label: missings-simple
#| echo: false
```
:::
::: {.column}
```{r}
#| label: count-gaps
count_gaps(harvest, .full=TRUE)
```
One missing in each level, although it is a different year.
Notice how `tsibble` handles this summary so neatly.
:::
::::
## Checking and filling missings [(3/4)]{.smallest}
:::: {.columns}
::: {.column}
```{r}
#| label: missings-simple
#| echo: false
```
:::
::: {.column}
Make the implicit missing values [explicit]{.monash-blue2}.
```{r}
#| label: fill-gaps
harvest <- fill_gaps(harvest,
.full=TRUE)
harvest
```
:::
::::
## Checking and filling missings [(4/4)]{.smallest}
:::: {.columns}
::: {.column}
```{r}
#| label: missings-simple
#| echo: false
```
:::
::: {.column}
We have already seen `na_ma()` function, that imputes using a moving average. Alternatively, `na_interpolation()` uses the previous and next values to impute.
```{r}
#| label: impute-gaps
harvest_nomiss <- harvest |>
group_by(fruit) |>
mutate(kilo =
na_interpolation(kilo)) |>
ungroup()
harvest_nomiss
```
:::
::::
## Case study: pedestrian sensors [(1/3)]{.smallest}
```{r}
#| fig-width: 15
#| fig-height: 10
#| out-width: 100%
#| code-fold: true
ped_QV <- pedestrian |>
mutate(year = year(Date),
month = month(Date)) |>
dplyr::filter(Sensor == "QV Market-Elizabeth St (West)",
year == 2016, month %in% c(2:5))
ggplot(ped_QV) +
geom_line(aes(x=Time, y=Count)) +
facet_calendar(~Date) +
theme(axis.title = element_blank(),
axis.text.y = element_blank(),
aspect.ratio = 0.5)
```
## Case study: pedestrian sensors [(2/3)]{.smallest}
```{r echo=TRUE}
#| code-fold: true
has_gaps(pedestrian, .full = TRUE)
ped_gaps <- pedestrian |>
dplyr::filter(Sensor == "QV Market-Elizabeth St (West)") |>
count_gaps(.full = TRUE)
ped_gaps
ped_qvm_filled <- pedestrian |>
dplyr::filter(Sensor == "QV Market-Elizabeth St (West)") |>
fill_gaps(.full = TRUE) |>
mutate(Date = as.Date(Date_Time, tz="Australia/Melbourne"),
Time = hour(Date_Time)) |>
mutate(year = year(Date),
month = month(Date))
```
There is a block of missing values on the last day of 2015. The other two missings correspond to an hour in April, when the clocks change from summer to regular time.
## Case study: pedestrian sensors [(3/3)]{.smallest}
```{r}
#| fig-width: 15
#| fig-height: 3
#| out-width: 100%
#| code-fold: true
# Which days are holidays
hol <- holiday_aus(2015:2016, state = "VIC")
ped_qvm_filled <- ped_qvm_filled |>
mutate(hol = is.weekend(Date)) |>
mutate(hol = if_else(Date %in% hol, TRUE, hol)) |>
mutate(Time = factor(Time))
# Fit linear model to inpute missings
ped_qvm_lm <- lm(Count~Time*hol, data=ped_qvm_filled)
ped_qvm_filled <- ped_qvm_filled |>
mutate(pCount = predict(ped_qvm_lm, ped_qvm_filled))
# Select subset and plot results
ped_qvm_sub <- ped_qvm_filled |>
filter(Date > ymd("2015-12-17"), Date < ymd("2016-01-01"))
ggplot(ped_qvm_sub) +
geom_line(aes(x=Date_Time, y=Count)) +
geom_line(data=filter(ped_qvm_sub, is.na(Count)),
aes(x=Date_Time,
y=pCount),
colour="seagreen3") +
scale_x_datetime("", date_breaks = "1 day",
date_labels = "%a %d") +
theme(aspect.ratio = 0.15)
```
```{r}
#| fig-width: 5
#| fig-height: 5
#| out-width: 50%
#| echo: false
#| eval: false
ggplot() +
geom_line(data=filter(ped_qvm_sub),
aes(x=Time, y=Count, group=Date)) +
geom_line(data=filter(ped_qvm_sub, is.na(Count)),
aes(x=Time,
y=pCount,
group=Date),
colour="seagreen3",
linewidth=2) +
facet_wrap(~hol) +
theme(aspect.ratio = 1)
```
## {.information-slide .center style="text-align: center;"}
The resources at [Forecasting: Principles and Practice (3e)](https://otexts.com/fpp3/) explains how to effectively do forecasting.
## Resources
- [tidymodels.org](https://www.tidymodels.org/)
- [Tidy Modeling with R](https://www.tmwr.org/)
- [Tidy Text](https://www.tidytextmining.com/)
- [Interpretable Machine Learning](https://christophm.github.io/interpretable-ml-book/)
- [Forecasting: Principles and Practice (3e)](https://otexts.com/fpp3/)
::: bottom
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
:::