Nhap du lieu

dirty.data <- read.csv("C:/Users/LENOVO/Desktop/BTL XSTK/Dirty_data/water_potability.csv")
View(dirty.data)

Doc du lieu

str(dirty.data)

## 'data.frame':    3276 obs. of  10 variables:
##  $ ph             : num  NA 3.72 8.1 8.32 9.09 ...
##  $ Hardness       : num  205 129 224 214 181 ...
##  $ Solids         : num  20791 18630 19910 22018 17979 ...
##  $ Chloramines    : num  7.3 6.64 9.28 8.06 6.55 ...
##  $ Sulfate        : num  369 NA NA 357 310 ...
##  $ Conductivity   : num  564 593 419 363 398 ...
##  $ Organic_carbon : num  10.4 15.2 16.9 18.4 11.6 ...
##  $ Trihalomethanes: num  87 56.3 66.4 100.3 32 ...
##  $ Turbidity      : num  2.96 4.5 3.06 4.63 4.08 ...
##  $ Potability     : int  0 0 0 0 0 0 0 0 0 0 ...

Dat lai ten bien

colnames(dirty.data) <- c("ph","har","TDS","chl","sul","EC","TOC","THMs","NTU","P")

Doc va kiem tra coi bien da duoc doi ten chua

str(dirty.data)

## 'data.frame':    3276 obs. of  10 variables:
##  $ ph  : num  NA 3.72 8.1 8.32 9.09 ...
##  $ har : num  205 129 224 214 181 ...
##  $ TDS : num  20791 18630 19910 22018 17979 ...
##  $ chl : num  7.3 6.64 9.28 8.06 6.55 ...
##  $ sul : num  369 NA NA 357 310 ...
##  $ EC  : num  564 593 419 363 398 ...
##  $ TOC : num  10.4 15.2 16.9 18.4 11.6 ...
##  $ THMs: num  87 56.3 66.4 100.3 32 ...
##  $ NTU : num  2.96 4.5 3.06 4.63 4.08 ...
##  $ P   : int  0 0 0 0 0 0 0 0 0 0 ...

Kiem tra so luong NA tung cot

apply(is.na(dirty.data),2,sum)

##   ph  har  TDS  chl  sul   EC  TOC THMs  NTU    P 
##  491    0    0    0  781    0    0  162    0    0

Xu ly du lieu khuyet

listMissingColumns <- colnames(dirty.data)[apply(dirty.data, 2, anyNA)]
medianMissing <- apply(dirty.data[,colnames(dirty.data) %in% listMissingColumns],
                2, median, na.rm = TRUE)
medianMissing

##         ph        sul       THMs 
##   7.036752 333.073546  66.622485

Xu ly NA bang median

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

clean.data <- dirty.data %>% mutate(    
              ph = ifelse(is.na(ph), medianMissing[1], ph),
              sul = ifelse(is.na(sul), medianMissing[2], sul),
              THMs = ifelse(is.na(THMs), medianMissing[3], THMs))

Kiem tra lai xem da xu ly dc NA chua

apply(is.na(clean.data),2,sum)

##   ph  har  TDS  chl  sul   EC  TOC THMs  NTU    P 
##    0    0    0    0    0    0    0    0    0    0

In clean.data

View(clean.data)

Dat bien phan loai va tom tat du lieu

clean.data$P <- as.factor(clean.data$P)
summary(clean.data)

##        ph              har              TDS               chl        
##  Min.   : 0.000   Min.   : 47.43   Min.   :  320.9   Min.   : 0.352  
##  1st Qu.: 6.278   1st Qu.:176.85   1st Qu.:15666.7   1st Qu.: 6.127  
##  Median : 7.037   Median :196.97   Median :20927.8   Median : 7.130  
##  Mean   : 7.074   Mean   :196.37   Mean   :22014.1   Mean   : 7.122  
##  3rd Qu.: 7.870   3rd Qu.:216.67   3rd Qu.:27332.8   3rd Qu.: 8.115  
##  Max.   :14.000   Max.   :323.12   Max.   :61227.2   Max.   :13.127  
##       sul              EC             TOC             THMs        
##  Min.   :129.0   Min.   :181.5   Min.   : 2.20   Min.   :  0.738  
##  1st Qu.:317.1   1st Qu.:365.7   1st Qu.:12.07   1st Qu.: 56.648  
##  Median :333.1   Median :421.9   Median :14.22   Median : 66.622  
##  Mean   :333.6   Mean   :426.2   Mean   :14.28   Mean   : 66.407  
##  3rd Qu.:350.4   3rd Qu.:481.8   3rd Qu.:16.56   3rd Qu.: 76.667  
##  Max.   :481.0   Max.   :753.3   Max.   :28.30   Max.   :124.000  
##       NTU        P       
##  Min.   :1.450   0:1998  
##  1st Qu.:3.440   1:1278  
##  Median :3.955           
##  Mean   :3.967           
##  3rd Qu.:4.500           
##  Max.   :6.739

Ve do thi thong ke Potability

bar.P <- barplot(table(clean.data$P),main = "Thống kê số mẫu theo khả năng uống được",xlab = "Khả năng", ylab = "Số mẫu", col = c(1,2), ylim = c(0,2500))

Ve Histogram phan phoi cac bien

hist(clean.data$ph, main = "Histogram của pH", xlab = "Độ pH", ylab = "Tần số", labels = T, col = 5, xlim = c(min(clean.data$ph),max(clean.data$ph)), ylim = c(0,1300))

hist(clean.data$har, main = "Histogram của độ cứng", xlab = "Độ cứng, mg/L", ylab = "Tần số", labels = T, col = 5, xlim = c(min(clean.data$har),max(clean.data$har)), ylim = c(0,900))

hist(clean.data$TDS, main = "Histogram của Solids (TDS)", xlab = "Solids, ppm", ylab = "Tần số", labels = T, col = 5, ylim = c(0,800))

hist(clean.data$chl, main = "Histogram của Chloramines", xlab = "Chloramines, ppm", ylab = "Tần số", labels = T, col = 5, xlim = c(min(clean.data$chl),max(clean.data$chl)),ylim = c(0,900))

hist(clean.data$sul, main = "Histogram của Sulfate", xlab = "Sulfate, mg/L", ylab = "Tần số", labels = T, col = 5,xlim = c(min(clean.data$sul),max(clean.data$sul)), ylim = c(0,1400))

hist(clean.data$EC, main = "Histogram của Tính dẫn điện", xlab = "Tính dẫn điện, micro S/cm", ylab = "Tần số", labels = T, col = 5,xlim = c(min(clean.data$EC),max(clean.data$EC)), ylim = c(0,800))

hist(clean.data$TOC, main = "Histogram của TOC", xlab = "TOC, ppm", ylab = "Tần số", labels = T, col = 5,xlim = c(min(clean.data$TOC),max(clean.data$TOC)), ylim = c(0,800))

hist(clean.data$THMs, main = "Histogram của THMs", xlab = "THMs, micro g/L", ylab = "Tần số", labels = T, col = 5,xlim = c(min(clean.data$THMs),max(clean.data$THMs)), ylim = c(0,1000))

hist(clean.data$NTU, main = "Histogram của Độ đục", xlab = "Độ đục, NTU", ylab = "Tần số", labels = T, col = 5,xlim = c(min(clean.data$NTU),max(clean.data$NTU)), ylim = c(0,900))

Ve do thi boxplot giua cac bien lien tuc va P

boxplot(clean.data$ph ~ clean.data$P, main = "Boxplot khả năng uống được theo pH", xlab = "Khả năng", ylab = "Độ pH", col = c(3,4))

boxplot(clean.data$har ~ clean.data$P, main = "Boxplot khả năng uống được theo độ cứng", xlab = "Khả năng", ylab = "Độ cứng, mg/L", col = c(3,4))

boxplot(clean.data$TDS ~ clean.data$P, main = "Boxplot khả năng uống được theo TDS", xlab = "Khả năng", ylab = "TDS, ppm", col = c(3,4))

boxplot(clean.data$chl ~ clean.data$P, main = "Boxplot khả năng uống được theo Chloramines", xlab = "Khả năng", ylab = "Chloramines, ppm", col = c(3,4))

boxplot(clean.data$sul ~ clean.data$P, main = "Boxplot khả năng uống được theo Sulfate", xlab = "Khả năng", ylab = "Sulfate, mg/L", col = c(3,4))

boxplot(clean.data$EC ~ clean.data$P, main = "Boxplot khả năng uống được theo độ dẫn điện", xlab = "Khả năng", ylab = "Độ dẫn điện, micro S/cm", col = c(3,4))

boxplot(clean.data$TOC ~ clean.data$P, main = "Boxplot khả năng uống được theo TOC", xlab = "Khả năng", ylab = "TOC, ppm", col = c(3,4))

boxplot(clean.data$THMs ~ clean.data$P, main = "Boxplot khả năng uống được theo THMs", xlab = "Khả năng", ylab = "THMs, micro g/L", col = c(3,4))

boxplot(clean.data$NTU ~ clean.data$P, main = "Boxplot khả năng uống được theo Độ đục", xlab = "Khả năng", ylab = "Độ đục, NTU", col = c(3,4))

Ma tran va do thi he so tuong quan

p > 0 thì x tăng, y tăng p < 0 thì x tăng, y giảm (và ngược lại) p = 0 thì không tương quan>

round(cor(clean.data[,1:9]),3)

##          ph    har    TDS    chl    sul     EC    TOC   THMs    NTU
## ph    1.000  0.076 -0.082 -0.032  0.014  0.017  0.040  0.003 -0.036
## har   0.076  1.000 -0.047 -0.030 -0.093 -0.024  0.004 -0.013 -0.014
## TDS  -0.082 -0.047  1.000 -0.070 -0.150  0.014  0.010 -0.009  0.020
## chl  -0.032 -0.030 -0.070  1.000  0.024 -0.020 -0.013  0.017  0.002
## sul   0.014 -0.093 -0.150  0.024  1.000 -0.014  0.027 -0.026 -0.010
## EC    0.017 -0.024  0.014 -0.020 -0.014  1.000  0.021  0.001  0.006
## TOC   0.040  0.004  0.010 -0.013  0.027  0.021  1.000 -0.013 -0.027
## THMs  0.003 -0.013 -0.009  0.017 -0.026  0.001 -0.013  1.000 -0.021
## NTU  -0.036 -0.014  0.020  0.002 -0.010  0.006 -0.027 -0.021  1.000

library(corrplot)

## corrplot 0.92 loaded

corrplot(cor(clean.data[,1:9]), method = "number")

Logistic

Tach du lieu (de xay dung va kiem dinh)

library(caTools)
set.seed(6)
samp = sample.split(clean.data$P, SplitRatio = 0.65)
xd = subset(clean.data, samp == TRUE)
kd = subset(clean.data, samp == FALSE)

Xay dung mo hinh

build <- glm(P ~., data = xd, family = binomial)
mh <- step(build)

## Start:  AIC=2855.73
## P ~ ph + har + TDS + chl + sul + EC + TOC + THMs + NTU
## 
##        Df Deviance    AIC
## - ph    1   2835.7 2853.7
## - NTU   1   2835.9 2853.9
## - sul   1   2835.9 2853.9
## - THMs  1   2836.0 2854.0
## - har   1   2836.3 2854.3
## - EC    1   2836.6 2854.6
## <none>      2835.7 2855.7
## - TOC   1   2839.1 2857.1
## - chl   1   2839.1 2857.1
## - TDS   1   2840.6 2858.6
## 
## Step:  AIC=2853.73
## P ~ har + TDS + chl + sul + EC + TOC + THMs + NTU
## 
##        Df Deviance    AIC
## - NTU   1   2835.9 2851.9
## - sul   1   2835.9 2851.9
## - THMs  1   2836.0 2852.0
## - har   1   2836.3 2852.3
## - EC    1   2836.6 2852.6
## <none>      2835.7 2853.7
## - TOC   1   2839.1 2855.1
## - chl   1   2839.1 2855.1
## - TDS   1   2840.6 2856.6
## 
## Step:  AIC=2851.88
## P ~ har + TDS + chl + sul + EC + TOC + THMs
## 
##        Df Deviance    AIC
## - sul   1   2836.0 2850.0
## - THMs  1   2836.2 2850.2
## - har   1   2836.5 2850.5
## - EC    1   2836.7 2850.7
## <none>      2835.9 2851.9
## - chl   1   2839.2 2853.2
## - TOC   1   2839.2 2853.2
## - TDS   1   2840.8 2854.8
## 
## Step:  AIC=2850.03
## P ~ har + TDS + chl + EC + TOC + THMs
## 
##        Df Deviance    AIC
## - THMs  1   2836.3 2848.3
## - har   1   2836.6 2848.6
## - EC    1   2836.9 2848.9
## <none>      2836.0 2850.0
## - chl   1   2839.3 2851.3
## - TOC   1   2839.4 2851.4
## - TDS   1   2841.3 2853.3
## 
## Step:  AIC=2848.32
## P ~ har + TDS + chl + EC + TOC
## 
##        Df Deviance    AIC
## - har   1   2836.9 2846.9
## - EC    1   2837.2 2847.2
## <none>      2836.3 2848.3
## - chl   1   2839.6 2849.6
## - TOC   1   2839.7 2849.7
## - TDS   1   2841.7 2851.7
## 
## Step:  AIC=2846.86
## P ~ TDS + chl + EC + TOC
## 
##        Df Deviance    AIC
## - EC    1   2837.8 2845.8
## <none>      2836.9 2846.9
## - chl   1   2840.2 2848.2
## - TOC   1   2840.2 2848.2
## - TDS   1   2842.3 2850.3
## 
## Step:  AIC=2845.76
## P ~ TDS + chl + TOC
## 
##        Df Deviance    AIC
## <none>      2837.8 2845.8
## - chl   1   2841.0 2847.0
## - TOC   1   2841.1 2847.1
## - TDS   1   2843.2 2849.2

Tom tat mo hinh xay dung => Mo hinh log = $\beta_0 + \beta_1.TDS + \beta_2.chl + \beta_3.TOC$ với $\beta_0 = -7.156*10^-1, \beta_1 = 9.83*10^-6, \beta_2 = 5.889*10^-2, \beta_3 = -2.59*10^-2$ Kiểm định hệ số $\beta_i$ $H_0: \beta_i = 0$ với i = 0, 1, 2, 3 $H_1: \beta_i$ khác 0 với i = 0, 1, 2, 3 => Giá trị P-value (xem Pr>|z|) đều < 0.05 nên bác bỏ H = 0 nên các $\beta_i$ đều có ý nghĩa thống kê

summary(mh)

## 
## Call:
## glm(formula = P ~ TDS + chl + TOC, family = binomial, data = xd)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -7.298e-01  3.073e-01  -2.375   0.0176 *
## TDS          1.190e-05  5.122e-06   2.322   0.0202 *
## chl          5.119e-02  2.837e-02   1.804   0.0712 .
## TOC         -2.412e-02  1.323e-02  -1.823   0.0682 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2849.1  on 2129  degrees of freedom
## Residual deviance: 2837.8  on 2126  degrees of freedom
## AIC: 2845.8
## 
## Number of Fisher Scoring iterations: 4

Tim khoang tin cay 95% của $\beta_i$

confint(mh, level = 0.95)

## Waiting for profiling to be done...

##                     2.5 %        97.5 %
## (Intercept) -1.334026e+00 -0.1287788356
## TDS          1.856755e-06  0.0000219461
## chl         -4.348774e-03  0.1069342458
## TOC         -5.010180e-02  0.0017729512

Tim OR và khoảng tin cậy 95% của OR

exp(cbind(OR = coef(mh), confint(mh, level = 0.95)))

## Waiting for profiling to be done...

##                    OR     2.5 %    97.5 %
## (Intercept) 0.4820069 0.2634147 0.8791684
## TDS         1.0000119 1.0000019 1.0000219
## chl         1.0525225 0.9956607 1.1128611
## TOC         0.9761690 0.9511326 1.0017745

Kiểm tra tính chính xác của mô hình

pred <- predict(mh, newdata = kd, type = "response",)
kd$pred <- round(pred)
head(kd, 10)

##           ph      har      TDS       chl      sul       EC       TOC     THMs
## 2   3.716080 129.4229 18630.06  6.635246 333.0735 592.8854 15.180013 56.32908
## 5   9.092223 181.1015 17978.99  6.546600 310.1357 398.4108 11.558279 31.99799
## 6   5.584087 188.3133 28748.69  7.544869 326.6784 280.4679  8.399735 54.91786
## 7  10.223862 248.0717 28749.72  7.513408 393.6634 283.6516 13.789695 84.60356
## 8   8.635849 203.3615 13672.09  4.563009 303.3098 474.6076 12.363817 62.79831
## 12  7.974522 218.6933 18767.66  8.110385 333.0735 364.0982 14.525746 76.48591
## 15  7.496232 205.3450 28388.00  5.072558 333.0735 444.6454 13.228311 70.30021
## 18  9.181560 273.8138 24041.33  6.904990 398.3505 477.9746 13.387341 71.45736
## 20  7.371050 214.4966 25630.32  4.432669 335.7544 469.9146 12.509164 62.79728
## 21  7.036752 227.4350 22305.57 10.333918 333.0735 554.8201 16.331693 45.38282
##         NTU P pred
## 2  4.500656 0    0
## 5  4.075075 0    0
## 6  2.559708 0    0
## 7  2.672989 0    0
## 8  4.401425 0    0
## 12 4.011718 0    0
## 15 4.777382 0    0
## 18 4.503661 0    0
## 20 2.560299 0    0
## 21 4.133423 0    0

Tạo bảng so sánh dữ liệu thực/dự đoán

quansat <- table(kd$P)
dudoan <- table(kd$pred)
rbind(quansat, dudoan)

##            0   1
## quansat  699 447
## dudoan  1139   7

Ma tran nham lan

table(kd$P, kd$pred)

##    
##       0   1
##   0 696   3
##   1 443   4

Tinh toan do chinh xac, do nhay, do dac hieu

dcx <- (table(kd$P, kd$pred)[1,1]+table(kd$P, kd$pred)[2,2])/nrow(kd)
dcx

## [1] 0.6108202

dn <- table(kd$P, kd$pred)[2,2]/(table(kd$P, kd$pred)[2,1]+table(kd$P, kd$pred)[2,2])
dn

## [1] 0.008948546

ddh <- table(kd$P, kd$pred)[1,1]/(table(kd$P, kd$pred)[1,1]+table(kd$P, kd$pred)[1,2])
ddh

## [1] 0.9957082

Ve ROC

library(ROCR)
ROCpred <- prediction(pred, kd$P)
ROCperf <- performance(ROCpred, "tpr", "fpr")
plot(ROCperf)

Tinh AUC

library(Metrics)
auc(kd$P, kd$pred)

## [1] 0.5023284

Anova 1 nhan to (Phan tich phg sai TDS, chl, TOC anh huong toi P nhu the nao)

Theo TDS

Phan tich phuong sai < Doc code: Df (degrees of freedom) là bậc tự do; Sum Sq là tổng bình phương (sum of squares), Mean Sq là trung bình bình phương (mean square); F value là giá trị F; và Pr(>F) là trị số P liên quan đến kiểm định F. Dòng clean.data$P trong kết quả trên có nghĩa là bình phương giữa các nhóm (betweengroups) và residual là bình phương trong mỗi nhóm (within-group). >

anv.TDS <- aov(clean.data$TDS ~ clean.data$P)

Chi tiet phan tich phuong sai

summary(anv.TDS)

##                Df    Sum Sq   Mean Sq F value Pr(>F)  
## clean.data$P    1 2.867e+08 286711019   3.732 0.0535 .
## Residuals    3274 2.515e+11  76823747                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theo chl

Phan tich phuong sai

anv.chl <- aov(clean.data$chl ~ clean.data$P)

Chi tiet phan tich phuong sai

summary(anv.chl)

##                Df Sum Sq Mean Sq F value Pr(>F)
## clean.data$P    1      5   4.641   1.852  0.174
## Residuals    3274   8203   2.506

Theo TOC

Phan tich phuong sai

anv.TOC <- aov(clean.data$TOC ~ clean.data$P)

Chi tiet phan tich phuong sai

summary(anv.TOC)

##                Df Sum Sq Mean Sq F value Pr(>F)  
## clean.data$P    1     32   32.26    2.95  0.086 .
## Residuals    3274  35809   10.94                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova 2 nhan to (theo dieu kien la P va pH)

Bien pH thanh bien dieu kien

clean.data$ph <- as.factor(clean.data$ph)

Theo TDS

Phan tich

anv2w.TDS <- aov(clean.data$TDS ~ clean.data$P + clean.data$ph)

Chi tiet

summary(anv2w.TDS)

##                 Df    Sum Sq   Mean Sq F value Pr(>F)  
## clean.data$P     1 2.867e+08 286711019   3.506 0.0617 .
## clean.data$ph 2784 2.115e+11  75952877   0.929 0.8623  
## Residuals      490 4.007e+10  81771710                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theo chl

Phan tich

anv2w.chl <- aov(clean.data$chl ~ clean.data$P + clean.data$ph)

Chi tiet

summary(anv2w.chl)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## clean.data$P     1      5   4.641   1.946  0.164
## clean.data$ph 2784   7034   2.527   1.059  0.209
## Residuals      490   1169   2.385

Theo TOC

Phan tich

anv2w.TOC <- aov(clean.data$TOC ~ clean.data$P + clean.data$ph)

Chi tiet

summary(anv2w.TOC)

##                 Df Sum Sq Mean Sq F value Pr(>F)  
## clean.data$P     1     32   32.26   2.916 0.0883 .
## clean.data$ph 2784  30389   10.92   0.987 0.5825  
## Residuals      490   5420   11.06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

XSTK

Truong Nguyen Trung Tin

2024-04-05

Nhap du lieu

Doc du lieu

Dat lai ten bien

Doc va kiem tra coi bien da duoc doi ten chua

Kiem tra so luong NA tung cot

Xu ly du lieu khuyet

Xu ly NA bang median

Kiem tra lai xem da xu ly dc NA chua

In clean.data

Dat bien phan loai va tom tat du lieu

Ve do thi thong ke Potability

Ve Histogram phan phoi cac bien

Ve do thi boxplot giua cac bien lien tuc va P

Ma tran va do thi he so tuong quan

Logistic

Anova 1 nhan to (Phan tich phg sai TDS, chl, TOC anh huong toi P nhu the nao)

Theo TDS

Theo chl

Theo TOC

Anova 2 nhan to (theo dieu kien la P va pH)

Bien pH thanh bien dieu kien

Theo TDS

Theo chl

Theo TOC