Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
library(tidyverse)
regress <- function(x_bar, y_bar, s_x, s_y, rho, n, alpha = 0.05) {
beta_hat <- rho * s_y / s_x
alpha_hat <- y_bar - beta_hat * x_bar
ssr <- (1 - rho^2) * s_y^2 * (n - 1)
sigma_sq_hat <- ssr / (n - 2)
sxx <- s_x^2 * (n - 1)
se_beta_hat <- sqrt(sigma_sq_hat / sxx)
se_alpha_hat <- sqrt(sigma_sq_hat * (1 / n + x_bar^2 / sxx))
ci_quantiles <- c(alpha/2, 1 - alpha/2)
alpha_ci <- alpha_hat + qt(ci_quantiles, df = n - 2) * se_alpha_hat
beta_ci <- beta_hat + qt(ci_quantiles, df = n - 2) * se_beta_hat
table <- tibble(
term = c("intercept", "slope"),
estimate = c(alpha_hat, beta_hat),
std.error = c(se_alpha_hat, se_beta_hat),
conf.low = c(alpha_ci[1], beta_ci[1]),
conf.high = c(alpha_ci[2], beta_ci[2])
)
slr <- list(
table = table,
alpha_hat = alpha_hat,
beta_hat = beta_hat
)
class(slr) <- "slr"
slr
}
print.slr <- function(x, ...) {
print(x$table)
}
predict.slr <- function(object, x, ...) {
object$alpha_hat + object$beta_hat * x
}
model <- regress(
x_bar = 101.70,
y_bar = 123.76,
s_x = 11.63,
s_y = 14.19,
rho = 0.28,
n = 927
)
model
data <- tibble(
iq = seq(0, 150)
) %>%
mutate(
.fitted = predict(model, iq),
overconfidence = .fitted - iq
)
ggplot(data) +
aes(iq, .fitted) +
geom_line(color = "steelblue") +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
theme_minimal() +
expand_limits(x = 0, y = 0) +
labs(
x = "IQ",
y = "SAIQ",
title = "High IQ individuals overestimate their IQ by less than low IQ individuals",
caption = "Dashed line is identity, blue line is SAIQ ~ IQ simple regression "
)
data %>%
ggplot(aes(iq, overconfidence)) +
geom_line(color = "steelblue") +
theme_minimal() +
labs(
x = "IQ",
y = "Overconfidence on self-assessed IQ",
title = "Overconfidence decreases on average with IQ"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment