07-supervised_learning.Rmd

---
output:
  html_document:
    toc: yes
  pdf_document:
    toc: yes
  html_notebook: default
---
```{r echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE, purl=FALSE}
library(knitr)
options(scipen = 999)
#This code automatically tidies code so that it does not reach over the page
opts_chunk$set(tidy.opts=list(width.cutoff=50),tidy=TRUE, rownames.print = FALSE, rows.print = 10)
opts_chunk$set(cache=T)
```

```{r message=FALSE, warning=FALSE, echo=F, eval=TRUE,paged.print = FALSE}
options(digits = 8)
```

# Supervised learning



## Linear regression

### Correlation


Before we start with regression analysis, we will review the basic concept of correlation first. Correlation helps us to determine the degree to which the variation in one variable, X, is related to the variation in another variable, Y. 

#### Correlation coefficient

The correlation coefficient summarizes the strength of the linear relationship between two metric (interval or ratio scaled) variables. Let's consider a simple example. Say you conduct a survey to investigate the relationship between the attitude towards a shop and the duration of being of its customer. The "Attitude" variable can take values between 1 (very unfavorable) and 12 (very favorable), and the "Duration" is measured in months. Let's further assume for this example that the attitude measurement represents an interval scale (although it is usually not realistic to assume that the scale points on an itemized rating scale have the same distance). To keep it simple, let's further assume that you only asked 12 people. We can create a short data set like this:    

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,paged.print = FALSE}
library(psych)
attitude <- c(6,9,8,3,10,4,5,2,11,9,10,2)
duration <- c(10,12,12,4,12,6,8,2,18,9,17,2)
att_data <- data.frame(attitude, duration)
att_data <- att_data[order(-attitude), ]
att_data$respodentID <- c(1:12)
str(att_data)
psych::describe(att_data[, c("attitude","duration")])
att_data
```

Let's look at the data first. The following graph shows the individual data points for the "duration" variable, where the y-axis shows the duration of residency in years and the x-axis shows the respondent ID. The blue horizontal line represents the mean of the variable (`r round(mean(att_data$duration),2)`) and the vertical lines show the distance of the individual data points from the mean.

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Scores for duration variable"}
library(ggplot2)
h <- round(mean(att_data$duration), 2)
ggplot(att_data, aes(x = respodentID, y = duration)) + 
  geom_point(size = 3, color = "deepskyblue4") + 
  scale_x_continuous(breaks = 1:12) + 
  geom_hline(data = att_data, aes(yintercept = mean(duration)), color ="deepskyblue4") +
  labs(x = "Observations",y = "Duration", size = 11) +
  coord_cartesian(ylim = c(0, 18)) +
  geom_segment(aes(x = respodentID,y = duration, xend = respodentID, 
                   yend = mean(duration)), color = "deepskyblue4", size = 1) + 
  theme(axis.title = element_text(size = 16),
        axis.text  = element_text(size=16),
        strip.text.x = element_text(size = 16),
        legend.position="none") + 
  theme_bw()
```

You can see that there are some respondents that have been the store's customers longer than average and some - shorter than average. Let's do the same for the second variable ("Attitude"). Again, the y-axis shows the observed scores for this variable and the x-axis shows the respondent ID.  

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Scores for attitude variable"}
ggplot(att_data, aes(x = respodentID, y = attitude)) + 
  geom_point(size = 3, color = "#f9756d") + 
  scale_x_continuous(breaks = 1:12) + 
  geom_hline(data = att_data, aes(yintercept = mean(attitude)), color = "#f9756d") +
  labs(x = "Observations",y = "Attitude", size = 11) +
  coord_cartesian(ylim = c(0,18)) +
  geom_segment(aes(x = respodentID, y = attitude, xend = respodentID, 
                   yend = mean(attitude)), color = "#f9756d", size = 1) + 
  theme_bw()
```

Again, we can see that some respondents have an above average attitude towards the store (more favorable) and some respondents have a below average attitude. Let's combine both variables in one graph now to see if there is some co-movement: 

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE,fig.align="center", fig.cap = "Scores for attitude and duration variables"}
ggplot(att_data) + 
  geom_point(size = 3, aes(respodentID, attitude), color = "#f9756d") +
  geom_point(size = 3, aes(respodentID, duration), color = "deepskyblue4")  + 
  scale_x_continuous(breaks = 1:12) + 
  geom_hline(data = att_data, aes(yintercept = mean(duration)), color = "deepskyblue4") +
  geom_hline(data = att_data, aes(yintercept = mean(attitude)), color = "#f9756d") +
  labs(x = "Observations", y = "Duration/Attitude", size = 11) +
  coord_cartesian(ylim = c(0, 18)) +
  scale_color_manual(values = c("#f9756d", "deepskyblue4")) +
  theme_bw()
```

We can see that there is indeed some co-movement here. The variables <b>covary</b> because respondents who have an above (below) average attitude towards the store also appear to have been its customers for an above (below) average amount of time and vice versa. Correlation helps us to quantify this relationship. Before you proceed to compute the correlation coefficient, you should first look at the data. We usually use a scatterplot to visualize the relationship between two metric variables:

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Scatterplot for durationand attitute variables"}
ggplot(att_data) + 
  geom_point(size = 3, aes(duration, attitude)) + 
  labs(x = "Duration", y = "Attitude", size = 11) +
  theme_bw()
```

How can we compute the correlation coefficient? Remember that the variance measures the average deviation from the mean of a variable:

\begin{equation} 
\begin{split}
s_x^2&=\frac{\sum_{i=1}^{N} (X_i-\overline{X})^2}{N-1} \\
     &= \frac{\sum_{i=1}^{N} (X_i-\overline{X})*(X_i-\overline{X})}{N-1}
\end{split}
(\#eq:variance)
\end{equation} 

When we consider two variables, we multiply the deviation for one variable by the respective deviation for the second variable: 

<p style="text-align:center;">
$(X_i-\overline{X})*(Y_i-\overline{Y})$
</p>

This is called the cross-product deviation. Then we sum the cross-product deviations:

<p style="text-align:center;">
$\sum_{i=1}^{N}(X_i-\overline{X})*(Y_i-\overline{Y})$
</p>

... and compute the average of the sum of all cross-product deviations to get the <b>covariance</b>:

\begin{equation} 
Cov(x, y) =\frac{\sum_{i=1}^{N}(X_i-\overline{X})*(Y_i-\overline{Y})}{N-1}
(\#eq:covariance)
\end{equation} 

You can easily compute the covariance manually as follows

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
x <- att_data$duration
x_bar <- mean(att_data$duration)
y <- att_data$attitude
y_bar <- mean(att_data$attitude)
N <- nrow(att_data)
cov <- (sum((x - x_bar)*(y - y_bar))) / (N - 1)
cov
```

Or you simply use the built-in ```cov()``` function:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
cov(att_data$duration, att_data$attitude)          # apply the cov function 
```

A positive covariance indicates that as one variable deviates from the mean, the other variable deviates in the same direction. A negative covariance indicates that as one variable deviates from the mean (e.g., increases), the other variable deviates in the opposite direction (e.g., decreases).

However, the size of the covariance depends on the scale of measurement. Larger scale units will lead to larger covariance. To overcome the problem of dependence on measurement scale, we need to convert the covariance to a standard set of units through standardization by dividing the covariance by the standard deviation (similar to how we compute z-scores).

With two variables, there are two standard deviations. We simply multiply the two standard deviations. We then divide the covariance by the product of the two standard deviations to get the standardized covariance, which is known as a correlation coefficient r:

\begin{equation} 
r=\frac{Cov_{xy}}{s_x*s_y}
(\#eq:corcoeff)
\end{equation} 

This is known as the product moment correlation (r) and it is straight-forward to compute:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
x_sd <- sd(att_data$duration)
y_sd <- sd(att_data$attitude)
r <- cov/(x_sd*y_sd)
r
```

Or you could just use the ```cor()``` function:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
cor(att_data[, c("attitude", "duration")], method = "pearson", use = "complete")
```

The properties of the correlation coefficient ('r') are:

* ranges from -1 to + 1
* +1 indicates perfect linear relationship
* -1 indicates perfect negative relationship
* 0 indicates no linear relationship
* ± .1 represents small effect
* ± .3 represents medium effect
* ± .5 represents large effect

#### Significance testing

How can we determine if our two variables are significantly related? To test this, we denote the population moment correlation *&rho;*. Then we test the null of no relationship between variables:

$$H_0:\rho=0$$
$$H_1:\rho\ne0$$

The test statistic is: 

\begin{equation} 
t=\frac{r*\sqrt{N-2}}{\sqrt{1-r^2}}
(\#eq:cortest)
\end{equation} 

It has a t distribution with n - 2 degrees of freedom. You can simply use the ```cor.test()``` function, which also produces the 95% confidence interval:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
cor.test(att_data$attitude, att_data$duration, alternative = "two.sided", method = "pearson", conf.level = 0.95)
```

To determine the linear relationship between variables, the data only needs to be measured using interval scales. If you want to test the significance of the association, the sampling distribution needs to be normally distributed (we usually assume this when our data are normally distributed or when N is large). If parametric assumptions are violated, you should use non-parametric tests:

* Spearman's correlation coefficient: requires ordinal data and ranks the data before applying Pearson's equation.
* Kendall's tau: use when N is small or the number of tied ranks is large.

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
cor.test(att_data$attitude, att_data$duration, alternative = "two.sided", method = "spearman", conf.level = 0.95)
cor.test(att_data$attitude, att_data$duration, alternative = "two.sided", method = "kendall", conf.level = 0.95)
```

Report the results:

A Pearson product-moment correlation coefficient was computed to assess the relationship between the duration of being a customer of a store and the attitude toward this store. There was a positive correlation between the two variables, r = 0.936, n = 12, p < 0.05. A scatterplot summarizes the results (Figure XY).

**A note on the interpretation of correlation coefficients:**

As we have already seen in Chapter 1, correlation coefficients give no indication of the direction of causality. In our example, we can conclude that the attitude toward the store is more positive as the months of being a customer increase. However, we cannot say that the duration causes the attitudes to be more positive. There are two main reasons for caution when interpreting correlations:

* Third-variable problem: there may be other unobserved factors that affect both the 'attitude' and the 'duration' variables
* Direction of causality: Correlations say nothing about which variable causes the other to change (reverse causality: attitudes may just as well cause the duration variable).


### Regression analysis


Correlations measure relationships between variables (i.e., how much two variables covary). Using regression analysis we can predict the outcome of a dependent variable (Y) from one or more independent variables (X). For example, we could be interested in how many products will we will sell if we increase the advertising expenditures by 1000 Euros? In regression analysis, we fit a model to our data and use it to predict the values of the dependent variable from one predictor variable (bivariate regression) or several predictor variables (multiple regression). The following table shows a comparison of correlation and regression analysis:

<br>

&nbsp; | Correlation	 | Regression	
-------------|--------------------------  | -------------------------- 
Estimated coefficient  | Coefficient of correlation (bounded between -1 and +1) | Regression coefficient (not bounded a priori)
Interpretation  | Linear association between two variables; Association is bidirectional | (Linear) relation between one or more independent variables and dependent variable; Relation is directional
Role of theory | Theory neither required nor testable  | Theory required and testable

<br>

#### Simple linear regression

In simple linear regression, we assess the relationship between one dependent (regressand) and one independent (regressor) variable. The goal is to fit a line through a scatterplot of observations in order to find the line that best describes the data (scatterplot).

Suppose you are a marketing research analyst at a big retail group and your task is to suggest, on the basis of historical data, a marketing plan for the next year that will maximize product sales. The product in question is beer brand Budweiser. The data set that is available to you includes information on the sales of Budweiser *move_ounce* (in ounces, FYI: 1 oz = 29,57 ml), prices *price_ounce* (in dollars per ounce), and several other variables: bonus buy - a price reduction if customers buy a certain quantity of a product (*sale_B*), price reduction in % (*sale_S*), and others. Let's load and inspect the data first: 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
regression <- read.table("https://raw.githubusercontent.com/dariayudaeva/RMA2024/main/data/bud_store102.csv", 
                          sep = ",", 
                          header = TRUE) # read in data
str(regression)
regression$store <- as.factor(regression$store) #convert grouping variable to factor
regression$brand_id <- as.factor(regression$brand_id) #convert grouping variable to factor
head(regression)
```

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE, paged.print = FALSE}
psych::describe(regression) #descriptive statistics using psych
```

As stated above, regression analysis may be used to relate a quantitative response ("dependent variable") to one or more predictor variables ("independent variables"). In a simple linear regression, we have one dependent and one independent variable and we regress the dependent variable on the independent variable.  

Here are a few important questions that we might seek to address based on the data:

* Is there a relationship between prices and sales? 
* How strong is the relationship between prices and sales?
* Which other variables contribute to sales?
* How accurately can we estimate the effect of each variable on sales?
* How accurately can we predict future sales?
* Is the relationship linear?
* Is there synergy among the advertising activities?

We may use linear regression to answer these questions. We will see later that the interpretation of the results strongly depends on the goal of the analysis - whether you would like to simply predict an outcome variable or you would like to explain the causal effect of the independent variable on the dependent variable (see Chapter 1). Let's start with the first question and investigate the relationship between advertising and sales. 

##### Estimating the coefficients

A simple linear regression model only has one predictor and can be written as:

\begin{equation} 
Y=\beta_0+\beta_1X+\epsilon
(\#eq:regequ)
\end{equation} 

In our specific context, let's consider only the influence of prices on sales for now:

\begin{equation} 
Sales=\beta_0+\beta_1*price+\epsilon
(\#eq:regequadv)
\end{equation} 

The word "price" represents data on advertising expenditures that we have observed and &beta;<sub>1</sub> (the "slope"") represents the unknown relationship between prices and sales. It tells you by how much sales will increase or decrease for an additional dollar added to price. &beta;<sub>0</sub> (the "intercept") is the number of sales we would expect if the price is set to 0. Note that the last assumption is highly theoretical: in the majority of real world scenarios, we never have such variables as prices, advertising expenditures, kilometers to the nearest store set to 0. Hence, it is incorrect to interpret the intercept like this. Together, &beta;<sub>0</sub> and &beta;<sub>1</sub> represent the model coefficients or *parameters*. The error term (&epsilon;) captures everything that we miss by using our model, including, (1) misspecifications (the true relationship might not be linear), (2) omitted variables (other variables might drive sales), and (3) measurement error (our measurement of the variables might be imperfect).

Once we have used our training data to produce estimates for the model coefficients, we can predict future sales on the basis of a particular value of price by computing:

\begin{equation} 
\hat{Sales}=\hat{\beta_0}+\hat{\beta_1}*price
(\#eq:predreg)
\end{equation} 

We use the hat symbol, <sup>^</sup>, to denote the estimated value for an unknown parameter or coefficient, or to denote the predicted value of the response (sales). In practice, &beta;<sub>0</sub> and &beta;<sub>1</sub> are unknown and must be estimated from the data to make predictions. In the case of our pricing example, the data set consists of the prices and product sales for 220 weeks (n = 220). Our goal is to obtain coefficient estimates such that the linear model fits the available data well. In other words, we fit a line through the scatterplot of observations and try to find the line that best describes the data. The following graph shows the scatterplot for our data, where the black line shows the regression line. The grey vertical lines shows the difference between the predicted values (the regression line) and the observed values. This difference is referred to as the residuals ("e").

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Ordinary least squares (OLS)"}
library(dplyr)
options(scipen = 999)
regression1 <- read.table("https://raw.githubusercontent.com/dariayudaeva/RMA2024/main/data/bud_store102.csv", 
                          sep = ",", 
                          header = TRUE)
lm <- lm(move_ounce ~ price_ounce, data = regression1)
regression1$yhat <- predict(lm)

ggplot(regression1, aes(x = price_ounce,y = move_ounce)) + 
  geom_point(size = 2, color = "deepskyblue4") + 
  labs(x = "Advertising expenditure (in Euros)", y = "Sales", size = 11) +
  geom_segment(aes(x = price_ounce, y = move_ounce, xend = price_ounce, yend = yhat), 
               color = "grey",size = 0.5, linetype = "solid", alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme(axis.title = element_text(size = 16),
        axis.text  = element_text(size = 16),
        strip.text.x = element_text(size = 16),
        legend.position="none") + 
  theme_minimal()
```

The estimation of the regression function is based on the idea of the method of least squares (OLS = ordinary least squares). The first step is to calculate the residuals by subtracting the observed values from the predicted values.

<p style="text-align:center;">
$e_i = Y_i-(\hat{\beta_0}+\hat{\beta_1}X_i)$
</p>

This difference is then minimized by minimizing the sum of the squared residuals:

\begin{equation} 
\sum_{i=1}^{N} e_i^2= \sum_{i=1}^{N} [Y_i-(\hat{\beta_0}+\hat{\beta_1X_i)}]^2\rightarrow min!
(\#eq:rss)
\end{equation} 

e<sub>i</sub>: Residuals (i = 1,2,...,N)<br>
Y<sub>i</sub>: Values of the dependent variable (i = 1,2,...,N) <br>
&\hat{beta;<sub>0</sub>}: Intercept<br>
&\hat{beta;<sub>1</sub>}: Regression coefficient / slope parameters<br>
X<sub>ni</sub>: Values of the nth independent variables and the i*th* observation<br>
N: Number of observations<br>

This is also referred to as the <b>residual sum of squares (RSS)</b>. Now we need to choose the values for &beta;<sub>0</sub> and &beta;<sub>1</sub> that minimize RSS. So how can we derive these values for the regression coefficient? The equation for &beta;<sub>1</sub> is given by:

\begin{equation} 
\hat{\beta_1}=\frac{COV_{XY}}{s_x^2}
(\#eq:slope)
\end{equation} 

The exact mathematical derivation of this formula is beyond the scope of this script, but the intuition is to calculate the first derivative of the squared residuals with respect to &beta;<sub>1</sub> and set it to zero, thereby finding the &beta;<sub>1</sub> that minimizes the term. Using the above formula, you can easily compute &beta;<sub>1</sub> using the following code:

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
cov_y_x <- cov(regression$price_ounce, regression$move_ounce)
cov_y_x
var_x <- var(regression$price_ounce)
var_x
beta_1 <- cov_y_x/var_x
beta_1
```

The interpretation of &beta;<sub>1</sub> is as follows: 

For every extra dollar increase of the price, sales can be expected to decrease by `r round(beta_1,3)` oz, which is around 270 liters.

Using the estimated coefficient for &beta;<sub>1</sub>, it is easy to compute &beta;<sub>0</sub> (the intercept) as follows:

\begin{equation} 
\hat{\beta_0}=\overline{Y}-\hat{\beta_1}\overline{X}
(\#eq:intercept)
\end{equation} 

The R code for this is:

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
beta_0 <- mean(regression$move_ounce) - beta_1*mean(regression$price_ounce)
beta_0
```

You may also verify this based on a scatterplot of the data. The following plot shows the scatterplot including the regression line, which is estimated using OLS.  

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE, fig.align="center", fig.cap = "Scatterplot"}
ggplot(regression, mapping = aes(price_ounce, move_ounce)) + 
  geom_point(shape = 1) +
  geom_smooth(method = "lm", fill = "blue", alpha = 0.1) + 
  labs(x = "Price ($ per oz)", y = "Sales (oz)") + 
  theme_bw()
```

The slope coefficient (&beta;<sub>1</sub>) tells you by how much sales (on the y-axis) would decrease if the price (on the x-axis) is increased by one unit ($).   

##### Significance testing

In a next step, we assess if the effect of prices on sales is statistically significant. This means that we test the null hypothesis H<sub>0</sub>: "There is no relationship between prices and sales" versus the alternative hypothesis H<sub>1</sub>: "The is some relationship between prices and sales". Or, to state this formally:

$$H_0:\beta_1=0$$
$$H_1:\beta_1\ne0$$

How can we test if the effect is statistically significant? Recall the generalized equation to derive a test statistic:

\begin{equation} 
test\ statistic = \frac{effect}{error}
(\#eq:teststatgeneral)
\end{equation} 

The effect is given by the &beta;<sub>1</sub> coefficient in this case. To compute the test statistic, we need to come up with a measure of uncertainty around this estimate (the error). This is because we use information from a sample to estimate the least squares line to make inferences regarding the regression line in the entire population. Since we only have access to one sample, the regression line will be slightly different every time we take a different sample from the population. This is sampling variation and it is perfectly normal! It just means that we need to take into account the uncertainty around the estimate, which is achieved by the standard error. Thus, the test statistic for our hypothesis is given by:

\begin{equation} 
t = \frac{\hat{\beta_1}}{SE(\hat{\beta_1})}
(\#eq:teststatreg)
\end{equation} 

After calculating the test statistic, we compare its value to the values that we would expect to find if there was no effect based on the t-distribution. In a regression context, the degrees of freedom are given by ```N - p - 1``` where N is the sample size and p is the number of predictors. In our case, we have 220 observations and one predictor. Thus, the degrees of freedom is 220 - 1 - 1 = 218. In the regression output below, R provides the exact probability of observing a t value of this magnitude (or larger) if the null hypothesis was true. This probability  is the p-value. A small p-value indicates that it is unlikely to observe such a substantial association between the predictor and the outcome variable due to chance in the absence of any real association between the predictor and the outcome.

To estimate the regression model in R, you can use the ```lm()``` function. Within the function, you first specify the dependent variable ("move_ounce") and independent variable ("price_ounce") separated by a ```~``` (tilde). As mentioned previously, this is known as _formula notation_ in R. The ```data = regression``` argument specifies that the variables come from the data frame named "regression". Strictly speaking, you use the ```lm()``` function to create an object called "sales_reg," which holds the regression output. You can then view the results using the ```summary()``` function: 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
sales_reg <- lm(move_ounce ~ price_ounce, data = regression) #estimate linear model
summary(sales_reg) #summary of results
```

Note that the estimated coefficients for &beta;<sub>0</sub> (`r round(summary(sales_reg)$coefficients[1],3)`) and &beta;<sub>1</sub> (`r round(summary(sales_reg)$coefficients[2],3)`) correspond to the results of our manual computation above. The associated t-values and p-values are given in the output. The t-values are larger than the critical t-values for the 95% confidence level, since the associated p-values are smaller than 0.05. In case of the coefficient for &beta;<sub>1</sub>, this means that the probability of an association between the prices and sales of the observed magnitude (or larger) is smaller than 0.05, if the value of &beta;<sub>1</sub> was, in fact, 0. This finding leads us to reject the null hypothesis of no association between prices and sales. 

The coefficients associated with the respective variables represent <b>point estimates</b>. To obtain a better understanding of the range of values that the coefficients could take, it is helpful to compute <b>confidence intervals</b>. A 95% confidence interval is defined as a range of values such that with a 95% probability, the range will contain the true unknown value of the parameter. For example, for &beta;<sub>1</sub>, the confidence interval can be computed as.

\begin{equation} 
CI = \hat{\beta_1}\pm(t_{1-\frac{\alpha}{2}}*SE(\beta_1))
(\#eq:regCI)
\end{equation} 

It is easy to compute confidence intervals in R using the ```confint()``` function. You just have to provide the name of you estimated model as an argument:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
confint(sales_reg)
```

For our model, the 95% confidence interval for &beta;<sub>0</sub> is [`r round(confint(sales_reg)[1,1],2)`,`r round(confint(sales_reg)[1,2],2)`], and the 95% confidence interval for &beta;<sub>1</sub> is [`r round(confint(sales_reg)[2,1],2)`,`r round(confint(sales_reg)[2,2],2)`]. Thus, we can conclude that when we do not spend any money on advertising, sales will be somewhere between `r round(confint(sales_reg)[1,1],0)` and `r round(confint(sales_reg)[1,2],0)` units on average. In addition, for each increase in advertising expenditures by one Euro, there will be an average increase in sales of between `r round(confint(sales_reg)[2,1],2)` and `r round(confint(sales_reg)[2,2],2)`. If you revisit the graphic depiction of the regression model above, the uncertainty regarding the intercept and slope parameters can be seen in the confidence bounds (blue area) around the regression line. 

##### Assessing model fit

Once we have rejected the null hypothesis in favor of the alternative hypothesis, the next step is to investigate how well the model represents ("fits") the data. How can we assess the model fit?

* First, we calculate the fit of the most basic model (i.e., the mean)
* Then, we calculate the fit of the best model (i.e., the regression model)
* A good model should fit the data significantly better than the basic model
* R<sup>2</sup>: Represents the percentage of the variation in the outcome that can be explained by the model
* The F-ratio measures how much the model has improved the prediction of the outcome compared to the level of inaccuracy in the model

Similar to ANOVA, the calculation of model fit statistics relies on estimating the different sum of squares values. SS<sub>T</sub> is the difference between the observed data and the mean value of Y (aka. total variation). In the absence of any other information, the mean value of Y ($\overline{Y}$) represents the best guess on where a particular observation $Y_{i}$ at a given level of advertising will fall:

\begin{equation} 
SS_T= \sum_{i=1}^{N} (Y_i-\overline{Y})^2
(\#eq:regSST)
\end{equation} 

The following graph shows the total sum of squares:

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Total sum of squares"}
library(dplyr)
options(scipen = 999)

ggplot(regression1, aes(x = price_ounce, y = move_ounce)) + 
  geom_point(size = 2, color = "deepskyblue4") + 
  labs(x = "Price ($ per oz)", y = "Sales (oz)", size = 11) +
  geom_segment(aes(x = price_ounce, y = move_ounce, xend = price_ounce, 
                   yend = mean(move_ounce)), color = "grey", 
               size = 0.5, linetype = "solid", alpha = 0.8) + 
  geom_hline(data = regression, aes(yintercept = mean(move_ounce)), color = "black", size = 1) +
  theme(axis.title = element_text(size = 16),
        axis.text  = element_text(size = 16),
        strip.text.x = element_text(size = 16),
        legend.position="none") + 
  theme_minimal()
```

Based on our linear model, the best guess about the sales level at a given level of prices is the predicted value $\hat{Y}_i$. The model sum of squares (SS<sub>M</sub>) therefore has the mathematical representation:

\begin{equation} 
SS_M= \sum_{i=1}^{N}  (\hat{Y}_i-\overline{Y})^2
(\#eq:regSSM)
\end{equation} 

The model sum of squares represents the improvement in prediction resulting from using the regression model rather than the mean of the data. The following graph shows the model sum of squares for our example:

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Ordinary least squares (OLS)"}
ggplot(regression1, aes(x = price_ounce, y = move_ounce)) + 
  geom_point(size = 2, color = "deepskyblue4") + 
  labs(x = "Price ($ per oz)",y = "Sales (oz)", size = 11) +
  geom_segment(aes(x = price_ounce, y = yhat, xend = price_ounce, yend = mean(move_ounce)), 
               color = "grey", size = 0.5, linetype = "solid", alpha = 0.8) + 
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  geom_hline(data = regression, aes(yintercept = mean(move_ounce)), color = "black", size = 1) +
  theme(axis.title = element_text(size = 16),
        axis.text  = element_text(size = 16),
        strip.text.x = element_text(size = 16),
        legend.position = "none") + 
  theme_minimal()
```

The residual sum of squares (SS<sub>R</sub>) is the difference between the observed data points ($Y_{i}$) and the predicted values along the regression line ($\hat{Y}_{i}$), i.e., the variation *not* explained by the model.

\begin{equation} 
SS_R= \sum_{i=1}^{N} ({Y}_{i}-\hat{Y}_{i})^2
(\#eq:regSSR)
\end{equation} 

The following graph shows the residual sum of squares for our example:

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.cap = "Ordinary least squares (OLS)"}
ggplot(regression1, aes(x = price_ounce, y = move_ounce)) + 
  geom_point(size = 2, color = "deepskyblue4") + 
  labs(x = "Price ($ per oz)",y = "Sales (oz)", size = 11) +
  geom_segment(aes(x = price_ounce, y = move_ounce, xend = price_ounce, yend = yhat),
               color = "grey", size = 0.5, linetype = "solid", alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, color = "black") +
  theme(axis.title = element_text(size = 16),
        axis.text  = element_text(size = 16),
        strip.text.x = element_text(size = 16),
        legend.position = "none") + 
  theme_minimal()
```

Based on these statistics, we can determine how well the model fits the data as we will see next. 

###### R-squared {-}

The R<sup>2</sup> statistic represents the proportion of variance that is explained by the model and is computed as:

\begin{equation} 
R^2= \frac{SS_M}{SS_T}
(\#eq:regSSR)
\end{equation} 

It takes values between 0 (very bad fit) and 1 (very good fit). Note that when the goal of your model is to *predict* future outcomes, a "too good" model fit can pose severe challenges. The reason is that the model might fit your specific sample so well, that it will only predict well within the sample but not generalize to other samples. This is called **overfitting** and it shows that there is a trade-off between model fit and out-of-sample predictive ability of the model, if the goal is to predict beyond the sample. We will come back to this point later in this chapter. 

You can get a first impression of the fit of the model by inspecting the scatter plot as can be seen in the plot below. If the observations are highly dispersed around the regression line (left plot), the fit will be lower compared to a data set where the values are less dispersed (right plot).

```{r message=FALSE, warning=FALSE, echo=FALSE, eval=TRUE, fig.align="center", fig.height = 4, fig.width = 10, fig.cap="Good vs. bad model fit"}
library(cowplot)
library(gridExtra)
library(grid)
set.seed(44)
#x3 <- rlnorm(250, log(1), log(0.6))
options(scipen = 999)
options(digits = 2)
x1 <- rnorm(200,614.41,485)
x <- as.data.frame(subset(x1,x1>0))
names(x)<-c('adspend')
error <- rnorm(nrow(x))
sales <- round(134 + 0.094*x$adspend + 50*error)
sales_data <- data.frame(sales,x$adspend) 
names(sales_data)<-c('sales','adspend')
examplereg <- subset(sales_data,sales>0 & adspend>0)
#summary(examplereg)
lm <- lm(sales ~ adspend, data = examplereg)
#summary(lm)
examplereg$yhat <- predict(lm)

scatter_plot1 <- ggplot(examplereg,aes(adspend,sales)) +  
  geom_point(size=2,shape=1) +    # Use hollow circles
  geom_smooth(method="lm") + # Add linear examplereg line (by default includes 95% confidence region);
  scale_x_continuous(name="advertising expenditures", limits=c(0, 1800)) +
  scale_y_continuous(name="sales", limits=c(0, 400)) +
  theme_bw() +
  labs(title = paste0("R-squared: ",round(summary(lm)$r.squared,2))) 

#x3 <- rlnorm(250, log(1), log(0.6))
options(scipen = 999)
options(digits = 2)
x1 <- rnorm(200,614.41,485)
x <- as.data.frame(subset(x1,x1>0))
names(x)<-c('adspend')
error <- rnorm(nrow(x))
sales <- round(134 + 0.094*x$adspend + 20*error)
sales_data <- data.frame(sales,x$adspend) 
names(sales_data)<-c('sales','adspend')
#summary(sales_data)
examplereg <- subset(sales_data,sales>0 & adspend>0)
#summary(examplereg)

lm <- lm(sales ~ adspend, data = examplereg)
examplereg$yhat <- predict(lm)

scatter_plot2 <- ggplot(examplereg,aes(adspend,sales)) +  
  geom_point(size=2,shape=1) +    # Use hollow circles
  geom_smooth(method="lm") + # Add linear examplereg line (by default includes 95% confidence region);
  scale_x_continuous(name="advertising expenditures", limits=c(0, 1800)) +
  scale_y_continuous(name="sales", limits=c(0, 400)) +
  theme_bw() +
  labs(title = paste0("R-squared: ",round(summary(lm)$r.squared,2))) 

#p <- plot_grid(plot1, plot2, ncol = 2)
p <- plot_grid(scatter_plot1,scatter_plot2, ncol = 2)
print(p)
# now add the title
#title <- ggdraw() + draw_label("", fontface='bold')
#plot_full <- plot_grid(title, p, ncol=1, rel_heights=c(0.1, 1)) # rel_heights values control title margins
#print(plot_full)
```

The R<sup>2</sup> statistic is reported in the regression output, so you don't need to compute it manually.


###### Adjusted R-squared {-}

Due to the way the R<sup>2</sup> statistic is calculated, it will never decrease if a new explanatory variable is introduced into the model. This means that every new independent variable either doesn't change the R<sup>2</sup> or increases it, even if there is no real relationship between the new variable and the dependent variable. Hence, one could be tempted to just add as many variables as possible to increase the R<sup>2</sup> and thus obtain a "better" model. However, this actually only leads to more noise and therefore a worse model. 

To account for this, there exists a test statistic closely related to the R<sup>2</sup>, the **adjusted R<sup>2</sup>**. It can be calculated as follows:

\begin{equation} 
\overline{R^2} = 1 - (1 - R^2)\frac{n-1}{n - k - 1}
(\#eq:adjustedR2)
\end{equation} 

where ```n``` is the total number of observations and ```k``` is the total number of explanatory variables. The adjusted R<sup>2</sup> is equal to or less than the regular R<sup>2</sup> and can be negative. It will only increase if the added variable adds more explanatory power than one would expect by pure chance. Essentially, it contains a "penalty" for including unnecessary variables and therefore favors more parsimonious models. As such, it is a measure of suitability, good for comparing different models and is very useful in the model selection stage of a project. In R, the standard ```lm()``` function automatically also reports the adjusted R<sup>2</sup> as you can see above.

###### F-test {-}

Similar to the ANOVA, another significance test is the F-test, which tests the null hypothesis:

$$H_0:R^2=0$$

<br>

Or, to state it slightly differently: 

$$H_0:\beta_1=\beta_2=\beta_3=\beta_k=0$$
<br>
This means that we test whether any of the included independent variables has a significant effect on the dependent variable. So far, we have only included one independent variable, but we will extend the set of predictor variables below.   

The F-test statistic is calculated as follows:

\begin{equation} 
F=\frac{\frac{SS_M}{k}}{\frac{SS_R}{(n-k-1)}}=\frac{MS_M}{MS_R}
(\#eq:regSSR)
\end{equation} 

which has a F distribution with k number of predictors and n degrees of freedom. In other words, you divide the systematic ("explained") variation due to the predictor variables by the unsystematic ("unexplained") variation. 

The result of the F-test is provided in the regression output as well. However, you might manually compute the F-test using the ANOVA results from the model:  

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE,paged.print = FALSE}
f_calc <- anova(sales_reg)$'Mean Sq'[1]/anova(sales_reg)$'Mean Sq'[2] #compute F
f_calc
f_crit <- qf(.95, df1 = 1, df2 = 100) #critical value
f_crit
f_calc > f_crit #test if calculated test statistic is larger than critical value
```

##### Using the model

After fitting the model, we can use the estimated coefficients to predict sales of Budweiser for different values of prices. Suppose the store plans to set the price per ounce to 2 dollars. How much will it sell? You can easily compute this either by hand:

$$\hat{sales}=63950.6 + (-9060.5)*2=45,829.6$$

<br>

... or by extracting the estimated coefficients from the model summary:

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
prediction <- summary(sales_reg)$coefficients[1,1] + 
              summary(sales_reg)$coefficients[2,1]*2 # the slope * 2 EUR
prediction
```

The predicted value of the dependent variable is 45,829.6 oz, i.e., the store will sell around 45,829.6 oz (~1,355 liters) of Budweiser.


#### Log-Log transformation

*(For more details about log-log transformation, a.k.a. multiplicative model, see chapter 6.1.3.3).* Have a look at the plots above again. You might notice some data specific pattern, making the data points look odd: they are pulled to lower edge of the scatterplot. In this particular case, we're dealing with different measurement scales of our independent and dependent variables. Moreover, you could also notice how odd the interpretation of the regression coefficients sounds. 

It is very rare that in the retailing context, the predictions are made as we did before. The concept that is used instead is familiar to you from the microeconomics course - elasticity is a measure that is used by retail managers and researchers much more often than mere unit changes.

Let's have a look at the plot again:  

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
ggplot(regression, mapping = aes(price_ounce, move_ounce)) + 
  geom_point(shape = 1) +
  geom_smooth(method = "lm", color = "lavenderblush4", fill = "red", alpha = 0.1) + 
  labs(x = "Price (ounce)", y = "Sales (ounce)") + 
  theme_minimal()
```

The way of obtaining a more reasonable view and interpretation in this case is called "multiplicative modeling", or log-log transformation (you can find additional details about log-log transformations below with a slightly different motivation and example). 

The multiplicative model has the following formal representation: 

\begin{equation} 
Y =\beta_0 *X_1^{\beta_1}*X_2^{\beta_2}*...*X_J^{\beta_J}*\epsilon
(\#eq:multiplicative)
\end{equation} 

This functional form can be linearized by taking the logarithm of both sides of the equation:

\begin{equation} 
log(Y) =log(\beta_0) + \beta_1*log(X_1) + \beta_2*log(X_2) + ...+ \beta_J*log(X_J) + log(\epsilon)
(\#eq:multiplicativetransformed)
\end{equation}

This means that taking logarithms of both sides of the equation makes linear estimation possible. The above transformation follows from two logarithm rules that we apply here: 

1. the product rule states that $log(xy)=log(x)+log(y)$; thus, when taking the logarithm of the right hand side of the multiplicative model, we can write $log(X_1) + log(X_2)... log(X_J)$ instead of $log(X_1*X_2*...X_J)$, and
2. the power rule states that $log(x^y) = ylog(x)$; thus, we can write $\beta*log(X)$ instead of $X^{\beta}$

Let's test how the scatterplot would look like if we use the logarithm of our variables using the ```log()``` function instead of the original values.  


```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
ggplot(regression, mapping = aes(log(price_ounce), log(move_ounce))) + 
  geom_point(shape = 1) +
  geom_smooth(method = "lm", color = "lavenderblush4", fill = "red", alpha = 0.1) + 
  labs(x = "Price", y = "Sales") + 
  theme_minimal()
```

You can see how the scales changed, and how the observations got more normally distributed. Hence, we can log-transform our variables and estimate the following equation: 

\begin{equation} 
log(sales) = log(\beta_0) + \beta_1*log(price) + log(\epsilon)
(\#eq:multiplicativetransformed1)
\end{equation}

Now, let's estimate a new regression by applying ```log()``` function to both sales and prices:

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE,paged.print = FALSE}
sales_reg2 <- lm(log(move_ounce) ~ log(price_ounce), data = regression)
summary(sales_reg2) #remember that now the interpretation changed
```

In this example, you would interpret the coefficient as follows: **A 1% increase in price leads to a 2.23% decrease in sales**. Hence, the interpretation is in proportional terms and no longer in units. This means that the coefficients in a log-log model can be directly interpreted as elasticities, which also makes communication easier. We can generally also inspect the R<sup>2</sup> statistic to see that the model fit has increased compared to the linear specification (i.e., R<sup>2</sup> has increased to 0.08 from 0.06). However, please note that the variables are now measured on a different scale, which means that the model fit in theory is not directly comparable.



#### Multiple linear regression


Multiple linear regression is a statistical technique that simultaneously tests the relationships between two or more independent variables and an interval-scaled dependent variable. The general form of the equation is given by:

\begin{equation} 
Y=(\beta_0+\beta_1*X_1+\beta_2*X_2+\beta_n*X_n)+\epsilon
(\#eq:regequ)
\end{equation} 

Again, we aim to find the combination of predictors that correlate maximally with the outcome variable. Note that if you change the composition of predictors, the partial regression coefficient of an independent variable will be different from that of the bivariate regression coefficient. This is because the regressors are usually correlated, and any variation in Y that was shared by X1 and X2 was attributed to X1. The interpretation of the partial regression coefficients is the expected change in Y when X is changed by one unit and all other predictors are held constant. 

Let's extend the previous example. Say, in addition to the influence of price itself, you are interested in estimating the influence of sales promotion on the amount of Budweiser. The corresponding equation, including bonus buy, would then be given by:

$$ Sales=\beta_0+\beta_1*price+\beta_2*bonus\_buy+\epsilon$$


&beta;<sub>1</sub> and &beta;<sub>2</sub> represent the unknown relationship between sales and independent variables (price and bonus buy, respectively). The corresponding coefficients tell you by how much sales will change for an additional dollar increase of price (when the other IVs are held constant) and by how much sales will change for an additional unit of price reduction (when price is held constant), etc. Thus, we can make predictions about sales using all these variables.

With several predictors, the partitioning of sum of squares is the same as in the bivariate model, except that the model is no longer a 2-D straight line. With two predictors, the regression line becomes a 3-D regression plane. While multiple regression models that have more than two predictors are not as easy to visualize, you may apply the same principles when interpreting the model outcome:

* Total sum of squares (SS<sub>T</sub>) is still the difference between the observed data and the mean value of Y (total variation)
* Residual sum of squares (SS<sub>R</sub>) is still the difference between the observed data and the values predicted by the model (unexplained variation)
* Model sum of squares (SS<sub>M</sub>) is still the difference between the values predicted by the model and the mean value of Y (explained variation)
* R measures the multiple correlation between the predictors and the outcome
* R<sup>2</sup> is the amount of variation in the outcome variable explained by the model

Estimating multiple regression models is straightforward using the ```lm()``` function. You just need to separate the individual predictors on the right hand side of the equation using the ```+``` symbol. In addition, as discussed before, we would need to use log-log transformation for our use case, which can be done in multiple regression context as well. Hence, we would specify the model as follows (note that bonus buy is already percentage in our data set, i.e., 0.2 value of price reduction is translated as 20% price decrease; we won't take a logarithm of it but rather interpret the results differently):

$$ log(Sales) =log(\beta_0) + \beta_1*log(Price) + \beta_2*bonus\_buy + log(\epsilon) $$


This regression could be estimated as follows: 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
multiple_sales_reg <- lm(log(move_ounce) ~ log(price_ounce) + sale_B, data = regression) # estimate the model
summary(multiple_sales_reg) #summary of results
```

The interpretation of the coefficients is as follows: 

* price (&beta;<sub>1</sub>): when price increases by 1%, sales will change by `r round(summary(multiple_sales_reg)$coefficients[2],3)`%
* bonus buy (&beta;<sub>2</sub>): when bonus buy increases by 1% (which is one step on a scale from 0 to 100 for a discount, i.e., this is equal to 1-unit increase), sales will change by `r round(summary(multiple_sales_reg)$coefficients[3],3)*100`%


The associated t-values and p-values are also given in the output. You can see that the p-values are smaller than 0.05 for price, while bonus sale is insignificant. Moreover, the p-value for F-test is smaller than 0.05. This means that if the null hypothesis was true (i.e., there was no effect between the variables and sales), the probability of observing associations of the estimated magnitudes (or larger) is very small (e.g., smaller than 0.05).     

Again, to get a better feeling for the range of values that the coefficients could take, it is helpful to compute <b>confidence intervals</b>. 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
confint(multiple_sales_reg)
```

What does this tell you? Recall that a 95% confidence interval is defined as a range of values such that with a 95% probability, the range will contain the true unknown value of the parameter. For example, for &beta;<sub>1</sub>, the confidence interval is [`r confint(multiple_sales_reg)[2,1]`,`r confint(multiple_sales_reg)[2,2]`]. Thus, although we have computed a point estimate of `r round(summary(multiple_sales_reg)$coefficients[2],3)` for the effect of price on sales based on our sample, the effect might actually just as well take any other value within this range, considering the sample size and the variability in our data. You could also visualize the output from your regression model including the confidence intervals using the `ggstatsplot` package as follows: 

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE,fig.align="center", fig.cap = "Confidence intervals for regression model"}
library(ggstatsplot)
ggcoefstats(x = multiple_sales_reg,
            title = "Sales predicted by price, bonus buy, and price reduction")
```

The output also tells us that `r summary(multiple_sales_reg)$r.squared*100`% of the variation can be explained by our model. You may also visually inspect the fit of the model by plotting the predicted values against the observed values. We can extract the predicted values using the ```predict()``` function. So let's create a new variable ```yhat```, which contains those predicted values.  

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE}
regression$logmove_ounce_hat <- fitted(multiple_sales_reg)
```

We can now use this variable to plot the predicted values against the observed values. In the following plot, the model fit would be perfect if all points would fall on the diagonal line. The larger the distance between the points and the line, the worse the model fit. In other words, if all points would fall exactly on the diagonal line, the model would perfectly predict the observed values. 

```{r message=FALSE, warning=FALSE, eval=TRUE, echo=TRUE, fig.align="center", fig.cap = "Model fit"}
ggplot(data = regression, aes(week, log(move_ounce))) +
  geom_vline(xintercept = regression$promoweek, colour = "lightgrey") +
  geom_line(aes(y = log(move_ounce), colour = "logsales"), size = 0.5) +
  geom_line(aes(y = (logmove_ounce_hat), colour = "logsales (predicted)"), size = 0.5) +
  scale_color_manual(values = c("black", "gold")) +
  theme_minimal()
```

**Partial plots**

In the context of a simple linear regression (i.e., with a single independent variable), a scatter plot of the dependent variable against the independent variable provides a good indication of the nature of the relationship. If there is more than one independent variable, however, things become more complicated. The reason is that although the scatter plot still show the relationship between the two variables, it does not take into account the effect of the other independent variables in the model. Partial regression plot show the effect of adding another variable to a model that already controls for the remaining variables in the model. In other words, it is a scatterplot of the residuals of the outcome variable and each predictor when both variables are regressed separately on the remaining predictors. In our example, the partial plot would show the effect of adding price as an explanatory variables while controlling for the variation that is explained by sales promotions in both variables (sales and price). Think of it as the purified relationship between price and sales that remains after controlling for other factors. The partial plots can easily be created using the ```avPlots()``` function from the ```car``` package:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Partial plots"}
library(car)
avPlots(multiple_sales_reg)
```



### Categorical predictors

#### Two categories

```{r}
categories <- read.table("https://raw.githubusercontent.com/WU-RDS/RMA2024/main/data/beer_categorical", 
                         sep = ",", 
                         header = TRUE) 
categories$store <- as.factor(categories$store)
categories$brand <- as.factor(categories$brand)
str(categories)
```


We will use a slightly different data set to explore additional opportunities for regression analysis. Suppose, you wish to investigate the effect of the variable "store" on sales, which is a categorical variable that can only take two levels (i.e., 98 = store with ID 98, and 100 = store with ID 100). Categorical variables with two levels are also called binary predictors; in our example, however, they are not decoded into typical binary view (i.e., they are not 0 and 1). It is straightforward to include these variables in your model as "dummy" variables. Dummy variables are factor variables that can only take two values. For our "store" variable, we can create a new predictor variable that takes the form:

\begin{equation} 
x_4 =
  \begin{cases}
    0       & \quad \text{if } i \text{th observation comes from store 98}\\
    1  & \quad \text{if } i \text{th observation comes from store 100}
  \end{cases}
(\#eq:dummycoding)
\end{equation} 

This new variable is then added to our regression equation from before, so that the equation becomes 

\begin{align}
Sales =\beta_0 &+\beta_1*price\\
      &+\beta_2*bonus\_buy\\
      &+\beta_3*store+\epsilon
\end{align}

where "store" represents the new dummy variable and is the coefficient associated with this variable. Estimating the model is straightforward - you just need to include the variable as an additional predictor variable. Note that the variable needs to be specified as a factor variable before including it in your model. If you haven't converted it to a factor variable before, you could also use the wrapper function ```as.factor()``` within the equation. 

First, let's reestimate the regression we had before (note that the result slightly changes because we are using a different data set - you can recall that with different samples, the estimation of true value changes). For the sake of easier interpretation, we use regular regression specification (i.e., not log-log transformed).

```{r}
multiple_regression_new <- lm(move_ounce ~ price_ounce + sale_B, data = categories) 
summary(multiple_regression_new)
```

Now, let's add the store variable:

```{r}
multiple_regression_store <- lm(move_ounce ~ price_ounce + sale_B + store, data = categories) 
summary(multiple_regression_store)
```


You can see that we now have an additional coefficient in the regression output, which tells us the effect of the dummy predictor. The dummy variable can generally be interpreted as the average difference in the dependent variable between the two groups, conditional on the other variables you have included in your model. In this case, the coefficient tells you the difference in sales between store 98 and 100 artists, and whether this difference is significant. Specifically, it means that sales in store 100 are on average 1,724.92 oz higher than in store 98, and this difference is significant (i.e., p < 0.05).  


#### More than two categories

Predictors with more than two categories, like our "brand"" variable, can also be included in your model. However, in this case one dummy variable cannot represent all possible values, since there are many brands (i.e., 1 = Amstel, 2 = Budweiser, 3 = Corona, 4 = Fosters, 5 = Heineken, 6 = Old Milwaukee). Thus, we need to create additional dummy variables. For example, for our "brand" variable, we create five dummy variables as follows:

\begin{equation} 
x_5 =
  \begin{cases}
    1       & \quad \text{if } i \text{th  product is Budweiser}\\
    0  & \quad \text{if } i \text{th product is Amstel}
  \end{cases}
(\#eq:dummycoding1)
\end{equation} 

\begin{equation} 
x_6 =
  \begin{cases}
    1       & \quad \text{if } i \text{th  product is Corona}\\
    0  & \quad \text{if } i \text{th product is Amstel}
  \end{cases}
(\#eq:dummycoding2)
\end{equation} 

and so on. 


We would then add these variables as additional predictors in the regression equation and obtain the following model

\begin{align}
Sales =\beta_0 &+\beta_1*price\\
      &+\beta_2*bonus\_buy\\
      &+\beta_3*store\\
      &+\beta_4*Budweiser\\
      &+\beta_5*Corona\\
      &+\beta_6*Fosters\\
      &+\beta_7*Heineken\\
      &+\beta_8*Old\_Milwaukee+\epsilon
\end{align}

where "Budweiser", "Corona", "Fosters", "Heineken", and "Old Milwaukee" represent our new dummy variables, and refer to the associated regression coefficients. You don't have to create the dummy variables manually as R will do this automatically when you add the variable to your equation.

The interpretation of the coefficients is as follows: $\beta_5$ is the difference in average sales between the brands "Amstel" and "Budweiser", $\beta_6$ is the difference in average sales between the brands "Amstel" and "Corona", and so on. Note that the level for which no dummy variable is created is also referred to as the *baseline*. In our case, "Amstel" would be the baseline brand. This means that there will always be one fewer dummy variable than the number of levels.

```{r}
multiple_regression_ext <- lm(move_ounce ~ price_ounce + sale_B + store + brand, data = categories) 
summary(multiple_regression_ext) 
```

How can we interpret the coefficients? It is estimated based on our model that products from the "Budweiser" brand will on average sell 22,298.63 oz more than products from the "Amstel" brand, and that products from the "Corona" brand will sell on average 2,288.88 oz more than the products from the "Amstel" brand, etc. The p-value of both these and some other brand-variables is smaller than 0.05, suggesting that there is statistical evidence for a real difference in sales between the brands

The level of the baseline category is arbitrary. As you have seen, R simply selects the first level as the baseline. If you would like to use a different baseline category, you can use the ```relevel()``` function and set the reference category using the ```ref``` argument. The following would estimate the same model using the second category as the baseline:

```{r}
multiple_regression_ext <- lm(move_ounce ~ price_ounce + sale_B + store + relevel(brand, ref = 2), data = categories) 
summary(multiple_regression_ext) 
```
  
Note that while your choice of the baseline category impacts the coefficients and the significance level, the prediction for each group will be the same regardless of this choice.


#### <span style="color: purple;">Non-linear relationships (extra)</span>

##### Multiplicative model 

In many practical applications, linear relationship might not be the case. Let's review the implications of a linear specification again:

* Constant marginal returns (e.g., an increase in ad-spend from 10€ to 11€ yields the same increase in sales as an increase from 100,000€ to 100,001€)
* Elasticities increase with X (e.g., advertising becomes relatively more effective; i.e., a relatively smaller change in advertising expenditure will yield the same return)

In many marketing contexts, these might not be reasonable assumptions. Consider the case of advertising. It is unlikely that the return on advertising will not depend on the level of advertising expenditures. It is rather likely that saturation occurs at some level, meaning that the return from an additional Euro spend on advertising is decreasing with the level of advertising expenditures (i.e., decreasing marginal returns). In other words, at some point the advertising campaign has achieved a certain level of penetration and an additional Euro spend on advertising won't yield the same return as in the beginning.

Let's use an example data set, containing the advertising expenditures of a company and the sales (in thousand units).

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
non_linear_reg <- read.table("https://raw.githubusercontent.com/IMSMWU/Teaching/master/MRDA2017/non_linear.dat", 
                          sep = "\t", 
                          header = TRUE) #read in data
head(non_linear_reg)
```

Now we inspect if a linear specification is appropriate by looking at the scatterplot:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Non-linear relationship"}
ggplot(data = non_linear_reg, aes(x = advertising, y = sales)) +
  geom_point(shape=1) + 
  geom_smooth(method = "lm", fill = "blue", alpha=0.1) + 
  theme_bw()
```

It appears that a linear model might **not** represent the data well. It rather appears that the effect of an additional Euro spend on advertising is decreasing with increasing levels of advertising expenditures. Thus, we have decreasing marginal returns. We could put this to a test and estimate a linear model:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
linear_reg <- lm(sales ~ advertising, data = non_linear_reg)
summary(linear_reg)
```

Advertising appears to be positively related to sales with an additional Euro that is spent on advertising resulting in 0.0005 additional sales. The R<sup>2</sup> statistic suggests that approximately 51% of the total variation can be explained by the model

To test if the linear specification is appropriate, let's inspect some of the plots that are generated by R. We start by inspecting the residuals plot. 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Residuals vs. Fitted"}
plot(linear_reg,1)
```

The plot suggests that the assumption of homoscedasticity is violated (i.e., the spread of values on the y-axis is different for different levels of the fitted values). In addition, the red line deviates from the dashed grey line, suggesting that the relationship might not be linear. Finally, the Q-Q plot of the residuals suggests that the residuals are not normally distributed. 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Q-Q plot"}
plot(linear_reg,2)
```

To sum up, a linear specification might not be the best model for this data set. 

In this case, a multiplicative model might be a better representation of the data. The multiplicative model has the following formal representation: 

\begin{equation} 
Y =\beta_0 *X_1^{\beta_1}*X_2^{\beta_2}*...*X_J^{\beta_J}*\epsilon
(\#eq:multiplicative)
\end{equation} 

This functional form can be linearized by taking the logarithm of both sides of the equation:

\begin{equation} 
log(Y) =log(\beta_0) + \beta_1*log(X_1) + \beta_2*log(X_2) + ...+ \beta_J*log(X_J) + log(\epsilon)
(\#eq:multiplicativetransformed)
\end{equation}

This means that taking logarithms of both sides of the equation makes linear estimation possible. The above transformation follows from two logarithm rules that we apply here: 

1. the product rule states that $log(xy)=log(x)+log(y)$; thus, when taking the logarithm of the right hand side of the multiplicative model, we can write $log(X_1) + log(X_2)... log(X_J)$ instead of $log(X_1*X_2*...X_J)$, and
2. the power rule states that $log(x^y) = ylog(x)$; thus, we can write $\beta*log(X)$ instead of $X^{\beta}$

Let's test how the scatterplot would look like if we use the logarithm of our variables using the ```log()``` function instead of the original values.  

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Linearized effect"}
ggplot(data = non_linear_reg, aes(x = log(advertising), y = log(sales))) + 
  geom_point(shape=1) + 
  geom_smooth(method = "lm", fill = "blue", alpha=0.1) +
  theme_bw()
```

It appears that now, with the log-transformed variables, a linear specification is a much better representation of the data. Hence, we can log-transform our variables and estimate the following equation: 

\begin{equation} 
log(sales) = log(\beta_0) + \beta_1*log(advertising) + log(\epsilon)
(\#eq:multiplicativetransformed1)
\end{equation}

This can be easily implemented in R by transforming the variables using the ```log()``` function:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE}
log_reg <- lm(log(sales) ~ log(advertising), data = non_linear_reg)
summary(log_reg)
```

Note that this specification implies decreasing marginal returns (i.e., the returns of advertising are decreasing with the level of advertising), which appear to be more consistent with the data. The specification is also consistent with proportional changes in advertising being associated with proportional changes in sales (i.e., advertising does not become more effective with increasing levels). This has important implications on the interpretation of the coefficients. In our example, you would interpret the coefficient as follows: **A 1% increase in advertising leads to a 0.3% increase in sales**. Hence, the interpretation is in proportional terms and no longer in units. This means that the coefficients in a log-log model can be directly interpreted as elasticities, which also makes communication easier. We can generally also inspect the R<sup>2</sup> statistic to see that the model fit has increased compared to the linear specification (i.e., R<sup>2</sup> has increased to 0.681 from 0.509). However, please note that the variables are now measured on a different scale, which means that the model fit in theory is not directly comparable. Also, we could use the residuals plot to confirm that the revised specification is more appropriate:

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE, fig.align="center", fig.cap = c("Residuals plot","Q-Q plot")}
plot(log_reg,1)
plot(log_reg,2)
```

Finally, we can plot the predicted values against the observed values to see that the results from the log-log model (red) provide a better prediction than the results from the linear model (blue). 

```{r message=FALSE, warning=FALSE, echo=TRUE, eval=TRUE,fig.align="center", fig.cap = "Comparison if model fit"}
non_linear_reg$pred_lin_reg <- predict(linear_reg)
non_linear_reg$pred_log_reg <- predict(log_reg)
ggplot(data = non_linear_reg) +
  geom_point(aes(x = advertising, y = sales),shape=1) + 
  geom_line(data = non_linear_reg,aes(x=advertising,y=pred_lin_reg),color="blue", size=1.05) + 
  geom_line(data = non_linear_reg,aes(x=advertising,y=exp(pred_log_reg)),color="red", size=1.05) + theme_bw()
```



## <span style="color: purple;">Logistic regression (extra)</span>


### Motivation and intuition

In the last section we saw how to predict continuous outcomes (e.g., sales) via linear regression models. Another interesting case is that of binary outcomes, i.e. when the variable we want to model can only take two values (yes or no, group 1 or group 2, dead or alive, etc.). To this end we would like to estimate how our predictor variables change the *probability* of a value being 0 or 1. In this case we can technically still use a linear model (e.g. OLS). However, its predictions will most likely not be particularly useful. A more useful method is the logistic regression. In particular we are going to have a look at the logit model. In the following dataset we are trying to predict whether a customer will churn (i.e., stop being our customer) any time soon. In the first step we are going to use only the "cash-back amount" index as a predictor. Later we are going to add more independent variables. 

```{r, echo=TRUE, warning=FALSE, message=FALSE, fig.align='center', fig.cap="Simulated binary outcome data"}
library(ggplot2)
library(gridExtra)

churn_data <- read.csv("https://raw.githubusercontent.com/WU-RDS/RMA2024/main/data/e_com_data.csv", 
                         sep = ",", header = T) 
head(churn_data)
str(churn_data)
# Variable "Churn" is 1 if a customer left and 0 else

# correct the variables types
churn_data$CustomerID <- as.factor(churn_data$CustomerID)
churn_data$Gender <- as.factor(churn_data$Gender)
churn_data$Tenure <- as.factor(churn_data$Tenure)
churn_data$PreferredLoginDevice <- as.factor(churn_data$PreferredLoginDevice)
churn_data$CityTier <- as.factor(churn_data$CityTier)
```

Below are two attempts to model the data. The left assumes a linear probability model (calculated with the same methods that we used in the last chapter), while the right model is a __logistic regression model__. As you can see, the linear probability model is able of producing probabilities that are above 1 and below 0 (see the lower right corner), which are not valid probabilities, while the logistic model stays between 0 and 1. Notice that songs with a higher cash-back value (on the right of the x-axis) seem to cluster more at $0$ and those with a lower more at $1$ so we expect a negative influence of cash-back value on the probability of churn. 

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.align='center', fig.cap="The same binary data explained by two models; A linear probability model (on the left) and a logistic regression model (on the right)"}
#Scatterplot showing the association between two variables using a linear model
plt.lin <- ggplot(churn_data, aes(CashbackAmount, Churn)) +  
  geom_point(shape = 1) +
  geom_smooth(method = "lm") +
  theme_minimal()

#Scatterplot showing the association between two variables using a glm
plt.logit <- ggplot(churn_data, aes(CashbackAmount, Churn)) +  
  geom_point(shape = 1) +
  geom_smooth(method = "glm", 
              method.args = list(family = "binomial"), 
              se = FALSE) +
  theme_minimal()

grid.arrange(plt.lin, plt.logit, ncol = 2)
```

A key insight at this point is that the connection between $\mathbf{X}$ and $Y$ is __non-linear__ in the logistic regression model. 

### Technical details of the model

As the name suggests, the logistic function is an important component of the logistic regression model. It has the following form:

$$
f(\mathbf{X}) = \frac{1}{1 + e^{-\mathbf{X}}}
$$
This function transforms all real numbers into the range between 0 and 1. We need this to model probabilities, as probabilities can only be between 0 and 1. 

```{r, echo = FALSE}
library(latex2exp)
x <- seq(-10, 10, length.out = 1000)
fx <- 1/(1+exp(-x))
df <- data.frame(x = x, fx = fx)
ggplot(df, aes(x = x, y = fx)) + 
  geom_line()+
  labs(y = TeX("$\\frac{1}{1+e^{-\\mathbf{X}}}$"), x = TeX("$\\mathbf{X}$"))+
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5)) + 
  theme_bw()
```



The logistic function on its own is not very useful yet, as we want to be able to determine how predictors influence the probability of a value to be equal to 1. To this end we replace the $\mathbf{X}$ in the function above with our familiar linear specification, i.e.

$$
\mathbf{X} = \beta_0 + \beta_1 * x_{1,i} + \beta_2 * x_{2,i} + ... +\beta_m * x_{m,i}\\
f(\mathbf{X}) = P(y_i = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 * x_{1,i} + \beta_2 * x_{2,i} + ... +\beta_m * x_{m,i})}}
$$

In our case we only have $\beta_0$ and $\beta_1$, the coefficient associated with cash-back. 

In general we now have a mathematical relationship between our predictor variables $(x_1, ..., x_m)$ and the probability of $y_i$ being equal to one. The last step is to estimate the parameters of this model $(\beta_0, \beta_1, ..., \beta_m)$ to determine the magnitude of the effects.  

### Estimation in R

We are now going to show how to perform logistic regression in R. Instead of ```lm()``` we now use ```glm(Y~X, family=binomial(link = 'logit'))``` to use the logit model. We can still use the ```summary()``` command to inspect the output of the model. 

```{r}
#Run the glm
logit_model <- glm(Churn ~ CashbackAmount, family = binomial(link = 'logit'), data = churn_data)
#Inspect model summary
summary(logit_model)
```

Noticeably this output does not include an $R^2$ value to asses model fit. Multiple "Pseudo $R^2$s", similar to the one used in OLS, have been developed. There are packages that return the $R^2$ given a logit model:

```{r}
library(DescTools)
PseudoR2(logit_model, which = "CoxSnell") # you can also use "McFadden", "McFaddenAdj", "Nagelkerke", "AldrichNelson", "VeallZimmermann", "Efron", "McKelveyZavoina", "Tjur", "all"
```

The coefficients of the model give the change in the [log odds](https://en.wikipedia.org/wiki/Odds#Statistical_usage) of the dependent variable due to a unit change in the regressor. This makes the exact interpretation of the coefficients difficult, but we can still interpret the signs and the p-values which will tell us if a variable has a significant positive or negative impact on the probability of the dependent variable being $1$. In order to get the odds ratios we can simply take the exponent of the coefficients. 

```{r}
exp(coef(logit_model))
```

This now gives the effect on the dependent variable: an additional dollar paid back as cash-back, on average, makes it $0.99$ time more likely (= by a constant factor of 0.99 = 1% less) for a customer to churn. 

Note: In some cases, depending on the scales, the coefficient can be extremely large (e.g., due to the fact that the variable is constrained to values between $0$ and $1$ and the coefficients are for a unit change). We can make the "unit-change" interpretation more meaningful by multiplying the variable by $100$, thus changing its scale. This linear transformation does not affect the model fit or the p-values.

We observe that cash-back negatively affects the likelihood of churning. To get the confidence intervals for the coefficients we can use the same function as with OLS.

```{r message=FALSE, warning=FALSE}
confint(logit_model)
```


To get the effect of an additional point at a specific value, we can calculate the odds ratio by predicting the probability at a value and at the value $+1$. For example, if we are interested in how much more (or, in our case, less) likely a customer with cash-back value of 201 compared to 200 is to churn, we can simply calculate the following:

```{r}
#Probability of churn with a cashback amount of 200
prob_200 <- exp(-(-summary(logit_model)$coefficients[1,1]-summary(logit_model)$coefficients[2,1]*200))
prob_200

#Probability of churn with a cashback amount of 201
prob_201 <- exp(-(-summary(logit_model)$coefficients[1,1]-summary(logit_model)$coefficients[2,1]*201))
prob_201

#Odds ratio
prob_201/prob_200
```

This is essentially what we got in the regression output earlier. So the odds are 1% lower at 201 than at 200 cash-back. 

#### Logistic model with multiple predictors

Of course we can also use multiple predictors in logistic regression as shown in the formula above. We might want to add order amount from last year, days since last order, warehouse-to-home distance, and order count.

Again, the familiar formula interface can be used with the ```glm()``` function. All the model summaries shown above still work with multiple predictors.

```{r}
multiple_logit_model <- glm(Churn ~ OrderAmountHikeFromlastYear + DaySinceLastOrder +
                              WarehouseToHome + OrderCount + CashbackAmount, 
                            family = binomial(link = 'logit'), data = churn_data)
summary(multiple_logit_model)
PseudoR2(multiple_logit_model, which = "CoxSnell")
```

Again, to properly interpret the coefficient, we extract odds ratio:
```{r}
exp(coef(multiple_logit_model))
confint(multiple_logit_model)
```

We can see that such variables as warehouse-to-home distance and order count increase the probability of churn, while days since last order and cash-back value decrease it.


#### Model selection

The question remains, whether a variable *should* be added to the model. We will present two methods for model selection for logistic regression. The first is based on the _Akaike Information Criterium_ (AIC). It is reported with the summary output for logit models. The value of the AIC is __relative__, meaning that it has no interpretation by itself. However, it can be used to compare and select models. The model with the lowest AIC value is the one that should be chosen. Note that the AIC does not indicate how well the model fits the data, but is merely used to compare models. 

For example, consider the following model, where we exclude the ```order count``` predictor. Seeing as it was able to contribute significantly to the explanatory power of the model, the AIC increases, indicating that the model including ```order count``` is better suited to explain the data. We always want the *lowest* possible AIC. 

```{r, message = FALSE, warning = FALSE}
multiple_logit_model2 <- glm(Churn ~ OrderAmountHikeFromlastYear + DaySinceLastOrder +
                              WarehouseToHome + CashbackAmount, 
                            family = binomial(link = 'logit'), data = churn_data)
summary(multiple_logit_model2)
```

As a second measure for variable selection, you can use the pseudo $R^2$s as shown above. The fit is worse according to all three values presented here, when excluding the order count. 

```{r}
PseudoR2(multiple_logit_model2, which = "CoxSnell")
```


#### Predictions

We can predict the probability given an observation using the ```predict(my_logit, newdata = ..., type = "response")``` function. Replace ```...``` with the observed values for which you would like to predict the outcome variable.


```{r, message = FALSE, warning = FALSE}
# Prediction for one observation
predict(multiple_logit_model, 
        newdata = data.frame(OrderAmountHikeFromlastYear = 5, DaySinceLastOrder = 30,
                             WarehouseToHome = 10, OrderCount = 10, CashbackAmount = 300), 
        type = "response")
```

The prediction indicates that a customer who made $5$ orders last year, $10$ orders in total, last order happened $30$ days ago, warehouse-to-home distance is 10 km, and cash-back amount is 300 dollars, has a tiny ($0.2\%$ chance of churning. 



## Learning check {-}

**(LC4.1) What is a correlation coefficient?**

- [ ] It describes the difference in means of two variables
- [ ] It describes the causal relation between two variables
- [ ] It is the standardized covariance
- [ ] It describes the degree to which the variation in one variable is related to the variation in another variable
- [ ] None of the above 

**(LC4.2) Which line through a scatterplot produces the best fit in a linear regression model?**

- [ ] The line associated with the steepest slope parameter
- [ ] The line that minimizes the sum of the squared deviations of the predicted values (regression line) from the observed values
- [ ] The line that minimizes the sum of the squared residuals
- [ ] The line that maximizes the sum of the squared residuals
- [ ] None of the above 

**(LC4.3) Which of the following statements about the adjusted R-squared is TRUE?**

- [ ] It is always larger than the regular $R^{2}$
- [ ] It increases with every additional variable
- [ ] It increases only with additional variables that add more explanatory power than pure chance
- [ ] It contains a “penalty” for including unnecessary variables
- [ ] None of the above 

**(LC4.4) When do you use a logistic regression model?**

- [ ] When the dependent variable is continuous
- [ ] When the independent and dependent variables are binary
- [ ] When the dependent variable is binary
- [ ] None of the above 

**(LC4.5) What is the correct way to implement a linear regression model in R? (x = independent variable, y = dependent variable)?**

- [ ] `lm(y~x, data=data)`
- [ ] `lm(x~y + error, data=data)`
- [ ] `lm(x~y, data=data)`
- [ ] `lm(y~x + error, data=data)`
- [ ] None of the above 

**(LC4.6) Consider the output from a bivariate correlation below**

- [ ] …lower prices cause higher sales.
- [ ] …lower prices are associated with higher sales and vice versa.
- [ ] None of the above


```{r,echo=FALSE,out.width = '75%',fig.align='center',fig.cap = ""}
knitr::include_graphics("./images/cor_table11.png")
```
**(LC4.7) When interpreting the statistical significance of a regression coefficient, the p-value…**

- [ ] …of 0.05 means that, if the null hypothesis is true (i.e., if the independent variable would NOT affect the outcome), the odds are 19 in 20 of getting a regression coefficient as large or larger than the estimated coefficient
- [ ] …of 0.05 means that, if the null hypothesis is true (i.e., if the independent variable would NOT affect the outcome), the odds are 1 in 20 of getting a regression coefficient as large or larger than the estimated coefficient
- [ ] ...of 0.05 means that the effect is statistically significant at the 5% level. 
- [ ] ...does not tell you anything about the importance of the effect
- [ ] …will get smaller, the larger the calculated value of the test statistic (t-value).

**(LC4.8) In which setting(s) would a regression coefficient be interpreted as "statistically significant"?**

- [ ] When the absolute value of the calculated test-statistic (e.g., t-value) exceeds the critical value of the test statistic at your specified significance level (e.g., 0.05)
- [ ] When the test-statistic (e.g., t-value) is lower than the critical value of the test statistic at your specified significance level (e.g., 0.05)
- [ ] When the confidence interval associated with the test does not contain zero
- [ ] When the p-value is smaller than your specified significance level (e.g., 0.05)

**(LC4.9) When interpreting the significance of the coefficients in a regression model, what is the relationship between the test statistic (e.g., t-value) and the p-value?**

- [ ] The lower the absolute value of the test statistic, the lower the p-value
- [ ] The higher the absolute value of the test statistic, the higher the p-value
- [ ] There is no connection between the test statistic and the p-value
- [ ] The higher the absolute value of the test statistic, the lower the p-value

**(LC4.10) What does the term overfitting refer to?**

- [ ] A regression model that fits to a specific data set so poorly, that it will not generalize to other samples
- [ ] A regression model that fits to a specific data set so well, that it will generalize to other samples particularly well
- [ ] A regression model that has too many predictor variables
- [ ] A regression model that fits to a specific data set so well, that it will only predict well within the sample but not generalize to other samples

## References {-}

* Field, A., Miles J., & Field, Z. (2012): Discovering Statistics Using R. Sage Publications (**chapters 6, 7, 8**).
* James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013): An Introduction to Statistical Learning with Applications in R, Springer (**chapter 3**)