Introduction to Linear Regression

pkgs <- c("tidyverse", "dplyr", "haven", "foreign", "lme4", "plyr", "nlme", "lsr", "emmeans", "afex", "knitr", "kableExtra", "QuantPsyc", "car", "readxl", "pastecs")

packages <- rownames(installed.packages())
p_to_install <- pkgs[!(pkgs %in% packages)]

if(length(p_to_install) > 0){
  install.packages(p_to_install)
}

lapply(pkgs, library, character.only = TRUE)

album_df <- read.delim("Album Sales 2.dat")

---

A note on assumptions

Here are the 3 assumptions that apply to all parametric tests:

1. Normally distributed data
   • Shapiro-Wilk test: caution–it’s easy to get a significant result if your sample is large enough
   • Q-Q plot: each value is compared to the expected value. If normally distributed, you get a diagonal line.
   • Plots()
     • REGRESSION
       • Assumption of linearity
     • ANOVA
       • Assumption that distributions within groups are normally distributed
2. Homogeneity of variance: The variances should be the same throughout the data.
   • REGRESSION
     • Homoscedasticity: inspect plot comparing residuals vs predicted values
   • ANOVA
     • Levene’s or Bartlett’s test (compares variances between groups)
     • Mauchly’s sphericity test (repeated measures ANOVA)
3. Independence
   • REGRESSION
     • Durbin-Watson: residuals are correlated.
     • Multicollinearity: predictors are highly correlated. Test with VIF and/or correlation matrix

---

Linear Regression

Does money spent on advertising predict album sales?

Here is the dataset we’ll be using in addition to Madalina’s dataset.

A record company wants to know if they spend $100,000 on advertising, how many albums could they expect to sell?

• adverts: 1 unit = $1,000
• sales: albums sold

head(album_df)

##    adverts sales airplay attract
## 1   10.256   330      43      10
## 2  985.685   120      28       7
## 3 1445.563   360      35       7
## 4 1188.193   270      33       7
## 5  574.513   220      44       5
## 6  568.954   170      19       5

# Create the figure
fig <- ggplot(data = album_df, aes(x = adverts, y = sales)) + 
  geom_point(alpha = 0.3, color = "#C06C84") +
  geom_smooth(method = "lm", formula = y ~ x, color = "#7D0552", alpha = .3, size = 1) +
  labs(x = "Advertisements", y = "Sales") +
  theme_bw() +
  labs(title = "Albums Sold Predicted by Money Spent on Advertisements")

# Print the figure
print(fig)

# Save figure
# ggsave("figure.jpg", dpi = 300, width = 5, height = 4, units = "in", type = "cairo")

# Sometimes I find scientific notation super confusing. 
options(scipen = 999)

# Run a linear regression (OLS) using adverts to predict sales
model_lm <- lm(sales ~ adverts, data = album_df)
summary(model_lm)

## 
## Call:
## lm(formula = sales ~ adverts, data = album_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -152.949  -43.796   -0.393   37.040  211.866 
## 
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 134.139938   7.536575  17.799 <0.0000000000000002 ***
## adverts       0.096124   0.009632   9.979 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.99 on 198 degrees of freedom
## Multiple R-squared:  0.3346, Adjusted R-squared:  0.3313 
## F-statistic: 99.59 on 1 and 198 DF,  p-value: < 0.00000000000000022

How do we interpret the output?

1. How do we interpret the estimates? How about the intercept?

2. What does it mean if R-squared is .33? What does it mean for the remaining .67?

Because there is only one predictor, this value represents the square of the simple correlation between advertising and album sales – we can find the square root of R2 by running:

sqrt(0.3346)

## [1] 0.5784462

So, a 1-unit increase in advertising predicts an increase in .096 albums sold. In other words, for every $1,000 we spend on advertising, we predict we’d sell an extra 96 albums.

---

What if we spend $100,000?

Remember that the equation for linear regression is:

\[Y_i = b_0 + b_1X_i + e_i\]

If we apply it to our data, the equation becomes:

\[album\ sales_i = b_0 + b_1advertising\ budget_i\] \[album\ sales_i = 134.14 + (.096*advertising\ budget)\] \[= 134.14 + (.096*100)\] \[= 143.74\]

In other words, if they spend $100,000, we predict they’d sell around 144,000 albums

---

Multiple regression

Let’s go over those assumptions in more detail:

1. Multicollinearity: there should be no perfect linear relationship between two or more predictors. This is an issue because if the predictors are related, then we can’t tease apart their unique contributions
   1. Otherwise, you may have:
      1. Untrustworthy coefficients: once the variance is accounted for by the first predictor, the second predictor accounts for little unique variance.
2. Homoscedasticity: at each level of the predictor, the variance of the residuals should be constant
3. Independent errors: for any two observations, the residual terms should be uncorrelated (use Durbin-Watson test)

---

Here, we’re going to include the variable airplay, which is how much the songs have been played on the radio and attract, which is the attractiveness of the band.

# Run a multiple linear regression (OLS)
model_mult <- lm(sales ~ adverts + airplay + attract, data = album_df)

# OR

model_mult <- update(model_lm, .~. + airplay + attract) # update prior models instead of writing them out again.

Let’s compare the simple regression to the multiple regression

summary(model_lm)

## 
## Call:
## lm(formula = sales ~ adverts, data = album_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -152.949  -43.796   -0.393   37.040  211.866 
## 
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 134.139938   7.536575  17.799 <0.0000000000000002 ***
## adverts       0.096124   0.009632   9.979 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.99 on 198 degrees of freedom
## Multiple R-squared:  0.3346, Adjusted R-squared:  0.3313 
## F-statistic: 99.59 on 1 and 198 DF,  p-value: < 0.00000000000000022

summary(model_mult)

## 
## Call:
## lm(formula = sales ~ adverts + airplay + attract, data = album_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -121.324  -28.336   -0.451   28.967  144.132 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -26.612958  17.350001  -1.534                0.127    
## adverts       0.084885   0.006923  12.261 < 0.0000000000000002 ***
## airplay       3.367425   0.277771  12.123 < 0.0000000000000002 ***
## attract      11.086335   2.437849   4.548           0.00000949 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.09 on 196 degrees of freedom
## Multiple R-squared:  0.6647, Adjusted R-squared:  0.6595 
## F-statistic: 129.5 on 3 and 196 DF,  p-value: < 0.00000000000000022

R-squared: how much of the variance in the DV is accounted for by the regression model

Adjusted r-squared: how much of the variance in the DV would be accounted for if the model had been taken from the population (conservative version)

1. What do you think of the model fit?
2. How would you interpret the estimates for each predictor?

---

Which model is better?

The big problem with R2 is that when you add more variables to the model, it will always go up. If you are deciding which of two models fits the data better, the model with more predictor variables in will always fit better. The Akaike information criterion (AIC) is a measure of fit which penalizes the model for having more variables – a little like adjusted R2.

AIC(model_lm, model_mult)

##            df      AIC
## model_lm    3 2247.375
## model_mult  5 2114.337

anova(model_lm, model_mult)

## Analysis of Variance Table
## 
## Model 1: sales ~ adverts
## Model 2: sales ~ adverts + airplay + attract
##   Res.Df    RSS Df Sum of Sq      F                Pr(>F)    
## 1    198 862264                                              
## 2    196 434575  2    427690 96.447 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Standardized betas

These estimates tell us the number of SDs by which the outcome will change as a result of one SD change in the predictor. This makes it easier to compare each effect to each other because they are not dependent on the units of measurements of the variables (e.g., 1 unit adverts = $1k).

For reference, here are the standard deviations of our predictors:

sd_table <- function(v1,v2,v3) {
  sd1 <- sd(v1)
  sd2 <- sd(v2)
  sd3 <- sd(v3)
  
  print(c(sd1, sd2, sd3))}

sd_table(album_df$adverts, album_df$airplay, album_df$attract)

## [1] 485.65521  12.26958   1.39529

# standardizes the betas. do NOT use this with interaction terms.
lm.beta(model_mult)

##   adverts   airplay   attract 
## 0.5108462 0.5119881 0.1916834

• adverts: as the advertising budget increases by 1 SD ($485,655), album sales increase by 0.511 SD.
• airplay: as the number of radio plays increases by 1 SD (12.27), album sales increase by 0.512 SD.
• attract: a band rated 1 SD (1.40 units) higher on attractiveness can expect an additioinal 0.192 SD in sales.

Diagnostics

Durbin-Watson: testing the assumption of independence

durbinWatsonTest(model_mult)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.0026951      1.949819   0.702
##  Alternative hypothesis: rho != 0

VIF: testing the assumption of no multicollinearity

VIF is variance inflation factor and used to measure multicollinearity.

tolerance = 1/VIF

• If the largest VIF is greater 10, then there’s cause for concern
• If the average VIF is substantially greater than 1, the model may be biased
• Tolerance below 0.1 means there’s a serious problem
• Tolerance below 0.2 means there’s a potential problem

Otherwise, there is no collinearity

# vif
vif(model_mult)

##  adverts  airplay  attract 
## 1.014593 1.042504 1.038455

# average vif
mean(vif(model_mult))

## [1] 1.03185

# tolerance
1/vif(model_mult)

##   adverts   airplay   attract 
## 0.9856172 0.9592287 0.9629695

Using plots for diagnostics

Residuals vs Fitted: is used to check the assumptions of linearity. If the residuals are spread equally around a horizontal line without distinct patterns (red line is approximately horizontal at zero), that is a good indication of having a linear relationship.

Normal Q-Q: is used to check the normality of residuals assumption. If the majority of the residuals follow the straight dashed line, then the assumption is fulfilled.

Scale-Location: is used to check the homoscedasticity of residuals (equal variance of residuals). If the residuals are spread randomly and the see a horizontal line with equally (randomly) spread points, then the assumption is fulfilled.

Residuals vs Leverage: is used to identify any influential value in our dataset. Influential values are extreme values that might influence the regression results when included or excluded from the analysis. Look for cases outside of a dashed line (if there is a concerning point, you’ll see red lines pop up).

plot(model_mult)

hist(model_mult$residuals)

shapiro.test(album_df$sales)

## 
##  Shapiro-Wilk normality test
## 
## data:  album_df$sales
## W = 0.98479, p-value = 0.02965

Added variable plots

These plots tell us whether the variables we’re including in the model are relevant to the model. If the slope is not close to 0 then we conclude that the variable is relevant to the model.

The x-axis represents the partial residuals of the predictor variable, and the y-axis represents the partial residuals of the response variable. The line in the plot represents the fitted values from a linear regression model of the partial residuals of the response variable on the partial residuals of the predictor variable.

avPlots(model_mult)

model_cor <- cor(album_df %>% 
      dplyr::select(adverts, airplay, attract))

corrplot::corrplot(model_cor)

What if the assumptions are violated?

• Violation of linearity: incorporate nonlinear transformation (e.g., log transformation); polynomial terms
• Violation of independence: use multilevel modeling
• Violation of homoscedasticity: transform the DV
• Violation of normality: transform the DV

Transformations

What taking the log of a number does is transform it, essentially shrinking it and making growth at higher numbers non-linear. When you take the log of a series of numbers, it can make them linear if they were formerly non-linear.

Transforming the data won’t change the relationships between variables (the relative differences between people for a given variable stay the same), but it does change the differences between different variables (because it changes the units of measurement).

*Best for data where the 0’s are naturally occurring and not due to error or recording negative values as zero.

Problem	Transformation	Considerations
positive skew, unequal variances	log, sqrt	for log, add constant if working with 0’s or negative numbers
negative skew	reverse score	remember to reverse your reverse when interpreting

skew_df <- read.delim("/Users/kareenadelrosario/Downloads/DownloadFestival(No Outlier).dat")

Day 1

hist.day1 <- ggplot(skew_df, aes(day1)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black", fill = "white", bins = 40) + 
  labs(x = "Hygiene score on day 1", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = mean(skew_df$day1, na.rm = TRUE), sd = sd(skew_df$day1, na.rm = TRUE)), colour = "black", size = 1)

print(hist.day1)

ggplot(skew_df, aes(sample = day1)) +
  stat_qq() +
  stat_qq_line(colour = "red") + # Add a reference line in red
  theme_minimal() + # Use a minimal theme for aesthetics
  ggtitle("Q-Q Plot")

Day 2

hist.day2 <- ggplot(skew_df, aes(day2)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black", fill = "white", bins = 40) + 
  labs(x = "Hygiene score on day 2", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = mean(skew_df$day2, na.rm = TRUE), sd = sd(skew_df$day2, na.rm = TRUE)), colour = "black", size = 1)

print(hist.day2)

ggplot(skew_df, aes(sample = day2)) +
  stat_qq() +
  stat_qq_line(colour = "red") + # Add a reference line in red
  theme_minimal() + # Use a minimal theme for aesthetics
  ggtitle("Q-Q Plot")

Skew & kurtosis should be near 0.

• SKEW: positive score = right skew, negative values = left skew
• KURTOSIS: positive score = fat tail, negative values = flat tail

pastecs::stat.desc(skew_df[, c("day1", "day2")], basic = FALSE, norm = TRUE)

##                     day1               day2
## median        1.79000000 0.7900000000000000
## mean          1.77113580 0.9609090909090909
## SE.mean       0.02436847 0.0443609457146470
## CI.mean.0.95  0.04783289 0.0873478099781075
## var           0.48099624 0.5195238852402352
## std.dev       0.69353892 0.7207800533035270
## coef.var      0.39157862 0.7501022314417026
## skewness     -0.00442835 1.0828112084309722
## skew.2SE     -0.02577395 3.6115738298309901
## kurtosis     -0.42159405 0.7554615276645142
## kurt.2SE     -1.22838457 1.2645083284699845
## normtest.W    0.99591522 0.9083191225801341
## normtest.p    0.03198482 0.0000000000128163

# log transform (add constant b/c we have zeros)
skew_df$day1_log <- log(skew_df$day1 + 1)
skew_df$day2_log <- log(skew_df$day2 + 1)

# sqrt transformations
skew_df$day1_sqrt <- sqrt(skew_df$day1)
skew_df$day2_sqrt <- sqrt(skew_df$day2)

hist.day2_log <- ggplot(skew_df, aes(day2_log)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black", fill = "white", bins = 40) + 
  labs(x = "Hygiene score on day 2", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = mean(skew_df$day2_log, na.rm = TRUE), sd = sd(skew_df$day2_log, na.rm = TRUE)), colour = "black", size = 1)

print(hist.day2_log)

hist.day2_sqrt <- ggplot(skew_df, aes(day2_sqrt)) +
  geom_histogram(aes(y = ..density..), 
                 colour = "black", fill = "white", bins = 40) + 
  labs(x = "Hygiene score on day 2", y = "Density") +
  stat_function(fun = dnorm, args = list(mean = mean(skew_df$day2_sqrt, na.rm = TRUE), sd = sd(skew_df$day2_sqrt, na.rm = TRUE)), colour = "black", size = 1)

print(hist.day2_sqrt)

Considerations

Sometimes transformations create more problems than they solve (not shown, but day 1 scores were normal, but have a slight skew after the transformations). What are your other options for non-normal data?

• Non-parametric tests
• Bootstrapping
• For zero-inflated count data, we can use poisson or negative binomial models

---

Replicating Madalina’s Python code

df <- read_xlsx("/Users/kareenadelrosario/Downloads/data (1).xlsx")

# Create a simple regression plot
ggplot(df, aes(x = Age, y = BELIEFcc)) +
  geom_point(color = "#C06C84", alpha = 0.3) + # Scatter plot
  geom_smooth(method = "lm", color = "#7D0552", size = 1) + # Regression line
  labs(x = "Age", y = "Belief in science") + # Label the x and y axis
  theme_minimal() + # Use a minimal theme
  theme(legend.position = "none") # Remove legend

# Save figure
# ggsave("figure.tif", width = 5, height = 4, dpi = 300, type = "cairo-png")

md <- lm(BELIEFcc ~ Age, data = df)

summary(md)

## 
## Call:
## lm(formula = BELIEFcc ~ Age, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.718 -14.050   6.061  20.728  28.975 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 76.41025    5.68958  13.430 <0.0000000000000002 ***
## Age         -0.07277    0.10968  -0.664               0.508    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.68 on 219 degrees of freedom
## Multiple R-squared:  0.002006,   Adjusted R-squared:  -0.002551 
## F-statistic: 0.4402 on 1 and 219 DF,  p-value: 0.5077

df_plot <- df %>% 
  mutate(Gender.f = case_when(
    Gender == 1 ~ "Female",
    Gender == 2 ~ "Male")) %>%
  mutate(Gender.f = factor(Gender.f, levels = c("Female", "Male"))) %>% 
  filter(!is.na(Gender.f))

ggplot(df_plot, aes(x=Gender.f, y=BELIEFcc, fill = Gender.f)) +
  geom_bar(stat="summary", position = "dodge") +
  ylim(0, 100) +
  theme_minimal() + # Optional: Adds a minimal theme for aesthetics
  labs(x="Party", y="Trust") # Labeling axes

Numeric vs Categorical Predictors

# Numeric
df$Gender <- as.numeric(df$Gender)
model <- lm(BELIEFcc ~ Age + Gender, data = df)
summary(model) 

## 
## Call:
## lm(formula = BELIEFcc ~ Age + Gender, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -74.366 -12.143   5.994  21.255  33.830 
## 
## Coefficients:
##             Estimate Std. Error t value          Pr(>|t|)    
## (Intercept) 60.79657    7.74482   7.850 0.000000000000186 ***
## Age         -0.07665    0.10786  -0.711           0.47808    
## Gender       9.70256    3.32806   2.915           0.00392 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.25 on 218 degrees of freedom
## Multiple R-squared:  0.03946,    Adjusted R-squared:  0.03064 
## F-statistic: 4.477 on 2 and 218 DF,  p-value: 0.01243

# Factor
df$Gender <- factor(df$Gender)
model <- lm(BELIEFcc ~ Age + Gender, data = df)
summary(model) 

## 
## Call:
## lm(formula = BELIEFcc ~ Age + Gender, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -74.85 -12.13   6.14  20.93  34.36 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 70.62086    5.96097  11.847 < 0.0000000000000002 ***
## Age         -0.08698    0.10909  -0.797              0.42615    
## Gender2     10.45065    3.51412   2.974              0.00327 ** 
## Gender4     13.28167   25.63469   0.518              0.60491    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.28 on 217 degrees of freedom
## Multiple R-squared:  0.04144,    Adjusted R-squared:  0.02819 
## F-statistic: 3.127 on 3 and 217 DF,  p-value: 0.02668

What if I want to change the reference group? How can I do that?