Mixed Effects Model - Day 2
Nguyen Van Tuan
30/12/2019
Làm quen với những package sau đây: lattice, ggplot2, gridExtra
Task 1: Biểu đồ phân bố với hist()
Đọc dữ liệu PISA Data Vietnam 2015.csv và gọi là đối tượng pisa. Thử nghiệm với hist()
pisa= read.csv("C:/Users/ADMIN/OneDrive/Statistical courses/Ton Duc Thang_Jan2020/TDTU Datasets for 2020 Workshop/PISA Data Vietnam 2015.csv")
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")
Nếu muốn để dấu tiếng Việt trong các labels cần phải dùng mã Java theo đường link. Ví dụ:
hist(pisa$Science, col="blue", border="white", xlab="\u0110i\u1EC3m khoa h\u1ECDc", ylab="T\u1EA7n s\u1ED1", main="Ph\u00E2n b\u1ED1 \u0111i\u1EC3m khoa h\u1ECDc")
Bảng mã màu của R có thể tải về theo link này
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ữ
p1 = hist(pisa$Science[pisa$Gender=="Boys"], plot=F)
p2 = hist(pisa$Science[pisa$Gender=="Girls"], plot=F)
plot(p1, col="skyblue", border="white")
plot(p2, add=T, col=scales::alpha("green", 0.5), border="white")
Bạn có nhận xét gì về biểu đồ?
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|Gender, 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"))
p2 = densityplot(~Math, groups=Area, data=pisa)
p3 = densityplot(~Read, groups=Area, data=pisa)
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))
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 Math và Science
plot(pisa$Science~ pisa$Math, col="blue")
plot(pisa$Science ~ pisa$Math, col="blue")
abline(lm(pisa$Science ~ pisa$Math), 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
ob= read.csv("C:/Users/ADMIN/OneDrive/Statistical courses/Ton Duc Thang_Jan2020/TDTU Datasets for 2020 Workshop/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: 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"
ob$OB = factor(ob$OB, levels=c("Underweight", "Normal", "Overweight", "Obese"))
table (ob$OB)##
## Underweight Normal Overweight Obese
## 107 857 230 15Tì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 + theme(legend.position="none")
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")
p = p + theme(plot.title=element_text(hjust=0.5))
p
- 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()
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)
Có thể thử nghiệm thêm với type ("density", "violin", "boxplot")
- Biểu đồ tương quan đa biến (package
GGally) Chọn biến quan tâm
library(GGally)
dat = ob[, c("gender", "age", "bmi", "weight", "height", "pcfat")]
ggpairs(dat)
Thêm màu theo nhóm
ggpairs(data=dat, mapping = aes(color = gender))
so sánh với:
ggpairs(data=ob, mapping = aes(color = gender), columns = c("age", "weight", "bmi", "pcfat"))
Task 7: Vẽ biểu đồ phân bố (histogram)
Đọc dữ liệu PISA Data Vietnam 2015.csv và gọi là đối tượng pisa
pisa = read.csv ("C:/Users/ADMIN/OneDrive/Statistical courses/Ton Duc Thang_Jan2020/TDTU Datasets for 2020 Workshop/PISA Data Vietnam 2015.csv")- 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)
- 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
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= Math, fill= factor(PARED), color= factor(PARED)))
p2 = p + geom_boxplot(col="black") + geom_jitter(alpha=0.05)
grid.arrange(p1, p2, ncol=2)