Làm quen biểu đồ
Khac Thuan
8/5/2019
Nạp và hợp nhất 2 dữ liệu PISA
#Nạp 2 dữ liệu
students = read.csv('C:/Users/Dell/Desktop/Dataset/PISA VN 2015.csv')
schools = read.csv('C:/Users/Dell/Desktop/Dataset/PISA VN SCHOOLS 2015.csv')
#Merge 2 dữ liệu
pisa = merge(students, schools, by = "CNTSCHID")
#Biến trong dữ liệu khá dài, có thể đổi tên bằng hàm label
#labels(pisa$PV1MATH) = "Math"
#labels(pisa$PV1SCIE) = "Science"
#labels(pisa$PV1READ) = "Reading"
#Summary
summary(pisa)## CNTSCHID AGE Gender PARED
## Min. :70400001 Min. :15.33 Boys :2786 Min. : 3.000
## 1st Qu.:70400052 1st Qu.:15.50 Girls:3040 1st Qu.: 9.000
## Median :70400096 Median :15.75 Median : 9.000
## Mean :70400097 Mean :15.78 Mean : 9.374
## 3rd Qu.:70400143 3rd Qu.:16.00 3rd Qu.:12.000
## Max. :70400188 Max. :16.25 Max. :17.000
## NA's :14
## HEDRES MISCED FISCED HISCED
## Min. :-4.3706 Min. :0.000 Min. :0.000 Min. :0.00
## 1st Qu.:-1.5169 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:2.00
## Median :-1.1119 Median :2.000 Median :2.000 Median :2.00
## Mean :-1.0470 Mean :2.069 Mean :2.296 Mean :2.58
## 3rd Qu.:-0.7026 3rd Qu.:3.000 3rd Qu.:4.000 3rd Qu.:4.00
## Max. : 1.1767 Max. :6.000 Max. :6.000 Max. :6.00
## NA's :20 NA's :34 NA's :88 NA's :14
## WEALTH ESCS INSTSCIE SCIEEFF
## Min. :-7.635 Min. :-5.657 Min. :-1.9301 Min. :-3.7565
## 1st Qu.:-2.829 1st Qu.:-2.539 1st Qu.: 0.0125 1st Qu.:-0.8585
## Median :-2.163 Median :-1.982 Median : 0.3708 Median :-0.3473
## Mean :-2.219 Mean :-1.822 Mean : 0.4835 Mean :-0.2662
## 3rd Qu.:-1.504 3rd Qu.:-1.221 3rd Qu.: 1.0218 3rd Qu.: 0.3155
## Max. : 3.211 Max. : 1.950 Max. : 1.7359 Max. : 3.2775
## NA's :15 NA's :1 NA's :17 NA's :19
## JOYSCIE ICTRES HOMEPOS HEDRES.1
## Min. :-2.1154 Min. :-3.508 Min. :-8.955 Min. :-4.3706
## 1st Qu.: 0.5094 1st Qu.:-2.587 1st Qu.:-2.669 1st Qu.:-1.5169
## Median : 0.5094 Median :-1.855 Median :-2.047 Median :-1.1119
## Mean : 0.6448 Mean :-1.795 Mean :-2.042 Mean :-1.0470
## 3rd Qu.: 1.1049 3rd Qu.:-1.117 3rd Qu.:-1.354 3rd Qu.:-0.7026
## Max. : 2.1635 Max. : 3.497 Max. : 2.770 Max. : 1.1767
## NA's :19 NA's :34 NA's :2 NA's :20
## CULTPOSS PV1MATH PV1READ PV1SCIE
## Min. :-1.8413 Min. :201.7 Min. :107.1 Min. :292.7
## 1st Qu.:-0.8310 1st Qu.:440.0 1st Qu.:442.5 1st Qu.:470.9
## Median :-0.4113 Median :493.4 Median :489.5 Median :523.9
## Mean :-0.4444 Mean :496.1 Mean :489.9 Mean :524.8
## 3rd Qu.: 0.1157 3rd Qu.:551.5 3rd Qu.:537.6 3rd Qu.:574.8
## Max. : 2.1655 Max. :820.1 Max. :744.1 Max. :807.3
## NA's :52
## STRATUM SCHSIZE CLSIZE STRATIO
## VNM0313:989 Min. : 113 Min. :13.00 Min. : 4.314
## VNM0208:884 1st Qu.: 650 1st Qu.:38.00 1st Qu.:14.024
## VNM0101:806 Median :1090 Median :38.00 Median :16.627
## VNM0207:790 Mean :1082 Mean :40.57 Mean :16.497
## VNM0102:764 3rd Qu.:1419 3rd Qu.:43.00 3rd Qu.:18.983
## VNM0314:679 Max. :4016 Max. :53.00 Max. :38.651
## (Other):914 NA's :34
## SCHLTYPE Region Area
## Min. :1.000 CENTRAL:2006 REMOTE: 410
## 1st Qu.:3.000 NORTH :1958 RURAL :2368
## Median :3.000 SOUTH :1862 URBAN :3048
## Mean :2.849
## 3rd Qu.:3.000
## Max. :3.000
## NA's :35Một số thử nghiệm với plot_frq trong sjPlot : vẽ Barchart
#Barchart by region, kèm đổi tên tittle
library(sjPlot) ## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!plot_frq(pisa$Region, title = "Tần số theo vùng miền")#Đây là 1 hàm trong sjPlot(x,..) library dùng để phân tích rồi vẽ biểu đồ về số lượng và tần số của 1 biến số
plot_frq(pisa$Area, axis.title = "Tỉ lệ theo khu vực", axis.labels = c("Vùng xa", "Nông thôn", "Thành thị"))
Phân tích từng thành phần con của biến phân loại
sjp.xtab(pisa$Area, pisa$Region) #sjp.xtab là hàm thuộc sjPlot, ở đây pisa$Area là chính, pisa$region là con, được dùng phân các thành phần trong pisa$Area thành các phần nhỏ hơn
sjp.xtab(pisa$Area, pisa$Region, coord.flip=T)#mặc định sjp.xtab là đặt dọc, nay ta thêm coord.flip = T , True tức sẽ thực hiện flip(sưap) xoay biểu đồ 90 độ -> biểu đồ ngang 
#Biểu đồ đẹp hơn, không còn chồng số, chữ cái
2. Task 2
Vẽ biểu đồ Histogram
#using function hist() to create a histogram (Biểu đồ tần số) with R built-in base function
hist(pisa$PV1SCIE, col = "lightblue")
#Một chút sửa đổi, customize
hist(pisa$PV1SCIE, col = "lightblue", border = "white", main ="Phân bố điểm Khoa học", xlab = "Điểm khoa học", ylab = "Số lượng")
#more graphical
hist(pisa$PV1SCIE, col = "lightblue", border = "white", xlab = "Science Score", ylab = "Frequency", main = "Distribution of Science Score", breaks = 20)
hist(pisa$PV1SCIE, col = "lightblue", border = "white", xlab = "Science Score", ylab = "Frequency", main = "Distribution of Science Score", breaks = 20, prob = 1)
Adding line representatives of density
hist(pisa$PV1SCIE, col = "lightblue", border = "white", prob = T)
lines(density(pisa$PV1SCIE), col = "red")
Biểu đồ phân bố histogram cho từng gender nam và nữ được ghép lại với nhau
#Distribution for both sex
p1 = hist(pisa$PV1SCIE[pisa$Gender=="Boys"], plot =F) #plot = F thì trả về sẽ chỉ là các dữ liệu số, không vẽ ra hình liền khi chạy trên workpace, plot = T thì sẽ vẽ hình liền
p2 = hist(pisa$PV1SCIE[pisa$Gender=="Girls"], plot=F)
plot(p1, col="yellow", border="white")
plot(p2,add = T, col=scales::alpha("red", 0.6), border = 'white') # scales::alpha('color', number) là hàm tự chỉnh màu theo ý thích
#Rất khó nhìn do 2 biểu đồ chồng lên nhau, vùng nào chồng thì 2 màu sắc sẽ pha lại, ở đây pha thành màu cam
Task 3 : Biểu đồ phân bố mật độ với nhiều đối tượng, với package lattice
Nếu ở base package, ta có thể vẽ density dạng đường, nhưng chỉ với dữ liệu một biến số liên tục, còn khi muốn so sánh/vẽ nhiều biến số trên một biểu đồ, ta phải dùng packages khác
library(lattice)
densityplot(data = pisa,~PV1SCIE, groups = Gender)#Gender có 2 phân nhóm
densityplot(~PV1SCIE, data = pisa, groups = Area)#Are có 3 phân nhóm
densityplot(data = pisa, ~PV1SCIE, groups = Area, auto.key = list(space = 'right'))
#auto.key là hàm base, chức năng hiển thị plot của x và y, ở đây chỉ có biến x, và x có 3 phân nhóm nên sẽ hiện chú thích 3 nhóm đó
#list là hàm xây dựng, liệt kê
#space = 'ở đây có thể chọn top, bottom, left, right,..' để xác định khoảng không gian chú thích( legend) xuất hiện tương ứng với biểu đồtable(pisa$Area)##
## REMOTE RURAL URBAN
## 410 2368 3048DIsplay various graph in a plot (use gridExtra package)
p1 = densityplot(~PV1SCIE, groups = Area, data = pisa)#densityplot là hàm trong lattice library
p2 = densityplot(~PV1MATH, groups = Area, data = pisa)
p3 = densityplot(~PV1READ, groups = Area, data = pisa)
library(gridExtra)
grid.arrange(p1, p2, p3, ncol=3, bottom = "Các biểu đồ điểm khoa học, toán và đọc")#Xếp theo cột
grid.arrange(p1, p2, p3, nrow = 3)#Xếp theo hàng
Biểu đồ hộp ( box plot ): biểu diễn phân bố biến liên tục
#Biểu đồ box bằng thư viện sjPlot với hàm sjp.grapfrq(,type ='box')
sjp.grpfrq(pisa$PV1SCIE, pisa$Region, type = "box")#sjp.grapfrq bắt buộc phải có 2 đối số (argument)## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
Task 4: Thử với hàm boxplot mặc định
par(mfrow=c(1, 3))#chia màn hình thành 3 cửa sổ ( 1 dòng, 3 cột )
boxplot(pisa$PV1SCIE, col="purple")
boxplot(pisa$PV1SCIE ~ pisa$Gender, col=c("blue", "red"))
boxplot(pisa$PV1SCIE ~ pisa$Region, col=c("blue", "red", "yellow"))
Task 5: Scatter plot ( biểu đồ tương quan )
plot(pisa$PV1SCIE ~ pisa$PARED, pch = 16, col="blue", xlab = "PARED", ylab = "Điểm khoa học", main = 'Mỗi liên quan giữa Pared và điểm khoa học')
#Hàm plot mặc định, dùng vẽ scatter plot, xlab là tiêu đề trục hoành, ylab là tiêu đề trục tung, main là tiêu đề biểu đồ
abline(lm(pisa$PV1SCIE ~ pisa$PARED), col="red", lwd = 2) 
#Hàm vẽ hồi quy tuyến tính mặc định, nhưng không nhiều thông tin và đẹp
#lwd là số chỉ độ dày của đường tuyến tính#Dùng boxplot trong sjPlot
library(sjPlot)
sjp.grpfrq(pisa$PV1SCIE, pisa$PARED, type="box")## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
Hai cách thể hiện biểu đồ khác nhau trên cùng đối tượng, biểu đồ nào trực quan hơn.Task 6 Biểu đồ Likert
big5 = read.csv('C:/Users/Dell/Desktop/Dataset/Big Five Personality Data.csv')
dim(big5)## [1] 19719 57head(big5)## race age engnat gender hand source country E1 E2 E3 E4 E5 E6 E7 E8 E9
## 1 3 53 1 1 1 1 US 4 2 5 2 5 1 4 3 5
## 2 13 46 1 2 1 1 US 2 2 3 3 3 3 1 5 1
## 3 1 14 2 2 1 1 PK 5 1 1 4 5 1 1 5 5
## 4 3 19 2 2 1 1 RO 2 5 2 4 3 4 3 4 4
## 5 11 25 2 2 1 2 US 3 1 3 3 3 1 3 1 3
## 6 13 31 1 2 1 2 US 1 5 2 4 1 3 2 4 1
## E10 N1 N2 N3 N4 N5 N6 N7 N8 N9 N10 A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 C1 C2
## 1 1 1 5 2 5 1 1 1 1 1 1 1 5 1 5 2 3 1 5 4 5 4 1
## 2 5 2 3 4 2 3 4 3 2 2 4 1 3 3 4 4 4 2 3 4 3 4 1
## 3 1 5 1 5 5 5 5 5 5 5 5 5 1 5 5 1 5 1 5 5 5 4 1
## 4 5 5 4 4 2 4 5 5 5 4 5 2 5 4 4 3 5 3 4 4 3 3 3
## 5 5 3 3 3 4 3 3 3 3 3 4 5 5 3 5 1 5 1 5 5 5 3 1
## 6 5 1 5 4 5 1 4 4 1 5 2 2 2 3 4 3 4 3 5 5 3 2 5
## C3 C4 C5 C6 C7 C8 C9 C10 O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
## 1 5 1 5 1 4 1 4 5 4 1 3 1 5 1 4 2 5 5
## 2 3 2 3 1 5 1 4 4 3 3 3 3 2 3 3 1 3 2
## 3 5 1 5 1 5 1 5 5 4 5 5 1 5 1 5 5 5 5
## 4 4 5 1 4 5 4 2 3 4 3 5 2 4 2 5 2 5 5
## 5 5 3 3 1 1 3 3 3 3 1 1 1 3 1 3 1 5 3
## 6 4 3 3 4 5 3 5 3 4 2 1 3 3 5 5 4 5 3#Ta thấy các biến E1, E2…E10, N1…10 dữ liệu là các con số (numreic)
test = big5[, c("gender", "E1", "E2", "E3", "E4", "E5")]
test$E1 = as.factor(test$E1) #Chuyển numeric sang character, từ đó có thể đếm số lượng 'số 1,2..'
test$E2 = as.factor(test$E2)
test$E3 = as.factor(test$E3)
test$E4 = as.factor(test$E4)
test$E5 = as.factor(test$E5)
library(sjPlot) #Gọi thư viện sjPlot
plot_likert(test) #sjp.likert(test) là code cũ nên không đúng## Warning: Detected uneven category count in items. Dropping last category.## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
## Warning in freq[valid] <- counts: number of items to replace is not a
## multiple of replacement length
Biểu đồ chất lượng cao ( dùng thư viện ggplot2 )
Task 7: Biểu đồ tương quan
ob = read.csv('C:/Users/Dell/Desktop/Dataset/Obesity data.csv')
#Tạo ra biến mới gọi là "OB" dựa vào biến bmi: nếu bmi < 18.5 thì OB = "Underweight"; bmi từ 18.5 đến 24.9 thì OB="Normal"; bmi từ 25.0 đến 29.9 thì OB="Overweight"; bmi từ 30.0 trở lên, OB="Obese"
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"
#Nếu bên công nghệ thông tin thay vì tốn công tận 4 dòng như trên ta dùng hàm ifelse - là hàm base
#ob$OB = ifelse(ob$bmi < 18.5, 'Underweight',ifelse(ob$bmi >= 18.5 & ob$bmi < 24.9,"Normal",ifelse(ob$bmi >= 25.0 & ob$bmi < 29.9, "Overweight", "Obese")))
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 + xlab("Obesity group") + ylab("Frequency")
p = p + ggtitle("Số lượng của từng nhóm cân nặng")
p + theme(legend.position = "right") #legend.position = 'top, bottom, right, left, none'
Task 8: Vẽ biểu đồ tương quan giữa weight và pcfat dùng ggplot2
#Biểu đồ đơn giản
p = ggplot(data=ob, aes(x=weight, y=pcfat))
p = p +geom_point() +geom_smooth()
p## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Ta thấy chưa đẹp lắm, nên chỉnh lại chút về các trục, tên biểu đồ
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat")
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 đồ trên gộp cả nam và nữ ( tức che dấu ), biết đâu giữa nam và nữ lại có khác biệt
#Tách theo từng giới: nam và nữ
p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
p## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Đôi lúc, ta thích một themes nào đó ( ví dụ themes/chủ đề của tờ new york time, washington post, economist,..) có một package giúp ta thực hiện mong muốn
library(ggthemes)
p + theme_economist()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#Vẽ theo style Economist#Vẽ theo các theme khác
p + theme_wsj()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
p + theme_tufte()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
p + theme_clean()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
p + theme_hc()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Tương quan đa chiều
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 Normalvars = ob[,c("gender", "age", "weight", "bmi", "pcfat")]
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(data=vars, 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`.
#Biểu đồ tương quan + biểu đồ phân bố dùng package ggExtra
p = ggplot(data=ob, aes(x=weight, y=pcfat, fill=gender, col=gender))
p = p + geom_point() + geom_smooth()
library(ggExtra)
p = 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
students = read.csv("C:/Users/Dell/Desktop/Dataset/PISA VN 2015.csv")
schools = read.csv("C:/Users/Dell/Desktop/Dataset/PISA VN SCHOOLS 2015.csv")
pisa = merge(students,schools, by="CNTSCHID")
head(pisa,3)## CNTSCHID AGE Gender PARED HEDRES MISCED FISCED HISCED WEALTH ESCS
## 1 70400001 15.58 Boys 9 -1.0418 1 2 2 -2.0697 -1.7899
## 2 70400001 15.92 Boys 12 -2.3041 1 4 4 -1.7903 -1.5423
## 3 70400001 15.42 Girls 9 -0.7218 1 2 2 -2.1942 -2.0475
## INSTSCIE SCIEEFF JOYSCIE ICTRES HOMEPOS HEDRES.1 CULTPOSS PV1MATH
## 1 0.9798 0.1421 2.1635 -1.5244 -2.0537 -1.0418 -0.7273 439.923
## 2 1.7359 -0.8432 2.1635 -1.9305 -2.2627 -2.3041 -0.2031 406.251
## 3 -0.2063 -0.1824 -0.1808 -1.6093 -1.9675 -0.7218 -0.2220 414.369
## PV1READ PV1SCIE STRATUM SCHSIZE CLSIZE STRATIO SCHLTYPE Region Area
## 1 412.290 475.612 VNM0313 883 18 22.075 3 SOUTH URBAN
## 2 409.598 450.320 VNM0313 883 18 22.075 3 SOUTH URBAN
## 3 384.307 405.787 VNM0313 883 18 22.075 3 SOUTH URBANp= ggplot(data=pisa, aes(x=Region, y=PV1SCIE, fill=Region))
p= p+geom_boxplot()
p
#Điểm có phụ thuộc thành phần kinh tế, trình độ học vấn cha mẹ
p= ggplot(data=pisa, aes(x=factor(PARED), y=PV1SCIE, fill=factor(PARED)))
p = p + geom_boxplot() + geom_jitter(alpha=0.05)
p
#WEALTH có ảnh hưởng?
p= ggplot(data=pisa, aes(x=WEALTH, y=PV1SCIE, fill=Region))
p = p + geom_point() + geom_smooth()
p## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Removed 15 rows containing non-finite values (stat_smooth).## Warning: Removed 15 rows containing missing values (geom_point).
p= ggplot(data=pisa, aes(x=WEALTH, y=PV1SCIE, col=Area, fill=Area))
p = p + geom_point() + geom_smooth()
p## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Removed 15 rows containing non-finite values (stat_smooth).## Warning: Removed 15 rows containing missing values (geom_point).