PISA 2022 資料分析

required_packages <- c("haven", "ggplot2", "dplyr")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace,
                                               logical(1), quietly = TRUE)]
if (length(missing_packages) > 0) {
  install.packages(missing_packages, repos = "https://cloud.r-project.org")
}

library(haven)
library(ggplot2)
library(dplyr)

data_file <- "PISA_tawian2022_trimmed_lab.sav"

資料讀取與清理

本分析使用以下變項：

• WORKPAY：學生在上學前或放學後從事有薪工作的頻率
• Gender：學生性別
• PV1MATH：數學成就分數（第一個 plausible value）

WORKPAY 的有效工作頻率是 0–10。原始資料中的 95 是 Valid Skip， 97 是 Not Applicable，98 是 Invalid，99 是 No Response；這些是問卷代碼， 不是工作頻率，因此不納入工作頻率分析。

pisa <- read_sav(data_file) %>%
  transmute(
    WORKPAY = as.numeric(WORKPAY),
    Gender = as.numeric(Gender),
    PV1MATH = as.numeric(PV1MATH)
  ) %>%
  mutate(
    WORKPAY = if_else(WORKPAY %in% 0:10, WORKPAY, NA_real_),
    Gender = recode(Gender, `1` = "Female", `2` = "Male", .default = NA_character_),
    WORKPAY_label = factor(
      WORKPAY,
      levels = 0:10,
      labels = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10+")
    )
  )

cat("資料筆數：", nrow(pisa), "\n")

## 資料筆數： 5599

cat("有效 PV1MATH：", sum(!is.na(pisa$PV1MATH)), "\n")

## 有效 PV1MATH： 5599

WORKPAY 分布

使用比例長條圖呈現不同工作頻率的學生比例。

workpay_counts <- pisa %>%
  filter(!is.na(WORKPAY_label)) %>%
  count(WORKPAY_label) %>%
  mutate(percent = n / sum(n) * 100)

ggplot(workpay_counts, aes(WORKPAY_label, percent)) +
  geom_col(fill = "#D97757", width = 0.75) +
  geom_text(aes(label = sprintf("%.1f%%", percent)), vjust = -0.35, size = 3) +
  scale_y_continuous(
    labels = function(x) paste0(x, "%"),
    expand = expansion(mult = c(0, 0.12))
  ) +
  labs(
    title = "Working for pay before or after school",
    x = "Times per week",
    y = "Students (%)"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

Gender 圓餅圖

gender_counts <- pisa %>%
  filter(!is.na(Gender)) %>%
  count(Gender) %>%
  mutate(percent = n / sum(n) * 100)

ggplot(gender_counts, aes(x = "", y = percent, fill = Gender)) +
  geom_col(width = 1, color = "white", show.legend = FALSE) +
  coord_polar(theta = "y") +
  geom_text(
    aes(label = sprintf("%s\n%d (%.1f%%)", Gender, n, percent)),
    position = position_stack(vjust = 0.5),
    size = 4,
    color = "white"
  ) +
  scale_fill_manual(values = c(Female = "#3B82A0", Male = "#E09F3E")) +
  labs(title = "Student gender", x = NULL, y = NULL) +
  theme_void() +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.background = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    plot.margin = margin(10, 10, 10, 10)
  )

PV1MATH 分布

使用直方圖與密度曲線觀察數學成就分數的分布，虛線代表平均數。

ggplot(pisa, aes(PV1MATH)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 30,
    fill = "#3B82A0",
    color = "white"
  ) +
  geom_density(color = "#B7472A", linewidth = 1) +
  geom_vline(
    xintercept = mean(pisa$PV1MATH, na.rm = TRUE),
    linetype = "dashed",
    color = "#B7472A"
  ) +
  labs(
    title = "Distribution of mathematics achievement",
    subtitle = "PV1MATH with density curve and mean (dashed line)",
    x = "PV1MATH score",
    y = "Density"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))

摘要統計

pisa %>%
  summarise(
    n = sum(!is.na(PV1MATH)),
    mean_PV1MATH = mean(PV1MATH, na.rm = TRUE),
    median_PV1MATH = median(PV1MATH, na.rm = TRUE),
    sd_PV1MATH = sd(PV1MATH, na.rm = TRUE),
    min_PV1MATH = min(PV1MATH, na.rm = TRUE),
    max_PV1MATH = max(PV1MATH, na.rm = TRUE)
  )

## # A tibble: 1 × 6
##       n mean_PV1MATH median_PV1MATH sd_PV1MATH min_PV1MATH max_PV1MATH
##   <int>        <dbl>          <dbl>      <dbl>       <dbl>       <dbl>
## 1  5599         538.           542.       110.        219.        917.