Show figure code
df <- bayestestR::simulate_correlation(n=500, r=0.7)
df |>
ggplot(aes(x=V1, y=V2)) +
geom_point() +
geom_smooth(method="lm", formula = 'y ~ x', se=FALSE, linewidth=2, color="red") +
theme_bw()Bayesian Statistics
2. On Regression Parameters in a Frequentist Framework
Recap
Solution
Quizz
- What is wrong with this Correlation plot?
- It is more a “regression” plot (scatter plot + regression line)
Solution
- A correlation is a measure of the strength of association [-1, 1]
- Better represented by an elipsis
- The line is a regression line
- The angle of the line is the regression coefficient (can be > 1)
Show figure code
df |>
ggplot(aes(x=V1, y=V2)) +
geom_point() +
geom_smooth(method="lm", formula = 'y ~ x', se=FALSE, linewidth=2, color="red") +
stat_ellipse(type = "norm", color="blue", linewidth=2) +
theme_bw()
Linear Model
- Linear function: \(y = ax + b\)
- Regression formula: \(\hat{y} = Intercept + \beta_{1}x\)
- \(\hat{y}\) (Y-hat): predicted response/outcome/dependent variable
- \(Intercept\) (\(\beta_{0}\)): “starting point”. Value of \(\hat{y}\) when \(x=0\)
- \(\beta_{1}\) (beta) : slope/“effect”. Change of \(\hat{y}\) from the intercept when \(x\) increases by 1
- \(x\): predictor/independent variable
- What are the parameters of the regression line below? Intercept = 1 and beta = 3
Code
ggplot(df, aes(x=V1, y=V2)) +
geom_point(alpha=0) +
geom_vline(xintercept = 0, linetype="dotted") +
geom_abline(intercept = 1, slope = 3, color="red", linewidth=2) +
geom_segment(aes(x = 0, y = 1, xend = 1, yend = 1), linewidth=1,
color="green", linetype="dashed") +
geom_segment(aes(x = 1, y = 1, xend = 1, yend = 4), linewidth=1,
color="green", linetype="solid", arrow=arrow(length = unit(0.1, "inches"))) +
geom_segment(aes(x = 0, y = 0, xend = 0, yend = 1), linewidth=1,
color="blue", linetype="solid", arrow=arrow(length = unit(0.1, "inches"))) +
geom_point(aes(x = 0, y = 0), color="purple", size=8, shape="+") +
labs(x="x", y="y") +
theme_bw() +
coord_cartesian(xlim = c(-1, 1.5), ylim = c(-3, 4))
Impact of Parameters
- What happens when we change the slope and the intercept?
Ordinary Least Squares (OLS)
- The Ordinary Least Squares (OLS) line is the line that minimizes the sum of squared residuals (Residual Sum of Squares - RSS)
- \(RSS = \sum(residuals^{2})\) (the square is to make them positive)
- The residuals are the differences between observed and predicted values
- \(RSS = \sum_{i=1}^{n_{obs}}(y_{i} - \hat{y}_{i})^{2}\)
- Squared residuals correspond to Errors (epsilon \(\varepsilon\))
Code
p <- mtcars |>
ggplot(aes(x=qsec, y=mpg)) +
geom_point(alpha=0.7, size=6) +
geom_smooth(method="lm", formula = 'y ~ x', se=FALSE, linewidth=2) +
geom_segment(aes(x = qsec, y = mpg, xend = qsec, yend = predict(lm(mpg ~ qsec, data=mtcars))),
color="red", linetype="dotted", linewidth=1) +
theme_bw() +
labs(x="x", y="y", title="The residuals are the vertical distances between each point and the line.")
p
Maximum Likelihood Estimation (MLE)
- The OLS “line” (i.e., the OLS coefficients) can be computed using various algorithms, but one is particularly useful
- This approach is called Maximum Likelihood Estimation (MLE)
- In a typical linear regression model, the assumption is made that the errors are normally distributed1
- Why? This “trick” allows to get OLS-equivalent coefficients and further enable inferential statistics
- MLE tries to maximize the likelihood of the data given the model2, by finding the best distribution parameters to “fit” the residuals
Ronald Fisher (1890 – 1962), “a genius who almost single-handedly created the foundations for modern statistical science”, introduced MLE in 1922
Linear Model = “Gaussian” Model
- Solving the linear model equation means finding the best fitting normal distribution around the residuals
- In other words, a typical “linear model” is actually a “Gaussian” model
Code
p <- mtcars |>
ggplot(aes(x=qsec, y=mpg)) +
geom_point(alpha=0.7, size=6) +
geom_smooth(method="lm", formula = 'y ~ x', se=FALSE, linewidth=2) +
geom_segment(aes(x = qsec, y = mpg, xend = qsec, yend = predict(lm(mpg ~ qsec, data=mtcars))),
color="red", linetype="dotted", linewidth=1) +
theme_bw() +
labs(x="x", y="y", title="The residuals are the vertical distances between each point and the line.")
p2 <- data.frame(Error=insight::get_residuals(lm(mpg ~ qsec, data=mtcars))) |>
ggplot(aes(x=Error)) +
geom_histogram(bins=10, fill="grey", color="black") +
geom_vline(xintercept = 0, linetype="dashed") +
geom_density(data=data.frame(Error=bayestestR::distribution_normal(n=100, sd = 2)),
aes(y=after_stat(density)*40), color="#F44336", linewidth=1, adjust=1) +
geom_density(data=data.frame(Error=bayestestR::distribution_normal(n=100, sd = 3)),
aes(y=after_stat(density)*50), color="#FF5722", linewidth=1, adjust=1) +
geom_density(data=data.frame(Error=bayestestR::distribution_normal(n=100, sd = 4)),
aes(y=after_stat(density)*60), color="#FF9800", linewidth=1, adjust=1) +
geom_point(aes(y=0), size=10, shape=16, alpha=0.3) +
theme_bw() +
coord_flip() +
labs(y = "Density")
p + theme(plot.title = element_blank()) | p2
Homoscedasticity
- The line that goes through the Normal distribution’s locations \(\mu\) also minimizes the sum of squared residuals (OLS)
Code
model <- lm(mpg ~ qsec, data=mtcars)
p <- mtcars |>
ggplot(aes(x=qsec, y=mpg)) +
geom_smooth(method="lm", formula = 'y ~ x', se=FALSE, linewidth=2) +
theme_bw()
# Function to add normal distribution curves
add_normals <- function(p, model) {
sigma <- summary(model)$sigma # Standard deviation of residuals
n <- 100 # Number of points for each curve
for(i in 1:nrow(mtcars)) {
x_val <- mtcars$qsec[i]
y_pred <- predict(model, newdata = data.frame(qsec = x_val))
# Create a sequence of y values for the normal curve
y_seq <- seq(y_pred - 3*sigma, y_pred + 3*sigma, length.out = n)
density <- dnorm(y_seq, mean = y_pred, sd = sigma)
# Adjust density to match the scale of the plot
max_width <- 1 # Max width of areas
density_scaled <- (density / max(density)) * max_width
# Create a dataframe for each path
path_df <- data.frame(x = x_val + density_scaled, y = y_seq)
path_dfv <- data.frame(x=path_df$x[1], ymin=min(path_df$y), ymax=max(path_df$y))
# Add the path to the plot
p <- p +
geom_segment(data = path_dfv, aes(x = x, xend=x, y = ymin, yend=ymax),
color = "#FF9800", linewidth = 0.7, alpha=0.8, linetype="dotted") +
geom_path(data = path_df, aes(x = x, y = y),
color = "#FF9800", linewidth = 1, alpha=0.8)
}
p
}
# Create the final plot
p <- add_normals(p, model) +
geom_segment(aes(x = qsec, y = mpg, xend = qsec, yend = predict(lm(mpg ~ qsec, data=mtcars))),
color="red", linetype="solid", linewidth=1) +
geom_point(alpha=0.8, size=6)
p
Linear Models - Deviation
- Once we have estimated the distribution of errors:
- We can obtain the regression line (the slope of the distribution’s locations \(\mu\))
- We can obtain the spread of the distribution, i.e., the standard deviation of the errors \(\sigma\) (sigma)1
- The standard deviation of the errors, also called the residual standard error, and is computed as:
- \(\sigma = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y_i})^2}{n_{obs} - k}}\)
- \(\sigma = \sqrt{\frac{RSS}{n_{obs} - k}}\)
- \(k\) is the number of predictors in the model (excluding the intercept)
- N observations - k predictors = degrees of freedom (df)
- \(\sigma = \sqrt{\frac{RSS}{df}}\)
The t-value
- The t-value shows how likely our coefficient compared to the assumed distributions of possible coefficient under the null hypothesis (i.e., if there is an absence of effect)
Code
# Plot t-distribution
x <- seq(-5, 5, length.out = 1000)
y <- dt(x, df=30)
df <- data.frame(x = x, y = y)
t_value <- insight::get_statistic(model)[2, 2]
df |>
ggplot(aes(x = x, y = y)) +
geom_area(color="black", fill="grey") +
geom_segment(aes(x = t_value, y = 0, xend = t_value, yend = dt(t_value, df=30)),
color="red", linetype="solid", linewidth=1) +
geom_point(aes(x = t_value, y = dt(t_value, df=30)), color="red", size=3) +
theme_minimal() +
labs(x = "\nt-values - standardized coefficient under the null hypothesis", y = "Probability",
title = paste0("t-Distribution (df=30); ", "t-value = ", round(t_value, 2), "\n"))
Is our coefficient precisely different from 0?
- Is our coefficient likely to be observed under the null hypothesis?
- What is the probability to obtain a coefficient at least as precise as ours (in both directions) under the null hypothesis?
- \(Prob > |t|\)
Code
df$Probability <- ifelse(df$x < -t_value, "< -t", "Smaller")
df$Probability <- ifelse(df$x > t_value, "> +t", df$Probability)
df |>
ggplot(aes(x = x, y = y)) +
geom_line() +
geom_area(aes(x = x, y = y, fill = Probability), alpha = 0.5) +
geom_segment(aes(x = t_value, y = 0, xend = t_value, yend = dt(t_value, df=30)),
color="red", linetype="solid", linewidth=1) +
geom_point(aes(x = t_value, y = dt(t_value, df=30)), color="red", size=3) +
theme_minimal() +
scale_fill_manual(values = c("red", "red", "grey")) +
labs(x = "\nt-values - standardized coefficient under the null hypothesis", y = "Probability",
title = paste0("t-Distribution (df=30); ", "t-value = ", round(t_value, 2), "\n"))
Confidence Intervals (CI)
- Confidence Intervals are often seen as a “range” of likely values for the coefficient, and the more we would get close to the coefficient the more “likely” the values would be. It is incorrect
- The confidence interval (CI) is the range of values that contains, e.g., with 95% of probability, the value of the coefficient if the t-distribution was centered around the t-value (as if the null effect was the estimated coefficient).
- The interpretation is, again, fairly convoluted
Code
# Get 95% ETI of the t distribution
ci <- qt(c(0.025, 0.975), df=30, lower.tail = TRUE) + t_value
# ci * parameters::parameters(model)$SE[2] # Indeed shows the right CI
df |>
ggplot(aes(x = x, y = y)) +
geom_area(fill="grey", alpha=0.5) +
geom_area(aes(x = x + t_value, y = y), color="red", fill="red", alpha=0.5) +
geom_segment(aes(x = t_value, y = 0, xend = t_value, yend = dt(t_value, df=30)),
color="red", linetype="solid", linewidth=1) +
geom_point(aes(x = t_value, y = dt(t_value, df=30)), color="red", size=3) +
# Horizontal segment indicating the CI range
geom_segment(aes(x = ci[1], y = 0.05, xend = ci[2], yend = 0.05),
color="orange", linetype="solid", linewidth=1, arrow = arrow(ends='both')) +
geom_vline(xintercept = ci[1], color="orange", linetype="dotted", linewidth=0.5) +
geom_vline(xintercept = ci[2], color="orange", linetype="dotted", linewidth=0.5) +
theme_minimal() +
labs(x = "\nt-values - standardized coefficient under the null hypothesis", y = "Probability",
title = paste0("Shifted t-Distribution (df=30); ", "location = ", round(t_value, 2), "(t-value)\n"))
The End (for now)
Thank you!
