library(nycflights13)
library(autoReg)
library(naniar)
library(patchwork)
library(tidyverse)
library(skimr)
library(showtext)
font_add("simhei", "C:/Users/21565/AppData/Local/Microsoft/Windows/Fonts/simhei.ttf") 
showtext_auto()

flights <- nycflights13::flights #主要表
glimpse(flights)

# 筛选出含NA的变量 
NA_variables <- flights %>% 
  summarise(across(everything(), ~ any(is.na(.))))  %>%  # 检查每个变量是否含NA 
  pivot_longer(everything(), names_to = "variable", values_to = "has_na") %>%  # 转换为长格式 
  filter(has_na) %>%  
  select(variable)   
NA_variables <- flights %>%
  select(all_of(NA_variables$variable))  # 选择含 NA 的变量
summary(NA_variables)

# arr_delay ~ dep_delay
set.seed(1042)
flights_sample <- slice_sample(flights, n = 5000) # 随机抽取n行 
# arr_delay ~ dep_delay
ggplot(flights_sample, aes(x = dep_delay, y = arr_delay)) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dotdash", linewidth = 0.8) +  # 添加 y = x的参考线
  geom_point(alpha = 0.2, color = "steelblue") +  # 设置点的透明度和颜色  
  labs(
    title = "到达延误 (arr_delay) ≈ 起飞延误 (dep_delay)", y = "到达延误 (arr_delay)",x = "起飞延误 (dep_delay)",
    subtitle = "y = x"
  ) +
  theme_minimal() +
  theme(
    plot.title  = element_text(face = "bold", size = 30),
    axis.title  = element_text(face = "bold", size = 20),
    axis.text.x  = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y  = element_text(size = 20, face = "bold.italic"),
    plot.subtitle = element_text(hjust = 0.8, vjust = -20, size = 20, color = "red")  # 坐标轴标题加粗 
  )

# 计算偏离程度 
flights <- flights %>%
  mutate(arr_delay_dep_delay_deviation = arr_delay - dep_delay)

# 统计偏离程度 
mode <- function(x) {  #构造函数计算众数
  as.numeric(names(which.max(table(x))) )
}
deviation_stats <- flights %>%
  summarise(
    mean_deviation = mean(arr_delay_dep_delay_deviation, na.rm  = TRUE),
    sd_deviation = sd(arr_delay_dep_delay_deviation, na.rm  = TRUE),
    median_deviation = median(arr_delay_dep_delay_deviation, na.rm  = TRUE),
    mode_deviation = mode(arr_delay_dep_delay_deviation)
  )
deviation_stats 

# 计算平均速度 
flights <- flights %>%
  mutate(speed = (distance / air_time) * 60)  # 速度单位：mph 
# speed ~ distance---------------------------------------------------------------------
flights_sample <- slice_sample(flights, n = 5000) # 随机抽取n行 
# speed ~ distance
ggplot(flights_sample, aes(x = distance, y = speed)) +
  geom_point(alpha = 0.3, color = "steelblue") +  # 设置点的透明度和颜色 
  stat_smooth(method = "lm", formula = y ~ log(x), color = "red3", fill = "red") +  # 对数回归曲线 
  labs(
    title = " 平均速度 (mph) ~ 飞行距离 (miles)", 
    y = "平均速度 (mph)",
    x = "飞行距离 (miles)"
  ) +
  theme_minimal() +
  theme(
    plot.title  = element_text(face = "bold", size = 30),
    axis.title  = element_text(face = "bold", size = 20),
    axis.text.x  = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y  = element_text(size = 20, face = "bold.italic"),
  )

planes <- nycflights13::planes %>% # 重命名 planes 表中的 year 变量为 manufacture_year 
  rename(manufacture_year = year)
plane_sample <- slice_sample(planes, n = 5000) # 随机抽取n行 
vis_miss(plane_sample) # 可视化检查缺失值

# model_mean_speed ~ model---------------------------------------------------------------------
# 按 model 分组计算平均速度 
model_speed <- flights_planes %>%
  group_by(model) %>%
  summarise(model_mean_speed = mean(speed, na.rm = TRUE))
# model_mean_speed ~ model
ggplot(model_speed, aes(x = reorder(model, model_mean_speed), y = model_mean_speed, fill = model)) +
  geom_bar(width = 0.6, stat = "identity") +
  coord_flip() +  # 水平柱状图 
  labs(
    title = "各机型平均速度分布图",
    y = "平均速度 (mph)",
    x = "机型 (model)"
  ) +
  theme_minimal() +
  theme(
    plot.title  = element_text(face = "bold", size = 30),
    axis.title  = element_text(face = "bold", size = 20),
    axis.text.x  = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y  = element_text(size = 8, face = "bold.italic"),
    legend.position  = "none",
    panel.grid.major.y = element_blank(),  # 去除主要水平网格线
    panel.grid.minor.y = element_blank()   # 去除次要水平网格线
  )

# 合并 flights 与 weather （通过 origin, time_hour）
weather <- nycflights13::weather
weather_sample <- slice_sample(weather, n = 5000) # 随机抽取n行 
vis_miss(weather_sample) # 可视化检查缺失值

# 遍历 weather 表中的所有变量
weather_vars <- names(weather) %>% 
  setdiff(c("origin", "time_hour", "year", "month", "day", "hour")) # 排除不需要的列 
# 创建一个空列表，用于存储所有热力图 
speed_weather_plot_list <- list() #创建列表
# 生成所有图形并存储到列表中 
for (var in weather_vars) {
  p <- ggplot(flights_weather_sample, aes(x = .data[[var]], y = speed)) +
    stat_binhex(bins = 30) +  
    scale_fill_gradient(low = "lightblue", high = "red") +
    labs(
      title = paste("平均速度 ~", var),
      y = "平均速度 (mph)",
      x = paste(var, "(单位根据变量而定)")
    ) +
    theme_minimal()
  
  speed_weather_plot_list[[var]] <- p +
    theme(plot.background  = element_rect(fill = "#CCCCCC", color = "white", size = 1),
          legend.key.width  = unit(0.1, "cm"))
}
# 使用 patchwork 将所有图形组合并显示 
wrap_plots(speed_weather_plot_list, ncol = 3)  # 每行显示n个图形 

# 遍历所有天气变量 
for (var in weather_vars) {
  # --- 动态分箱策略 --- 
  unique_values <- unique(flights_weather[[var]])
  if (length(unique_values) <= 5) {
    breaks <- unique_values  # 离散变量直接分箱（如雾：0/1）
  } else {
    breaks <- unique(quantile(flights_weather[[var]], probs = seq(0, 1, 0.2), na.rm  = TRUE))  # 连续变量分位数分箱 
  }
  
  # --- 分层抽样 --- 
  flights_weather_sample <- flights_weather %>%
    mutate(var_bin = cut(.data[[var]], breaks, include.lowest  = TRUE)) %>%
    group_by(var_bin) %>%
    sample_n(100, replace = TRUE) %>%  # 允许重复抽样（小样本时必需）
    ungroup()
  
  # --- 生成热力图（速度分析） --- 
  p <- ggplot(flights_weather_sample, aes(x = .data[[var]], y = speed)) +
    stat_binhex(bins = 30) +  # 六边形热力图 
    scale_fill_gradient(low = "lightblue", high = "red") +  # 颜色渐变 
    geom_smooth(method = "loess", color = "black", linewidth = 0.5, se = FALSE) +  # 添加拟合线 
    labs(
      title = paste("平均速度 ~", var),  # 动态标题 
      y = "平均速度 (mph)",  # Y轴标签 
      x = paste(var, "(单位根据变量而定)")  # X轴标签 
    ) +
    theme_minimal() +  # 简约主题 
    theme(
      plot.background  = element_rect(fill = "#CCCCCC", color = "white", size = 1),  # 灰色背景+白色边框 
      legend.key.width  = unit(0.1, "cm")  # 调整图例宽度 
    )
  
  # 存储图形 
  speed_weather_plot_list[[var]] <- p
}

# 使用 patchwork 将所有图形组合并显示 
wrap_plots(speed_weather_plot_list, ncol = 3)  # 每行显示3个图形 

#2.用众数填补 model 含有缺失值的数据(?)-------------------------------------------
mode_model <- names(sort(table(flights_planes$model), decreasing = TRUE))[1]
flights_planes <- flights_planes %>%
  mutate(model = ifelse(is.na(model),  mode_model, model))

# 处理dep_delay异常值------------------------------------------
# dep_delay分布箱线图 
dep_delay_boxplot_plot <- ggplot(flights, aes(y = dep_delay)) +
  geom_boxplot(fill = "steelblue", color = "black") +
  labs(title = "起飞延误时间箱线图", y = "延误时间（分钟）") +
  theme_minimal() +
  theme(
    axis.title.x  = element_blank(),  # 去除 x 轴标题 
    axis.text.x  = element_blank(),   # 去除 x 轴刻度标签 
    axis.ticks.x  = element_blank(),   # 去除 x 轴刻度线 
    axis.text.y  = element_text(size = 10, face = "bold.italic"),
    plot.title  = element_text(face = "bold", size = 15),
    axis.title  = element_text(face = "bold", size = 15)
  )
# dep_delay分布直方图
dep_delay_histogram_plot <- ggplot(flights, aes(x = dep_delay)) +
  geom_histogram(binwidth = 20, fill = "steelblue") +
  labs(title = "起飞延误时间分布图", x = "延误时间（分钟）", y = "频数") +
  theme_minimal() +
  theme(
    plot.title  = element_text(face = "bold", size = 15),
    axis.title  = element_text(face = "bold", size = 15),
    axis.text.x  = element_text(hjust = 0.5, size = 10, face = "bold"),
    axis.text.y  = element_text(size = 10, face = "bold.italic")
  )

dep_delay_plot <- dep_delay_boxplot_plot + dep_delay_histogram_plot 
dep_delay_plot 

# dep_delay ~ origin--------------------------
flights %>%
  group_by(origin) %>%
  summarise(mean_dep_delay = mean(dep_delay, na.rm  = TRUE))

# dep_delay ~ hour-------------------------
# 计算每小时平均延误 
hour_dep_delay <- flights %>%
  mutate(hour = sched_dep_time %/% 100) %>%  # 提取计划起飞小时 
  group_by(hour) %>%
  summarise(avg_delay = mean(dep_delay, na.rm  = TRUE))

# 绘制峰曲线图 
ggplot(hour_dep_delay, aes(x = hour, y = avg_delay)) +
  geom_line(color = "steelblue", linewidth = 0.8) +  
  geom_point(aes(fill = avg_delay), shape = 21, color = "steelblue", size = 3, alpha = 0.95) + 
  scale_x_continuous(breaks = 0:23) + 
  scale_fill_gradient(low = "blue4", high = "red") +
  labs(
    title = "不同计划起飞小时的平均起飞延误",
    x = "计划起飞小时", 
    y = "平均起飞延误 (dep_delay)"
  ) +
  theme_minimal() +
  theme(
    plot.title  = element_text(face = "bold", size = 30),
    axis.title  = element_text(face = "bold", size = 20),
    axis.text.x  = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y  = element_text(size = 20, face = "bold.italic"), 
    plot.subtitle  = element_text(hjust = 0.8, vjust = -20, size = 20, color = "red"),
    panel.grid.minor  = element_blank()   # 去除次网格线 
  )

# dep_delay ~ weather----------------------------------------------------------
# 随机抽取 4000 行数据 
flights_weather_sample <- slice_sample(flights_weather, n = 4000)

# 遍历 weather 表中的所有变量，排除不需要的列 
weather_vars <- names(weather) %>% 
  setdiff(c("origin", "time_hour", "year", "month", "day", "hour"))

# 创建一个空列表，用于存储所有热力图 
dep_delay_weather_plot_list <- list()

# 生成所有图形并存储到列表中 
for (var in weather_vars) {
  p <- ggplot(flights_weather_sample, aes(x = .data[[var]], y = dep_delay)) +
    stat_binhex(bins = 30) +  # 使用热力图 
    scale_fill_gradient(low = "lightyellow", high = "#154415") +  # 设置颜色渐变 
    labs(
      title = paste("起飞延误 ~", var),  # 动态设置标题 
      y = "起飞延误 (分钟)",
      x = paste(var, "(单位根据变量而定)")
    ) +
    theme_minimal()
  
  dep_delay_weather_plot_list[[var]] <- p +
    theme(plot.background  = element_rect(fill = "lightblue", color = "white", size = 1))
}

# 使用 patchwork 将所有图形组合并显示 
wrap_plots(dep_delay_weather_plot_list, ncol = 3)

# 遍历所有天气变量
for (var in weather_vars) {
  # --- 第一步：动态分箱策略 --- 
  # 检查当前变量的唯一值数量（例如：降水可能只有0和5两种值）
  unique_values <- unique(flights_weather[[var]])
  
  
  if (length(unique_values) <= 5) {
    breaks <- unique_values # 如果变量值种类≤5种，直接按实际值分箱 
  } else {
    # 变量值种类多时，按20%分位数切分5个区间 
    # 例如：温度可能分为[低温, 中低温, 中温, 中高温, 高温]
    # unique()用于防止分位数重复（如80%分位和100%分位值相同）
    breaks <- unique(quantile(flights_weather[[var]], 
                              probs = seq(0, 1, 0.2), 
                              na.rm  = TRUE))  # 忽略缺失值 
  }
  
  # --- 第二步：分层抽样 --- 
  flights_weather_sample <- flights_weather %>%
    # 创建分箱标签列（var_bin）
    mutate(
      var_bin = cut(
        .data[[var]], 
        breaks = breaks,  # 使用上一步计算的分箱边界 
        include.lowest  = TRUE  # 包含最小值所在的区间 
      )
    ) %>%
    # 按分箱结果分组，每组抽取100条数据 
    group_by(var_bin) %>%
    sample_n(500, replace = TRUE) %>%  # replace = TRUE 可重复抽样
    ungroup()  # 取消分组状态 
  
  # --- 第三步：生成热力图 --- 
  p <- ggplot(flights_weather_sample, aes(x = .data[[var]], y = dep_delay)) +
    stat_binhex(bins = 30) +  # 六边形热力图
    scale_fill_gradient(low = "lightyellow", high = "#154415") + 
    geom_smooth(method = "loess", color = "red", linewidth = 0.5, se = FALSE) +  # 添加拟合线 
    labs(
      title = paste("起飞延误 ~", var), 
      y = "起飞延误 (dep_delay)",  
      x = paste(var, "(单位根据变量而定)") 
    ) +
    theme_minimal() + 
    theme(plot.background  = element_rect(fill = "lightblue", color = "white", size = 1))
  
  # 存储图形到列表（按变量名作为键）
  dep_delay_weather_plot_list[[var]] <- p 
}

# 使用 patchwork 将所有图形组合并显示 
wrap_plots(dep_delay_weather_plot_list, ncol = 3)

# model ~ carrier---------------------------------------------------------------
#1.汇总数据：统计每个 carrier 对应的 model 数量 
carrier_model_count <- flights_planes %>%
  group_by(carrier, model) %>%
  summarise(count = n(), .groups = 'drop')
#2.计算每个 carrier 最多使用的 model 数量及其比例 
carrier_total <- carrier_model_count %>% #提取出各机型在各航空公司的数量
  group_by(carrier) %>%
  summarise(total = sum(count), .groups = 'drop')
carrier_most_model <- carrier_model_count %>% #提取出前三频率的
  group_by(carrier) %>%
  slice_max(order_by = count, n = 3) %>% 
  ungroup()
carrier_most_model <- carrier_most_model %>% #合并
  left_join(carrier_total, by = "carrier") %>%
  mutate(carrier_model_ratio = count / total)
#3. model ~ carrier
ggplot(carrier_most_model, aes(x = carrier, y = carrier_model_ratio, fill = model)) +
  geom_bar(stat = "identity", position = "stack", color = "#777777") +  # 绘制堆叠柱状图 
  geom_text(aes(label = model), 
            position = position_stack(vjust = 0.5), size = 4, angle = 30, color = "black", alpha = 0.8) +  # 仅添加机型名称 
  labs(
    title = "航空公司 (carrier) 和其前三频率使用 机型 (model) 的关系图",
    x = "航空公司 (carrier)",
    y = "占对应公司使用机型比例",
  ) +
  theme_minimal() +  # 使用简洁主题 
  theme(
    plot.title  = element_text(face = "bold", size = 20),  # 标题加粗 
    axis.title  = element_text(face = "bold", size = 20),  # 坐标轴标题加粗 
    axis.text.x  = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y  = element_text(size = 20, face = "bold.italic"),
    legend.position  = "none"  # 移除图例 
  ) +
  scale_y_continuous(labels = scales::percent)  # 将 y 轴标签显示为百分比 

# carrier_mean_dep_delay ~ carrier carrier----------------------------------------------------------------
carrier_counts <- flights %>%
  count(carrier, name = "count")
carrier_counts

# dep_delay ~ month----------------------------------------------------------------
# 按 month 分组计算平均起飞延误 
month_dep_delay <- flights %>%
  group_by(month) %>%
  summarise(month_mean_dep_delay = mean(dep_delay, na.rm  = TRUE),
            month_sd_dep_delay = sd(dep_delay, na.rm  = TRUE))

# mean_dep_delay ~ month 
ggplot(month_dep_delay, aes(x = month, y = month_mean_dep_delay, fill = factor(month))) +
  geom_bar(stat = "identity", color = "black") +  # 绘制柱状图 
  scale_x_continuous(breaks = 1:12, labels = month.abb)  +  # 将月份显示为缩写
  labs(
    title = "各月飞机平均延误时间",
    y = "平均延误 (minute)",
    x = "月份"
  ) +
  theme_minimal() +
  theme(
    plot.title   = element_text(face = "bold", size = 30),  # 标题加粗 
    axis.title   = element_text(face = "bold", size = 20),  # 坐标轴标题加粗 
    axis.text.x   = element_text(hjust = 0.5, vjust = 6, size = 20, face = "bold"),
    axis.text.y   = element_text(size = 20, face = "bold.italic"), 
    legend.position  = "none"  # 移除图例 
  )

# carrier ~ month----------------------------------------------------------------
month_carrier_count <- flights %>%
  group_by(month, carrier) %>%
  summarise(count = n(), .groups = 'drop')
# 提取出各月份的总航班数 
month_total <- month_carrier_count %>%
  group_by(month) %>%
  summarise(total = sum(count), .groups = 'drop')
# 合并数据，计算比例 
month_carrier <- month_carrier_count %>%
  left_join(month_total, by = "month") %>%
  mutate(month_carrier_ratio = count / total)
# carrier ~ month
ggplot(month_carrier, aes(x = month, y = month_carrier_ratio, fill = carrier)) +
  geom_bar(stat = "identity", position = "stack", color = "#777777") +  # 绘制堆叠柱状图 
  geom_text(aes(label = carrier), 
            position = position_stack(vjust = 0.5), size = 5, color = "black", alpha = 0.8) +  # 添加航空公司名称 
  labs(
    title = "各月份 (month) 航空公司 (carrier) 所占航班比例的关系图",
    x = "月份 (month)",
    y = "占对应月份使用航空公司比例",
  ) +
  theme_minimal() +  # 使用简洁主题 
  theme(
    plot.title   = element_text(face = "bold", size = 20),  # 标题加粗 
    axis.title   = element_text(face = "bold", size = 20),  # 坐标轴标题加粗 
    axis.text.x   = element_text(hjust = 0.5, size = 20, face = "bold"),
    axis.text.y   = element_text(size = 20, face = "bold.italic"), 
    legend.position   = "none"  # 移除图例 
  ) +
  scale_y_continuous(labels = scales::percent) +  # 将 y 轴标签显示为百分比 
  scale_x_continuous(breaks = 1:12, labels = month.abb)   # 将月份显示为缩写（Jan, Feb 等）