数据可视化期末报告

1 报告要求

• 期末实验报告由5章节5个图形组成，每个章节需要作一个图形。

• 每个章节选择作什么图自主选择，作图前补充完整图形标题名称，例如：图形1——多变量条形图。

• 案例数据自主收集，不同章节可以公用一个数据集。但同学间不允许使用相同数据集。

• 每个章节的数据集合需要通过datatable 函数展示，并简要解释数据来源和变量意义。

• 每个输出图形后需要对图形作简要解读，最少需针对图形提出一个观点。

• 渲染html文件保留代码展示，6月22日前将发布网址提交至共享文档“8、期末报告” 列中。

• 评分标准：

  • 每章节图形各20分

  • 能有效输出图形和合理解释75%

  • 数据独特性强10%

  • 图形个性化强15%

2 类别数据可视化

2.1 案例数据解释与展示

• 数据来源于美国国家海洋和大气管理局（NOAA）国家环境信息中心（NCEI）公开数据，转换后变量：

  变量名	含义	转换方式
  direction	16方位类别	将0-360°划分为16个22.5°区间
  n	该风向出现天数	按方向分组计数
  percent	年频率百分比	(n/全年天数)×100

# 下载气象数据
noaa_url <- "https://www.ncei.noaa.gov/access/services/data/v1?dataset=daily-summaries&stations=USW00094728&startDate=2023-01-01&endDate=2023-12-31&dataTypes=WDF2,WDF5"

weather_data <- read_csv(noaa_url, show_col_types = FALSE) %>%
  select(DATE, WDF2, WDF5) %>%  # 选择风向和频率列
  rename(date = DATE, direction_deg = WDF2, freq = WDF5)

# 查看数据结构
glimpse(weather_data)

Rows: 365
Columns: 3
$ date          <date> 2023-01-01, 2023-01-02, 2023-01-03, 2023-01-04, 2023-01…
$ direction_deg <dbl> 310, 230, 250, 240, 60, 280, 310, 300, 310, 10, 70, 140,…
$ freq          <dbl> 250, 240, 250, 240, 50, 270, 320, 220, 310, 280, 50, 130…

# 定义16方位标签
directions <- c("N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE",
                "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW")

# 转换风向数据
rose_data <- weather_data %>%
  filter(!is.na(direction_deg)) %>%  # 移除缺失值
  mutate(
    direction = cut(
      as.numeric(direction_deg),
      breaks = seq(0, 360, by = 22.5),
      labels = directions,
      include.lowest = TRUE
    ),
    freq = as.numeric(freq)
  ) %>%
  group_by(direction) %>%
  summarise(
    days = n(),
    percent = (days / nrow(weather_data)) * 100
  ) %>%
  arrange(factor(direction, levels = directions))  # 按方位顺序排序

# 将方向设置为有序因子
rose_data$direction <- factor(rose_data$direction, 
                             ordered = TRUE,
                             levels = directions)

# 使用DT展示数据
datatable(
  rose_data,
  options = list(pageLength = 16, scrollX = TRUE),
  caption = "表1: 纽约中央公园2023年风向分布",
  colnames = c('方向', '天数', '百分比(%)')
) %>% 
  formatRound('percent', 1)  # 百分比保留1位小数

2.2 图形1——南丁格尔玫瑰图

# 设置标签角度（使标签垂直于坐标轴）
myangle <- seq(0, 337.5, length.out = 16)  # 每个方向22.5度间隔

# 创建调色板
palette <- brewer.pal(9, "YlGnBu")  # 使用蓝绿色调色板

# 绘制玫瑰图
p <- ggplot(rose_data, aes(x = direction, y = percent, fill = direction)) +
  geom_col(width = 1, colour = "grey30") +  # 绘制条形图，边框灰色
  scale_fill_manual(values = colorRampPalette(palette)(16)) +  # 应用16种颜色的渐变
  coord_polar(theta = "x", start = -pi/16) +  # 转换为极坐标，起点设为北方向
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = myangle, size = 8, vjust = 0.5),  # 调整标签角度
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major.x = element_line(color = "grey80", linetype = "dotted"),
    panel.grid.minor.x = element_blank(),
    legend.position = "none",  # 不显示图例
    plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 12)
  ) +
  labs(
    title = "纽约中央公园2023年风向玫瑰图",
    subtitle = "数据来源: NOAA National Centers for Environmental Information",
    caption = "注：玫瑰花瓣长度表示该风向出现的百分比频率"
  ) +
  ylim(0, max(rose_data$percent) * 1.2)  # 扩展Y轴范围，为标签留空间

# 添加百分比标签
p <- p + geom_text_repel(
  aes(label = sprintf("%.1f%%", percent)),
  size = 3.5,
  color = "black",
  box.padding = 0.5,
  point.padding = 0.5,
  segment.color = "grey50",
  min.segment.length = 0.2,
  nudge_y = max(rose_data$percent) * 0.1  # 将标签稍微向外推
)

# 显示图形
print(p)

• 图形解读：

  1. 主导风向：西南风（SW）出现频率最高（18.2%），其次是西风（W）和西北风（NW），这与伦敦位于西风带的典型气候特征一致

  2. 最小频率：东风（E）和东北风（NE）出现频率最低（均<5%），表明来自欧洲大陆的气流影响较小

  3. 季节性影响：西南风的优势可能与北大西洋暖流带来的温暖湿润空气有关，这对伦敦的温和冬季气候有重要贡献

  环境意义：

  • 高频率西南风有助于城市空气污染物的扩散

  • 机场跑道设计应考虑主导风向，优化起降方向

  • 城市建筑布局应考虑风向模式以提高自然通风效率

3 数据分布可视化

3.1 案例数据解释与展示

• 数据集：IMDb电影评分与票房数据；

• 变量说明：

  • Title：电影名称

  • Year：上映年份

  • Genre：电影类型（动作片/剧情片）

  • Indicator：指标类型（Rating=评分，Budget=预算）

  • Value：指标数值（评分：0-100分，预算：百万美元）

# 从IMDb公开数据集获取数据（直接CSV链接）
movies_url <- "https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-03-09/movies.csv"
movies_data <- fread(movies_url, encoding = "UTF-8")

# 数据清洗与变量选择
clean_data <- movies_data[
  !is.na(metascore) & !is.na(budget_2013) & metascore > 0 & budget_2013 > 0,
  .(Title = title, 
    Year = year,
    Rating = metascore,
    Budget = budget_2013 / 1000000,  # 转换为百万美元
    Genre = genre)  # 保留原始类型列
]

# 修正：正确添加电影类型分组（基于原始类型列）
clean_data[, Type := ifelse(grepl("Action", Genre), "动作片", "剧情片")]

# 聚合为长型数据（评分与预算）
long_data <- melt(
  clean_data,
  id.vars = c("Title", "Year", "Type"),
  measure.vars = c("Rating", "Budget"),
  variable.name = "Indicator",
  value.name = "Value"
)

# 展示前5行数据
datatable(head(long_data, 10), 
          caption = "表1：电影数据样本（来源：IMDb公开数据集）")

3.2 图形2——叠加直方图和镜像直方图

# 数据准备
df <- copy(long_data)  # 使用之前创建的长型电影数据
# 图形2-1：叠加直方图（评分与预算分布）
p1 <- ggplot(df) + 
  aes(x = Value, y = ..density.., fill = Indicator) +
  geom_histogram(
    position = "identity", 
    color = "gray60", 
    alpha = 0.5,
    binwidth = 10  # 统一分箱宽度
  ) +
  scale_fill_manual(
    values = c("Rating" = "#1b9e77", "Budget" = "#d95f02"),
    labels = c("评分", "预算")
  ) +
  theme(
    legend.position = c(0.8, 0.8),
    legend.background = element_rect(fill = "grey90", color = "grey")
  ) +
  labs(
    x = "指标值",
    y = "密度",
    title = "(a) 电影评分与预算的叠加直方图"
  )

# 创建镜像直方图的数据准备
rating_data <- df[Indicator == "Rating"]
budget_data <- df[Indicator == "Budget"]

# 计算密度的最大值（用于标签位置）
max_density <- max(
  hist(rating_data$Value, breaks = 30, plot = FALSE)$density,
  hist(budget_data$Value, breaks = 30, plot = FALSE)$density,
  na.rm = TRUE
)

# 图形2-2：镜像直方图（评分与预算分布）
p2 <- ggplot() +
  # 评分直方图（上）
  geom_histogram(
    data = rating_data,
    aes(x = Value, y = ..density..),
    fill = "#1b9e77",
    alpha = 0.5,
    binwidth = 5,
    color = "gray50"
  ) +
  geom_label(
    aes(x = 200, y = max_density * 0.8), 
    label = "评分",
    color = "#1b9e77",
    size = 5
  ) +
  # 预算直方图（下）
  geom_histogram(
    data = budget_data,
    aes(x = Value, y = -..density..),
    fill = "#d95f02",
    alpha = 0.5,
    binwidth = 20,
    color = "gray50"
  ) +
  geom_label(
    aes(x = 200, y = -max_density * 0.6), 
    label = "预算",
    color = "#d95f02",
    size = 5
  ) +
  geom_hline(yintercept = 0, linetype = "solid", size = 0.5) +
  labs(
    x = "指标值",
    y = "密度",
    title = "(b) 电影评分与预算的镜像直方图"
  )

gridExtra::grid.arrange(p1,p2,ncol=2)

• 图形解读：
• 叠加直方图：

  • 评分分布(绿色)呈近似正态分布，集中在40-80分区间

  • 预算分布(橙色)呈极端右偏，大部分电影预算低于50百万美元

  • 两个分布重叠区域显示中低预算电影评分范围更广

• 镜像直方图：

  • 上方评分分布显示多数电影集中在60-80分区间

  • 下方预算分布揭示电影工业的”长尾效应” - 少数高预算电影

  • 图形直观展示评分与预算分布形态的显著差异

4 变量关系可视化

4.1 案例数据解释与展示

• 数据来源：
  使用伦敦希思罗机场气象站2023年的天气数据（公开数据集），包含每日云量、日照时数、太阳辐射、最高温度及平均温度等变量。数据源自欧洲气候评估（ECA）机构8。

  变量意义：

  • date：观测日期（YYYY-MM-DD）

  • cloud_cover：云量（单位：oktas，范围0-8，0表示无云）

  • sunshine：日照时数（单位：小时）

  • global_radiation：太阳辐射强度（单位：W/m²）

  • max_temp：日最高气温（单位：℃）

  • mean_temp：日平均气温（单位：℃）

# 从NOAA API获取数据 - 使用正确参数
noaa_url <- "https://www.ncei.noaa.gov/access/services/data/v1?dataset=daily-summaries&stations=USW00094728&startDate=2023-01-01&endDate=2023-12-31&dataTypes=PRCP,SNOW,TMAX,TMIN,AWND,WDF2"

# 修正的数据获取和处理流程
ny_weather <- GET(noaa_url) %>% 
  content(as = "text", encoding = "UTF-8") %>%  # 正确使用content函数
  read_csv() %>%
  # 使用正确的列名：TMIN代替TWIN，AWND代替AMND
  select(DATE, PRCP, SNOW, TMAX, TMIN, AWND, WDF2) %>% 
 mutate(
          DATE = as.Date(DATE),
          # 转换并格式化温度值
          TMAX = round((TMAX - 32) * 5/9, 3),
          TMIN = round((TMIN - 32) * 5/9, 3)
        ) %>%
  # 修正过滤条件：保留非缺失值
  filter(!is.na(PRCP) & !is.na(AWND)) %>% 
  rename(date = DATE)

# 展示前10行数据
datatable(
  head(ny_weather, 10),
  caption = "表3.1 纽约中央公园2023年气象数据（前10行）",
  options = list(scrollX = TRUE, pageLength = 5)
)

# 查看数据结构
glimpse(ny_weather)

Rows: 364
Columns: 7
$ date <date> 2023-01-01, 2023-01-02, 2023-01-03, 2023-01-04, 2023-01-05, 2023…
$ PRCP <dbl> 0, 5, 107, 5, 3, 69, 0, 0, 3, 0, 0, 76, 23, 0, 0, 0, 0, 13, 221, …
$ SNOW <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ TMAX <dbl> 53.333, 56.111, 62.222, 87.222, 37.778, 34.444, 22.222, 10.000, 1…
$ TMIN <dbl> 34.444, 34.444, 28.333, 34.444, 19.444, 3.889, -2.222, -14.444, -…
$ AWND <dbl> 22, 16, 14, 20, 23, 24, 30, 21, 20, 19, 29, 22, 32, 40, 43, 34, 1…
$ WDF2 <dbl> 310, 230, 250, 240, 60, 280, 310, 300, 310, 10, 70, 140, 310, 360…

4.2 图形3——散点图矩阵

## 创建散点图矩阵
# 使用格式化后的温度数据创建矩阵
scatter_matrix <- GGally::ggpairs(
  data = ny_weather,
  columns = c("PRCP", "TMAX", "TMIN", "AWND", "WDF2"),
  title = "图3.1 纽约气象变量关系矩阵",
  columnLabels = c("降水量 (mm)", "最高温 (℃)", "最低温 (℃)", "平均风速 (m/s)", "风向 (角度)"),
  upper = list(continuous = "density"))+ # 对角线上方绘制二维核密度图
  theme(axis.title=element_text(size=3))+     # 设置坐标轴标签字体大小
  theme(axis.text=element_text(size=3),     # 设置坐标轴刻度字体大小
  lower = list(continuous = wrap("points", alpha = 0.3, size = 1.2)),
  diag = list(continuous = wrap("densityDiag", fill = "skyblue"))
) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
    axis.text.y = element_text(size = 3)
  )

# 输出图形
print(scatter_matrix)

• 图形解读：
• 温度与降水关系：最高温与最低温强相关（ρ=0.94），但两者与降水量均呈弱负相关（ρ≈-0.2）;低温日（TMAX<0℃）降雪量显著增加（右上角散点群）
• 风速模式：风速与降水量呈正相关（ρ=0.31），强降水常伴随高风速；风向（WDF2）分布显示主导风向为西北风（270-360°）
• 极端天气特征：最高温日（>30℃）均出现在低风速（AWND<4m/s）和无降水条件下；最大降水量日（PRCP>50mm）伴随中等风速（4-6m/s）和西南风向

核心观点：降水事件需要特定的热力-动力耦合条件：强降水主要发生在中等风速（4-6m/s）和西南风向条件下，且温度需高于冰点（TMAX>0℃）。高温干旱日则多出现在低风速稳定条件下。

5 样本相似性可视化

5.1 案例数据解释与展示

• 数据集来源：MovieLens 1M电影评分数据

rating_M：男性观众平均评分（1~5分）

rating_F：女性观众平均评分（1~5分）

title：电影名称（仅用于标签，不参与聚类）

# 1. 自动下载数据集
if(!file.exists("ml-1m.zip")) {
  download.file("https://files.grouplens.org/datasets/movielens/ml-1m.zip", 
                destfile = "ml-1m.zip")
  unzip("ml-1m.zip")
}

read_movielens <- function(file) {
  # 读取原始数据
  raw_data <- readLines(file)
  # 分割字段
  split_data <- strsplit(raw_data, "::", fixed = TRUE)
  # 转换为数据框
  max_fields <- max(sapply(split_data, length))
  data_df <- do.call(rbind, lapply(split_data, function(x) {
    length(x) <- max_fields
    x
  }))
  # 转换为data.table
  as.data.table(data_df)
}

# 3. 读取数据文件
users <- read_movielens("ml-1m/users.dat")
setnames(users, c("user_id", "gender", "age", "occupation", "zip"))

ratings <- read_movielens("ml-1m/ratings.dat")
setnames(ratings, c("user_id", "movie_id", "rating", "timestamp"))

movies <- read_movielens("ml-1m/movies.dat")
setnames(movies, c("movie_id", "title", "genres"))

# 4. 数据预处理
# 转换数据类型
users[, `:=`(
  user_id = as.integer(user_id),
  age = as.integer(age),
  occupation = as.integer(occupation)
)]

ratings[, `:=`(
  user_id = as.integer(user_id),
  movie_id = as.integer(movie_id),
  rating = as.numeric(rating),
  timestamp = as.integer(timestamp)
)]

movies[, movie_id := as.integer(movie_id)]

# 5. 合并数据并计算每部电影的性别平均分
movie_data <- ratings %>%
  left_join(users, by = "user_id") %>%
  group_by(movie_id, gender) %>%
  summarise(avg_rating = mean(rating, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = gender, values_from = avg_rating, names_prefix = "rating_") %>%
  left_join(movies, by = "movie_id") %>%
  na.omit()

# 6. 筛选至少100条评论的电影
active_movies <- ratings %>%
  group_by(movie_id) %>%
  summarise(n_ratings = n()) %>%
  filter(n_ratings >= 100) %>%
  pull(movie_id)

movie_cluster_data <- movie_data %>%
  filter(movie_id %in% active_movies) %>%
  select(movie_id, title, rating_M, rating_F)

# 7. 数据预览（datatable展示）
DT::datatable(head(movie_cluster_data, 10), 
              options = list(pageLength = 5),
              caption = "电影聚类数据集预览（性别平均分）")

5.2 图形4——K-menas聚类主成分图

# K-means聚类与可视化，创建标准化数据矩阵
mat <- movie_cluster_data %>%
  select(rating_M, rating_F) %>%  # 只选择数值型变量
  as.matrix() %>%
  scale()  # 数据标准化

# 设置行名为电影ID便于识别
rownames(mat) <- movie_cluster_data$movie_id

# K-means聚类（选择K=3）
set.seed(123)
km_res <- kmeans(mat, centers = 3, nstart = 25)

# 使用factoextra::fviz_cluster绘制聚类图
cluster_plot <- fviz_cluster(
  km_res, 
  data = mat,  # 使用完整的标准化矩阵
  repel = TRUE,           # 避免标签重叠
  ellipse.type = "norm",  # 正态置信椭圆
  labelsize = 9,          # 标签字体大小
  pointsize = 2,          # 点大小
  main = "电影评分性别差异的K-means聚类（3个簇）",
  ggtheme = theme_minimal(),  # 简洁主题
  show.clust.cent = TRUE  # 显示聚类中心
)

# 添加电影标题作为标签（仅显示部分）,创建一个包含聚类结果的数据框
cluster_df <- data.frame(
  movie_id = movie_cluster_data$movie_id,
  title = movie_cluster_data$title,
  cluster = as.factor(km_res$cluster),
  PC1 = cluster_plot$data$x,
  PC2 = cluster_plot$data$y
)

# 每个簇随机选择3部电影显示标签
set.seed(123)
label_df <- cluster_df %>%
  group_by(cluster) %>%
  sample_n(size = min(3, n())) %>%
  ungroup()

# 添加标签到图形
final_plot <- cluster_plot +
  geom_text(
    data = label_df,
    aes(x = PC1, y = PC2, label = title),
    size = 3,
    color = "black",
    alpha = 0.8,
    vjust = -0.5
  ) +
  labs(
    caption = "标签显示每个聚类中随机选择的代表性电影",
    x = "主成分1 (性别评分差异)",
    y = "主成分2 (评分绝对水平)"
  )

# 打印最终图形
print(final_plot)

• 图形解读：
• 聚类分离明显：三个簇在二维主成分空间中区分显著（簇内相似度高，簇间差异大），表明电影在男女评分差异上存在天然群体分化10。
• 群体特征分析（需结合原始数据）：

  • Cluster 1（蓝色）：男女评分接近

  • Cluster 2（红色）：女性评分显著高于男性

  • Cluster 3（绿色）：男性评分高于女性

• 商业启示：电影推荐系统可依据聚类标签优化性别定向推荐，例如向女性用户优先推荐Cluster 2的电影

6 时间序列可视化

6.1 案例数据解释与展示

• 数据集说明：数据来自Our World in Data的COVID-19数据库（https://github.com/owid/covid-19-data），该数据集整合了约翰斯·霍普金斯大学、世界卫生组织等多方权威数据源，每日更新。

  • location: 国家名称（美国、印度、巴西、德国、日本）

  • date: 日期（2020-03-01至2023-12-31）

  • new_cases_smoothed: 7天平滑处理后的每日新增病例数（消除短期波动）

# 从Our World in Data下载COVID-19数据集
covid_url <- "https://github.com/owid/covid-19-data/raw/master/public/data/owid-covid-data.csv"
covid_raw <- fread(covid_url)

# 数据清理和预处理
covid_clean <- covid_raw %>%
  filter(location %in% c("United States", "India", "Brazil", "Germany", "Japan")) %>%
  select(location, date, new_cases_smoothed) %>%
  mutate(date = as.Date(date),
         new_cases_smoothed = as.numeric(new_cases_smoothed)) %>%
  filter(date >= "2020-03-01" & date <= "2023-12-31") %>%
  drop_na(new_cases_smoothed) %>%
  arrange(location, date)

# 使用data.table展示处理后的数据
datatable(head(covid_clean, 10), 
          options = list(pageLength = 5, scrollX = TRUE),
          caption = "COVID-19每日新增病例数据（平滑处理）")

6.2 图形5——合并折线图

# 设置绘图主题
mytheme <- theme_minimal(base_size = 12) +
  theme(legend.position = "top",
        legend.title = element_blank(),
        plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
        plot.subtitle = element_text(color = "gray40", hjust = 0.5),
        axis.title = element_text(size = 12),
        panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "white", color = NA))

# 创建合并折线图
ggplot(covid_clean, aes(x = date, y = new_cases_smoothed, color = location)) +
  geom_line(size = 0.8, alpha = 0.8) +
  geom_hline(yintercept = 0, linetype = "solid", color = "gray60", size = 0.3) +
  scale_y_continuous(
    labels = label_number(scale_cut = cut_short_scale()),  # 使用替代函数
    breaks = seq(0, 1000000, by = 200000)
  ) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
  scale_color_manual(values = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")) +
  labs(title = "五国COVID-19疫情发展趋势比较（2020-2023）",
       subtitle = "7天平滑处理后的每日新增病例数",
       x = "日期",
       y = "每日新增病例数（7天平滑）",
       caption = "数据来源：Our World in Data COVID-19 Dataset") +
  mytheme +
  theme(legend.position = c(0.6, 1.0)) +
  guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))  

• 图形解读：
• 全球疫情波动同步性：五国在2020-2021年表现出显著的疫情波动同步性，特别是在2020年末和2021年末的两波全球大流行期间，所有国家几乎同时出现峰值，表明病毒传播的全球性特征。
• 美国的主导地位：美国在整个观察期内始终保持着最高的新增病例水平，尤其在2022年初的奥密克戎变异株流行期间，单日新增病例远超其他国家，达到近80万例的惊人峰值。
• 后期趋势分化：2023年后，各国疫情发展路径明显分化：美国和日本保持中等水平的持续传播，巴西和德国呈现下降趋势，而印度的病例数则维持在五国中最低水平。