R语言基础：二维数据可视化与方法入门

# 连续型结果：明显非线性
n <- 300
x <- seq(0, 10, length.out = n)
y <- 2 + 1.8 * x - 0.18 * x^2 + rnorm(n, 0, 1.5)
df_nl <- tibble(x = x, y = y)

# 二元结果：适合逻辑回归
x2 <- rnorm(400)
p_true <- 1 / (1 + exp(-(-0.6 + 1.8 * x2)))
y_bin <- rbinom(400, 1, p_true)
df_logit <- tibble(x = x2, y = y_bin, p_true = p_true)

model_linear <- lm(y ~ x, data = df_nl)
model_quad   <- lm(y ~ x + I(x^2), data = df_nl)

pred_grid <- tibble(
  x = seq(min(df_nl$x), max(df_nl$x), length.out = 200)
)
pred_grid$linear <- predict(model_linear, newdata = pred_grid)
pred_grid$quadratic <- predict(model_quad, newdata = pred_grid)

p1 <- ggplot(df_nl, aes(x, y)) +
  geom_point(alpha = 0.5) +
  geom_line(data = pred_grid, aes(y = linear), linewidth = 1.2) +
  labs(
    title = "Linear Regression",
    subtitle = "A linear model can only fit a straight line",
    x = "X",
    y = "Y"
  )

p2 <- ggplot(df_nl, aes(x, y)) +
  geom_point(alpha = 0.5) +
  geom_line(data = pred_grid, aes(y = quadratic), linewidth = 1.2) +
  labs(
    title = "Nonlinear Regression with a Quadratic Term",
    subtitle = "The quadratic term allows the marginal effect to change with X",
    x = "X",
    y = "Y"
  )

p1 + p2

summary(model_quad)

logit_model <- glm(y ~ x, data = df_logit, family = binomial)

grid_logit <- tibble(
  x = seq(min(df_logit$x), max(df_logit$x), length.out = 300)
)
grid_logit$prob_hat <- predict(logit_model, newdata = grid_logit, type = "response")

ggplot(df_logit, aes(x, y)) +
  geom_jitter(height = 0.05, width = 0, alpha = 0.25) +
  geom_line(data = grid_logit, aes(y = prob_hat), linewidth = 1.5) +
  labs(
    title = "Logistic Regression",
    subtitle = "The fitted curve gives predicted probabilities between 0 and 1",
    x = "X",
    y = "Observed Y and Predicted Probability"
  )

summary(logit_model)

glm(y ~ x1 + x2, data = your_data, family = binomial)

set.seed(1)
df_tree <- tibble(
  x1 = rnorm(300),
  x2 = rnorm(300)
) |>
  mutate(y = ifelse(x1 + 0.8 * x2 + rnorm(300, 0, 0.6) > 0, 1, 0))

tree_model <- rpart(as.factor(y) ~ x1 + x2, data = df_tree, method = "class")
rpart.plot(tree_model, type = 2, extra = 104, under = TRUE, fallen.leaves = TRUE)
title("Decision Tree for Classification")

x_over <- seq(0, 10, length.out = 60)
y_over <- sin(x_over) + rnorm(60, 0, 0.25)
df_over <- tibble(x = x_over, y = y_over)

fit_simple <- lm(y ~ x, data = df_over)
fit_flex   <- lm(y ~ poly(x, 10), data = df_over)

grid_over <- tibble(x = seq(0, 10, length.out = 300))
grid_over$simple <- predict(fit_simple, newdata = grid_over)
grid_over$flexible <- predict(fit_flex, newdata = grid_over)

p_over1 <- ggplot(df_over, aes(x, y)) +
  geom_point() +
  geom_line(data = grid_over, aes(y = simple), linewidth = 1.2) +
  labs(
    title = "A Simpler Model",
    subtitle = "It may underfit, but it is often more stable",
    x = "X",
    y = "Y"
  )

p_over2 <- ggplot(df_over, aes(x, y)) +
  geom_point() +
  geom_line(data = grid_over, aes(y = flexible), linewidth = 1.2) +
  labs(
    title = "An Overly Flexible Model",
    subtitle = "It may learn noise as if it were signal: overfitting",
    x = "X",
    y = "Y"
  )

p_over1 + p_over2

for (i in 1:5) {
  print(paste("Current i =", i))
}

set.seed(123)
df_loop <- tibble(
  y  = rnorm(200),
  x1 = rnorm(200),
  x2 = rnorm(200),
  x3 = rnorm(200)
)

vars <- c("x1", "x2", "x3")
results <- c()

for (v in vars) {
  formula_now <- as.formula(paste("y ~", v))
  model_now <- lm(formula_now, data = df_loop)
  results[v] <- coef(model_now)[2]
}

results

square_num <- function(x) {
  return(x^2)
}

square_num(5)

get_slope <- function(data, varname) {
  formula_now <- as.formula(paste("y ~", varname))
  model_now <- lm(formula_now, data = data)
  return(coef(model_now)[2])
}

get_slope(df_loop, "x1")

slopes <- c()

for (v in vars) {
  slopes[v] <- get_slope(df_loop, v)
}

slopes

sapply(vars, function(v) get_slope(df_loop, v))

df_did <- expand_grid(
  treat = c(0, 1),
  post = c(0, 1),
  id = 1:100
) |>
  mutate(
    y = 50 +
      3 * treat +
      2 * post +
      6 * (treat * post) +
      rnorm(n(), 0, 3)
  )

summary_did <- df_did |>
  group_by(treat, post) |>
  summarise(mean_y = mean(y), .groups = "drop") |>
  mutate(
    group = ifelse(treat == 1, "Treatment", "Control"),
    time  = ifelse(post == 1, "After", "Before")
  )

ggplot(summary_did, aes(x = time, y = mean_y, group = group, color = group)) +
  geom_line(linewidth = 1.5) +
  geom_point(size = 3) +
  labs(
    title = "Difference-in-Differences",
    subtitle = "The extra increase in the treatment group after the policy is the DID effect",
    x = NULL,
    y = "Average Outcome",
    color = NULL
  ) +
  theme(legend.position = "top")

did_model <- lm(y ~ treat * post, data = df_did)
summary(did_model)

set.seed(100)
score <- seq(40, 80, length.out = 250)
cutoff <- 60
treat_rdd <- ifelse(score >= cutoff, 1, 0)
outcome <- 20 + 0.6 * score + 8 * treat_rdd + rnorm(250, 0, 2)

df_rdd <- tibble(score = score, treat = treat_rdd, outcome = outcome)

ggplot(df_rdd, aes(score, outcome)) +
  geom_point(alpha = 0.5) +
  geom_vline(xintercept = cutoff, linetype = 2, linewidth = 1) +
  geom_smooth(
    data = df_rdd |> filter(score < cutoff),
    method = "lm", se = FALSE, linewidth = 1.3
  ) +
  geom_smooth(
    data = df_rdd |> filter(score >= cutoff),
    method = "lm", se = FALSE, linewidth = 1.3
  ) +
  labs(
    title = "Regression Discontinuity Design",
    subtitle = "A visible jump at the cutoff suggests a treatment effect",
    x = "Running Variable",
    y = "Outcome"
  )

rdd_model <- lm(outcome ~ score + treat + score:treat, data = df_rdd)
summary(rdd_model)