Linear Regression

1 Linear Regression

Linear Regression is the workhorse of applied Data Science; it has long been the most commonly used method by scientists and can be applied to a wide variety of datasets and questions. Unlike more recently developed methods in Machine Learning, Linear Regression models can be used to predict new data points and help us understand the relative impact one variable has on another. For example, one well-designed regression model can answer:

• How, and to what extent, does advertising in print media effect the total sales of a product?
• What are the predicted total sales for a product, given the amount spent on print advertisement this month?

how to harness the malleability and explanatory power of Linear Regression models by following the four primary steps of statistical model building: confirming data assumptions, building a model on training data, assessing model fit, and analyzing model results. Using two real-life datasets, conversion and advertising, this lesson will also focus on the application of regression modeling for use in marketing data science.

conversion is an example of data that a scientist would receive from a large social media platform, like Facebook or Twitter, and can be used to answer questions about the performance of newsfeed advertisements. advertising is an example of data that would be recorded by the marketing department of a company; a marketing data scientist would often be given this data and expected to answer questions around of the return on investment for various media channels, like podcast or tv advertisements.

# load libraries and data
library(dplyr)
library(readr)
library(ggplot2)
# load data
conversion <- read.csv('conversion.csv', header= T)
advertising <- read.csv('advertising.csv', header= T)

1.1 Exploring the datasets

A good statistical workflow always involves and thorough understanding of the data available to model and a qualitative analysis of relevant variables. Use str() to write out the structure of the dataset and list of variables types. Which variables seem like possible predictors of purchase, or total_convert?

str(conversion)

'data.frame':   1138 obs. of  11 variables:
 $ ad_id           : int  708746 708749 708771 708815 708818 708820 708889 708895 708953 708958 ...
 $ campaign_id     : int  916 916 916 916 916 916 916 916 916 916 ...
 $ fb_campaign_id  : int  103916 103917 103920 103928 103928 103929 103940 103941 103951 103952 ...
 $ age             : chr  "30-34" "30-34" "30-34" "30-34" ...
 $ gender          : chr  "M" "M" "M" "M" ...
 $ interest        : int  15 16 20 28 28 29 15 16 27 28 ...
 $ impressions     : int  7350 17861 693 4259 4133 1915 15615 10951 2355 9502 ...
 $ clicks          : int  1 2 0 1 1 0 3 1 1 3 ...
 $ spent           : num  1.43 1.82 0 1.25 1.29 ...
 $ total_convert   : int  2 2 1 1 1 1 1 1 1 1 ...
 $ approved_convert: int  1 0 0 0 1 1 0 1 0 0 ...

str(advertising)

'data.frame':   199 obs. of  4 variables:
 $ tv       : num  44.5 17.2 151.5 180.8 8.7 ...
 $ podcast  : num  39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 5.8 ...
 $ newspaper: num  45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 24.2 ...
 $ sales    : num  10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 8.6 ...

Mean of sales:

mean(advertising$sales)

[1] 13.98191

ggplot(advertising, aes(x = tv, y = sales)) +
  geom_point()

ggplot(advertising, aes(x = podcast, y = sales)) +
  geom_point()

2 Simple Linear Regression

While the linear regression is perhaps the most widely applied method in Data Science, it relies on a strict set of assumptions about the relationship between predictor and outcome variables. The most obvious (but crucial!) assumption is a linear relationship between the predictor and outcome. Following from this assumption is one key observation about any variables we want to include in our model which must be tested before building a model:

The expected value of the outcome variable is be a straight-line function of exclusively the predictor variable. The best test for this relationship is quite straightforward–– we can just visualize the relationship between the predictor and outcome variables as a scatterplot. A linear relationship will resemble a straight line with a slope not equal to zero, like the relationship between spending on TV ads and the overall sales volume of the related product found in our advertising dataset.

ggplot(advertising, aes(x = tv, y = sales)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE)

We can also quantitatively test for a linear relationship by computing the correlation coefficient. The correlation coefficient is always between positive one and negative one. A coefficient close to 0 (roughly between -0.20 and 0.20) suggests a weak linear relationship between two variables. A coefficient closer to positive or negative one suggests a stronger linear relationship. In R, we can compute the correlation coefficient using the cor.test() method as follows:

coefficient <- cor.test(advertising$tv, advertising$sales)
coefficient$estimate

      cor 
0.7812404 

Use a combination of the base ggplot() function and geom_bar() to plot the distribution of the clicks variable, a measure of how many times a user clicked on an advertisement. Save the result to a variable called clicks_dist. Call clicks_dist.

# save viz to object
clicks_dist <- conversion %>% ggplot(aes(x = clicks)) + geom_bar()
clicks_dist

Take a closer look at the clicks_dist visualization. What is the approximate range of the clicks variable? What seems like the most common value (otherwise called the mode) of clicks? Set clicks_mode equal to approximate value of the clicks mode.

clicks_mode <- 1

Assign the result of calling cor.test(), with conversion$total_convert and conversion$clicks as input parameters, to a variable called correlation. Print out correlation$estimate. Does the coefficient value suggest that the variables have a linear relationship?

correlation <- cor.test(conversion$total_convert, conversion$clicks)
correlation$estimate

      cor 
0.6814551 

Yes, seems to have positive linear correlation.

2.1 Assumptions of Linear Regression (Outliers)

Our next step is to check for outlier data points. Linear regression models also assume that there are no extreme values in the data set that are not representative of the actual relationship between predictor and outcome variables. A box-and-whisker plot is a common method used to quickly determine whether a data set contains any outliers, or data points that differ significantly from other observations in a dataset. An outlier may be caused by variability in measurement, or it might be a sign of an error in the collection of data.

Regardless, ggplot’s geom_boxplot() method allows for the easy creation of box-and-whisker plots. To plot the distribution of a single variable, like advertising$sales––the total number of sales for a product in a month–– we pass in the same variable as both x and y in our call to geom:

plot <- advertising %>%
  ggplot(aes(sales, tv)) +
  geom_boxplot()
plot

If we have negative sales values in the dataset, this is not what we would expect given our understanding of the data; how could an entire market have negative average sales over an entire year? This seems like an error stemming from the collection of this data into a spreadsheet format. In this case, we will filter out these negative datapoints from our dataset using the filter() method. We can pass a boolean argument into filter() to exclude values that resolve to false.

advertising <- advertising %>% filter(sales > 0)

Let’s check our variables for any outliers. Use a combination of the base ggplot() function and geom_boxplot() to plot a boxplot of the clicks variable; assign the result to a variable called clicks_bx_plot. Don’t forget to pass clicks in as the y variable, and afterwards call clicks_bx_plot in order to view the plot!

# create viz object
clicks_bx_plot <- ggplot(conversion, aes(x = clicks, y = clicks)) +
geom_boxplot()

# print out object 
clicks_bx_plot

Do any data points look like outliers? Set threshold equal to the value of clicks above which all data points fall outside of the whiskers of our box plot.

# set threshold value
threshold <- 100

Use the filter() method to remove all outlying values from clicks. Save the resulting data frame to a variable called convert_clean.

# remove outliers
convert_clean <- conversion %>%
  filter(clicks < threshold)

Let’s check our work! Create another box plot of clicks, but this time, use the convert_clean dataset and save the plot to clean_bx_plot. Call clean_bx_plot and note any differences from clicks_bx_plot. Are the outliers gone?

# creat second box plot 
clean_bx_plot <- convert_clean %>%
  ggplot(aes(clicks, clicks)) +
  geom_boxplot()

clean_bx_plot

2.2 Building a Simple Linear Model

Simple linear regression is not a misnomer–– it is an uncomplicated technique for predicting a continuous outcome variable, Y, on the basis of just one predictor variable, X. As detailed in previous exercises, a number of assumptions are made so that we can model the relationship between X and Y as a linear function. Using our advertising dataset, we could model the relationship between the amount spent on podcast advertising in a month and the number of respective products eventually sold as follows:

$ Y= beta_0 + beta_1∗X + error$

Where…

Y: represents the dollar value of products sold

X: represents the amount spent on respective product podcast ads

Beta_0: is the intercept, or the number of products sold when no money has been spent on podcasts

Beta_1: is the coefficient, or the slope, of the line representing the relationship

Error: represents the random variation in the relationship between the two variables

To build this model in R, using the standard lm() package, we use the formula notation of Y ~ X:

model <- lm(sales ~ podcast, data = train)

But wait! Before building this model, we need to split our data into test and training sets For the development of this simple model, we’ll use a standard 60/40 split of our data; where 60% is used to train the model, and 40% is used to test the model’s accuracy and generalizability. We can randomly assign data points to test or training using base R’s sample() method and list indexing functionality

# specify 60/40 split
sample <- sample(c(TRUE, FALSE), nrow(advertising), replace = T, prob =     c(0.6,0.4))
# subset data points into train and test sets
train <- advertising[sample, ]
test <- advertising[!sample, ]

First, let’s split our conversion_clean dataset into 60/40 train/test subsets

Create a variable named data_sample by assigning the result of calling sample(), with c(TRUE, FALSE), nrow(conversion_clean), and prob = c(0.6,0.4) as parameters.

• Using list indexing, assign all data points in sample to a variable called train
• Using list indexing, assign all data points in not in sample to a variable called test

# set sampling seed
set.seed(123)
# specify 60/40 split
data_sample <- sample(c(TRUE, FALSE), nrow(convert_clean), replace = T, prob = c(0.6, 0.4))

# subset data points into train and test sets
train <- convert_clean[data_sample, ]
test <- convert_clean[!data_sample, ]

Let’s fit a linear model of the relationship between the number of products sold and the number of clicks on the respective product advertisement; this means that conversion’s total_convert value, the total number of product purchases by a single user, will be our Y variable, or outcome variable. clicks the total number of times a user clicks on a version of an ad, will be our X variable, or predictor variable.

Assign the result of calling lm(), using ~ formula notation to set a linear relationship between total_conversions and clicks, to a variable called model. Don’t forget to set the data parameter equal to train!

# build model
model <- lm(total_convert ~ clicks, data = train)
model

Call:
lm(formula = total_convert ~ clicks, data = train)

Coefficients:
(Intercept)       clicks  
    1.04046      0.05037  

2.3 Quantifying Model Fit

Once we have an understanding of the kind of relationship our model describes, we want to understand the extent to which this modeled relationship actually fits the data. This is typically referred to as the goodness-of-fit. In simple linear models, we can measure this quantitatively by assessing two things:

1. Residual standard error (RSE)

2. R squared (R^2) The RSE is an estimate of the standard deviation of the error of the model (error in our mathematical definition of linear regression). Roughly speaking, it is the average amount that the response will deviate from the true regression line. We get the RSE at the bottom of summary(model), we can also get it directly with

   sigma(model)

   #output 3.2

An RSE value of 3.2 means the actual sales in each market will deviate from the true regression line by approximately 3,200 units, on average. Is this too large of a deviation?

Well, that’s subjective, but when compared to the average value of sales over all markets the percentage error is 22%:

sigma(model)/mean(train$sales)

# output
[1] 0.2207373

The RSE provides an absolute measure of lack of fit of our model to the data. But since it is measured in the units of Y, it is not always clear what constitutes a good RSE.

The R^2 statistic provides an alternative measure of fit. It represents the proportion of variance explained, so it always takes on a value between 0 and 1, and is independent of the scale of Y, our outcome variable. Similar to RSE, the R^2 can be found at the bottom of summary(model) but we can also extract it directly by calling summary(model)$r.squared. The result below suggests that podcast advertising budget can explain 64% of the variability in the total sales value.

summary(model)$r.squared

# output
[1] 0.6372581

# compute avg_rse
avg_rse = sigma(model)/mean(train$total_convert)
#uncomment f-string below
sprintf("The percentage error of the model is %s . Any prediction drawn from this model could be %s percent off from the actual observed value.", avg_rse, avg_rse)

[1] "The percentage error of the model is 0.803572334607046 . Any prediction drawn from this model could be 0.803572334607046 percent off from the actual observed value."

Model fit is often quantified in comparison to other models, then used to determine which variation of a modeled relationship best fits the data. Let’s build a second model so that we can contextualize our fit metrics.

Assign the result of building a simple linear model of impressions, the total number of times a user views a version of an advertisement, on total_convert to the variable model_2.

model_2 <- lm(total_convert ~ impressions, data = train)

Let’s use a combination of R’s variable selection syntax, the $ character, and summary() to investigate the percent of variability explained by both model and model2.

Call extract the r-squared measure from model, and save the result to a variable called r_sq.

# compute r-squared
r_sq = summary(model)$r.squared

Call extract the r-squared measure from model, and save the result to a variable called r_sq_2.

r_sq_2 = summary(model_2)$r.squared

Print out both r-square variables. Which model better explains a user’s likelihood of purchasing a product they have been shown an advertisement for?

# print out r-squared values
r_sq

[1] 0.3716766

r_sq_2

[1] 0.5018523

# uncomment f-string below
 sprintf("Based on a pair of simple linear regression models, we have determined that %s percent of the variation in user purchase behavior can be explained by the number of times a user viewed on a relevant ad campaign; whereas only %s percent of this variation can be explained by the number of times a user clicked on a relevant ad.", r_sq_2, r_sq)

[1] "Based on a pair of simple linear regression models, we have determined that 0.501852326225873 percent of the variation in user purchase behavior can be explained by the number of times a user viewed on a relevant ad campaign; whereas only 0.371676623224315 percent of this variation can be explained by the number of times a user clicked on a relevant ad."

2.4 Checking Model Residuals

Great! We can build a model! But… how do we know if it’s any good? Also, if another data scientist builds a different model using a different independent variable, how can we tell which model is “best”? Even within Statistics, “best” can be a subjective qualifier. However, scientists who use regression models generally agree that the best model is the one that minimizes the distance between a data point and the estimation line drawn by a model. The vertical distance between a datapoint and the line estimated by a regression model is called a residual; residuals and their aggregations are the fundamental units of measures of regression model fit and accuracy.

Because residuals are based on cartesian distances, it often helps to visualize their values. For instance, consider the plot of a simple linear regression alongside its’ training data below. Note one point is 4 units above the regression estimate line; in this example, the residual for that point is 4. Meanwhile, another point is 2 units below the regression estimate line; the residual for that point is -2. A data point is best fit by the model which results in the smallest residual for that point.

When scientists make quantitative arguments for a best fit model, they rely on an aggregation, often the sum or average, of residual values across an entire dataset. While is it easy to be overwhelmed by the variety of measures used to argue that one model is better than the other, it is crucial to realize that all measures are grounded in the simple difference between regression estimate and observed data point. Let’s produce a visualization of our own model of clicks on total_convert to better understand our model residuals.

First, let’s pull out the points that make up the estimate line, and respective residual values, from our model. We can save them to columns back in our train dataset called estimate and residuals in our main train training dataset.

Call predict() on model and save the result to train$estimate.

train$estimate <- predict(model)

Call residuals() on model and save the result to train$residuals

train$residuals <- residuals(model)

lot the values of clicks and total_convert using a combination of ggplot() and geom_point(). Save the result to a variable called plot. Don’t forget to pass in clicks as as X variable, and total_convert as a Y variable. Make sure you’re using train as the dataset.

Call plot to view your visualization.

#create visualization
plot <- ggplot(train, aes(clicks, total_convert)) +
geom_point()

plot

Plot the observed data points in our train dataset by adding another call to geom_point(). Make sure to explicitly pass in estimate as a y value, and set the color parameter equal to “blue”.

plot <- ggplot(train, aes(clicks, total_convert)) +
  geom_point() +
  geom_point(aes(y = estimate), color = "blue") 

plot

Let’s explicitly plot the vertical distance between estimate values and their respective observed data point. Add a call to geom_segment(), passing in xend = clicks and yend = estimate as arguments to aes(). Don’t forget to set color = “gray”

plot <- ggplot(train, aes(clicks, total_convert)) +
  geom_point() +
  geom_point(aes(y = estimate), color = "blue") +
  geom_segment(aes(xend = clicks, yend = estimate), color = "gray")

plot

Finally, we should provide another way to observe the size of residuals. Update our first call to geom_point() by passing in size = abs(residuals) as an argument to aes(). As total_convert increases, how do the value of model residuals change?

plot <- ggplot(train, aes(clicks, total_convert)) +
  geom_point(aes(size = abs(residuals))) +
  geom_point(aes(y = estimate), color = "blue") +
  geom_segment(aes(xend = clicks, yend = estimate), color = "gray")

plot

2.5 Visualizing Model Fit

In addition to the quantitative measures that characterize our model accuracy, it is alway a best practice to produce visual summaries to assess our model. First, we should always visualize our model within our data. For simple linear regression this is quite simple; we can use geom_point() to plot our observed values, and geom_smooth(method = “lm”) to plot our model. In addition, we can include a second call to geom_smooth(), with parameters (se = FALSE, color = “red”). This combination of function calls allows us to compare the linearity of our model, visualized below as the blue line with the 95% confidence interval covering the shaded region, in comparison to a non-linear LOESS smoother visualized in red.

ad_sample <- sample(c(TRUE, FALSE), nrow(advertising), replace = T, prob = c(0.6, 0.4))

# subset data points into train and test sets
advertising[ad_sample, ] %>% 
ggplot(aes(podcast, sales)) +
geom_point() +
geom_smooth(method = "lm") +
geom_smooth(se = FALSE, color = "red") 

LOESS smoothers plot a line based on the weighted value of data points; the line produced by a LOESS smoother is similar to taking a moving average of data points as our x-axis variable increases. The smoother should not be used to predict new values, as it relies heavily on our training data, but it is a helpful tool for visualizing where our linear model diverges from our training data.

Considering the LOESS smoother remains within the confidence interval of our model, we can assume the linear trend fits the essence of this relationship. However, we should note that as the podcast advertising budget gets closer to 0 there is a stronger reduction in sales beyond what the linear trend follows; this means that our model might be less accurate in instances where the podcast budget is very low.

We’ve plotted clicks against total converts. Let’s add a LOESS smoother. Add two calls of geom_smooth() to plot. The first should use the parameter method = “lm”. The second should use the parameters se = FALSE and color = “red”.

# build plot of clicks on total_convert below
plot <- ggplot(train, aes(clicks, total_convert)) +
geom_point() +
geom_smooth(method = "lm") +
geom_smooth(se = FALSE, color = "red")

plot

How closely does the relationship between clicks and conversion follow a linear trend? Set the variable linear_relationship equal to either “a”, “b”, “c”, or “d” depending on the statement that best characterizes the relations:

A. The relationship is less linear when clicks approaches large values.

B. There is a clear divergence from a linear relationship when clicks approaches zero or when clicks approaches infinity.

C. The relationship between clicks and total_conversion is perfectly linear.

D. There is no linear relationship between clicks and total_conversion

linear_relationship = "a"

Let’s extend our linearity analysis to our model2, which describes the relationship between impressions and total_conversion. Add the two calls to geom_smooth() to plot_2 to make a comparison to a LOESS model.

# build plot of impressions on total_convert below
plot_2 <- ggplot(train, aes(impressions, total_convert)) +
geom_point() +
geom_smooth(method = "lm") +
geom_smooth(se = FALSE, color = "red")

plot_2

How closely does the relationship between impressions and conversion follow a linear trend? Set the variable linear_relationship_2 equal to “a”, “b”, “c”, or “d” depending on the statement that best characterizes the relations:

A. The relationship between impressions and total_conversion is perfectly linear.

B. There is a clear divergence from a linear relationship when impressions approaches zero and when impressions is around 500,000.

C. The relationship is less linear when impressions approaches very large values and when impressions is around 500,000.

D. There is no linear relationship between impressions and total_conversion

linear_relationship_2  = "c"

2.6 Reading Model Results

Ready for the real fun? We’ve done our due diligence and confirmed that our data fulfills the assumptions of simple linear regression models; we’ve split our data into test and training subsets, and properly built a model using y ~ x + b notation; we’ve even taken the time to assess the fit of our model using both quantitative and qualitative approaches. Now we can finally analyze the results of our model and discover the relationship between user advertisement clicks and the purchase rate of related products!

We can view the results of a linear regression model in R by calling summary() on the model variable to which we saved the results of our call to lm(). The summary() function will print out a lot of information about our model–– but don’t be overwhelmed! There are four primary sections of quantitative results that are crucial to interpreting regression models:

### Call

This section simple displays the call to lm() which created these model results. It’s a helpful reminder of which version of a outcome-predictor pair is currently under analysis.

2.6.1 Residuals

As covered in our earlier exercises, a residual is the difference between the value of an outcome variable predicted by the model and the actual observed value of the variable. The summary() output displays a set of numbers that summarize the distribution of residuals in our model, including minimum/maximum residual values, the values first/third quantiles, and the median residual value for the model. We’ve already analyzed our residual values by creating a plot in an earlier exercise, but these summary values are a helpful reminder of the overall spread of our model errors.

2.6.2 Coefficients

2.6.2.1 Estimate

Coefficients are most important results in the interpretation of regression models. The number you see in the Estimate column, (a value of 0.048939 for clicks) is called a regression coefficient. Looking back to formal definition of a linear regression model:

$ Y= beta_0 + beta_1∗X + error$

The regression coefficient is represented by the beta_1 variable. This linear regression equation tells us that the regression coefficient represents the expected change in the dependent variable (in our case total_convert) for a one-unit increase in the independent variable (clicks). In other words, for every additional click on an advertisement, the expected sales of a related product are estimated to increase by 0.049 dollars. In addition to the size of the coefficient, it is also important to note the sign of the coefficient. If our clicks coefficient was negative, our model would be estimating that the sales of a product actually decreases every time an advertisement is clicked.

2.6.2.2 Std. Error

The column adjacent to Estimate is called Std. Error; the standard error of each coefficient is the estimate of the standard deviation of the coefficient. It is crucial to note that the standard error is not a quantity of interest by itself, but depends on the value of our regression coefficient

2.6.2.3 T-value and Pr(>|t|)

The t value and Pr(>|t|) inherently answer the same question–– given the value of our variable’s regression coefficient and its’ standard error, does the variable explain a significant part of the change in our outcome variable? However, the Pr(>|t|) column purposely provides a more concise response to this question, using the asterisk notation that corresponds with the Signif. codes legend at the bottom of the Coefficients results section. In R model output, one asterisk means “p < .05”. Two asterisks mean “p < .01”; and three asterisks mean “p < .001”. These values are referred to as p-values in scientific literature. How can we use p-values to answer our question around model significance?

Asterisks in a regression table indicate the level of the statistical significance of a regression coefficient. Our understanding of statistical significance is based off of the idea of a random sample. When interpreting these asterisk values, we ask ourselves: if there truly is no relationship between clicks on an advertisement and product sales, then what are chances that, across many user clicks on an ad, we see behavior that suggests that there is no relationship?

For our clicks variable, with *** in the Pr(>|t|) column, the answer is very unlikely. The value of ***, or p < .001, means that random sample resulting in the regression coefficient and standard error that we observed for clicks-given that there was truly no difference relationship between clicks and product purchase—would occur in less than one time in a random draw of 100, on average. Given that we would so rarely observe the situation that suggests that there is no relationship between our click and total_convert, we can say that there is a statistically significant relationship between the two variables. Generally speaking, scientists accept that a variable coefficient with p-value less than 0.05 is statistically significant.

2.6.3 Measures of Model Fit

At the bottom of the output of summary() are a series of labeled metrics, like Residual Standard Error (RSE) and Multiple R-squared, which quantify the fit of our model. Our previous exercises have covered how to interpret and plot many of these measures, but it’s a helpful reminder for them to be summarized along with other model output.

2.7 Assessing Simple Linear Regression

Let’s practice our model interpretation skills! We know that for continuous independent variables, like podcasts, the regression coefficient represents the difference in the predicted value of sales for each one-dollar increase in podcasts. Given the output of calling summary(model) below, we can correctly say that for every one dollar increase in podcast advertisement spending, total sales of the related project increases by 1.742 dollars.

summary(model)

#output 
Call:
lm(formula = sales ~ radio, data = train)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.57927    0.91176  10.506  < 2e-16 ***
podcast      1.74240    0.03255   5.353 3.71e-07 ***
---    
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We could also extract the value of the podcast coefficient, the second coefficient returned by our model, using list indexing as follows:

podcast_coefficent <- model$coefficients[2] 

It is important to note that the interpretation of the intercept coefficient is slightly different from that of variable coefficients. The intercept coefficient represents the value we would predict for our outcome variable, sales, given that podcast spending is equal to zero. It is crucial to remember that the intercept coefficient is only interpretable if we can reasonably expect a zero value for all independent variables in a model. Assuming, just as our simple linear model does, that spending on podcasts is the only variable that explains changes in sales, it does not make sense for any sales to occur without podcast spending. Therefore, for this model, our intercept coefficient is not interpretable.

However, in many cases, intercept coefficients are interpretable! The analysis of any model results requires a thorough understanding of our data, the system that produces this data, and a critical approach to interpretation of coefficient values.

Use summary() to print out the results of model and model2.

summary(model)

Call:
lm(formula = total_convert ~ clicks, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0196 -0.4219 -0.0908 -0.0405 16.1315 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.040462   0.074873   13.90   <2e-16 ***
clicks      0.050369   0.002628   19.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.513 on 621 degrees of freedom
Multiple R-squared:  0.3717,    Adjusted R-squared:  0.3707 
F-statistic: 367.3 on 1 and 621 DF,  p-value: < 2.2e-16

summary(model_2)

Call:
lm(formula = total_convert ~ impressions, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0444 -0.4048  0.0140  0.0784 15.7470 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 9.154e-01  6.640e-02   13.79   <2e-16 ***
impressions 9.793e-06  3.915e-07   25.01   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.347 on 621 degrees of freedom
Multiple R-squared:  0.5019,    Adjusted R-squared:  0.5011 
F-statistic: 625.6 on 1 and 621 DF,  p-value: < 2.2e-16

Save the results of calling model$coefficients[2] to a variable called clicks_coefficient;

clicks_coefficient <- model$coefficients[2]

sprintf("Based on a simple regression of `total_convert` by `clicks`, we estimate that for every additional click, the number of product purchases increases by %s.", clicks_coefficient)

[1] "Based on a simple regression of `total_convert` by `clicks`, we estimate that for every additional click, the number of product purchases increases by 0.0503687412817466."

How might we interpret the coefficient estimate for our model intercept? Set the variable intercept_coefficient equal to the lowercase letter (ex. “a”) associated with the statement that correctly interprets the estimate value:

A. The intercept coefficient is not significantly different from zero, therefore the model suggests that when clicks equal zero, the number of total_converts is likely to be around zero. This is reasonable because if a user has not clicked on an ad, we do not expect them to purchase the related product.

B. The intercept coefficient is negative and significant; this means that when clicks equal zero, we expect the number of total_converts to be less than zero. This makes sense because the user did not click on an ad, so we expect zero or less products to be purchased.

C. The intercept coefficient is greater than zero, so when clicks equals zero we expect total_converts to be somewhen between 0.90 and 1.25. Even though a user does not click on an ad, they have been exposed to the product brand and might purchase the product by later searching for it through a search engine.

D. The The intercept coefficient is close to zero; this means that when when clicks equal zero, the number of total_converts is likely to be around zero. This is reasonable because if a user has not clicked on an ad, we do not expect them to purchase the related product.

intercept_coefficient <- "c"

2.8 Making Predictions

Data Scientists are often interested in building models to make predictions on new data. While the add_predictions() function from the modelr package makes it easy to predict new values from a technical standpoint, it is far more difficult to develop and assess accurate predicted values.

The most common metric used to compute the accuracy of predicted values is mean squared error on test data. Similar to residual squared error (RSE) and R-squared, MSE measures the average squared difference between predicted and observed values. When we are working with just one model, it is helpful to compare the difference between MSE on our training dataset, and MSE on test data. We can calculate training MSE for model using a combination of add_predictions() and summarise(). add_predictions() creates and adds predicted values from a model to a column called pred. summarise() then allows us to calculate the mean of the squared difference between our observed values (sales) and predicted values (pred).

train %>% 
  add_predictions(model) %>%
  summarise(MSE = mean((sales - pred)^2))

#output

       MSE
  31.60713

We can use the same combination of functions to calculate MSE for our test dataset, which results in a MSE of around 32.5. Testing MSE will almost always be higher than training MSE, as the model has been built off of training data; however, it is important to confirm that there is not a substantial difference between model training and test MSE. The value of using MSE to quantify prediction accuracy is more clear when comparing multiple models, as it allows us to determine which versions of a model best predicts an outcome variable. For instance, we could compute the MSE for a model of tv spending on sales.

model2 <- lm(sales ~ tv, data = train)

train %>% 
  add_predictions(model2) %>%
  summarise(MSE = mean((sales - pred)^2))

#output
    MSE
    27.28415

Comparing the train MSE for our tv-based model, at 27.28, to our train MSE for a podcast-based model, at 31.60, it is clear that the predictions from the tv-based model are more accurate, as the model’s MSE is lower. If a data scientist was trying to predict the expected volume of sales for a future business quarter, it would be a better idea for them to base their estimations off of a tv-based model.

Use add_predictions() and summarize() to calculate the test mean squared error of model (MSE), and save the result to a variable called mse_clicks.

Make sure to use the test dataset.

library(modelr)

mse_clicks <- test %>%
  add_predictions(model) %>%
  summarise(MSE = mean((total_convert - pred)^2))
mse_clicks

Print out both mse_clicks and mse_impressions. Which one is smaller? What does this tell us about the accuracy of our models on the test data?

mse_impressions <- test %>%
  add_predictions(model_2) %>%
  summarise(MSE = mean((total_convert - pred)^2))
mse_impressions

Let’s plot the predicted test values of the model with the smallest MSE against the observed test data.

First, use a combination of add_predictions(), ggplot() and geom_point()to plot the observed values from the test data.

When creating the ggplot() add an aes() where the x and y values correspond to the variables used when making the correct model.

Save the visualization to a variable called plot, then print out the variable.

plot <- test %>%
    add_predictions(model_2) %>%
    ggplot(aes(impressions, total_convert)) +
    geom_point() +
    geom_point(aes(y = pred), color = "blue")
plot

3 Multiple Linear Regression

linearvsmultiple

We’ve been able to really dig into the results of simple linear regression models and show how the results convey a substantial amount of information about the relationship between two variables. However, by this point you might be wondering–– what if I think variables other than podcast have contribute to the total sales of a product? You might remember that the primary assumption behind simple linear models is that the expected value of the outcome variable is a straight-line function of exclusively the predictor variable. This means that our simple linear models assume that all variation in the outcome variable is explained by the predictor variable. In the case of our sales dataset, we know this is almost certainly not true; oftentimes more money is spent on TV or newspaper ads than on podcasts, so this spending might have an even larger effect than podcast spend.

Thankfully, there are methods to include the effects of TV and newspaper in linear regression models. We can expand our model definition from a simple model of one predictor variable to a multiple model of, you guessed it, multiple predictor variables. The formal definition of multiple linear regression models is a direct extension of the formula for simple linear regression:

$ Y= beta_0 + beta_1∗X + beta_2∗X+error$

As in a simple linear model, Y represents the dollar value of products sold, X represents the amount spent on respective product podcast ads, and Beta_0 is the model intercept. Now, Beta_1, Beta_2, and Beta_3 represent each the coefficients of predictor variables. To build a similar model in R, using the standard lm() package, we still use the formula notation of Y ~ X:

model <- lm(sales ~ podcast + tv, data = train)

While building a multiple regression model is a straightforward extension of the code used to a build a simple model and the output of the model results below looks quite similar, a bit more effort goes into the interpretation of the results of this model. Remember that in a simple linear regression model, the regression coefficient represents the expected change in the dependent variable for a one-unit increase in the independent variable. In other words, the coefficient for podcast represents the expected increase in sales given a one dollar increase podcast advertisement spend. Because multiple linear regression includes more than one predictor variable, the coefficient estimates must be interpreted differently. In multiple linear regression, the regression coefficient represents the expected change in the dependent variable for a one-unit increase in the independent variable, holding all other variables in the model constant.

summary(model)

#output
Call:
lm(formula = sales ~ TV + podcast, data = train)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.583386   1.024616   4.473 1.65e-05 ***
TV          3.006340   1.004924   7.380 1.62e-11 ***
podcast     1.049249   1.027665   5.395 3.10e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

For example, a call to summary(model) shows that the coefficient for podcasts is equal to 1.04944. This means that, when one more dollar is spent on podcast advertising, about 1.049 more dollars of the related product is sold, given that there is no increase in the amount of money spent on tv advertisements. In this way, multiple linear regression models allow us to isolate the unique effect of one predictor variable on the outcome variable.

As this example shows, the selection of variables in a regression model can have wide-ranging impacts on the results and interpretation of our models! Let’s dive into one more exercise to practice building and interpreting multiple linear regression models.

3.1 Assessing Multiple Linear Regression

Time to pull it all together! The interpretation of coefficents in multiple linear regression is slightly different than that of coefficents in simple linear regression. Coefficent of independent continunous variables, like podcasts, represents the difference in the predicted value of sales for each one-dollar increase in podcasts, given that all other variables in the model, including tv, are held constant. Given the output of calling summary(model) below, we can correctly say that for every one dollar increase in podcast advertisement spending, while holding the amount spent on tv and newspaper constant, the total sales of the related product increases by 1.049 dollars.

summary(model)

#output
Call:
lm(formula = sales ~ TV + podcast + newspaper, data = train)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.583386   1.024616   4.473 1.65e-05 ***
TV          3.006340   1.004924   7.380 1.62e-11 ***
podcast     1.049249   1.027665   5.395 3.10e-07 ***
newspaper   1.006340   1.002924   6.380 1.12e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In addition, the interpretation of boolean categorical variables differs slightly from that of continous variables. The coefficent value associated with a boolean categorical variable represents the effect of changing from one category to another. for instance, the coefficient value of 1.006 for newspaper tell us that running print advertisements results in a 1.006 dollar increase in sales, holding the values of TV and podcast constant.

Data scientists often build many variations of a model with different combinations of independent variables before ultimately commiting to the model that best fits test data. Let’s practice building, interpreting, and selecting the best fit multi-linear model for our convert_clean dataset!

Build a multiple linear regression model which regresses impressions, clicks, and gender on total_convert, using our train dataset. Save the result to a variable called model; then call summary(model) to view the model results.

model <- lm(formula = total_convert ~ impressions + clicks + gender, data = train)

summary(model)

Call:
lm(formula = total_convert ~ impressions + clicks + gender, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1869 -0.3854 -0.1383  0.1819 15.7961 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.148e+00  8.821e-02  13.016  < 2e-16 ***
impressions  1.717e-05  1.158e-06  14.825  < 2e-16 ***
clicks      -4.672e-02  6.941e-03  -6.731 3.86e-11 ***
genderM     -3.356e-01  1.090e-01  -3.079  0.00217 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.301 on 619 degrees of freedom
Multiple R-squared:  0.5369,    Adjusted R-squared:  0.5346 
F-statistic: 239.2 on 3 and 619 DF,  p-value: < 2.2e-16

How might we interpret the coefficient estimate for gender? Set the variable gender_coefficient equal to the statement that most correctly interprets the estimate value — either “a”, “b”, or “c”:

A. The coefficient of the gender variable is not statistically significant, so we cannot come to any substantive conclusion from its’ value.

B. The coefficient of the gender variable is negative, which means that as total_convert, clicks, and impressions increases, men are less likely to purchase a advertised product.

C. The coefficient of the gender variable is negative. This means that a men are less likely than women––with the same value of clicks and impressions–– to purchase an advertised product.

gender_coefficient = "c"

Let’s build a second, simpler model so that we can confirm adding gender to our model increases its’ accuracy. Build a multiple linear regression model which regresses impressions, and clicks on total_convert. Save the result to a variable called model2

model_2 <- lm(formula = total_convert ~ impressions + clicks, data = train)

Compute the R-squared value for model and model2, and save the results to rsq_model and rsq_model2 respectively. Call both variables to view their values.

# compute r-squared below
rsq_model <- summary(model)$r.squared 
rsq_model

[1] 0.536877

rsq_model2 <- summary(model_2)$r.squared 
rsq_model2

[1] 0.5297831

Which model best fits our data? Set the variable best_fit equal to the larger r-squared value.

# define best_fit below
best_fit <- rsq_model

Set the variable gender_diff equal to the difference between rsq_model and rsq_model2.

# define gender_diff below
gender_diff <- rsq_model  - rsq_model2
sprintf("Based on the results of a series of multiple linear regressions on total_convert, we estimate that user gender accounts for approximately %s percent of the variation in product purchase rate.", gender_diff)

[1] "Based on the results of a series of multiple linear regressions on total_convert, we estimate that user gender accounts for approximately 0.00709394202260971 percent of the variation in product purchase rate."

4 Review

Whew, that’s a wrap! You’ve covered a lot of material related to linear regression and its implementation in R. Here are the main concepts we’ve covered:

• Statistical model building entails four main steps: 1) confirming data assumptions, 2) building a model on training data, 3) assessing model fit, and 4) analyzing model results.

• We can use a combination of qualitative methods, such as box-plots, and quantitative methods, like the correlation coefficient, to assess that data meets our assumptions

• We use the lm() method and Y ~ X notation to build a linear regression model. The Y variable is referred to as the outcome variable of the model, and any X variable is referred to as the predictor variable.

• We can use a similar set of qualitative and quantitative methods to evaluate the fit of our model, including a comparison of the plotted model to a LOESS smoother and statistics like mean squared error (MSE) and R-squared.

• MSE and R-squared statistics are summaries of the overall value of the model residuals. A residual is the difference between the value of a data point predicted by a model and its actual observed value.

• The results of a linear regression model include regression coefficients. These coefficients represent the effect their respective predictor variable has on the model’s outcome variable.

• In a simple linear regression, the regression coefficient represents the effect of a one-unit increase in the predictor variable.

T- he intercept coefficient represents the value of the outcome variable given that the predictor variable is equal to zero; this coefficient isn’t always meaningful and depends on the situation being modeled.

• The p-value associated with a regression coefficient helps us understand whether the effect of a variable is statistically significant.

• Multiple linear regression is similar to simple regression, except that it includes multiple predictor variables.

• In a multiple linear regression, the regression coefficient represents the effect of a one-unit increase in the respective predictor variable, given that all other predictor variables are held constant.

• In both simple and multiple linear regression, boolean categorical variables represent the total effect of switching from one category to another.