We are really proud of both editions of the book. This book can help an R user directly experience the data science style of working with data and machine learning techniques.
Get a signed copy off us! We will be giving away some e-copies and a few signed physical copies at various conferences and meet-ups
(for example at PyData LA 2019).
Regularization is a way of avoiding overfit by restricting the magnitude of model coefficients (or in deep learning, node weights). A simple example of regularization is the use of ridge or lasso regression to fit linear models in the presence of collinear variables or (quasi-)separation. The intuition is that smaller coefficients are less sensitive to idiosyncracies in the training data, and hence, less likely to overfit.
Cross-validation is a way to safely reuse training data in nested model situations. This includes both the case of setting hyperparameters before fitting a model, and the case of fitting models (let’s call them base learners) that are then used as variables in downstream models, as shown in Figure 1. In either situation, using the same data twice can lead to models that are overtuned to idiosyncracies in the training data, and more likely to overfit.
Figure 1 Properly nesting models with cross-validation
In general, if any stage of your modeling pipeline involves looking at the outcome (we’ll call that a y-aware stage), you cannot directly use the same data in the following stage of the pipeline. If you have enough data, you can use separate data in each stage of the modeling process (for example, one set of data to learn hyperparameters, another set of data to train the model that uses those hyperparameters). Otherwise, you should use cross-validation to reduce the nested model bias.
Cross-validation is relatively computationally expensive; regularization is relatively cheap. Can you mitigate nested model bias by using regularization techniques instead of cross-validation?
The short answer: no, you shouldn’t. But as, we’ve written before, demonstrating this is more memorable than simply saying “Don’t do that.”
A simple example
Suppose you have a system with two categorical variables. The variable x_s has 10 levels, and the variable x_n has 100 levels. The outcome y is a function of x_s, but not of x_n (but you, the analyst building the model, don’t know this). Here’s the head of the data.
With most modeling techniques, a categorical variable with K levels is equivalent to K or K-1 numerical (indicator or dummy) variables, so this system actually has around 110 variables. In real life situations where a data scientist is working with high-cardinality categorical variables, or with a lot of categorical variables, the number of actual variables can begin to swamp the size of training data, and/or bog down the machine learning algorithm.
One way to deal with these issues is to represent each categorical variable by a single variable model (or base learner), and then use the predictions of those base learners as the inputs to a bigger model. So instead of fitting a model with 110 indicator variables, you can fit a model with two numerical variables. This is a simple example of nested models.
Figure 2 Impact coding as an example of nested model
We refer to this procedure as “impact coding,” and it is one of the data treatments available in the vtreat package, specifically for dealing with high-cardinality categorical variables. But for now, let’s go back to the original problem.
The naive way
For this simple example, you might try representing each variable as the expected value of y - mean(y) in the training data, conditioned on the variable’s level. So the ith “coefficient” of the one-variable model would be given by:
vi = E[y|x = si] − E[y]
Where si is the ith level. Let’s show this with the variable x_s (the code for all the examples in this article is here):
In other words, whenever the value of x_s is s_01, the one variable model vs returns the value 0.8503282, and so on. If you do this for both variables, you get a training set that looks like this:
Now fit a linear model for y as a function of vs and vn.
model_raw = lm(y ~ vs + vn,
data=dtrain_treated)
summary(model_raw)
##
## Call:
## lm(formula = y ~ vs + vn, data = dtrain_treated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.33068 -0.57106 0.00342 0.52488 2.25472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05050 0.05597 -0.902 0.368
## vs 0.77259 0.05940 13.006 <2e-16 ***
## vn 0.61201 0.06906 8.862 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8761 on 242 degrees of freedom
## Multiple R-squared: 0.6382, Adjusted R-squared: 0.6352
## F-statistic: 213.5 on 2 and 242 DF, p-value: < 2.2e-16
Note that this model gives significant coefficients to both vs and vn, even though y is not a function of x_n (or vn). Because you used the same data to fit the one variable base learners and to fit the larger model, you have overfit.
The right way: cross-validation
The correct way to impact code (or to nest models in general) is to use cross-validation techniques. Impact coding with cross-validation is already implemented in vtreat; note the similarity between this diagram and Figure 1 above.
Figure 3 Cross-validated data preparation with vtreat
The training data is used both to fit the base learners (as we did above) and to also to create a data frame of cross-validated base learner predictions (called a cross-frame in vtreat). This cross-frame is used to train the overall model. Let’s fit the correct nested model, using vtreat.
library(vtreat)
library(wrapr)
xframeResults = mkCrossFrameNExperiment(dtrain,
qc(x_s, x_n), "y",
codeRestriction =qc(catN),
verbose =FALSE)
# the plan uses the one-variable models to treat data
treatmentPlan = xframeResults$treatments
# the cross-frame
dtrain_treated = xframeResults$crossFrame
head(dtrain_treated)
##
## Call:
## lm(formula = mk_formula("y", variables), data = dtrain_treated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2157 -0.7343 0.0225 0.7483 2.9639
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.04169 0.06745 -0.618 0.537
## x_s_catN 0.92968 0.06344 14.656 <2e-16 ***
## x_n_catN 0.10204 0.06654 1.533 0.126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.055 on 242 degrees of freedom
## Multiple R-squared: 0.4753, Adjusted R-squared: 0.471
## F-statistic: 109.6 on 2 and 242 DF, p-value: < 2.2e-16
This model correctly determines that x_n (and its one-variable model x_n_catN) do not affect the outcome. We can compare the performance of this model to the naive model on holdout data.
rmse
rsquared
ypred_naive
1.303778
0.2311538
ypred_crossval
1.093955
0.4587089
The correct model has a much smaller root-mean-squared error and a much larger R-squared than the naive model when applied to new data.
An attempted alternative: regularized models.
But cross-validation is so complicated. Can’t we just regularize? As we’ll show in the appendix of this article, for a one-variable model, L2-regularization is simply Laplace smoothing. Again, we’ll represent each “coefficient” of the one-variable model as the Laplace smoothed value minus the grand mean.
vi = ∑xj = siyi/(counti + λ) − E[yi]
Where counti is the frequency of si in the training data, and λ is the smoothing parameter (usually 1). If λ = 1 then the first term on the right is just adding one to the frequency of the level and then taking the “adjusted conditional mean” of y.
##
## Call:
## lm(formula = y ~ vs + vn, data = dtrain_treated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.30354 -0.57688 -0.02224 0.56799 2.25723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.06665 0.05637 -1.182 0.238
## vs 0.81142 0.06203 13.082 < 2e-16 ***
## vn 0.85393 0.09905 8.621 8.8e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8819 on 242 degrees of freedom
## Multiple R-squared: 0.6334, Adjusted R-squared: 0.6304
## F-statistic: 209.1 on 2 and 242 DF, p-value: < 2.2e-16
Again, both variables look significant. Even with regularization, the model is still overfit. Comparing the performance of the models on holdout data, you see that the regularized model does a little better than the naive model, but not as well as the correctly cross-validated model.
rmse
rsquared
ypred_naive
1.303778
0.2311538
ypred_crossval
1.093955
0.4587089
ypred_reg
1.267648
0.2731756
The Moral of the Story
Unfortunately, regularization is not enough to overcome nested model bias. Whenever you apply a y-aware process to your data, you have to use cross-validation methods (or a separate data set) at the next stage of your modeling pipeline.
Appendix: Derivation of Laplace Smoothing as L2-Regularization
Without regularization, the optimal one-variable model for y in terms of a categorical variable with K levels {sj} is a set of K coefficients v such that
is minimized (N is the number of data points). L2-regularization adds a penalty to the magnitude of v, so that the goal is to minimize
where λ is a known smoothing hyperparameter, usually set (in this case) to 1.
To minimize the above expression for a single coefficient vj, take the deriviative with respect to vj and set it to zero:
Where countj is the number of times the level sj appears in the training data. Now solve for vj:
This is Laplace smoothing. Note that it is also the one-variable equivalent of ridge regression.
rquery is a data wrangling system designed to express complex data manipulation as a series of simple data transforms. This is in the spirit of R’s base::transform(), or dplyr’s dplyr::mutate() and uses a pipe in the style popularized in R with magrittr. The operators themselves follow the selections in Codd’s relational algebra, with the addition of the traditional SQL “window functions.” More on the background and context of rquery can be found here.
The R/rquery version of this introduction is here, and the Python/data_algebra version of this introduction is here.
In transform formulations data manipulation is written as transformations that produce new data.frames, instead of as alterations of a primary data structure (as is the case with data.table). Transform system can use more space and time than in-place methods. However, in our opinion, transform systems have a number of pedagogical advantages.
In rquery’s case the primary set of data operators is as follows:
These operations break into a small number of themes:
Simple column operations (selecting and re-naming columns).
Simple row operations (selecting and re-ordering rows).
Creating new columns or replacing columns with new calculated values.
Aggregating or summarizing data.
Combining results between two data.frames.
General conversion of record layouts (supplied by the cdata package).
The point is: Codd worked out that a great number of data transformations can be decomposed into a small number of the above steps. rquery supplies a high performance implementation of these methods that scales from in-memory scale up through big data scale (to just about anything that supplies a sufficiently powerful SQL interface, such as PostgreSQL, Apache Spark, or Google BigQuery).
We will work through simple examples/demonstrations of the rquery data manipulation operators.
This section reflects an important design decision in the book: teach model evaluation first, and as a step separate from model construction.
It is funny, but it takes some effort to teach in this way. New data scientists want to dive into the details of model construction first, and statisticians are used to getting model diagnostics as a side-effect of model fitting. However, to compare different modeling approaches one really needs good model evaluation that is independent of the model construction techniques.
This teaching style has worked very well for us both in R and in Python (it is considered one of the merits of our LinkedIn AI Academy course design):
Nina Zumel finished new documentation on how vtreat‘s cross validation works, which I want to share here.
vtreat is a system that makes data preparation for machine learning a “one-liner” (available in R or available in Python). We have a set of starting off points here. These documents describe what vtreat does for you, you just find the one that matches your task and you should have a good start for solving data science problems in R or in Python.
The latest documentation is a bit about how vtreat works, and how to control some of the details of the work it is doing for you.
Nina Zumel finished some great new documentation showing how to use Pythonvtreat to prepare data for multinomial classification mode. And I have finally finished porting the documentation to Rvtreat. So we now have good introductions on how to use vtreat to prepare data for the common tasks of:
That is now 8 introductions to start with. To use vtreat you only have to work through one introduction (the one helping with the task you have at hand in the language you are using).
As I have said before:
vtreat helps with project blocking issues commonly seen in real world data: missing values, re-coding categorical variables, and dealing high cardinality categorical variables.
If you aren’t using a tool like vtreat in your data science projects: you are really missing out (and making more work for yourself).
When studying regression models, One of the first diagnostic plots most students learn is to plot residuals versus the model’s predictions (that is, with the predictions on the x-axis). Here’s a basic example.
# build an "ideal" linear process.
set.seed(34524)
N = 100
x1 = runif(N)
x2 = runif(N)
noise = 0.25*rnorm(N)
y = x1 + x2 + noise
df = data.frame(x1=x1, x2=x2, y=y)
# Fit a linear regression model
model = lm(y~x1+x2, data=df)
summary(model)
##
## Call:
## lm(formula = y ~ x1 + x2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73508 -0.16632 0.02228 0.19501 0.55190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.16706 0.07111 2.349 0.0208 *
## x1 0.90047 0.09435 9.544 1.30e-15 ***
## x2 0.81444 0.09288 8.769 6.07e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2662 on 97 degrees of freedom
## Multiple R-squared: 0.6248, Adjusted R-squared: 0.6171
## F-statistic: 80.78 on 2 and 97 DF, p-value: < 2.2e-16
In the above plot, we’re plotting the residuals as a function of model prediction, and comparing them to the line y = 0, using a smoothing curve through the residuals. The idea is that for a well-fit model, the smoothing curve should approximately lie on the line y = 0. This is true not only for linear models, but for any model that captures most of the explainable variance, and for which the unexplainable variance (the noise) is IID and zero mean.
If the residuals aren’t zero mean independently of the model’s predictions, then either you are missing some explanatory variables, or your model does not have the correct structure, or an appropriate inductive bias. A simple example of the second case is trying to fit a linear model to a process where the outcome is quadratically (or otherwise non-linearly) related to the outcome. To see this, let’s make an example quadratic system while deliberately failing to supply that structure to the model.
In this case, the smoothing line on the residuals doesn’t approximate the line y = 0; when the model predicts a value in the range 1 to about 2.3, it tends to be overpredicting; otherwise, it tends to underpredict. This is an instance of a pathology called “structure in the residuals.”
The Peril of Outcomes on the x-axis
What happens if you erroneously plot the residuals versus the true outcome, instead of the predictions? Let’s try this with the model for the linear process (which we know is a well-fit model):
# the wrong residual graph
ggplot(df, aes(x=y, y=residual)) +
geom_point(alpha=0.5) + geom_hline(yintercept=0, color="red") +
geom_smooth(se=FALSE) +
ggtitle("Incorrect residual plot",
subtitle = "linear model and process")
If you make this plot when you meant to make the other, you will give yourself a nasty shock. Plotting residuals versus the outcome will always look more or less like the above graph. You might think that for a good model, the outcome and the prediction are close to each other, so the residual graphs should look about the same no matter which quantity you plot on the x-axis, right? Why do they look so different?
Reversion to mediocrity (or the mean)
One reason that the proper residual graph (for a well fit model) should smooth out to the line y=0 is known as reversion to mediocrity, or regression to the mean.
Imagine that you have an ideal process that always produces a single value y. You don’t actually observe this “true value”; instead, what you observe is y plus (IID, zero mean) noise. You can build a “model” for this process that predicts the mean of the observations, in this case the value 0.1033149. Then you can calculate the residuals of your “model” in the usual way.
When you plot the residuals as a function of the prediction, all the datums fall at the same horizontal coordinate of the graph, centered around zero, and approximately equally distributed between positive and negative. The “smoothing line” through this graph is simply the point (0.1033149, 0) – that is, the graph is centered at zero.
On the other hand, if you plot the residuals as a function of the observed outcome, all the observations will be sorted so that the observations with positive noise are to the right of the observations with negative noise, and the smoothing line through the graph no longer looks like the line y = 0.
For a process that varies as a function of the input, you can think of the prediction corresponding to an input X as the mean of all the observations corresponding to X, and the idea is the same.
Incidentally, this regression to the mean is also why model predictions tend to have less range than the original training data.
Plotting observations versus predictions
Sometimes instead of plotting residuals versus the predictions, I plot observations versus predictions. In this case, you want to check that the predictions lie approximately on the line y = x. This isn’t a standard diagnostic plot, but it does give a better sense of the magnitude of the errors relative to the magnitudes of the outcomes. Again, the important thing to remember is that the predictions go on the x-axis.
# the "wrong" way
ggplot(df, aes(x=y, y=pred)) +
geom_point(alpha=0.5) + geom_abline(color="red") +
geom_smooth(se=FALSE) +
ggtitle("Incorrect prediction plot")
Notice how the wrong plot again seems to show pathological structure where none exists.
Conclusion
The above examples show why you should always take care to plot your model diagnostics as functions of the predictions and not of the observations. Most students have heard this already, but we feel that demonstrating why will be more memorable that simply saying “make it so.”
