thuc hanh
ly
2025-06-16
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-library(lessR)##
## lessR 4.4.3 feedback: gerbing@pdx.edu
## --------------------------------------------------------------
## > d <- Read("") Read data file, many formats available, e.g., Excel
## d is default data frame, data= in analysis routines optional
##
## Many examples of reading, writing, and manipulating data,
## graphics, testing means and proportions, regression, factor analysis,
## customization, forecasting, and aggregation from pivot tables
## Enter: browseVignettes("lessR")
##
## View lessR updates, now including time series forecasting
## Enter: news(package="lessR")
##
## Interactive data analysis
## Enter: interact()##
## Attaching package: 'lessR'## The following object is masked from 'package:table1':
##
## label## The following object is masked from 'package:base':
##
## sort_bylibrary(ggplot2)bw = read.csv("D:\\PhD 7May2025\\GS Nguyễn Văn Tuấn khóa học phương pháp nghiên cứu khoa học ppnckh\\16Jun2025 Phân tích DL với R\\birthwt.csv")
head(bw)## id low age lwt race smoke ptl ht ui ftv bwt
## 1 85 0 19 182 2 0 0 0 1 0 2523
## 2 86 0 33 155 3 0 0 0 0 3 2551
## 3 87 0 20 105 1 1 0 0 0 1 2557
## 4 88 0 21 108 1 1 0 0 1 2 2594
## 5 89 0 18 107 1 1 0 0 1 0 2600
## 6 91 0 21 124 3 0 0 0 0 0 2622tail(bw, 10)## id low age lwt race smoke ptl ht ui ftv bwt
## 180 71 1 17 120 2 0 0 0 0 2 2438
## 181 75 1 26 154 3 0 1 1 0 1 2442
## 182 76 1 20 105 3 0 0 0 0 3 2450
## 183 77 1 26 190 1 1 0 0 0 0 2466
## 184 78 1 14 101 3 1 1 0 0 0 2466
## 185 79 1 28 95 1 1 0 0 0 2 2466
## 186 81 1 14 100 3 0 0 0 0 2 2495
## 187 82 1 23 94 3 1 0 0 0 0 2495
## 188 83 1 17 142 2 0 0 1 0 0 2495
## 189 84 1 21 130 1 1 0 1 0 3 2495bw$ethnicity [bw$race==1] <- "White"
bw$ethnicity [bw$race==2] <- "Black"
bw$ethnicity [bw$race==3] <- "Others"
bw$ethnicity = factor(bw$ethnicity, levels = c("White", "Others", "Black"))
bw$smoking = ifelse(bw$smoke == 1, "Smoking", "Non smoking")
bw$low.bw = ifelse(bw$low==1, "Low BW", "Normal")
bw$mother.wt = bw$lwt * 0.45
head(bw)## id low age lwt race smoke ptl ht ui ftv bwt ethnicity smoking low.bw
## 1 85 0 19 182 2 0 0 0 1 0 2523 Black Non smoking Normal
## 2 86 0 33 155 3 0 0 0 0 3 2551 Others Non smoking Normal
## 3 87 0 20 105 1 1 0 0 0 1 2557 White Smoking Normal
## 4 88 0 21 108 1 1 0 0 1 2 2594 White Smoking Normal
## 5 89 0 18 107 1 1 0 0 1 0 2600 White Smoking Normal
## 6 91 0 21 124 3 0 0 0 0 0 2622 Others Non smoking Normal
## mother.wt
## 1 81.90
## 2 69.75
## 3 47.25
## 4 48.60
## 5 48.15
## 6 55.80library(table1)
bw = read.csv("D:\\PhD 7May2025\\GS Nguyễn Văn Tuấn khóa học phương pháp nghiên cứu khoa học ppnckh\\16Jun2025 Phân tích DL với R\\birthwt.csv")
head(bw)## id low age lwt race smoke ptl ht ui ftv bwt
## 1 85 0 19 182 2 0 0 0 1 0 2523
## 2 86 0 33 155 3 0 0 0 0 3 2551
## 3 87 0 20 105 1 1 0 0 0 1 2557
## 4 88 0 21 108 1 1 0 0 1 2 2594
## 5 89 0 18 107 1 1 0 0 1 0 2600
## 6 91 0 21 124 3 0 0 0 0 0 2622table1(~age + lwt + bwt, data=bw)| Overall (N=189) |
|
|---|---|
| age | |
| Mean (SD) | 23.2 (5.30) |
| Median [Min, Max] | 23.0 [14.0, 45.0] |
| lwt | |
| Mean (SD) | 130 (30.6) |
| Median [Min, Max] | 121 [80.0, 250] |
| bwt | |
| Mean (SD) | 2940 (729) |
| Median [Min, Max] | 2980 [709, 4990] |
table1(~age + lwt + bwt | low, data=bw)| 0 (N=130) |
1 (N=59) |
Overall (N=189) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 23.7 (5.58) | 22.3 (4.51) | 23.2 (5.30) |
| Median [Min, Max] | 23.0 [14.0, 45.0] | 22.0 [14.0, 34.0] | 23.0 [14.0, 45.0] |
| lwt | |||
| Mean (SD) | 133 (31.7) | 122 (26.6) | 130 (30.6) |
| Median [Min, Max] | 124 [85.0, 250] | 120 [80.0, 200] | 121 [80.0, 250] |
| bwt | |||
| Mean (SD) | 3330 (478) | 2100 (391) | 2940 (729) |
| Median [Min, Max] | 3270 [2520, 4990] | 2210 [709, 2500] | 2980 [709, 4990] |
table1(~low + smoke + race, data=bw)| Overall (N=189) |
|
|---|---|
| low | |
| Mean (SD) | 0.312 (0.465) |
| Median [Min, Max] | 0 [0, 1.00] |
| smoke | |
| Mean (SD) | 0.392 (0.489) |
| Median [Min, Max] | 0 [0, 1.00] |
| race | |
| Mean (SD) | 1.85 (0.918) |
| Median [Min, Max] | 1.00 [1.00, 3.00] |
table1(~factor(smoke) + factor(race) | factor(low), data=bw)| 0 (N=130) |
1 (N=59) |
Overall (N=189) |
|
|---|---|---|---|
| factor(smoke) | |||
| 0 | 86 (66.2%) | 29 (49.2%) | 115 (60.8%) |
| 1 | 44 (33.8%) | 30 (50.8%) | 74 (39.2%) |
| factor(race) | |||
| 1 | 73 (56.2%) | 23 (39.0%) | 96 (50.8%) |
| 2 | 15 (11.5%) | 11 (18.6%) | 26 (13.8%) |
| 3 | 42 (32.3%) | 25 (42.4%) | 67 (35.4%) |
# slide 238
head(bw)## id low age lwt race smoke ptl ht ui ftv bwt
## 1 85 0 19 182 2 0 0 0 1 0 2523
## 2 86 0 33 155 3 0 0 0 0 3 2551
## 3 87 0 20 105 1 1 0 0 0 1 2557
## 4 88 0 21 108 1 1 0 0 1 2 2594
## 5 89 0 18 107 1 1 0 0 1 0 2600
## 6 91 0 21 124 3 0 0 0 0 0 2622#trang slide 242
library(lessR)
Histogram(bwt, fill="blue", xlab="Birthweight (g)", ylab="Frequency", data=bw)
## >>> Suggestions
## bin_width: set the width of each bin
## bin_start: set the start of the first bin
## bin_end: set the end of the last bin
## Histogram(bwt, density=TRUE) # smoothed curve + histogram
## Plot(bwt) # Violin/Box/Scatterplot (VBS) plot
##
## --- bwt ---
##
## n miss mean sd min mdn max
## 189 0 2944.59 729.21 709.00 2977.00 4990.00
##
##
##
## --- Outliers --- from the box plot: 1
##
## Small Large
## ----- -----
## 709.0
##
##
## Bin Width: 500
## Number of Bins: 9
##
## Bin Midpnt Count Prop Cumul.c Cumul.p
## -----------------------------------------------------
## 500 > 1000 750 1 0.01 1 0.01
## 1000 > 1500 1250 4 0.02 5 0.03
## 1500 > 2000 1750 14 0.07 19 0.10
## 2000 > 2500 2250 40 0.21 59 0.31
## 2500 > 3000 2750 38 0.20 97 0.51
## 3000 > 3500 3250 45 0.24 142 0.75
## 3500 > 4000 3750 38 0.20 180 0.95
## 4000 > 4500 4250 7 0.04 187 0.99
## 4500 > 5000 4750 2 0.01 189 1.00
## BarChart(race, data=bw)
## >>> Suggestions
## BarChart(race, horiz=TRUE) # horizontal bar chart
## BarChart(race, fill="reds") # red bars of varying lightness
## PieChart(race) # doughnut (ring) chart
## Plot(race) # bubble plot
## Plot(race, stat="count") # lollipop plot
##
## --- race ---
##
## Missing Values: 0
##
## 1 2 3 Total
## Frequencies: 96 26 67 189
## Proportions: 0.508 0.138 0.354 1.000
##
## Chi-squared test of null hypothesis of equal probabilities
## Chisq = 39.270, df = 2, p-value = 0.000Plot(lwt, bwt, data=bw)
##
## >>> Suggestions or enter: style(suggest=FALSE)
## Plot(lwt, bwt, enhance=TRUE) # many options
## Plot(lwt, bwt, fill="skyblue") # interior fill color of points
## Plot(lwt, bwt, fit="lm", fit_se=c(.90,.99)) # fit line, stnd errors
## Plot(lwt, bwt, MD_cut=6) # Mahalanobis distance from center > 6 is an outlier
##
##
## >>> Pearson's product-moment correlation
##
## Number of paired values with neither missing, n = 189
## Sample Correlation of lwt and bwt: r = 0.186
##
## Hypothesis Test of 0 Correlation: t = 2.585, df = 187, p-value = 0.011
## 95% Confidence Interval for Correlation: 0.044 to 0.320
## library(lessR)
BarChart(race, data=bw)
## >>> Suggestions
## BarChart(race, horiz=TRUE) # horizontal bar chart
## BarChart(race, fill="reds") # red bars of varying lightness
## PieChart(race) # doughnut (ring) chart
## Plot(race) # bubble plot
## Plot(race, stat="count") # lollipop plot
##
## --- race ---
##
## Missing Values: 0
##
## 1 2 3 Total
## Frequencies: 96 26 67 189
## Proportions: 0.508 0.138 0.354 1.000
##
## Chi-squared test of null hypothesis of equal probabilities
## Chisq = 39.270, df = 2, p-value = 0.000#trang slide 242
library(lessR)
Histogram(bwt, fill="blue", xlab="Birthweight (g)", ylab="Frequency", data=bw)
## >>> Suggestions
## bin_width: set the width of each bin
## bin_start: set the start of the first bin
## bin_end: set the end of the last bin
## Histogram(bwt, density=TRUE) # smoothed curve + histogram
## Plot(bwt) # Violin/Box/Scatterplot (VBS) plot
##
## --- bwt ---
##
## n miss mean sd min mdn max
## 189 0 2944.59 729.21 709.00 2977.00 4990.00
##
##
##
## --- Outliers --- from the box plot: 1
##
## Small Large
## ----- -----
## 709.0
##
##
## Bin Width: 500
## Number of Bins: 9
##
## Bin Midpnt Count Prop Cumul.c Cumul.p
## -----------------------------------------------------
## 500 > 1000 750 1 0.01 1 0.01
## 1000 > 1500 1250 4 0.02 5 0.03
## 1500 > 2000 1750 14 0.07 19 0.10
## 2000 > 2500 2250 40 0.21 59 0.31
## 2500 > 3000 2750 38 0.20 97 0.51
## 3000 > 3500 3250 45 0.24 142 0.75
## 3500 > 4000 3750 38 0.20 180 0.95
## 4000 > 4500 4250 7 0.04 187 0.99
## 4500 > 5000 4750 2 0.01 189 1.00
## # Đọc file CSV từ máy (thay đường dẫn bằng file bạn đã lưu)
data <- read.csv("D:/PhD 7May2025/GS Nguyễn Văn Tuấn khóa học phương pháp nghiên cứu khoa học ppnckh/16Jun2025 Phân tích DL với R/CHNS data full.csv")
library(ggplot2)
summary(data$income)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0 1525 5773 15830 19600 2400000 224ggplot(data, aes(x = income)) +
geom_histogram(binwidth = 1000, fill = "skyblue", color = "blue") +
labs(title = "Biểu đồ Histogram của Thu nhập",
x = "Thu nhập",
y = "Tần số") +
theme_minimal()
library(ggplot2)
ggplot(data, aes(x = log(income))) +
geom_histogram(binwidth = 0.2, fill = "blue", color = "white") +
labs(
title = "Phân bố thu nhập (logarithm)",
x = "log(Thu nhập)",
y = "Số người"
) +
theme_minimal()
Nếu chưa cài ggplot2 thì chạy lệnh này một lần:
install.packages(“ggplot2”)
library(ggplot2)
ggplot(data, aes(x = income)) +
geom_histogram(binwidth = 1000, fill = "white", color = "blue") +
labs(
title = "Phân bố thu nhập",
x = "Thu nhập",
y = "Số người"
) +
theme_minimal()
library(ggplot2)
# Chuyển gender thành factor
data$gender <- factor(data$gender, labels = c("Nam", "Nữ"))
# Lọc income > 0 để tránh lỗi log
filtered_data <- subset(data, income > 0)
# Vẽ biểu đồ
ggplot(filtered_data, aes(x = log(income), fill = gender)) +
geom_histogram(binwidth = 0.2, position = "identity", alpha = 0.6, color = "white") +
scale_fill_manual(values = c("Nam" = "red", "Nữ" = "blue")) +
labs(
title = "Phân bố log(Thu nhập) theo giới tính",
x = "log(Thu nhập)",
y = "Số người",
fill = "Giới tính"
) +
theme_minimal()
# Chuyển 'edu' thành factor với nhãn rõ ràng
data$edu <- factor(data$edu, levels = c(1, 2, 3),
labels = c("Primary", "Secondary", "Tertiary"))
library(ggplot2)
# Loại bỏ giá trị NA và income <= 0
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(edu))
# Vẽ biểu đồ
ggplot(filtered_data, aes(x = edu, y = income, fill = edu)) +
stat_summary(fun = median, geom = "bar", color = "black") +
scale_fill_manual(values = c("Primary" = "blue",
"Secondary" = "red",
"Tertiary" = "orange")) +
labs(
title = "Trung vị thu nhập theo trình độ học vấn",
x = "Trình độ học vấn",
y = "Thu nhập (trung vị)",
fill = "Trình độ"
) +
theme_minimal()
library(ggplot2)
# Chuẩn bị dữ liệu: loại bỏ NA, income <= 0, định dạng edu
data$edu <- factor(data$edu, levels = c(1, 2, 3),
labels = c("Primary", "Secondary", "Tertiary"))
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(edu))
ggplot(filtered_data, aes(x = edu, y = income, fill = edu)) +
# Vẽ cột trung vị
stat_summary(fun = median, geom = "bar", color = "black") +
# Thêm nhãn số (giá trị trung vị)
stat_summary(fun = median, geom = "text", aes(label = round(..y.., 0)),
vjust = -0.5, size = 5) +
scale_fill_manual(values = c("Primary" = "blue",
"Secondary" = "red",
"Tertiary" = "orange")) +
labs(
title = "Trung vị thu nhập theo trình độ học vấn",
x = "Trình độ học vấn",
y = "Thu nhập (trung vị)",
fill = "Trình độ"
) +
theme_minimal()
library(ggplot2)
# Dữ liệu gốc có gender = 1 (Nam), 2 (Nữ)
# Chuyển về factor nhưng giữ nguyên giá trị gốc
data$gender <- factor(data$gender, levels = c(1, 2), labels = c("Nam", "Nữ"))
# Lọc dữ liệu hợp lệ: income > 0, age không NA
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(age) & !is.na(gender))
# Vẽ biểu đồ scatter + đường LOESS riêng cho từng giới
ggplot(filtered_data, aes(x = age, y = log(income), color = gender)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = TRUE) +
scale_color_manual(values = c("Nam" = "red", "Nữ" = "blue")) +
labs(
title = "Mối liên quan giữa tuổi và log(Thu nhập) theo giới tính",
x = "Tuổi",
y = "log(Thu nhập)",
color = "Giới tính"
) +
theme_minimal()
# Đảm bảo 'gender' là factor có nhãn rõ ràng
data$gender <- factor(data$gender, levels = c(1, 2), labels = c("Male", "Female"))
# Lọc dữ liệu có income hợp lệ
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(gender))
library(ggplot2)
ggplot(filtered_data, aes(x = gender, y = log(income), fill = gender)) +
geom_boxplot(outlier.color = "black", outlier.shape = 16, outlier.size = 2, alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.5, color = "darkgray") +
scale_fill_manual(values = c("Male" = "red", "Female" = "blue")) +
labs(
title = "Phân bố log(Thu nhập) theo giới tính",
x = "Giới tính",
y = "log(Thu nhập)"
) +
theme_minimal()
# Chuyển gender thành factor có nhãn đúng
data$gender <- factor(data$gender, levels = c(1, 2), labels = c("Male", "Female"))
library(ggplot2)
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(gender))
ggplot(filtered_data, aes(x = gender, y = log(income), fill = gender)) +
geom_boxplot(outlier.color = "black", outlier.shape = 16, outlier.size = 2, alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.5, color = "darkgray") +
scale_fill_manual(values = c("Male" = "red", "Female" = "blue")) +
labs(
title = "Phân bố log(Thu nhập) theo giới tính",
x = "Giới tính",
y = "log(Thu nhập)"
) +
theme_minimal()
str(data$gender)## Factor w/ 2 levels "Male","Female": NA NA NA NA NA NA NA NA NA NA ...table(data$gender)##
## Male Female
## 0 0# Đảm bảo gender là factor với nhãn "Male" và "Female"
data$gender <- factor(data$gender, levels = c(1, 2), labels = c("Male", "Female")) # hoặc c("1", "2") nếu là ký tự
library(ggplot2)
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(gender))
ggplot(filtered_data, aes(x = gender, y = log(income), fill = gender)) +
geom_boxplot(outlier.color = "black", outlier.shape = 16, outlier.size = 2, alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.5, color = "darkgray") +
scale_fill_manual(values = c("1" = "red", "2" = "blue")) +
labs(
title = "Phân bố log(Thu nhập) theo giới tính",
x = "Giới tính",
y = "log(Thu nhập)"
) +
theme_minimal()
# Chuyển gender về factor nếu chưa đúng
data$gender <- factor(data$gender,
levels = c(1, 2),
labels = c("Male", "Female"))
# Lọc dữ liệu hợp lệ
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(age) & !is.na(gender))
library(ggplot2)
ggplot(filtered_data, aes(x = age, y = log(income), color = gender)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = TRUE) +
scale_color_manual(values = c("Male" = "red", "Female" = "blue")) +
labs(
title = "Mối liên quan giữa tuổi và log(Thu nhập) theo giới tính",
x = "Tuổi",
y = "log(Thu nhập)",
color = "Giới tính"
) +
theme_minimal()
library(ggplot2)
# Dữ liệu gốc có gender = 1 (Nam), 2 (Nữ)
# Chuyển về factor nhưng giữ nguyên giá trị gốc
data$gender <- factor(data$gender, levels = c(1, 2), labels = c("Nam", "Nữ"))
# Lọc dữ liệu hợp lệ: income > 0, age không NA
filtered_data <- subset(data, !is.na(income) & income > 0 & !is.na(age) & !is.na(gender))
# Vẽ biểu đồ scatter + đường LOESS riêng cho từng giới
ggplot(filtered_data, aes(x = age, y = log(income), color = gender)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = TRUE) +
scale_color_manual(values = c("NamNam" = "red", "Nữ" = "blue")) +
labs(
title = "Mối liên quan giữa tuổi và log(Thu nhập) theo giới tính",
x = "Tuổi",
y = "log(income)",
color = "Giới tính"
) +
theme_minimal()
ob = read.csv("D:\\PhD 7May2025\\GS Nguyễn Văn Tuấn khóa học phương pháp nghiên cứu khoa học ppnckh\\16Jun2025 Phân tích DL với R\\Obesity data.csv")
summary(lm(pcfat ~ bmi, data=ob))##
## Call:
## lm(formula = pcfat ~ bmi, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.612 -4.181 1.392 4.690 18.241
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.39889 1.36777 6.141 0.00000000111 ***
## bmi 1.03619 0.06051 17.123 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.45 on 1215 degrees of freedom
## Multiple R-squared: 0.1944, Adjusted R-squared: 0.1937
## F-statistic: 293.2 on 1 and 1215 DF, p-value: < 0.00000000000000022library(ggplot2)
p = ggplot(data=ob, aes(x=bmi, y=pcfat, col=gender))
p + geom_point() + geom_smooth(method="lm")## `geom_smooth()` using formula = 'y ~ x'
# slide 266