What do we already know?

• test association between two categorical variables (chi-square)
• find relationships between categorical and continuous variables (ANOVA, t-test)
• test relationships between two continuous variables (correlations)

Why do we need regression?

• it is an umbrella method for all of the above
• more than one predictor variable is allowed
• relationship can be expressed as a number (e.g., every extra year of work experience increases the expected salary by £50)

Omitted variable effect (1)

Omitted variable effect (2)

Regression modelling

• describes relationships between variables with the following formula: \[y_i = a + b*x_i + e_i\]
  • a - intercept. This is the value of Y when X = 0.
  • b - slope. The is a regression coefficient showing the change in Y when X increases by 1.

Toy data example

There are variables X and Y, n = 50, normally distributed

##        y                x          
##  Min.   :-8.146   Min.   :-2.7091  
##  1st Qu.:-3.083   1st Qu.:-0.9292  
##  Median :-1.367   Median :-0.3983  
##  Mean   :-1.138   Mean   :-0.3787  
##  3rd Qu.: 1.265   3rd Qu.: 0.2604  
##  Max.   : 7.922   Max.   : 2.4979

Visualise their relationship

plot(x, y, pch = 16)

Use correlation to estimate the strength of their relationship

cor.test(x, y, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 21.127, df = 48, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9134637 0.9715847
## sample estimates:
##       cor 
## 0.9502104

Estimate the model

model1 <- lm(y ~ x)
summary(model1)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.39265 -0.72426  0.04457  0.70943  1.95334 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.009131   0.153084   -0.06    0.953    
## x            2.980894   0.141096   21.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.014 on 48 degrees of freedom
## Multiple R-squared:  0.9029, Adjusted R-squared:  0.9009 
## F-statistic: 446.3 on 1 and 48 DF,  p-value: < 2.2e-16

What are the conclusions?

• The value of Y when X = 0, e.g. expected salary when work experience is 0 years.
• How much will Y change when X increases by 1?

Plot

Coefficient of determination R²

• shows the relation of variance explained by the model to total variance of the outcome variable \[R^2 = \frac{SS_{reg}}{SS_{tot}}\]
• varies from 0 to 1 (or between 0-100 per cent)
• R² above 0.1 demonstrates a satisfactory predictive power (but that is not a hard fact)
• as the number of predictors grows, R² grows as well. Therefore,
• adjusted R²: \[R^2_{adj} = \frac{SS_{reg} / (n – k))}{(SS_{tot} / (n – 1))}\]
• where n is the sample size, and k is the number of predictors

Other ways to evaluate model fit:

• AIC — Akaike information criterion \[AIC = \frac{2k}{n} + ln \frac{SS_{res}}{n}\]
• BIC — Bayesian information criterion \[BIC = \frac{k*ln(n)}{n} + ln \frac{SS_{res}}{n}\]
• the smaller the value (the closer it is to 0), the better
• suitable for model comparison only

Several predictors

• The main advantage of regression modelling is that it can digest several predictors in one model

Dummy-variables

• dichotomised
• show the differences in intercepts between groups (there is a ‘reference group’)
• if class(var) = factor in R, it is recoded into dummies by default

Example: salary of males and females

200 males and females

##       age            salary           sex           
##  Min.   :17.00   Min.   : 49.63   Length:200        
##  1st Qu.:26.75   1st Qu.: 84.08   Class :character  
##  Median :30.00   Median : 94.42   Mode  :character  
##  Mean   :30.09   Mean   : 95.50                     
##  3rd Qu.:34.25   3rd Qu.:106.66                     
##  Max.   :46.00   Max.   :146.75

Example: exploratory plots

Example: the model

model2 <- lm(salary ~ age + sex, data = genderdata)

tab_model(model2, show.ci = F)

Example: plot

What if there are more than 2 categories? (e.g. educational attainment)

• create n-1 variables where n is the original number of categories
• the first category is now the ‘reference’
• What does the coding look like?

##        educ dummy1 dummy2
## 1  tertiary      0      0
## 2 secondary      1      0
## 3   primary      0      1

Another cool example

Artwork by @allison_horst

Interaction effect

• the model in the previous example implies that the difference between males and females is constant
• What if this gender gap grows with age?

Model Comparison

model3 <- lm(salary ~ age + sex, data = genderdata2)
model4 <- lm(salary ~ age * sex, data = genderdata2)
tab_model(model3, model4, show.ci = F)

What the Interaction Effect Model Looks Like

What happens if we ignore the interaction effect?

Can we identify a possible interaction effect visually?

Sometimes yes. Heteroskedasticity may indicate a moderation/ interaction. Let’s recall the assumptions first.

Linear Regression Assumptions

• linearity
• normally distributed residuals of the outcome
• predictors have an unlimited range
• independent observations (no correlation between their residuals)
• independent variables (no multicollinearity)
• no outliers
• homoskedasticity (residual variances are the same along the values of Y)

Homoskedasticity/Heteroskedasticity

x <- 1:100
y1 <- rnorm(n = 100, mean = x, sd = 10)
y2 <- rnorm(n = 100, mean = x, sd = 0.4*x)
par(mfrow = c(1, 2))
plot(x, y1, pch = 16); abline(lm(y1 ~ x), col = "red")
plot(x, y2, pch = 16); abline(lm(y2 ~ x), col = "red")

Pane on the right-hand side with data points ‘fanning out’ shows there is a third variable which comes into play at high values of X

Anscombe’s Quartet

Summary: nuts and bolts of regression modelling

• Input
  • outcome - 1 pc.
  • predictors - >= 1 pc.
• Output
  • regression coefficients
  • R²

What’s next?

• more on interaction effects
• model diagnostics (next year)
• logistic regression (after that)