---
title: "polypoly package overview"
author: "Tristan Mahr"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: TRUE
vignette: >
  %\VignetteIndexEntry{polypoly package overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  comment = "#>",
  collapse = TRUE
)
```


## Background

The `poly()` function in the `stats` package creates a matrix of (orthogonal)
polynomials over a set of values. The code below shows some examples of these
matrices.

```{r}
# orthogonal polynomials
m <- poly(1:6, degree = 3, simple = TRUE)
m

# the terms are uncorrelated. that's why they are "orthogonal".
zapsmall(cor(m))

# raw polynomials
m <- poly(1:6, degree = 3, simple = TRUE, raw = TRUE)
m

# raw polynomials are highly correlated.
round(cor(m), 2)
```

This package provides some helpful functions for working with these matrices.

## Examples

### Tidying

This package provides a tidying function `poly_melt()`. 

```{r}
library(polypoly)
xs <- 1:40
poly_mat <- poly(xs, degree = 5)
poly_melt(poly_mat)
```

The returned dataframe has one row per cell of the original matrix. Essentialy, 
the columns of the matrix are stacked on top of each other to create a long 
dataframe. The `observation` and `degree` columns record each values' original
row number and column name, respectively.

### Plotting

Plot a matrix with `poly_plot()`.

```{r example, fig.width = 5, fig.height = 3}
poly_plot(poly_mat)
```

We can also plot raw polynomials, but that display is less useful because the
_x_-axis corresponds to the row number of polynomial matrix.

```{r raw-example, fig.width = 5, fig.height = 3}
poly_raw_mat <- poly(-10:10, degree = 3, raw = TRUE)
poly_plot(poly_raw_mat)
```

We can make the units clearer by using `by_observation = FALSE` so that the
_x_-axis corresponds to the first column of the polynomial matrix.

```{r raw-by-degree1,  fig.width = 5, fig.height = 3}
poly_plot(poly_raw_mat, by_observation = FALSE)
```

`poly_plot()` returns a plain ggplot2 plot, so we can further customize the
output. For example, we can use ggplot2 to compute the sum of the individual
polynomials and re-theme the plot.

```{r sum, fig.width = 5, fig.height = 3}
library(ggplot2)
poly_plot(poly_mat) + 
  stat_summary(
    aes(color = "sum"), 
    fun = "sum", 
    geom = "line", 
    size = 1
  ) + 
  theme_minimal()
```

For total customization, `poly_plot_data()` will return the dataframe that
_would_ have been plotted by `poly_plot()`.

```{r}
poly_plot_data(poly_mat, by_observation = FALSE)
```


### Rescaling a matrix

The ranges of the terms created by `poly()` are sensitive to repeated values.

```{r}
# For each column in a matrix, return difference between max and min values
col_range <- function(matrix) {
  apply(matrix, 2, function(xs) max(xs) - min(xs))  
}

p1 <- poly(0:9, degree = 2)
p2 <- poly(rep(0:9, 18), degree = 2)

col_range(p1)
col_range(p2)
```

Thus, two models fit with `y ~ poly(x, 3)` will not have comparable coefficients
when the number of rows changes, even if the unique values of `x` did not
change!

`poly_rescale()` adjusts the values in the polynomial matrix so that the linear
component has a specified range. The other terms are scaled by the same factor. 

```{r rescaled,  fig.width = 5, fig.height = 3}
col_range(poly_rescale(p1, scale_width = 1))
col_range(poly_rescale(p2, scale_width = 1))

poly_plot(poly_rescale(p2, scale_width = 1), by_observation = FALSE)
```



### Adding columns to a dataframe

`poly_add_columns()` adds orthogonal polynomial transformations of a predictor
variable to a dataframe.

Here's how we could add polynomials to the `sleepstudy` dataset.

```{r}
df <- tibble::as_tibble(lme4::sleepstudy)
print(df)

poly_add_columns(df, Days, degree = 3)
```

We can optionally customize the column names and rescale the polynomial terms.

```{r}
poly_add_columns(df, Days, degree = 3, prefix = "poly_", scale_width = 1)
```

We can confirm that the added columns are orthogonal.

```{r sleepstudy, fig.width = 5, fig.height = 3}
df <- poly_add_columns(df, Days, degree = 3, scale_width = 1)
zapsmall(cor(df[c("Days1", "Days2", "Days3")]))
```


### Experimental

#### Splines

This package also (accidentally) works on splines. Splines are not officially
supported, but they could be an avenue for future development.

```{r splines,  fig.width = 5, fig.height = 3}
poly_plot(splines::bs(1:100, 10, intercept = TRUE))
poly_plot(splines::ns(1:100, 10, intercept = FALSE))
```

## Longer example: Visualizing growth curve contributions

This section illustrates a use case that may or may not be included in the 
package someday: Visualizing the weighting of polynomial terms from a linear
model. For now, here's how to do that task with this package.

Suppose we want to model some change over time using a cubic polynomial. For
example, the growth of trees.

```{r trees1, fig.width = 5, fig.height = 3}
library(lme4)
df <- tibble::as_tibble(Orange)
df$Tree <- as.character(df$Tree)
df

ggplot(df) + 
  aes(x = age, y = circumference, color = Tree) + 
  geom_line()
```

We can bind the polynomial terms onto the data and fit a model.

```{r}
df <- poly_add_columns(Orange, age, 3, scale_width = 1)

model <- lmer(
  scale(circumference) ~ age1 + age2 + age3 + (age1 + age2 + age3 | Tree), 
  data = df
)
summary(model)
```

How do we understand the contribution of each of these terms? We can recreate
the model matrix by attaching the intercept term to a polynomial matrix.

```{r}
poly_mat <- poly_rescale(poly(df$age, degree = 3), 1)

# Keep only seven rows because there are 7 observations per tree
poly_mat <- poly_mat[1:7, ]

pred_mat <- cbind(constant = 1, poly_mat)
pred_mat
```

Weight the predictors using the model fixed effects.

```{r}
weighted <- pred_mat %*% diag(fixef(model))
colnames(weighted) <- colnames(pred_mat)
```

And do some tidying to plot the two sets of predictors.

```{r trees-comparison, fig.width = 7, fig.height = 3}
df_raw <- poly_melt(pred_mat)
df_raw$predictors <- "raw"

df_weighted <- poly_melt(weighted)
df_weighted$predictors <- "weighted"

df_both <- rbind(df_raw, df_weighted)

# Only need the first 7 observations because that is one tree
ggplot(df_both[df_both$observation <= 7, ]) + 
  aes(x = observation, y = value, color = degree) + 
  geom_line() + 
  facet_grid(. ~ predictors) + 
  labs(color = "term")
```

The linear trend drives the growth curve. The quadratic and cubic terms make 
tiny contributions. We can see that the intercept term does nothing (because we 
used `scale()` in the model). 

Hmmm... perhaps we need to find a better example dataset for this example.




## Resources

If you searched for help on `poly()`, see also:

* [What does the R function `poly()` really do?](https://stackoverflow.com/questions/19484053/what-does-the-r-function-poly-really-do)
* [Source code for the poly function](https://github.com/wch/r-source/blob/af7f52f70101960861e5d995d3a4bec010bc89e6/src/library/stats/R/contr.poly.R#L85)