Đọ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).