Bài tập 3

title: “Bài tập 3” output: html_document —

library(readxl)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(plm)

Attaching package: 'plm'

The following objects are masked from 'package:dplyr':

    between, lag, lead

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.5.2

library(reshape2)
library(ggpubr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ stringr   1.5.2
✔ lubridate 1.9.4     ✔ tibble    3.3.0
✔ purrr     1.2.0     ✔ tidyr     1.3.1
✔ readr     2.1.5     

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ plm::between()  masks dplyr::between()
✖ dplyr::filter() masks stats::filter()
✖ plm::lag()      masks dplyr::lag(), stats::lag()
✖ plm::lead()     masks dplyr::lead()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

file_path <- "~/DATA MIDTERM PTĐL (1).xlsx"
sheet_names <- excel_sheets(file_path)

# Tạo vector năm theo chính xác năm có trong Excel
years <- c(2008,2010,2012,2014,2016,2018,2019,2020,2021,2022,2023,2024)

data_list <- lapply(seq_along(sheet_names), function(i) {
  df <- read_excel(file_path, sheet = sheet_names[i])
  df$year <- years[i]   # map chính xác năm từ vector
  return(df)
})

data <- bind_rows(data_list)

data <- data %>%
  mutate(
    ln_inc = log(inc + 1),
    ln_lab = log(lab + 1),
    ln_gdp = log(grdp + 1)
  )

panel_data <- pdata.frame(data, index = c("province","year"))

summary_compare <- data %>%
  summarise(
    mean_inc = mean(inc, na.rm=TRUE), sd_inc = sd(inc, na.rm=TRUE),
    mean_ln_inc = mean(ln_inc, na.rm=TRUE), sd_ln_inc = sd(ln_inc, na.rm=TRUE),
    
    mean_lab = mean(lab, na.rm=TRUE), sd_lab = sd(lab, na.rm=TRUE),
    mean_ln_lab = mean(ln_lab, na.rm=TRUE), sd_ln_lab = sd(ln_lab, na.rm=TRUE),
    
    mean_grdp = mean(grdp, na.rm=TRUE), sd_grdp = sd(grdp, na.rm=TRUE),
    mean_ln_gdp = mean(ln_gdp, na.rm=TRUE), sd_ln_gdp = sd(ln_gdp, na.rm=TRUE)
  )

summary_compare

# A tibble: 1 × 12
  mean_inc sd_inc mean_ln_inc sd_ln_inc mean_lab sd_lab mean_ln_lab sd_ln_lab
     <dbl>  <dbl>       <dbl>     <dbl>    <dbl>  <dbl>       <dbl>     <dbl>
1     3.21   1.48        1.37     0.364     20.0   8.50        2.97     0.390
# ℹ 4 more variables: mean_grdp <dbl>, sd_grdp <dbl>, mean_ln_gdp <dbl>,
#   sd_ln_gdp <dbl>

library(readxl)
library(dplyr)

file_path <- "~/DATA MIDTERM PTĐL (1).xlsx"
sheet_names <- excel_sheets(file_path)
sheet_names   # kiểm tra tên sheet, đặc biệt sheet 2024

 [1] "2010" "2012" "2014" "2016" "2018" "2019" "2020" "2021" "2022" "2023"
[11] "2024"

# Gán năm dựa trên tên sheet thực
# Giả sử sheet_names = c("2008","2010",...,"2024")
years <- as.integer(sheet_names)  # tự lấy năm từ tên sheet

data_list <- lapply(seq_along(sheet_names), function(i) {
  df <- read_excel(file_path, sheet = sheet_names[i])
  df$year <- years[i]   # map đúng năm theo sheet
  return(df)
})

data <- bind_rows(data_list)

# Kiểm tra xem có đủ năm không
unique(data$year)

 [1] 2010 2012 2014 2016 2018 2019 2020 2021 2022 2023 2024

years <- as.integer(sheet_names)

# Tạo biến log
data <- data %>%
  mutate(
    ln_inc = log(inc + 1),
    ln_lab = log(lab + 1),
    ln_grdp = log(grdp + 1)
  )

# Phân vùng 6 vùng Việt Nam theo province
data <- data %>%
  mutate(
    region = case_when(
      province %in% c("Hà Nội","Vĩnh Phúc","Bắc Ninh","Quảng Ninh","Hải Dương",
                      "Hải Phòng","Hưng Yên","Thái Bình","Hà Nam","Nam Định",
                      "Ninh Bình") ~ "Đồng bằng sông Hồng",
      
      province %in% c("Hà Giang","Cao Bằng","Bắc Kạn","Tuyên Quang","Lào Cai",
                      "Yên Bái","Thái  Nguyên","Lạng Sơn","Bắc Giang","Phú Thọ") ~ "Đông Bắc",

      province %in% c("Điện Biên","Lai Châu","Sơn La","Hòa Bình") ~ "Tây Bắc",

      province %in% c("Thanh Hóa","Nghệ An","Hà Tĩnh","Quảng Bình","Quảng Trị",
                      "Thừa Thiên-Huế") ~ "Bắc Trung Bộ",

      province %in% c("Đà Nẵng","Quảng  Nam","Quảng  Ngãi","Bình Định","Phú Yên",
                      "Khánh  Hòa","Ninh  Thuận","Bình Thuận") ~ "Nam Trung Bộ",

      province %in% c("Kon Tum","Gia Lai","Đắk Lắk","Đắk Nông","Lâm Đồng") ~ "Tây Nguyên",

      province %in% c("Bình Phước","Tây Ninh","Bình  Dương","Đồng Nai","Bà Rịa - Vũng Tàu",
                      "TP. Hồ Chí Minh","Long An","Tiền Giang","Bến Tre","Trà Vinh",
                      "Vĩnh Long","Đồng Tháp","An Giang","Kiên  Giang","Cần Thơ",
                      "Hậu Giang","Sóc Trăng","Bạc Liêu","Cà Mau") ~ "Đồng bằng sông Cửu Long",

      TRUE ~ "Khác"
    )
  )

# Tạo panel data
panel_data <- pdata.frame(data, index = c("province","year"))

# Phân bố GRDP (2010-2024)
ggplot(data, aes(x = grdp)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "pink", alpha = 0.8) +
  geom_density(color = "skyblue", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố GRDP (2010–2024)", x = "GRDP", y = "Mật độ xác suất") +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.

# Phân bố INC (2010-2024)
ggplot(data, aes(x = inc)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "pink", alpha = 0.8) +
  geom_density(color = "purple", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố INC (2010–2024)", x = "INC", y = "Mật độ xác suất") +
  theme_minimal()

# Phân bố LAB (2010-2024)
ggplot(data, aes(x = lab)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "pink", alpha = 0.8) +
  geom_density(color = "red", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố LAB (2010–2024)", x = "LAB", y = "Mật độ xác suất") +
  theme_minimal()

# Phân bố ln(GRDP) (2010-2024)
ggplot(data, aes(x = ln_grdp)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "steelblue3", alpha = 0.8) +
  geom_density(color = "darkblue", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố ln(GRDP) (2010–2024)", x = "ln(GRDP)", y = "Mật độ xác suất") +
  theme_minimal()

# Phân bố ln(INC) (2010-2024)
ggplot(data, aes(x = ln_inc)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "steelblue3", alpha = 0.8) +
  geom_density(color = "Violet", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố ln(INC) (2010–2024)", x = "ln(INC)", y = "Mật độ xác suất") +
  theme_minimal()

# Phân bố ln(lab)) (2010-2024)
ggplot(data, aes(x = ln_lab)) +
  geom_histogram(aes(y = ..density..), bins = 25, fill = "steelblue3", alpha = 0.8) +
  geom_density(color = "darkred", size = 1) +
  facet_wrap(~year, scales = "free") +
  labs(title = "Phân bố ln(LAB) (2010–2024)", x = "ln(LAB)", y = "Mật độ xác suất") +
  theme_minimal()

# Histogram & density INC vs ln_INC
ggplot(data) +
  geom_density(aes(x = inc), fill = "blue", alpha = 0.35) +
  geom_density(aes(x = ln_inc), fill = "purple", alpha = 0.35) +
  labs(title = "So sánh phân phối thu nhập trước và sau log", x = "Giá trị", y = "Mật độ") +
  theme_minimal()

# Boxplot và violin theo vùng
ggplot(data, aes(x = region, y = ln_inc, fill = region)) +
  geom_violin(trim = FALSE) +
  geom_boxplot(width = 0.2, color = "black") +
  coord_flip() +
  labs(title = "Phân phối thu nhập theo vùng", x = "Vùng", y = "ln(Thu nhập)") +
  theme_minimal()

# Line chart xu hướng thu nhập theo vùng
trend <- data %>%
  group_by(region, year) %>%
  summarise(mean_inc = mean(ln_inc, na.rm = TRUE))

`summarise()` has grouped output by 'region'. You can override using the
`.groups` argument.

ggplot(trend, aes(x = year, y = mean_inc, color = region)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  labs(title = "Xu hướng thu nhập trung bình theo vùng", x = "Năm", y = "ln(Thu nhập)") +
  theme_minimal()

library(dplyr)
library(tidyr)
library(ggplot2)

# Long format GDP
data_long_gdp <- data %>%
  select(region, grdp, ln_grdp) %>%
  pivot_longer(cols = c(grdp, ln_grdp), names_to = "variable", values_to = "value")

# Long format INC
data_long_inc <- data %>%
  select(region, inc, ln_inc) %>%
  pivot_longer(cols = c(inc, ln_inc), names_to = "variable", values_to = "value")

# Boxplot INC trước/sau log theo vùng
ggplot(data_long_inc, aes(x = region, y = value, fill = variable)) +
  geom_boxplot(position = position_dodge(0.8)) +
  coord_flip() +
  labs(title = "So sánh phân phối INC và ln(INC) theo vùng", x = "Vùng", y = "Giá trị") +
  theme_minimal()

# Line chart ln(GRDP) trung bình theo vùng theo năm
trend_gdp <- data %>%
  group_by(region, year) %>%
  summarise(mean_ln_gdp = mean(ln_grdp, na.rm = TRUE))

`summarise()` has grouped output by 'region'. You can override using the
`.groups` argument.

ggplot(trend_gdp, aes(x = year, y = mean_ln_gdp, color = region, group = region)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  labs(title = "Xu hướng ln(GRDP) trung bình theo vùng (2010–2024)",
       x = "Năm", y = "Ln(GRDP)") +
  theme_minimal()

# Quan hệ giữa ln(inc) và Gini
ggplot(data, aes(x = ln_inc, y = gini, color = as.factor(year))) +
  geom_point(size = 3, alpha = 0.7) +
  geom_smooth(method = "loess", se = FALSE, color = "red", size = 1.2) +
  labs(
    title = "Quan hệ giữa ln(Thu nhập) và Gini (có đường xu hướng)",
    x = "ln(Thu nhập)",
    y = "Gini",
    color = "Năm"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

# Phân bổ hệ số Gini theo năm
ggplot(data, aes(x = as.factor(year), y = gini, fill = as.factor(year))) +
  geom_boxplot(alpha = 0.85, color = "black") +
  labs(
    title = "Phân bố hệ số Gini theo năm",
    x = "Năm",
    y = "Hệ số Gini",
    fill = "Năm"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

mean_ln <- data %>%
  group_by(year) %>%
  summarise(mean_ln_inc = mean(ln_inc, na.rm = TRUE))

ggplot(mean_ln, aes(x = year, y = mean_ln_inc, color = as.factor(year))) +
  geom_line(size = 1.3) +
  geom_point(size = 4) +
  labs(
    title = "Sự thay đổi ln(Thu nhập bình quân) theo năm",
    x = "Năm",
    y = "ln(Thu nhập bình quân)",
    color = "Năm"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

data <- data %>% rename(labor = lab)

# Bubble chart: Gini - Gini - Lao động
library (ggplot2)
ggplot(data, aes(x = grdp, y = gini, size = labor, color = as.factor(year))) +
  geom_point(alpha = 0.7) +
  scale_size(range = c(2, 10)) +
  scale_color_brewer(palette = "Set1") +  # màu nổi bật
  labs(
    title = "Bubble Chart: GRDP - GINI - Lao động",
    x = "GRDP bình quân đầu người",
    y = "Chỉ số GINI",
    size = "Lao động",
    color = "Năm"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.title = element_text(face = "bold"),
    legend.title = element_text(face = "bold")
  )

Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set1 is 9
Returning the palette you asked for with that many colors

Warning: Removed 126 rows containing missing values or values outside the scale range
(`geom_point()`).

# Thu nhập bình quân theo năm và vùng
ggplot(data, aes(x = as.factor(year), y = ln_inc, fill = region)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(
    title = "Thu nhập bình quân theo năm và theo vùng",
    x = "Năm",
    y = "Giá trị trung bình ln(thu nhập)",
    fill = "Vùng"
  ) +
  theme_minimal()

# Phân bổ thu 
ggplot(data, aes(x = as.factor(year), y = ln_inc, fill = region)) +
  geom_bar(stat = "summary", fun = "mean") +
  labs(
    title = "Phân bổ thu nhập theo vùng qua các năm",
    x = "Năm",
    y = "Tổng ln(thu nhập)",
    fill = "Vùng"
  ) +
  theme_minimal()

ggplot(data, aes(x = gini, y = ln_grdp, color = as.factor(year))) +
  geom_point(alpha = 0.7, size = 3) +               # biểu đồ điểm
  geom_smooth(method = "loess", se = TRUE, color = "red", size = 1) +  # đường xu hướng
  labs(
    title = "Ảnh hưởng của Gini đến ln(GRDP) theo năm",
    x = "Hệ số Gini",
    y = "ln(GRDP)",
    color = "Năm"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.title = element_text(face = "bold"),
    legend.title = element_text(face = "bold")
  )

`geom_smooth()` using formula = 'y ~ x'