---
title: "Multiple Regression"
subtitle: "Chapter 6.1-6.2"
author: "Qi Yi"
date: "10.02.2025"
format: 
  revealjs:
    code-line-numbers: true
    transition: slide
    transition-speed: fast
    slide-number: true
    scrollable: true
    incremental: true
    footer: "Quarto presentation is DOPE 🤓"
    logo: regression logo.png
    background-image: content_background1.png
    background-size: contain
    title-slide-attributes: 
      data-background-image: front page background.png
    theme: [simple,theme.scss] 
    linestretch: 1.7
    embed-resources: true
presentation: slide
editor: visual
---

## In today's session, we will cover

1.  what is multiple regression

2.  how to build a multiple regression model in R

3.  how to plot a multiple regression model

------------------------------------------------------------------------

### Why we need multiple regression? {background-image="content_background1.png"}

::: incremental
-   In *multiple regression,* we attempt to predict a dependent or response variable y on the basis of an assumed linear relationshop with several independent variables $x_1$, $x_2$, ..., $x_k$.
:::

------------------------------------------------------------------------

### The multiple linear regression model can be expressed as {background-image="content background1.png"}

-   $$
    y=β_0 +β_1x_1 +β_2x_2 + ... +β_kx_k + e
    $$
-   Y is the *dependent variable*, $x_1$, $x_2$, ..., $x_k$ are *independent variables*
-   k is the number of independent variables
-   $\beta_0$, $\beta_1$ ..., $\beta_k$ are *regression coefficients*

------------------------------------------------------------------------

### Example with data

We will use the data from the paper "Which words are most iconic? Iconicity in English sensory words" (Winter et al, 2017).

-   Whether word forms resemble their meaning??

    This paper investigated iconicity of English words through the relationship between systematicity and the sensory properties of words. 3001 English words were rated by English native speakers.

-   Iconicity as a function of sensory experience

------------------------------------------------------------------------

### Load the data

```{r}
#| echo: true
#| slide-fragment: true
#| code-line-numbers: True 
library(tidyverse)
library(broom)

rawdata<- read_csv("/Users/qiyi/Github_repository/R Group/Ch6.1-6.2/winter_2017_iconicity.csv")
rawdata %>% print(n=10, width=Inf) # width = Inf is to see all the columns of the tibble
```

:::: fragment
::: smaller
-   SER - sensory experience, Syst - systematicity, Freq - word frequency, CoreseImag - imageability, Conc - concreteness

-   Sensory experiences were rated from 1 to 7. Five common senses were considered (taste, touch, sound, sight and smell)

-   Iconicity was rated from -5 to 5
:::
::::

------------------------------------------------------------------------

### Read the data {background-image="content background1.png" style="color: #708993;"}

```{r}
#| echo: true
range(rawdata$Iconicity, na.rm=TRUE)
#check missing data
not_missing_imag<- sum(is.na(rawdata$CorteseImag))
not_missing_SER <- sum(is.na(rawdata$SER))
not_missing_Conc <- sum(is.na(rawdata$Conc))
not_missing_syst <- sum(is.na(rawdata$Syst))
not_missing_freq <- sum(is.na(rawdata$Freq))
not_missing_iconicity <- sum(is.na(rawdata$Iconicity))

not_missing_imag
not_missing_SER
not_missing_Conc
not_missing_syst
not_missing_freq
not_missing_iconicity

```

------------------------------------------------------------------------

### Read the data 2 {background-image="content background1.png" style="color: #708993;"}

```{r}
#| echo: true
most_iconic<- rawdata %>% mutate(Iconicity = round(Iconicity,3)) %>% 
  filter(Iconicity==max(Iconicity))
least_iconic <-rawdata %>% mutate(Iconicity = round(Iconicity,3)) %>% 
  filter(Iconicity==min(Iconicity))

```

<br>

```{r}
#| echo: true
most_iconic
least_iconic
```

------------------------------------------------------------------------

### Pre-processing the Data

Plot the distribution of the raw frequency data

::: fragment
```{r}
#| echo: true
#| output-location: fragment
hist(rawdata$Freq,main = "Histogram of Raw Frequency",xlab = "Frequency",
              ylab = "Count",col = "#A376A2",border = "#DDC3C3")

```
:::

-   The raw frequency data is extremely skewed, therefore, we need to log transform the frequency data

------------------------------------------------------------------------

#### Log transform frequency

```{r}
#| slide-title: Histogram of Log Frequency
#| echo: true
#| code-line-numbers: true
#| output-location: fragment
rawdata <- mutate(rawdata, logfreq = log10(Freq))
hist(
  rawdata$logfreq,
  main = "Histogram of Log Frequency",
  xlab = "Log10(Freq)",
  ylab = "Count",
  col = "#A376A2",
  border = "#DDC3C3"
)

```

------------------------------------------------------------------------

### Build the model

model1 \<- lm (Iconicity\~ SER + CorteseImag + Syst + logfreq, data= rawdata)

```{r}

#| echo: true
model1 <- lm (Iconicity~ SER + CorteseImag + Syst + logfreq, data= rawdata)

```

<br>

```{r}
#| echo: true
glance(model1)

```

<br>

```{r}
#| echo: true
summary(model1)
```

------------------------------------------------------------------------

#### Check R square

The model explained about 21% of the variation in iconicity ratings.

```{r}
#| echo: true
glance(model1)$r.squared
```

------------------------------------------------------------------------

#### Check coefficients {auto-animate="true"}

We can use `round()` to round the coefficients to make them easier for explanation.

-   `round()` takes two arguments, 1) a vector of numbers, 2) a number indicating to how many decimals the vector should be rounded to.

::: fragment
```{r}
#| echo: true
#| output-location: fragment
tidy (model1) %>% select(term, estimate) %>% mutate (estimate= round(estimate, 2))
```
:::

------------------------------------------------------------------------

Based on the estimates, a predictive equation can be expressed as:

$$
iconicity = 1.54 + 0.50 * SER + (-0.26) * CorteseImag + 401.52 * Syst + (-0.25) * logfreq
$$

-   What is the the strongest predictor from this result?

-   What does a one-unit change mean for systematicity?

:::: fragment
::: callout-caution
One-unit change of systematicity leads to 401.5 change on the iconicity rating. 😱
:::
::::

<br>

::: fragment
```{r}
#| echo: true
range(rawdata$Syst, na.rm=TRUE)
```
:::

------------------------------------------------------------------------

## Standardization !!

To make the independent variables more comparable!

-   `scale( )` centering and z-scoring

::: fragment
```{r}
#| echo: true
#| output-location: fragment
st_data<- mutate(rawdata,
                 SER_st=scale(SER),
                 CorteseImag_st= scale(CorteseImag),
                 Syst_st=scale(Syst),
                 Freq_st=scale(logfreq))
```
:::

------------------------------------------------------------------------

### Re-build the model

```{r}
#| echo: true
#| output-location: fragment
model2<- lm(Iconicity~SER_st+CorteseImag_st+Syst_st+Freq_st, data=st_data)
summary(model2)
```

:::: fragment
::: callout-important
R-square value DID NOT change!

Standardization DID NOT change the underlying model!
:::

-   What does one-unit change of sensory experience mean?

-   It means "For each increase in sensory experience by 1 STD, iconicity rating increases by 0.5".
::::

------------------------------------------------------------------------

### Plots

::: smaller
-   In the first plot, we want to see the red line be as horizontal as possible, to indicate that the only variation left in the data is due to the unexplained errors.
-   The second plot is to check if the residuals are normally distributed
-   The third plot shows far individual values are from the regression line. Again we want to see a horizontal red line with no particular patterns
-   The fourth plot checks outliers
:::

```{r}
#| echo: true
plot(model2)

```

------------------------------------------------------------------------

### Partial Regression

Plot one dependent variable while controlling for others

```{r}
#| echo: true
library (car)
avPlots(model2)
```

::: smaller
-   The slope of the line in the partial regression plot is the coefficient of the corresponding variable in the multiple regression model.
:::

------------------------------------------------------------------------

#### Reference

::: callout-note
Data for *Statistics For Linguistics: An Introduction Using R*, can be downloaded here.

<https://osf.io/34mq9/files/osfstorage>
:::
