Lecture 2 - Market Response Model

1 Read.me

The structure of the notebook is as the following:

• Data Description: description of the data and mini case
• 5-Step Framework for Regression Analysis
• Model Comparison: a model without Brand Equity for comparison
• Risk Control: to check the normality and equal variance assumptions

For this notebook, you need to download and load a R package ggplot2.

library("ggplot2")

2 Data Description

You can download the “Marks Spencer.RData” file from Canvas. Please put the data in the same folder as the R notebook. After this, you can load the data into R.

load("Marks Spencer.RData")

The data includes the following variables:

• Week: the id’s of weeks in total 100 weeks
• Advertising: the advertising intensity
• Promotion: a factor of two-levels - “Yes” or “No”; the baseline is already set to “No”
• Price: the price index
• Brand_Equity: a measure of weekly brand equity of Marks & Spencer
• Sales: the overall sales

From the data, we know the sales from Week 1 to Week 75. Using the first 75 weeks of data, we want to predict the sales from Week 76 to Week 100. It is known ad hoc that the sales of Marks & Spencer decreased after Week 75. The question is: can we predict the sales decrease?

We are using two models, one model with Brand Equity and the other model without Brand Equity.

Run the code below to get a sneak peak of the data:

head(marks_spencer)

##   Week Advertising Promotion    Price Brand_Equity    Sales
## 1    1    6.643165        No 2.664169     99.01042 208.9463
## 2    2    6.551415       Yes 4.799289     99.31082 159.9861
## 3    3    3.606885        No 4.798553     99.44092 157.0999
## 4    4    4.323537        No 4.674879     99.52237 168.7835
## 5    5    5.859314        No 3.101983     99.53748 203.5160
## 6    6    6.366751        No 4.583886     99.67113 165.1949

3 Running Regression Analysis Step-by-Step

In this section, we will run a linear regression to predict sales by following the framework discussed in class step-by-step. For this analysis, we are using Sales as the DV and all the other variables as the IVs.

Note: The framework is presented to give a structure to follow as a novice. If you have become an expert, no need to follow it exactly.

We are using the first 75 weeks to estimate the model, and the remaining 25 weeks to do prediction. So, from the full data, we first obtain a train set.

train <- marks_spencer[1:75,]
str(train)

## 'data.frame':    75 obs. of  6 variables:
##  $ Week        : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ Advertising : num  6.64 6.55 3.61 4.32 5.86 ...
##  $ Promotion   : Factor w/ 2 levels "No","Yes": 1 2 1 1 1 1 1 1 2 1 ...
##  $ Price       : num  2.66 4.8 4.8 4.67 3.1 ...
##  $ Brand_Equity: num  99 99.3 99.4 99.5 99.5 ...
##  $ Sales       : num  209 160 157 169 204 ...

3.1 Examining the data

To get a good prediction, we want good IVs that are relevant (or correlated with the DV). In addition, we do NOT want IVs to be highly correlated with one another, which lead to information redundancy and eventually bad predictions.

3.1.1 Correlation between IVs and the DV

You can compute correlation coefficients as you learn from CMR. Alternatively, you can visualize the relationships between IVs and the DV with scatter plots. Here, we will only plot the relationship between Price and Sales as an example. Try other variables on your own.

ggplot(data = train, aes(x = Price, y = Sales)) + geom_point()

From the plot, it is clear that Price is negatively correlated with Sales.

3.1.2 Checking the VIFs

We also want to check if there’s any multi-collinearity problem. Essentially, we do not want IVs to be highly correlated. To this end, we calculate the VIFs. You are already given a function called vif along with the data. This function takes two inputs: the names of IVs that you want to check and a data frame containing those IVs, and outputs the VIF values of the variables (as a vector).

# getting the variable names of the IVs
vnames <- colnames(train)[2:5]
vif(vnames,train)

## PromotionYes  Advertising        Price Brand_Equity 
##     1.112713     1.047414     1.008407     1.068874

The criterion for the VIF values is if they are larger than 10. In the output, no VIFs are larger than 10. So, we do not have collinearity issue.

3.1.3 Formulating the model

We are using all the explanatory variables: Advertising, Promotion, Price and Brand_Equity to predict Sales. The regression equation is as the following: \[Sales = \beta_0+\beta_1Advertising + \beta_2Promotion + \beta_3Price + \beta_4BrandEquity+e\] ## Estimating The Model Next, we estimate our linear regression model with the lm function by translating the equation into an R formula. Please see Canvas for the formula language.

model_with_brand <- lm(Sales ~ 1 + Advertising + Promotion + Price + 
                        Brand_Equity, data = train)
summary(model_with_brand)

## 
## Call:
## lm(formula = Sales ~ 1 + Advertising + Promotion + Price + Brand_Equity, 
##     data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.493  -3.082   0.006   4.011  13.866 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  102.34877    7.33386  13.956   <2e-16 ***
## Advertising    1.74314    0.67864   2.569   0.0123 *  
## PromotionYes  -0.78282    1.68069  -0.466   0.6428    
## Price        -20.21582    0.85321 -23.694   <2e-16 ***
## Brand_Equity   1.46121    0.05805  25.170   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.363 on 70 degrees of freedom
## Multiple R-squared:  0.9485, Adjusted R-squared:  0.9455 
## F-statistic: 322.1 on 4 and 70 DF,  p-value: < 2.2e-16

3.2 Validating the model

From the results above, we first validate the model. We check a few things.

3.2.1 The overall significance

The idea of the overall significance is to compare the model we run to a null model with no IVs. If our model significantly outperforms the null model, it means it’s worth it. Or in statistical term, we test the null hypothesis that the coefficients (\(\beta\)’s) of all IVs are zero. Or, \[H_0:\beta_1=\beta_2=...=0\] For this model, the F-stat is 322.1, with a p-value < 2.2e-16 < 0.05. So, we can reject the null hypothesis, and the model is of predictive value.

3.2.2 The significance of coefficients

We have included multiple IVs in the model. However, it’s possible that an IV does not contribute to the prediction accuracy. To check this, we focus on an individual coefficient and test the following null hypothesis for a particular IV: \[H_0:\beta_k=0\] For example, if we check Brand Equity, the coefficient \(\beta_{brandequity}\) is 1.46, with a standard error of 0.06. The p-value < 2.2e-16 < 0.05. So, we reject the null hypothesis, and conclude that Brand Equity is of predictive value.

3.2.3 Model fit

The first index to look at is the residual standard errors, which form the foundation for other indexes such as R-squared. In our model, the value is 6.36. This is relatively low as the average Sales is 158.23.

mean(train$Sales)

## [1] 158.2321

To see how well the model fits with our data, we look at R-squared. The direct interpretation of R-squared is the percentage of variation explained by the model. Here, the R-squared of the model is 0.9485 or 94.85%. It means 94.85% of variation in the sales is captured by the model.

There is no clear cutoff value for R-squared. You need to consider the specific setting to judge whether the R-squared is good enough. For example, in sales prediction, we usually expect a big R-squared (e.g., 90%). This is because the retailers oftentimes are faced with a relatively stable environment where consumers show persistent habits of buying products. In our case, the value of 94.85% indicates a good fit.

3.3 Predicting with the model

After estimating the model, we can use it to do prediction. In the case, we want to predict the sales from Week 76 to Week 100. Let’s first create a prediction data set (or test set).

test <- marks_spencer[76:100,]
str(test)

## 'data.frame':    25 obs. of  6 variables:
##  $ Week        : num  76 77 78 79 80 81 82 83 84 85 ...
##  $ Advertising : num  6.57 6.61 6.55 6.94 5.43 ...
##  $ Promotion   : Factor w/ 2 levels "No","Yes": 1 2 1 1 1 1 2 2 2 1 ...
##  $ Price       : num  3.62 3.76 4.52 4.67 4.51 ...
##  $ Brand_Equity: num  53 50.8 45.8 44.9 44.5 ...
##  $ Sales       : num  NA NA NA NA NA NA NA NA NA NA ...

To predict the sales, we use the predict function. Please use ?predict for more information of this function. We can also obtain the 95% confidence interval of the prediction so we can prepare for the worst- and best-case scenario.

sales_with_brand <- predict(model_with_brand,
                            newdata = test,
                            interval = "confidence")
#coerce into a data frame
sales_with_brand <- as.data.frame(sales_with_brand)

Next, let’s create a line plot of the predicted sales to see the predicted trend.

# add week id for plotting
sales_with_brand$week <- 76:100
ggplot(data = sales_with_brand,
       aes(x = week, y = fit)) + geom_line()

From the plot, we can see a decreasing trend of sales.

4 Model Comparision

In this section, we are running a model without Brand Equity. The analyses are the same as in Section 3. We would not do things step-by-step. If you are interested, you can practice with this model by yourself.

The model specification is as the following (with Brand Equity not included): \[Sales = \beta_0+\beta_1Advertising + \beta_2Promotion + \beta_3Price +e\] Next, we estimate our linear regression model with the lm function by translating the equation into an R formula. Please see Canvas for the formula language.

model_without_brand <- lm(Sales ~ 1 + Advertising + Promotion + Price,
                          data = train)
summary(model_without_brand)

## 
## Call:
## lm(formula = Sales ~ 1 + Advertising + Promotion + Price, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.413 -17.598   2.055  17.177  33.894 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   225.630     17.182  13.131  < 2e-16 ***
## Advertising     2.916      2.131   1.368   0.1755    
## PromotionYes  -11.473      5.119  -2.241   0.0281 *  
## Price         -20.183      2.686  -7.515 1.32e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.03 on 71 degrees of freedom
## Multiple R-squared:  0.4822, Adjusted R-squared:  0.4603 
## F-statistic: 22.04 on 3 and 71 DF,  p-value: 3.416e-10

To compare this model with the model with Brand Equity, we use adjusted R-squared. This is because the two model have different no. of IVs. So, we must adjust for this difference. The adjusted R-squared of the model with vs. without Brand Equity is 0.9455 vs. 0.4603. Therefore, the model with Brand Equity fits the data better.

Next, let’s get the predicted sales of this model, using the same code as above.

sales_without_brand <- predict(model_without_brand,
                            newdata = test,
                            interval = "confidence")
#coerce into a data frame
sales_without_brand <- as.data.frame(sales_without_brand)

Again, we plot the predicted sales:

# add week id for plotting
sales_without_brand$week <- 76:100
ggplot(data = sales_without_brand,
       aes(x = week, y = fit)) + geom_line()

From this plot, we do not see a clear decreasing trend in sales. To conclude, to correctly predict the sales decrease, the company also needs to consider consumer metrics such as Brand Equity.

5 Risk Control

After running the models, we want to check if the statistical assumptions of linear regression are Okay. As discussed in class, the violation of these assumptions leads to biased estimation and bad prediction. Over the years, people have developed many tests for these assumptions. We will focus on two assumptions: 1) the normality assumption; and 2) the equal variance assumption.

We will use the full model (the model with brand equity or ) as an example. For other models, the procedure works the same.

5.1 The normality assumption

For the normality assumption, we are using KS test. KS test checks a variable against a normal distribution and the value of the KS test tells you whether the variable follows a normal distribution.

Note that the null hypothesis of the test is: \[H_0: The \: variable \: follows \:a \: normal \: distribution.\] For the normality assumption to hold, the null hypothesis should not be rejected. In other words, the KS test should NOT be significant if the normality assumption holds.

First, we obtain the residuals from the estimation results and visualize the distribution with a histogram.

# obtain the residuals using the residuals() function
error_term <- residuals(model_with_brand)

# build a histogram
# leave the data empty as we only plot a variable error_term
ggplot( , aes(x=error_term)) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

From the plot, it seems the error_term follows a normal distribution. We may further use KS test to do the formal testing. Note that by default, the KS test compares a variable to a standard normal distribution, or a normal distribution with mean 0 and standard deviation 1. Therefore, we need to first rescale the residuals to make it consistent with the default distribution.

# scale() function standardizes a variable use help(scale) to see more.
std_error_term <- scale(error_term)

# the ks test with ks.test(); 
#"pnorm" means testing against the standard normal distribution
ks.test(std_error_term,"pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  std_error_term
## D = 0.077301, p-value = 0.7318
## alternative hypothesis: two-sided

From the test results, p-value is 0.73 > 0.05. So, we cannot reject the null hypothesis. This implies the normality assumption is Okay.

5.2 The equal variance assumption

Another assumption of linear regression is all residuals follow the same normal distribution or the same variance. Although there is a formal test for this, it is beyond the scope of this course. Instead, we rely on plotting to eyeball this assumption.

To produce this plot, we need two values: the residuals and the DV (sales in this case). We do a scatter plot with the residuals as the Y-axis and the DV as the X-axis. In addition, sometimes it’s not so easy to see the patterns due to the scales of variables. So, we rescale the variables by standardization (with the scale() function).

ggplot( , aes(y = std_error_term, x = scale(train$Sales)))+ 
  geom_point() + geom_hline(yintercept = 0)

The trick is to check whether at different value of the DV (sales), the span or the range of the residuals are similar. You may add auxiliary lines to help you check.

Equal Variance Scatter Plot

From the plot, we an see that the range or span of residuals are similar. So, the equal variance assumption is likely to hold. Note that as we rely on eyeballing, it can be subjective. For the learning purpose, as long as you know how to use the plot to detect equal variance, your conclusion is not an issue.