library(haven)
library(lubridate)
##
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(dplyr)
##
## 다음의 패키지를 부착합니다: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(psych)
##
## 다음의 패키지를 부착합니다: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(knitr)
library(scales)
##
## 다음의 패키지를 부착합니다: 'scales'
## The following objects are masked from 'package:psych':
##
## alpha, rescale
library(tidyr)
# 데이터 불러오기
klips <- read.csv("C:/rhr/data/data/klips.csv")
# 성별 범주화 (1 = 남자, 2 = 여자)
klips$gender <- factor(klips$p130101, levels = c(1, 2), labels = c("남자", "여자"))
# 만나이 범주화
klips$age_group <- with(klips, ifelse(p130107 < 40, "1",
ifelse(p130107 >= 40 & p130107 < 50, "2",
ifelse(p130107 >= 50 & p130107 < 70, "3",
ifelse(p130107 >= 70, "2", NA)))))
# 학력 범주화
klips$education <- with(klips, ifelse(p130110 == -1, NA,
ifelse(p130110 >= 1 & p130110 < 5, "중졸 이하",
ifelse(p130110 == 5, "고졸",
ifelse(p130110 >= 6 & p130110 < 8, "대졸 이하",
ifelse(p130110 >= 8 & p130110 < 10, "대학원 졸", NA))))))
# 결과 출력
cat("성별 분포:\n")
## 성별 분포:
print(table(klips$gender, useNA = "ifany"))
##
## 남자 여자 <NA>
## 6763 7354 14119
cat("\n나이 범주 분포:\n")
##
## 나이 범주 분포:
print(table(klips$age_group, useNA = "ifany"))
##
## 1 2 3 <NA>
## 11338 8656 8240 2
cat("\n학력 분포:\n")
##
## 학력 분포:
print(table(klips$education, useNA = "ifany"))
##
## 고졸 대졸 이하 대학원 졸 중졸 이하 <NA>
## 4837 4570 400 4306 14123
# 1) 일자리 형태 범주화
klips$jobtype_cat <- factor(klips$jobtype, levels = c(1, 2), labels = c("임금근로자", "비임금근로자"))
# 2) 종사상 지위 재범주화 (정수형 범주)
klips$employment_status <- with(klips, ifelse(p130314 == -1, NA,
ifelse(p130314 == 1, 1,
ifelse(p130314 %in% c(2, 3), 2,
ifelse(p130314 == 4, 3,
ifelse(p130314 == 5, 4, NA))))))
# 범주 이름 붙이기
klips$employment_status <- factor(klips$employment_status,
levels = c(1, 2, 3, 4),
labels = c("상용직", "임시·일용직", "고용주/자영업자", "무급가족종사자"))
# 3) 정규직 여부
klips$regular_job <- with(klips, ifelse(p130317 == -1, NA,
ifelse(p130317 == 1, "정규직",
ifelse(p130317 == 2, "비정규직", NA))))
# 결과 출력
cat("1) 일자리 형태 분포:\n")
## 1) 일자리 형태 분포:
print(table(klips$jobtype_cat, useNA = "ifany"))
##
## 임금근로자 비임금근로자 <NA>
## 5289 2549 20398
cat("\n2) 종사상 지위 (재범주화) 분포:\n")
##
## 2) 종사상 지위 (재범주화) 분포:
print(table(klips$employment_status, useNA = "ifany"))
##
## 상용직 임시·일용직 고용주/자영업자 무급가족종사자 <NA>
## 3871 1418 2035 514 20398
cat("\n3) 정규직 여부 분포:\n")
##
## 3) 정규직 여부 분포:
print(table(klips$regular_job, useNA = "ifany"))
##
## 비정규직 정규직 <NA>
## 1902 3362 22972
# 국민연금
klips$national_pension <- ifelse(klips$p132101 == 1, "가입",
ifelse(klips$p132101 == 2, "미가입", NA))
# 특수직역연금
klips$special_pension <- ifelse(klips$p132102 == 1, "가입",
ifelse(klips$p132102 == 2, "미가입", NA))
# 건강보험
klips$health_insurance <- ifelse(klips$p132103 == 1, "가입",
ifelse(klips$p132103 == 2, "미가입", NA))
# 고용보험
klips$employment_insurance <- ifelse(klips$p132104 == 1, "가입",
ifelse(klips$p132104 == 2, "미가입", NA))
# 산재보험
klips$industrial_insurance <- ifelse(klips$p132105 == 1, "가입",
ifelse(klips$p132105 == 2, "미가입", NA))
# 기술통계 출력
cat("1) 국민연금 가입 여부:\n")
## 1) 국민연금 가입 여부:
print(table(klips$national_pension, useNA = "ifany"))
##
## 가입 미가입 <NA>
## 3176 2097 22963
cat("\n2) 특수직역연금 가입 여부:\n")
##
## 2) 특수직역연금 가입 여부:
print(table(klips$special_pension, useNA = "ifany"))
##
## 가입 미가입 <NA>
## 410 4859 22967
cat("\n3) 건강보험 가입 여부:\n")
##
## 3) 건강보험 가입 여부:
print(table(klips$health_insurance, useNA = "ifany"))
##
## 가입 미가입 <NA>
## 3715 1560 22961
cat("\n4) 고용보험 가입 여부:\n")
##
## 4) 고용보험 가입 여부:
print(table(klips$employment_insurance, useNA = "ifany"))
##
## 가입 미가입 <NA>
## 3258 1609 23369
cat("\n5) 산재보험 가입 여부:\n")
##
## 5) 산재보험 가입 여부:
print(table(klips$industrial_insurance, useNA = "ifany"))
##
## 가입 미가입 <NA>
## 3491 1727 23018
# 임금근로자 소득 (p131642)
# -1 결측 처리
klips$wage_income <- ifelse(klips$p131642 == -1, NA, klips$p131642)
# 임금근로자 평균 소득
wage_mean <- mean(klips$wage_income, na.rm = TRUE)
cat("1) 임금근로자 월평균 소득의 평균:", wage_mean, "\n")
## 1) 임금근로자 월평균 소득의 평균: 199.409
# 임금근로자 소득 구간화
klips$wage_group <- ifelse(klips$wage_income < 100, 1,
ifelse(klips$wage_income < 200, 2,
ifelse(klips$wage_income < 300, 3,
ifelse(klips$wage_income < 400, 4,
ifelse(klips$wage_income >= 400, 5, NA)))))
# 임금근로자 소득 구간별 비중
cat("\n1-2) 임금근로자 월평균 소득 구간별 비중:\n")
##
## 1-2) 임금근로자 월평균 소득 구간별 비중:
print(prop.table(table(klips$wage_group)))
##
## 1 2 3 4 5
## 0.16967166 0.39590055 0.24103245 0.10571266 0.08768267
# 비임금근로자 소득 (p131672)
# -1 결측 처리
klips$nonwage_income <- ifelse(klips$p131672 == -1, NA, klips$p131672)
# 비임금근로자 평균 소득
nonwage_mean <- mean(klips$nonwage_income, na.rm = TRUE)
cat("\n2) 비임금근로자 월평균 소득의 평균:", nonwage_mean, "\n")
##
## 2) 비임금근로자 월평균 소득의 평균: 235.2523
# 비임금근로자 소득 구간화
klips$nonwage_group <- ifelse(klips$nonwage_income < 100, 1,
ifelse(klips$nonwage_income < 200, 2,
ifelse(klips$nonwage_income < 300, 3,
ifelse(klips$nonwage_income < 400, 4,
ifelse(klips$nonwage_income >= 400, 5, NA)))))
# 비임금근로자 소득 구간별 비중
cat("\n2-2) 비임금근로자 월평균 소득 구간별 비중:\n")
##
## 2-2) 비임금근로자 월평균 소득 구간별 비중:
print(prop.table(table(klips$nonwage_group)))
##
## 1 2 3 4 5
## 0.2171923 0.2690743 0.2075280 0.1276704 0.1785351
income_labels <- c("1" = "<100", "2" = "100~199", "3" = "200~299",
"4" = "300~399", "5" = "400+")
# ▓ 임금근로자 소득 구간 비중 그래프
wage_plot_data <- as.data.frame(table(klips$wage_group))
colnames(wage_plot_data) <- c("income_group", "count")
wage_plot_data$proportion <- wage_plot_data$count / sum(wage_plot_data$count, na.rm = TRUE)
ggplot(wage_plot_data, aes(x = income_group, y = proportion, fill = income_group)) +
geom_bar(stat = "identity") +
scale_x_discrete(labels = income_labels) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(title = "임금근로자 소득 구간별 비중", x = "소득 구간 (만원)", y = "비중") +
theme_minimal() +
theme(legend.position = "none")
# ▓ 비임금근로자 소득 구간 비중 그래프
nonwage_plot_data <- as.data.frame(table(klips$nonwage_group))
colnames(nonwage_plot_data) <- c("income_group", "count")
nonwage_plot_data$proportion <- nonwage_plot_data$count / sum(nonwage_plot_data$count, na.rm = TRUE)
ggplot(nonwage_plot_data, aes(x = income_group, y = proportion, fill = income_group)) +
geom_bar(stat = "identity") +
scale_x_discrete(labels = income_labels) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
labs(title = "비임금근로자 소득 구간별 비중", x = "소득 구간 (만원)", y = "비중") +
theme_minimal() +
theme(legend.position = "none")
# 성별에 따라 차이가 있는지 교차표 + 카이검정
cat("\n\n===== 성별 vs 일자리 형태 =====\n")
##
##
## ===== 성별 vs 일자리 형태 =====
table_gender_jobtype <- table(klips$gender, klips$jobtype)
print(table_gender_jobtype)
##
## 1 2
## 남자 3161 1490
## 여자 2128 1059
print(chisq.test(table_gender_jobtype))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table_gender_jobtype
## X-squared = 1.172, df = 1, p-value = 0.279
cat("\n\n===== 성별 vs 정규직 여부 =====\n")
##
##
## ===== 성별 vs 정규직 여부 =====
table_gender_regular <- table(klips$gender, klips$p130317)
print(table_gender_regular)
##
## -1 1 2
## 남자 20 2242 899
## 여자 5 1120 1003
print(chisq.test(table_gender_regular))
##
## Pearson's Chi-squared test
##
## data: table_gender_regular
## X-squared = 194.81, df = 2, p-value < 2.2e-16
# 사회보험 가입 여부에 대해 반복문으로 처리
insurance_vars <- c("p132101", "p132102", "p132103", "p132104", "p132105")
insurance_names <- c("국민연금", "특수직역연금", "건강보험", "고용보험", "산재보험")
for (i in 1:length(insurance_vars)) {
cat("\n\n===== 성별 vs", insurance_names[i], "가입 여부 =====\n")
tbl <- table(klips$gender, klips[[insurance_vars[i]]])
print(tbl)
print(chisq.test(tbl))
}
##
##
## ===== 성별 vs 국민연금 가입 여부 =====
##
## 1 2 3
## 남자 2060 1090 11
## 여자 1116 1007 5
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 87.709, df = 2, p-value < 2.2e-16
##
##
##
## ===== 성별 vs 특수직역연금 가입 여부 =====
##
## 1 2 3
## 남자 254 2893 14
## 여자 156 1966 6
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 1.7894, df = 2, p-value = 0.4087
##
##
##
## ===== 성별 vs 건강보험 가입 여부 =====
##
## 1 2 3
## 남자 2417 735 9
## 여자 1298 825 5
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 147.25, df = 2, p-value < 2.2e-16
##
##
##
## ===== 성별 vs 고용보험 가입 여부 =====
##
## 1 2 3
## 남자 2131 767 16
## 여자 1127 842 13
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 140.89, df = 2, p-value < 2.2e-16
##
##
##
## ===== 성별 vs 산재보험 가입 여부 =====
##
## 1 2 3
## 남자 2305 821 35
## 여자 1186 906 36
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 167.51, df = 2, p-value < 2.2e-16
cat("\n\n===== 학력 vs 일자리 형태 =====\n")
##
##
## ===== 학력 vs 일자리 형태 =====
table_edu_jobtype <- table(klips$education, klips$jobtype)
print(table_edu_jobtype)
##
## 1 2
## 고졸 1795 933
## 대졸 이하 2295 589
## 대학원 졸 264 47
## 중졸 이하 934 980
print(chisq.test(table_edu_jobtype))
##
## Pearson's Chi-squared test
##
## data: table_edu_jobtype
## X-squared = 543.13, df = 3, p-value < 2.2e-16
cat("\n\n===== 학력 vs 정규직 여부 =====\n")
##
##
## ===== 학력 vs 정규직 여부 =====
table_edu_regular <- table(klips$education, klips$p130317)
print(table_edu_regular)
##
## -1 1 2
## 고졸 11 1027 757
## 대졸 이하 7 1846 442
## 대학원 졸 4 217 43
## 중졸 이하 3 272 659
print(chisq.test(table_edu_regular))
## Warning in chisq.test(table_edu_regular): 카이제곱 approximation은 정확하지
## 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: table_edu_regular
## X-squared = 848.18, df = 6, p-value < 2.2e-16
for (i in 1:length(insurance_vars)) {
cat("\n\n===== 학력 vs", insurance_names[i], "가입 여부 =====\n")
tbl <- table(klips$education, klips[[insurance_vars[i]]])
print(tbl)
print(chisq.test(tbl))
}
##
##
## ===== 학력 vs 국민연금 가입 여부 =====
##
## 1 2 3
## 고졸 1068 725 2
## 대졸 이하 1633 654 8
## 대학원 졸 161 103 0
## 중졸 이하 314 614 6
## Warning in chisq.test(tbl): 카이제곱 approximation은 정확하지 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 395.26, df = 6, p-value < 2.2e-16
##
##
##
## ===== 학력 vs 특수직역연금 가입 여부 =====
##
## 1 2 3
## 고졸 48 1744 3
## 대졸 이하 290 1995 10
## 대학원 졸 65 198 1
## 중졸 이하 7 921 6
## Warning in chisq.test(tbl): 카이제곱 approximation은 정확하지 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 314.67, df = 6, p-value < 2.2e-16
##
##
##
## ===== 학력 vs 건강보험 가입 여부 =====
##
## 1 2 3
## 고졸 1172 621 2
## 대졸 이하 1936 353 6
## 대학원 졸 227 36 1
## 중졸 이하 380 549 5
## Warning in chisq.test(tbl): 카이제곱 approximation은 정확하지 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 666.22, df = 6, p-value < 2.2e-16
##
##
##
## ===== 학력 vs 고용보험 가입 여부 =====
##
## 1 2 3
## 고졸 1096 648 3
## 대졸 이하 1643 364 14
## 대학원 졸 158 39 3
## 중졸 이하 361 557 9
## Warning in chisq.test(tbl): 카이제곱 approximation은 정확하지 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 556.47, df = 6, p-value < 2.2e-16
##
##
##
## ===== 학력 vs 산재보험 가입 여부 =====
##
## 1 2 3
## 고졸 1124 660 11
## 대졸 이하 1794 461 40
## 대학원 졸 203 51 10
## 중졸 이하 370 554 10
## Warning in chisq.test(tbl): 카이제곱 approximation은 정확하지 않을수도 있습니다
##
## Pearson's Chi-squared test
##
## data: tbl
## X-squared = 518.1, df = 6, p-value < 2.2e-16
# 4. 임금근로자의 월평균 임금 검정 (모집단 평균: 250만 원)
t.test(
klips$wage_group,
mu = 250,
alternative = "two.sided",
conf.level = 0.95
)
##
## One Sample t-test
##
## data: klips$wage_group
## t = -15604, df = 5268, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 250
## 95 percent confidence interval:
## 2.514746 2.576922
## sample estimates:
## mean of x
## 2.545834
# 5. 비임금 근로자의 월평균 임금 검정 (모집단 평균: 300만 원)
t.test(
klips$nonwage_group,
mu = 300,
alternative = "two.sided",
conf.level = 0.95
)
##
## One Sample t-test
##
## data: klips$nonwage_group
## t = -9479.3, df = 1965, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 300
## 95 percent confidence interval:
## 2.719790 2.842774
## sample estimates:
## mean of x
## 2.781282
# 성별에 따른 임금근로자의 월평균 임금 차이 검정 (wage_group 사용)
t.test(
klips$wage_group[klips$gender == "남자"],
klips$wage_group[klips$gender == "여자"],
alternative = "two.sided",
conf.level = 0.95
)
##
## Welch Two Sample t-test
##
## data: klips$wage_group[klips$gender == "남자"] and klips$wage_group[klips$gender == "여자"]
## t = 34.393, df = 5139.8, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.9128664 1.0232265
## sample estimates:
## mean of x mean of y
## 2.935515 1.967468
# 정규직 여부에 따른 임금근로자의 월평균 임금 차이 검정 (wage_group 사용)
t.test(
klips$wage_group[klips$regular_job == "정규직"],
klips$wage_group[klips$regular_job == "비정규직"],
alternative = "two.sided",
conf.level = 0.95
)
##
## Welch Two Sample t-test
##
## data: klips$wage_group[klips$regular_job == "정규직"] and klips$wage_group[klips$regular_job == "비정규직"]
## t = 41.601, df = 4790.8, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.074418 1.180690
## sample estimates:
## mean of x mean of y
## 2.948718 1.821164
# 평균값 계산
data <- read_sav("C:/rhr/data/data/dietstudy.sav")
summary_stats <- data %>%
group_by(gender) %>%
summarise(
tg0_mean = mean(tg0, na.rm = TRUE),
tg4_mean = mean(tg4, na.rm = TRUE),
wgt0_mean = mean(wgt0, na.rm = TRUE),
wgt4_mean = mean(wgt4, na.rm = TRUE)
)
print(summary_stats)
## # A tibble: 2 × 5
## gender tg0_mean tg4_mean wgt0_mean wgt4_mean
## <dbl+lbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 [Male] 147. 117. 224. 216.
## 2 1 [Female] 127 134. 166. 158.
# 시각화를 위한 데이터 재구성
plot_data <- data %>%
group_by(gender) %>%
summarise(
tg0 = mean(tg0, na.rm = TRUE),
tg4 = mean(tg4, na.rm = TRUE),
wgt0 = mean(wgt0, na.rm = TRUE),
wgt4 = mean(wgt4, na.rm = TRUE)
) %>%
tidyr::pivot_longer(
cols = -gender,
names_to = "variable",
values_to = "mean_value"
) %>%
mutate(
time = ifelse(grepl("0", variable), "Time0", "Time4"),
type = ifelse(grepl("tg", variable), "TG", "Weight"),
gender = factor(gender, labels = c("Male", "Female"))
)
# 중성지방 변화 그래프
ggplot(plot_data %>% filter(type == "TG"), aes(x = time, y = mean_value, group = gender, color = gender)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
labs(title = "중성지방 평균 변화", x = "시점", y = "중성지방 (TG)") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# 체중 변화 그래프
ggplot(plot_data %>% filter(type == "Weight"), aes(x = time, y = mean_value, group = gender, color = gender)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
labs(title = "체중 평균 변화", x = "시점", y = "체중 (Weight)") +
theme_minimal()
# 남성 데이터 필터링
male_data <- data %>% filter(gender == 0)
# 여성 데이터 필터링
female_data <- data %>% filter(gender == 1)
# 남성의 tg1 vs tg2 평균 차이 검정 (paired t-test)
t_test_male <- t.test(male_data$tg1, male_data$tg2, paired = TRUE, na.rm = TRUE)
print("남성 중성지방(tg1 vs tg2) 변화 검정 결과:")
## [1] "남성 중성지방(tg1 vs tg2) 변화 검정 결과:"
print(t_test_male)
##
## Paired t-test
##
## data: male_data$tg1 and male_data$tg2
## t = 0.26185, df = 8, p-value = 0.8001
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -30.35880 38.13658
## sample estimates:
## mean difference
## 3.888889
# 여성의 tg1 vs tg2 평균 차이 검정 (paired t-test)
t_test_female <- t.test(female_data$tg1, female_data$tg2, paired = TRUE, na.rm = TRUE)
print("여성 중성지방(tg1 vs tg2) 변화 검정 결과:")
## [1] "여성 중성지방(tg1 vs tg2) 변화 검정 결과:"
print(t_test_female)
##
## Paired t-test
##
## data: female_data$tg1 and female_data$tg2
## t = -0.32195, df = 6, p-value = 0.7584
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -39.31519 30.17233
## sample estimates:
## mean difference
## -4.571429