Unit 7 (Part Deux) - ANOVA & regression

PSYC 640 - Fall 2024

Dustin Haraden, PhD

Reminders

• Journal Entries
• Mid-course Feedback
• Presentation 2 Groups

Last Class

• ANOVA
  • Some theory, but more applied
  • easystats::report() makes things easier
  • Everyone getting better with R!

Today…

ANOVA

Regression!!

# File management
library(here)
# for dplyr, ggplot2
library(tidyverse)
#Loading data
library(rio)
# Estimating Marginal Means
library(emmeans)
# Pretty Tables
library(kableExtra)
# Pretty variable names
library(janitor)

#Remove Scientific Notation 
options(scipen=999)

Review of ANOVA

One-Way ANOVA

Goal: Inform of differences among the levels of our variable of interest (Omnibus Test)

Using the between and within group variance to create the \(F\)-statistic/ratio

Hypotheses:

\[ H_0: it\: is\: true\: that\: \mu_1 = \mu_2 = \mu_3 =\: ...\mu_k \\ H_1: it\: is\: \boldsymbol{not}\: true\: that\: \mu_1 = \mu_2 = \mu_3 =\: ...\mu_k \]

What is a Two-Way ANOVA?

Examines the impact of 2 nominal/categorical variables on a continuous outcome

We can now examine:

• The impact of variable 1 on the outcome (Main Effect)

• The impact of variable 2 on the outcome (Main Effect)

• The interaction of variable 1 & 2 on the outcome (Interaction Effect)

  • The effect of variable 1 depends on the level of variable 2

Main Effect & Interactions

Main Effect: Basically a one-way ANOVA

• The effect of variable 1 is the same across all levels of variable 2

Interaction:

• Able to examine the effect of variable 1 across different levels of variable 2

• Basically speaking, the effect of variable 1 on our outcome DEPENDS on the levels of variable 2

Ghost Data

Since there are 2 header rows, we need to only include the first one. To do that, we have to “skip” the first two and then give the names of the columns back to the data.

We also need to specify what the missing values are. Typically we have been working with NA which is more traditional. However, missing values in this dataset are DK/REF and a blank. This will need to be specified in the import function (used this website)

# get the names of your columns which is the first row
ghost_data_names <- read_csv(here("lectures", "data", "ghosts.csv")) %>% 
  names()

# import second time; skip row 2, and assign column names to argument col_names =
ghost_data <- read_csv(here("lectures", "data", "ghosts.csv"),
                       skip = 2,
                       col_names = ghost_data_names,
                       na = c("DK/REF", "", " ")
                      ) %>% 
  clean_names()

Running the Test

Let’s take a look at Income by Political Affiliation

aov1 <- aov(income ~ political_affiliation, 
            data = ghost_data)

summary(aov1)

                       Df        Sum Sq     Mean Sq F value Pr(>F)  
political_affiliation   2   44252069121 22126034561   4.007 0.0189 *
Residuals             393 2170016726333  5521671059                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
604 observations deleted due to missingness

Post-hoc tests

#get basic summary stats
ghost_data %>% 
  group_by(political_affiliation) %>% 
  summarise(mean = mean(income, na.rm = TRUE), 
            sd = sd(income, na.rm = TRUE))

# A tibble: 4 × 3
  political_affiliation    mean     sd
  <chr>                   <dbl>  <dbl>
1 Democrat               80431. 73147.
2 Independent           103919. 87848.
3 Republican             85623. 49064.
4 <NA>                   78233. 51022.

# conduct the comparisons
emmeans(aov1, pairwise ~ political_affiliation, 
        adjust = "none")

$emmeans
 political_affiliation emmean   SE  df lower.CL upper.CL
 Democrat               80431 6520 393    67618    93244
 Independent           103919 5870 393    92369   115468
 Republican             85623 7220 393    71433    99812

Confidence level used: 0.95 

$contrasts
 contrast                 estimate   SE  df t.ratio p.value
 Democrat - Independent     -23488 8770 393  -2.677  0.0077
 Democrat - Republican       -5192 9720 393  -0.534  0.5937
 Independent - Republican    18296 9310 393   1.966  0.0500

Visualizing

ggstatsplot https://indrajeetpatil.github.io/ggstatsplot/

#Welch's F-test
oneway.test(income ~ political_affiliation, 
            data = ghost_data)

    One-way analysis of means (not assuming equal variances)

data:  income and political_affiliation
F = 3.4289, num df = 2.00, denom df = 260.52, p-value = 0.03389

ggstatsplot::ggbetweenstats(ghost_data,
                            political_affiliation, 
                            income)

Regression

Today…

Regression

• Why use regression?

• One equation to rule them all

Today…

Regression

• Why use regression?

• One equation to rule them all

  • Ordinary Least Squares

  • Interpretation

#Don't know if I'm using all of these, but including theme here anyways
library(tidyverse)
library(rio)
library(broom)
library(psych)
library(gapminder)
library(psychTools)


#Remove Scientific Notation 
options(scipen=999)

Overview of Regression

Regression is a general data analytic system

• Lots of things fall under the umbrella of regression

• This system can handle a variety of forms of relations and types of variables

The output of regression includes both effect sizes and statistical significance

We can also incorporate multiple influences (IVs) and account for their intercorrelations

Uses for regression

• Adjustment: Take into account (control) known effects in a relationship

• Prediction: Develop a model based on what has happened previously to predict what will happen in the future

• Explanation: examining the influence of one or more variable on some outcome

Study Design & Collection

Design - When data are collected

• Retrospective/Prospective

• Longitudinal

• Cross-Sectional

Collection - How data are collected

• Experimental

• Field

• Observational

• Meta-analysis

• Neuroimaging/Psychophys

• Survey

• Quasi-Experimental

Reminder to Professor for Model Drawings

maybe?

Regression Equation

With regression, we are building a model that we think best represents the data at hand

At the most simple form we are drawing a line to characterize the linear relationship between the variables so that for any value of x we can have an estimate of y

\[ Y = mX + b \]

• Y = Outcome Variable (DV)

• m = Slope Term

• X = Predictor (IV)

• b = Intercept

Regression Equation

Overall, we are providing a model to give us a “best guess” on predicting

Let’s “science up” the equation a little bit:

\[ Y_i = b_0 + b_1X_i + e_i \]

This equation is capturing how we are able to calculate each observation ( \(Y_i\) )

\[ \hat{Y_i} = b_0 + b_1X_i \]

This one will give us the “best guess” or expected value of \(Y\) given \(X\)

Regression Equation

There are two ways to think about our regression equation. They’re similar to each other, but they produce different outputs.
\[Y_i = b_{0} + b_{1}X_i +e_i\]
\[\hat{Y_i} = b_{0} + b_{1}X_i\]
The model we are building by including new variables is to explain variance in our outcome

\(\hat{Y}\) signifies the fitted score – no error
The difference between the fitted and observed score is the residual (\(e_i\))
There is a different e value for each observation in the dataset

Expected vs. Actual

\[Y_i = b_{0} + b_{1}X_i + e_i\]

\[\hat{Y_i} = b_{0} + b_{1}X_i\]

\(\hat{Y}\) signifies that there is no error. Our line is predicting that exact value. We interpret it as being “on average”

Important to identify that that \(Y_i - \hat{Y_i} = e_i\).

\(\hat{Y}\) signifies the fitted score – no error

The difference between the fitted and observed score is the residual (\(e_i\))

There is a different e value for each observation in the dataset

OLS

• How do we find the regression estimates?

• Ordinary Least Squares (OLS) estimation

• Minimizes deviations

  • \[ min\sum(Y_{i} - \hat{Y} ) ^{2} \]

• Other estimation procedures possible (and necessary in some cases)

set.seed(142)
x.1 <- rnorm(10, 0, 1)
e.1 <- rnorm(10, 0, 2)
y.1 <- .5 + .55 * x.1 + e.1
d.1 <- data.frame(x.1,y.1)
m.1 <- lm(y.1 ~ x.1, data = d.1)
d1.f<- augment(m.1)


ggplot(d1.f , aes(x=x.1, y=y.1)) +
    geom_point(size = 2) +
  geom_smooth(method = lm, se = FALSE) +
  theme_bw(base_size = 20)

ggplot(d1.f , aes(x=x.1, y=y.1)) +
    geom_point(size = 2) +
  geom_point(aes(y = .fitted), shape = 1, size = 2) +
  theme_bw(base_size = 20)

ggplot(d1.f , aes(x=x.1, y=y.1)) +
    geom_point(size = 2) +
  geom_point(aes(y = .fitted), shape = 1, size = 2) +
  geom_segment(aes( xend = x.1, yend = .fitted))+
  theme_bw(base_size = 20)

ggplot(d1.f , aes(x=x.1, y=y.1)) +
    geom_point(size = 2) +
  geom_smooth(method = lm, se = FALSE) +
  geom_point(aes(y = .fitted), shape = 1, size = 2) +
  geom_segment(aes( xend = x.1, yend = .fitted))+
  theme_bw(base_size = 20)

compare to bad fit

OLS

The line that yields the smallest sum of squared deviations

\[\Sigma(Y_i - \hat{Y_i})^2\] \[= \Sigma(Y_i - (b_0+b_{1}X_i))^2\] \[= \Sigma(e_i)^2\]

In order to find the OLS solution, you could try many different coefficients \((b_0 \text{ and } b_{1})\) until you find the one with the smallest sum squared deviation. Luckily, there are simple calculations that will yield the OLS solution every time.

Regression coefficient, \(b_{1}\)

\[b_{1} = \frac{cov_{XY}}{s_{x}^{2}} = r_{xy} \frac{s_{y}}{s_{x}}\]

What units is the regression coefficient in?

The regression coefficient (slope) equals the estimated change in Y for a 1-unit change in X

\[\Large b_{1} = r_{xy} \frac{s_{y}}{s_{x}}\]

If the standard deviation of both X and Y is equal to 1:

\[\Large b_1 = r_{xy} \frac{s_{y}}{s_{x}} = r_{xy} \frac{1}{1} = r_{xy} = \beta_{yx} = b_{yx}^*\]

Standardized regression equation

\[\Large Z_{y_i} = b_{yx}^*Z_{x_i}+e_i\]

\[\Large b_{yx}^* = b_{yx}\frac{s_x}{s_y} = r_{xy}\]

According to this regression equation, when \(X = 0, Y = 0\). Our interpretation of the coefficient is that a one-standard deviation increase in X is associated with a \(b_{yx}^*\) standard deviation increase in Y. Our regression coefficient is equivalent to the correlation coefficient when we have only one predictor in our model.

Estimating the intercept, \(b_0\)

• intercept serves to adjust for differences in means between X and Y

\[\hat{Y_i} = \bar{Y} + r_{xy} \frac{s_{y}}{s_{x}}(X_i-\bar{X})\]

• if standardized, intercept drops out

• otherwise, intercept is where regression line crosses the y-axis at X = 0

• Also, notice that when \(X = \bar{X}\) the regression line goes through \(\bar{Y}\)

\[b_0 = \bar{Y} - b_1\bar{X}\]

The intercept adjusts the location of the regression line to ensure that it runs through the point \(\large (\bar{X}, \bar{Y}).\) We can calculate this value using the equation:

\[\Large b_0 = \bar{Y} - b_1\bar{X}\]

Visuals of OLS Regression

https://setosa.io/ev/ordinary-least-squares-regression/

https://observablehq.com/@yizhe-ang/interactive-visualization-of-linear-regression

Example

library(gapminder)
gapminder = gapminder %>% filter(year == 2007 & continent == "Asia") %>% 
  mutate(log_gdp = log(gdpPercap))
describe(gapminder[,c("log_gdp", "lifeExp")], fast = T)

        vars  n  mean   sd median   min   max range  skew kurtosis   se
log_gdp    1 33  8.74 1.24   8.41  6.85 10.76  3.91  0.21    -1.37 0.22
lifeExp    2 33 70.73 7.96  72.40 43.83 82.60 38.77 -1.07     1.79 1.39

cor(gapminder$log_gdp, gapminder$lifeExp)

[1] 0.8003474

If we regress lifeExp onto log_gdp:

r = cor(gapminder$log_gdp, gapminder$lifeExp)
m_log_gdp = mean(gapminder$log_gdp)
m_lifeExp = mean(gapminder$lifeExp)
s_log_gdp = sd(gapminder$log_gdp)
s_lifeExp = sd(gapminder$lifeExp)

b1 = r*(s_lifeExp/s_log_gdp)

[1] 5.157259

b0 = m_lifeExp - b1*m_log_gdp

[1] 25.65011

How will this change if we regress GDP onto life expectancy?

(b1 = r*(s_lifeExp/s_log_gdp))

[1] 5.157259

(b0 = m_lifeExp - b1*m_log_gdp)

[1] 25.65011

(b1 = r*(s_log_gdp/s_lifeExp))

[1] 0.1242047

(b0 = m_log_gdp - b1*m_lifeExp)

[1] -0.04405086

In R

fit.1 <- lm(lifeExp ~ log_gdp, data = gapminder)
summary(fit.1)

Call:
lm(formula = lifeExp ~ log_gdp, data = gapminder)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.314  -1.650  -0.040   3.428   8.370 

Coefficients:
            Estimate Std. Error t value     Pr(>|t|)    
(Intercept)  25.6501     6.1234   4.189     0.000216 ***
log_gdp       5.1573     0.6939   7.433 0.0000000226 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.851 on 31 degrees of freedom
Multiple R-squared:  0.6406,    Adjusted R-squared:  0.629 
F-statistic: 55.24 on 1 and 31 DF,  p-value: 0.00000002263

Things to discuss

• Coefficient estimates
• Statistical tests (covered in more detail soon)

Data, fitted, and residuals

library(broom)
model_info = augment(fit.1)
head(model_info)

# A tibble: 6 × 8
  lifeExp log_gdp .fitted .resid   .hat .sigma .cooksd .std.resid
    <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>
1    43.8    6.88    61.1 -17.3  0.101    3.63 0.796       -3.76 
2    75.6   10.3     78.8  -3.15 0.0802   4.89 0.0199      -0.676
3    64.1    7.24    63.0   1.08 0.0765   4.93 0.00224      0.233
4    59.7    7.45    64.1  -4.33 0.0646   4.86 0.0294      -0.923
5    73.0    8.51    69.5   3.43 0.0314   4.89 0.00836      0.718
6    82.2   10.6     80.3   1.94 0.100    4.92 0.00994      0.422

describe(model_info, fast = T)

           vars  n  mean   sd median    min   max range  skew kurtosis   se
lifeExp       1 33 70.73 7.96  72.40  43.83 82.60 38.77 -1.07     1.79 1.39
log_gdp       2 33  8.74 1.24   8.41   6.85 10.76  3.91  0.21    -1.37 0.22
.fitted       3 33 70.73 6.37  69.00  60.98 81.16 20.19  0.21    -1.37 1.11
.resid        4 33  0.00 4.77  -0.04 -17.31  8.37 25.68 -1.37     3.29 0.83
.hat          5 33  0.06 0.03   0.05   0.03  0.11  0.08  0.60    -0.98 0.00
.sigma        6 33  4.84 0.23   4.90   3.63  4.93  1.30 -4.53    20.90 0.04
.cooksd       7 33  0.04 0.14   0.01   0.00  0.80  0.80  5.08    25.12 0.02
.std.resid    8 33  0.00 1.02  -0.01  -3.76  1.77  5.53 -1.44     3.56 0.18

Point out the average of the residuals is 0, just like average deviation from the mean is 0.

The relationship between \(X_i\) and \(\hat{Y_i}\)

model_info %>% ggplot(aes(x = log_gdp, y = .fitted)) +
  geom_point() + geom_smooth(se = F, method = "lm") +
  scale_x_continuous("X") + scale_y_continuous(expression(hat(Y))) + theme_bw(base_size = 30)

The relationship between \(X_i\) and \(e_i\)

model_info %>% ggplot(aes(x = log_gdp, y = .resid)) +
  geom_point() + geom_smooth(se = F, method = "lm") + 
  scale_x_continuous("X") + scale_y_continuous("e") + theme_bw(base_size = 30)

The relationship between \(Y_i\) and \(\hat{Y_i}\)

model_info %>% ggplot(aes(x = lifeExp, y = .fitted)) +
  geom_point() + geom_smooth(se = F, method = "lm") + 
  scale_x_continuous("Y") + scale_y_continuous(expression(hat(Y))) + theme_bw(base_size = 30)

The relationship between \(Y_i\) and \(e_i\)

model_info %>% ggplot(aes(x = lifeExp, y = .resid)) +
  geom_point() + geom_smooth(se = F, method = "lm") + 
  scale_x_continuous("Y") + scale_y_continuous("e") + theme_bw(base_size = 25)

The relationship between \(\hat{Y_i}\) and \(e_i\)

model_info %>% ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() + geom_smooth(se = F, method = "lm") + 
  scale_y_continuous("e") + scale_x_continuous(expression(hat(Y))) + theme_bw(base_size = 30)

Using easystats

performance::check_model(fit.1)

Using sjPlot

sjPlot::tab_model(fit.1)

	life Exp
Predictors	Estimates	CI	p
(Intercept)	25.65	13.16 – 38.14	<0.001
log gdp	5.16	3.74 – 6.57	<0.001
Observations	33
R2 / R2 adjusted	0.641 / 0.629

Regression to the mean

An observation about heights was part of the motivation to develop the regression equation: If you selected a parent who was exceptionally tall (or short), their child was almost always not as tall (or as short).

library(psychTools)
library(tidyverse)
heights = psychTools::galton
mod = lm(child~parent, data = heights)
point = 902
heights = broom::augment(mod)


heights %>%
  ggplot(aes(x = parent, y = child)) +
  geom_jitter(alpha = .3) +
  geom_hline(aes(yintercept = 72), color = "red") +
  geom_vline(aes(xintercept = 72), color = "red") +
  theme_bw(base_size = 20)

Regression to the mean

This phenomenon is known as regression to the mean. This describes the phenomenon in which an random variable produces an extreme score on a first measurement, but a lower score on a second measurement.

Regression to the mean

This can be a threat to internal validity if interventions are applied based on first measurement scores.

Example in R

Try out a linear regression on the Ghosts Data!