--- 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)