x <- rnorm(1000, mean=50, sd=10)
MASS::fitdistr(x, "normal") mean sd
49.8828661 9.9825633
( 0.3156764) ( 0.2232169)Bayesian Statistics
2. On Regression Parameters in a Frequentist Framework
Recap
Solution
Fitting distributions to data
- Maximum Likelihood Estimation (MLE) is the act of estimating the parameters of a distribution most fitting the data.
- Typical algorithms that solve it proceed by iteratively adjusting the parameters of the distribution and computing the divergence between the observed and predicted probabilities, and finding the parameters that minimize this divergence (i.e., that maximizes the likelihood).
- They are not trivial algorithms
- We can use the
fitdistr()(“fit distribution”) function from theMASSpackage.
Show figure code
params <- MASS::fitdistr(x, "normal")
data.frame(x=x) |>
ggplot(aes(x=x)) +
geom_histogram(aes(y=after_stat(density)), bins=60, fill="grey", color="black") +
geom_line(
data = data.frame(
x=seq(0, 100, length.out = 500),
y=dnorm(seq(0, 100, length.out = 500), mean=params$estimate["mean"], sd=params$estimate["sd"])),
aes(x=x, y=y), color="red", linewidth=2) +
theme_bw() +
labs(x="x", y="Density")
x <- rgamma(1000, 2, 3)
MASS::fitdistr(x, "normal") mean sd
0.65923460 0.45948465
(0.01453018) (0.01027439)Show figure code
params <- MASS::fitdistr(x, "normal")
data.frame(x=x) |>
ggplot(aes(x=x)) +
geom_histogram(aes(y=after_stat(density)), bins=100, fill="grey", color="black") +
geom_line(
data = data.frame(
x=seq(-1.5, 3.5, length.out = 500),
y=dnorm(seq(-1.5, 3.5, length.out = 500), mean=params$estimate["mean"], sd=params$estimate["sd"])),
aes(x=x, y=y), color="red", linewidth=2) +
theme_bw() +
labs(x="x", y="Density")
Maximum Likelihood: Caveats
- MLE will always find some best fitting parameters: it doesn’t mean that the distribution is a good fit for the data.
x <- rbeta(2000, 0.5, 0.5)
MASS::fitdistr(x, "normal") mean sd
0.500284152 0.360421491
(0.008059270) (0.005698764)Show figure code
params <- MASS::fitdistr(x, "normal")
data.frame(x=x) |>
ggplot(aes(x=x)) +
geom_histogram(aes(y=after_stat(density)), bins=100, fill="grey", color="black") +
geom_line(
data = data.frame(
x=seq(-0.5, 1.5, length.out = 500),
y=dnorm(seq(-0.5, 1.5, length.out = 500), mean=params$estimate["mean"], sd=params$estimate["sd"])),
aes(x=x, y=y), color="red", linewidth=2) +
theme_bw() +
labs(x="x", y="Density")
Quizz
- Run the following:
df <- data.frame(y = rgamma(50000, 0.5, 8))Code
ggplot(df, aes(x=y)) +
geom_histogram(aes(y=after_stat(density)), bins=100, fill="grey", color="black") +
theme_bw() +
labs(x="y", y="Density")
- What are the parameters of the best fitting Normal distribution of
y?
MASS::fitdistr(df$y, "normal") mean sd
0.0626589570 0.0886953261
(0.0003966576) (0.0002804792)- Fit an intercept-only linear regression model to
yand look at its summary.
model <- lm(y ~ 1, data=df)
summary(model)
Call:
lm(formula = y ~ 1, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.06266 -0.05617 -0.03399 0.02038 1.30855
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0626590 0.0003967 158 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0887 on 49999 degrees of freedom- What do you notice?
Linear Model = “Gaussian” Model
- You might have heard that an “intercept only” model predicts the mean of the data
- That is true, but it is not the whole story…
- More generally, a linear model finds the parameters of the best fitting normal distribution
- The location \(\mu\) corresponds to the mean of the data
- The spread \(\sigma\) corresponds to the standard deviation of the “residuals” (i.e., the distance between the observed data and the predicted mean)
- In other words, a typical “linear model” is actually a “Gaussian” model (based on a Normal distribution)
Code
df <- data.frame(y = rnorm(5000, 10, 5))
ggplot(df, aes(y=y)) +
geom_jitter(aes(x=0), alpha=0.3) +
geom_linerange(data=data.frame(x=0, y=mean(df$y), ymin=mean(df$y)-0.5*sd(df$y), ymax=mean(df$y)+0.5*sd(df$y)),
aes(x=x, y=y, ymin=ymin, ymax=ymax), color="red", linewidth=2) +
geom_point(data=data.frame(x=0, y=mean(df$y)), aes(x=x, y=y), color="red", size=4) +
geom_hline(yintercept = 10, linetype="dashed", color="red", linewidth=1) +
coord_cartesian(xlim = c(-1, 0.4)) +
labs(x="", y="y") +
ggside::geom_ysidedensity(data=data.frame(y = rnorm(50000, 10, 5)), color = "red", linewidth = 2) +
theme_minimal() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank()) 
Homoscedasticity
- In typical linear models, we express \(\mu\) as a function of the predictors (e.g., \(\mu = \beta_{0} + \beta_{1}x\))
- \(\sigma\) is assumed to be constant across all observations (i.e., homoscedasticity)
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_ml <- 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_ml
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
From RSS to Sigma
- Under OLS, we can compute the regression line using a simple formula (which is the same as under MLE with a Gaussian distribution)
- Once we have the line, we can compute the residuals, from which we can work out the \(\sigma\) (sigma) of the best fitting Normal distribution
- 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}}\)
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
df <- bayestestR::simulate_correlation(n=500, r=0.7)
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?
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"))
