Đọc dữ liệu
t = "F:\\NCKH - PHAN TICH THONG KE\\BAI GIANG CAC KHOA PT DL\\Kho hoc - Du lieu da tang SG1.2020\\Data for practice\\PISA Data Vietnam 2015.csv"
pisa = read.csv(t, header = T)
head(pisa)
## School SchoolSize ClassSize STratio SchoolType Area Region Age
## 1 70400001 883 18 22.075 3 URBAN SOUTH 15.58
## 2 70400001 883 18 22.075 3 URBAN SOUTH 15.92
## 3 70400001 883 18 22.075 3 URBAN SOUTH 15.42
## 4 70400001 883 18 22.075 3 URBAN SOUTH 15.58
## 5 70400001 883 18 22.075 3 URBAN SOUTH 15.92
## 6 70400001 883 18 22.075 3 URBAN SOUTH 16.25
## Gender PARED HISCED WEALTH INSTSCIE JOYSCIE ICTRES Math Read
## 1 Boys 9 2 -2.0697 0.9798 2.1635 -1.5244 439.923 412.290
## 2 Boys 12 4 -1.7903 1.7359 2.1635 -1.9305 406.251 409.598
## 3 Girls 9 2 -2.1942 -0.2063 -0.1808 -1.6093 414.369 384.307
## 4 Girls 5 1 -2.0301 -0.3115 -0.4318 -1.6250 468.801 459.104
## 5 Girls 9 2 -1.0522 0.7648 1.3031 -0.5305 355.432 402.435
## 6 Girls 5 1 -3.0570 0.3708 0.5094 -2.5873 458.955 483.885
## Science
## 1 475.612
## 2 450.320
## 3 405.787
## 4 462.968
## 5 453.736
## 6 529.866
# Xem kich thuoc bang du lieu (so hang, cot)
dim(pisa)
## [1] 5826 18
# Tom tat du lieu
summary(pisa)
## School SchoolSize ClassSize STratio
## Min. :70400001 Min. : 113 Min. :13.00 Min. : 4.314
## 1st Qu.:70400052 1st Qu.: 650 1st Qu.:38.00 1st Qu.:14.024
## Median :70400096 Median :1090 Median :38.00 Median :16.627
## Mean :70400097 Mean :1082 Mean :40.57 Mean :16.497
## 3rd Qu.:70400143 3rd Qu.:1419 3rd Qu.:43.00 3rd Qu.:18.983
## Max. :70400188 Max. :4016 Max. :53.00 Max. :38.651
## NA's :34
## SchoolType Area Region Age Gender
## Min. :1.000 REMOTE: 410 CENTRAL:2006 Min. :15.33 Boys :2786
## 1st Qu.:3.000 RURAL :2368 NORTH :1958 1st Qu.:15.50 Girls:3040
## Median :3.000 URBAN :3048 SOUTH :1862 Median :15.75
## Mean :2.849 Mean :15.78
## 3rd Qu.:3.000 3rd Qu.:16.00
## Max. :3.000 Max. :16.25
## NA's :35
## PARED HISCED WEALTH INSTSCIE
## Min. : 3.000 Min. :0.00 Min. :-7.635 Min. :-1.9301
## 1st Qu.: 9.000 1st Qu.:2.00 1st Qu.:-2.829 1st Qu.: 0.0125
## Median : 9.000 Median :2.00 Median :-2.163 Median : 0.3708
## Mean : 9.374 Mean :2.58 Mean :-2.219 Mean : 0.4835
## 3rd Qu.:12.000 3rd Qu.:4.00 3rd Qu.:-1.504 3rd Qu.: 1.0218
## Max. :17.000 Max. :6.00 Max. : 3.211 Max. : 1.7359
## NA's :14 NA's :14 NA's :15 NA's :17
## JOYSCIE ICTRES Math Read
## Min. :-2.1154 Min. :-3.508 Min. :201.7 Min. :107.1
## 1st Qu.: 0.5094 1st Qu.:-2.587 1st Qu.:440.0 1st Qu.:442.5
## Median : 0.5094 Median :-1.855 Median :493.4 Median :489.5
## Mean : 0.6448 Mean :-1.795 Mean :496.1 Mean :489.9
## 3rd Qu.: 1.1049 3rd Qu.:-1.117 3rd Qu.:551.5 3rd Qu.:537.6
## Max. : 2.1635 Max. : 3.497 Max. :820.1 Max. :744.1
## NA's :19 NA's :34
## Science
## Min. :292.7
## 1st Qu.:470.9
## Median :523.9
## Mean :524.8
## 3rd Qu.:574.8
## Max. :807.3
##
str(pisa)
## 'data.frame': 5826 obs. of 18 variables:
## $ School : int 70400001 70400001 70400001 70400001 70400001 70400001 70400001 70400001 70400001 70400001 ...
## $ SchoolSize: int 883 883 883 883 883 883 883 883 883 883 ...
## $ ClassSize : int 18 18 18 18 18 18 18 18 18 18 ...
## $ STratio : num 22.1 22.1 22.1 22.1 22.1 ...
## $ SchoolType: int 3 3 3 3 3 3 3 3 3 3 ...
## $ Area : Factor w/ 3 levels "REMOTE","RURAL",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Region : Factor w/ 3 levels "CENTRAL","NORTH",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Age : num 15.6 15.9 15.4 15.6 15.9 ...
## $ Gender : Factor w/ 2 levels "Boys","Girls": 1 1 2 2 2 2 2 2 2 1 ...
## $ PARED : int 9 12 9 5 9 5 3 5 9 5 ...
## $ HISCED : int 2 4 2 1 2 1 0 1 2 1 ...
## $ WEALTH : num -2.07 -1.79 -2.19 -2.03 -1.05 ...
## $ INSTSCIE : num 0.98 1.736 -0.206 -0.311 0.765 ...
## $ JOYSCIE : num 2.163 2.163 -0.181 -0.432 1.303 ...
## $ ICTRES : num -1.52 -1.93 -1.61 -1.62 -0.53 ...
## $ Math : num 440 406 414 469 355 ...
## $ Read : num 412 410 384 459 402 ...
## $ Science : num 476 450 406 463 454 ...
Task 1: Biểu đồ phân bố với hist()
#Thử nghiệm với hist()
hist(pisa$Math, col="blue")

hist(pisa$Math, col="blue", border="white")

hist(pisa$Science, col="blue", border="white", xlab="Science Score", ylab="Frequency", main="Distribution of Science Scores")

#Thêm đường biểu diễn
hist(pisa$Science, col="blue", border="white", prob=T)
lines(density(pisa$Science), col="red")

#Vẽ 2 phân bố cho nam và nữ
##Gán bieur đồ cho p1 và p2 nhưng không vẻ (plot=F)
p1 = hist(pisa$Science[pisa$Gender=="Boys"], plot=F)
p2 = hist(pisa$Science[pisa$Gender=="Girls"], plot=F)
plot(p1, col="skyblue", border="white")
##Add thêm biểu đồ p2 với option add=T
plot(p2, add=T, col=scales::alpha("green", 0.5), border="white")

Task 2: Biểu đồ phân bố với lattice()
#Vẽ biểu đồ phân bố cho nhiều nhóm dùng hàm densityplot
library(lattice)
densityplot(~Science, groups=Gender, data=pisa)

densityplot(~Science, groups=Area, data=pisa)

densityplot(~Science, groups=Area, data=pisa, auto.key=list(space="top"))

#Chia 3 cửa sổ (dùng package lattice và gridExtra)
p1 = densityplot(~Science, groups=Area, data=pisa, auto.key=list(space="top"), ylim = c(0.000, 0.007))
p1

p2 = densityplot(~Math, groups=Area, data=pisa, auto.key=list(space="top"), ylim = c(0.000, 0.007))
p3 = densityplot(~Read, groups=Area, data=pisa, auto.key=list(space="top"), ylim = c(0.000, 0.007))
library(gridExtra)
grid.arrange(p1, p2, p3, ncol=3)

Task 3: Biểu đồ hộp với boxplot()
# chia màn hình ra 3 cửa sổ
par(mfrow=c(1, 3))
par(mfrow=c(1, 1))
boxplot(pisa$Science, col="purple")

boxplot(pisa$Science ~ pisa$Gender, col=c("blue", "red"))

boxplot(pisa$Science ~ pisa$Region, col=c("blue", "red", "purple"))

Task 4: Biểu đồ tương quan với plot()
#Tìm hiểu mối liên quan giữa PARED và điểm môn khoa học
plot(pisa$Science ~ pisa$PARED, col="blue")
abline(lm(pisa$Science ~ pisa$PARED), col="red")

#Tìm hiểu mối liên quan giữa điểm môn toán và điểm môn khoa học
plot(pisa$Math ~ pisa$Science, col="blue")
abline(lm(pisa$Math ~ pisa$Science), col="red")

Phân tích biểu đồ chất lượng cao (ggplot2)
Task 5: Biểu đồ tương quan
#Đọc dữ liệu "obesity data" vào R và gọi đối tượng là "ob"
t1 = "F:\\NCKH - PHAN TICH THONG KE\\BAI GIANG CAC KHOA PT DL\\Kho hoc - Du lieu da tang SG1.2020\\Data for practice\\Obesity data.csv"
ob = read.csv(t1, header = T)
head(ob)
## id gender height weight bmi age bmc bmd fat lean pcfat
## 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3
## 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8
## 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0
## 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8
## 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8
## 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2
#Tạo biến mới OB theo bmi
ob$OB [ob$bmi < 18.5] = "Underweight"
ob$OB [ob$bmi >= 18.5 & ob$bmi <= 24.9] = "Normal"
ob$OB [ob$bmi >= 25.0 & ob$bmi <= 29.9] = "Overweight"
ob$OB [ob$bmi >= 30] = "Obese"
#Đặt lại thứ tự nhãn của biến OB với hàm factor. Nếu không R sẽ mặc định theo thứ tự a,b,c
ob$OB = factor(ob$OB, levels=c("Underweight", "Normal", "Overweight", "Obese"))
#Tìm hiểu phân bố của OB
library(ggplot2)
p = ggplot(data=ob, aes(OB, fill=OB)) + geom_bar()
p

p = p + xlab("Obesity group") + ylab("Frequency")
p + theme(legend.position="none")

table(ob$OB)
##
## Underweight Normal Overweight Obese
## 107 865 230 15
Task 6: Vẽ biểu đồ tương quan giữa weight và pcfat dùng ggplot2
#Biểu đồ đơn giản
library(ggplot2)
library(ggthemes)
library(gridExtra)
p = ggplot(data=ob, aes(x=weight, y=pcfat))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat")
##Đổi tên biểu đồ sang tiếng việt
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Bi\u1EC3u \u0111\u1ED3 t\u01B0\u01A1ng quan gi\u1EEFa weight v\u00E0 pcfat")
p = p + theme(plot.title=element_text(hjust=0.5))
p
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

#Biểu đồ theo nhóm nam nữ
p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_economist()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Thử biểu đồ với các kiểu hình nền khác nhau gồm: theme_tufte(), theme_few(), theme_wsj(), theme_clean(), theme_hc()
p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_tufte()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_few()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_wsj()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_clean()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat for men and women separately")
p = p + theme(plot.title=element_text(hjust=0.5))
p + theme_hc()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#Biểu đồ tương quan + biểu đồ phân bố dùng package ggExtra
library(ggExtra)
p = ggplot(data=ob, aes(x=bmi, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
ggMarginal(p, type="histogram", groupColour=T, groupFill=T)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
##Thử nghiệm kết hợp với các loại biểu đồ khác với type ("density", "violin", "boxplot")
p = ggplot(data=ob, aes(x=bmi, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
ggMarginal(p, type="density", groupColour=T, groupFill=T)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
p = ggplot(data=ob, aes(x=bmi, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
ggMarginal(p, type="violin", groupColour=T, groupFill=T)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals
p = ggplot(data=ob, aes(x=bmi, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
ggMarginal(p, type="boxplot", groupColour=T, groupFill=T)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#Biểu đồ tương quan đa biến (package GGally)
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
head(ob)
## id gender height weight bmi age bmc bmd fat lean pcfat OB
## 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3 Normal
## 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8 Normal
## 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0 Normal
## 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8 Normal
## 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8 Normal
## 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2 Normal
dat = ob[, c("gender", "age", "bmi", "weight", "height", "pcfat")]
ggpairs(dat)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##Thêm màu theo nhóm
ggpairs(data=dat, mapping = aes(color = gender))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

##so sánh với
ggpairs(data=ob, mapping = aes(color = gender), columns = c("gender", "weight", "bmi", "pcfat"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Task 7: Vẽ biểu đồ phân bố (histogram)
#Vẽ biểu đồ phân bố điểm môn khoa học
library(gridExtra)
p = ggplot(data=pisa, aes(x=Science))
p1 = p + geom_histogram(color="white", fill="blue")
p = ggplot(data=pisa, aes(x=Science))
p = p + geom_histogram(aes(y=..density..), color="white", fill="blue")
p2 = p + geom_density(col="red")
grid.arrange(p1, p2, ncol=2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Phân bố điểm môn khoa học theo vùng (Area)
p = ggplot(data=pisa, aes(x=Science, fill=Area))
p1 = p + geom_histogram(position="dodge")
p2 = ggplot(data=pisa, aes(x=Science, fill=Area, color=Area)) + geom_density(alpha = 0.1)
grid.arrange(p1, p2, nrow=2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Histogram và xác suất tích lũy
p = ggplot(data=pisa, aes(HISCED))
p = p + stat_ecdf(color="red", lwd=1)
p = p + geom_bar(aes(y = (..count..)/sum(..count..)), fill="blue", colour='yellow')
p
## Warning: Removed 14 rows containing non-finite values (stat_ecdf).
## Warning: Removed 14 rows containing non-finite values (stat_count).

Task 8:Vẽ biểu đồ hộp (box plot)
#Biểu đồ hộp theo vùng
p = ggplot(data=pisa, aes(x=Area, y=Science, col=Area, fill=Area))
p1 = p + geom_boxplot(col="black")
p2 = p + geom_boxplot(col="black") + geom_jitter(alpha=0.05)
grid.arrange(p1, p2, ncol=2)

#Biểu đồ hộp theo kinh tế
p = ggplot(data=pisa, aes(x=PARED, y=Science, fill=PARED))
p1 = p + geom_boxplot(col="black") + geom_jitter(alpha=0.02)
p = ggplot(data=na.omit(pisa), aes(x=factor(PARED), y=Science, fill=factor(PARED), color=factor(PARED)))
p2 = p + geom_boxplot(col="black") + geom_jitter(alpha=0.05)
grid.arrange(p1, p2, ncol=2)
## Warning: Continuous x aesthetic -- did you forget aes(group=...)?
## Warning: Removed 14 rows containing missing values (stat_boxplot).
## Warning: Removed 14 rows containing missing values (geom_point).
