knitr::opts_chunk$set( comment = "#>", collapse = TRUE )
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.
# 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.
This package provides a tidying function poly_melt()
.
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.
Plot a matrix with poly_plot()
.
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.
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.
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.
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()
.
poly_plot_data(poly_mat, by_observation = FALSE)
The ranges of the terms created by poly()
are sensitive to repeated values.
# 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.
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)
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.
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.
poly_add_columns(df, Days, degree = 3, prefix = "poly_", scale_width = 1)
We can confirm that the added columns are orthogonal.
df <- poly_add_columns(df, Days, degree = 3, scale_width = 1) zapsmall(cor(df[c("Days1", "Days2", "Days3")]))
This package also (accidentally) works on splines. Splines are not officially supported, but they could be an avenue for future development.
poly_plot(splines::bs(1:100, 10, intercept = TRUE)) poly_plot(splines::ns(1:100, 10, intercept = FALSE))
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.
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.
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.
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.
weighted <- pred_mat %*% diag(fixef(model)) colnames(weighted) <- colnames(pred_mat)
And do some tidying to plot the two sets of predictors.
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.
If you searched for help on poly()
, see also:
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.