Thực tập buổi 2
Trần Ngọc Nguyên
8/5/2019
Làm quen với các packages lacttice ggplot2 gridExtra
Đọc dữ liệu “PISA VN 2015.csv” và gọi là đối tượng “students” Đọc dữ liệu “PISA VN SCHOOLS 2015.csv” và gọi là đối tượng “schools” Hợp nhất hai dữ liệu bằng biến CNTSCHID và gọi kết quả là “pisa”
students = read.csv("C:\\Users\\Nguyen\\Desktop\\Workshop NCKH 2019\\Dataset\\PISA VN 2015.csv")
schools = read.csv("C:\\Users\\Nguyen\\Desktop\\Workshop NCKH 2019\\Dataset\\PISA VN SCHOOLS 2015.csv")
pisa = merge(students, schools, by="CNTSCHID")
names(pisa)## [1] "CNTSCHID" "AGE" "Gender" "PARED" "HEDRES" "MISCED"
## [7] "FISCED" "HISCED" "WEALTH" "ESCS" "INSTSCIE" "SCIEEFF"
## [13] "JOYSCIE" "ICTRES" "HOMEPOS" "HEDRES.1" "CULTPOSS" "PV1MATH"
## [19] "PV1READ" "PV1SCIE" "STRATUM" "SCHSIZE" "CLSIZE" "STRATIO"
## [25] "SCHLTYPE" "Region" "Area"#Dán nhãn cho biến số
library(DescTools)
Label(pisa$PV1MATH) = "Math"
Label(pisa$PV1SCIE) = "Science"
Label(pisa$PV1READ) = "Reading"TASK 1
Biểu đồ phân bố với hist()
hist(pisa$PV1SCIE, col="blue") 
hist(pisa$PV1SCIE, col="blue", border="white") 
hist(pisa$PV1SCIE, col="blue", border="white", xlab="Science Score", ylab="Frequency", main="Distribution of Science Score")
Thêm đường biểu diễn mật độ
hist(pisa$PV1SCIE, col="blue", border="white", prob=T)
lines(density(pisa$PV1SCIE), col="red")
Vẽ hai phân bố cho nam và nữ
Bảng mã màu trong R Nếu dùng Windows, phải chuyển mã Java để chú thích biểu đồ bằng tiếng Việt. Hệ điều hành MacOs thì không cần. Chuyển mã Java
p1 = hist(pisa$PV1SCIE[pisa$Gender=="Boys"], plot=F)
p2 = hist(pisa$PV1SCIE[pisa$Gender=="Girls"], plot=F)
plot(p1, col="skyblue", border="white", xlab = '\u0110i\u1EC3m khoa h\u1ECDc', ylab = 'T\u1EA7n su\u1EA5t', main = 'Ph\u00E2n b\u1ED1 \u0111i\u1EC3m khoa h\u1ECDc')
plot(p2, add=T, col=scales::alpha("green", 0.5), border="white")
#Gán legend
legend('right', pch = c(16, 16), col = c('skyblue', 'green'), legend = c('Nam','N\u1EEF'), bty ='n', cex = 1)
TASK 2
Biểu đồ phân bố với densityplot() trong package lattice
library(lattice)
densityplot(~PV1SCIE, groups=Gender, data=pisa)
Mô tả điểm khoa học theo Area
densityplot(~PV1SCIE, groups=Area, data=pisa)
#Gán legend
densityplot(~PV1SCIE, groups=Area, data=pisa, auto.key=list(space="right"))
Chia 3 cửa sổ (dùng package gridExtra)
p1 = densityplot(~PV1SCIE, groups=Area, data=pisa, auto.key=list(space="top"))
p2 = densityplot(~PV1MATH, groups=Area, data=pisa)
p3 = densityplot(~PV1READ, groups=Area, data=pisa)
library(gridExtra)
grid.arrange(p1, p2, p3, ncol=3)
TASK 3
Biểu đồ thanh với sjPlot
Tìm hiểu bao nhiêu học sinh từng vùng và miền
library(sjPlot)
plot_frq(pisa$Area)
plot_frq(pisa$Region) 
Phân tích theo nhóm
sjp.xtab(pisa$Region, pisa$Area, margin="row", bar.pos="stack", show.summary=T)
sjp.xtab(pisa$Region, pisa$Area, margin="row", bar.pos="stack", show.summary=TRUE, coord.flip=T) 
Điểm môn khoa học theo vùng
sjp.grpfrq(pisa$PV1SCIE, pisa$Region, type="boxplot")## 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
Biểu đồ hộp với boxplot()
#Chia màn hình ra 3 cửa sổ
par(mfrow=c(1, 3))
boxplot(pisa$PV1SCIE, col="purple")
boxplot(pisa$PV1SCIE ~ pisa$Gender, col=c("blue", "red"))
boxplot(pisa$PV1SCIE ~ pisa$Region, col=c("blue", "red", "purple"))
TASK 5
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$PV1SCIE ~ pisa$PARED, col="blue")
abline(lm(pisa$PV1SCIE ~ pisa$PARED), col="red")
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
So sánh hai biểu đồ. Biểu đồ nào dễ cảm nhận hơn?!
TASK 6
Biểu đồ likert với dữ liệu “Big Five Personality Data.csv”
big5 = read.csv("C:\\Users\\Nguyen\\Desktop\\Workshop NCKH 2019\\Dataset\\Big Five Personality Data.csv")
test = big5[, c("gender", "E1", "E2", "E3", "E4", "E5")]
test$E1 = as.factor(test$E1)
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)
plot_likert(test)## 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
#Phân tích biểu đồ chất lượng cao (ggplot2) # TASK 7 ## Biểu đồ tương quan Đọc dữ liệu “obesity data” vào R và gọi đối tượng là “ob”. 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 = read.csv("C:\\Users\\Nguyen\\Desktop\\Workshop NCKH 2019\\Dataset\\Obesity data.csv")
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.2Tạo ra biến mới gọi là “OB” dựa vào biến 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"
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")
TASK 8
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## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
p = p + xlab("Weight") + ylab("Percent body fat") + ggtitle("Weight and Percent Body Fat")
p## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
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## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
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'
Có thể thử theme_tufte(), theme_few(), theme_wsj(), theme_clean(), theme_hc()
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'
Có thể thử nghiệm thêm với type (“density”, “violin”, “boxplot”):
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'
Biểu đồ tương quan đa biến (package GGally)
Chọn biến quan tâm
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2dat = 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`.
Sso sánh với:
ggpairs(data=ob, mapping = aes(color = gender), columns = c("age", "weight", "bmi", "pcfat"))
TASK 9
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=PV1SCIE))
p1 = p + geom_histogram(color="white", fill="blue")
p = ggplot(data=pisa, aes(x=PV1SCIE))
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=PV1SCIE, fill=Area))
p1 = p + geom_histogram(position="dodge")
p2 = ggplot(data=pisa, aes(x=PV1SCIE, 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 10
Vẽ biểu đồ hộp (box plot)
Biểu đồ hộp theo vùng
p = ggplot(data=pisa, aes(x=Area, y=PV1SCIE, 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=PV1SCIE, fill=PARED))
p1 = p + geom_boxplot(col="black") + geom_jitter(alpha=0.02)
p = ggplot(data=na.omit(pisa), aes(x=factor(PARED), y=PV1SCIE, 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).