XSTK

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.7)
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=3073.72
## P ~ ph + har + TDS + chl + sul + EC + TOC + THMs + NTU
## 
##        Df Deviance    AIC
## - ph    1   3053.7 3071.7
## - NTU   1   3053.7 3071.7
## - THMs  1   3054.0 3072.0
## - sul   1   3054.2 3072.2
## - EC    1   3054.4 3072.4
## - har   1   3054.9 3072.9
## <none>      3053.7 3073.7
## - TDS   1   3057.1 3075.1
## - TOC   1   3057.8 3075.8
## - chl   1   3058.4 3076.4
## 
## Step:  AIC=3071.73
## P ~ har + TDS + chl + sul + EC + TOC + THMs + NTU
## 
##        Df Deviance    AIC
## - NTU   1   3053.8 3069.8
## - THMs  1   3054.0 3070.0
## - sul   1   3054.2 3070.2
## - EC    1   3054.4 3070.4
## - har   1   3054.9 3070.9
## <none>      3053.7 3071.7
## - TDS   1   3057.1 3073.1
## - TOC   1   3057.8 3073.8
## - chl   1   3058.4 3074.4
## 
## Step:  AIC=3069.75
## P ~ har + TDS + chl + sul + EC + TOC + THMs
## 
##        Df Deviance    AIC
## - THMs  1   3054.0 3068.0
## - sul   1   3054.2 3068.2
## - EC    1   3054.4 3068.4
## - har   1   3054.9 3068.9
## <none>      3053.8 3069.8
## - TDS   1   3057.2 3071.2
## - TOC   1   3057.9 3071.9
## - chl   1   3058.5 3072.5
## 
## Step:  AIC=3068.04
## P ~ har + TDS + chl + sul + EC + TOC
## 
##        Df Deviance    AIC
## - sul   1   3054.5 3066.5
## - EC    1   3054.7 3066.7
## - har   1   3055.2 3067.2
## <none>      3054.0 3068.0
## - TDS   1   3057.5 3069.5
## - TOC   1   3058.1 3070.1
## - chl   1   3058.8 3070.8
## 
## Step:  AIC=3066.48
## P ~ har + TDS + chl + EC + TOC
## 
##        Df Deviance    AIC
## - EC    1   3055.1 3065.1
## - har   1   3055.5 3065.5
## <none>      3054.5 3066.5
## - TDS   1   3058.4 3068.4
## - TOC   1   3058.7 3068.7
## - chl   1   3059.2 3069.2
## 
## Step:  AIC=3065.14
## P ~ har + TDS + chl + TOC
## 
##        Df Deviance    AIC
## - har   1   3056.2 3064.2
## <none>      3055.1 3065.1
## - TDS   1   3059.0 3067.0
## - TOC   1   3059.3 3067.3
## - chl   1   3059.8 3067.8
## 
## Step:  AIC=3064.22
## P ~ TDS + chl + TOC
## 
##        Df Deviance    AIC
## <none>      3056.2 3064.2
## - TDS   1   3060.2 3066.2
## - TOC   1   3060.3 3066.3
## - chl   1   3060.9 3066.9

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.156e-01  2.974e-01  -2.406   0.0161 *
## TDS          9.830e-06  4.940e-06   1.990   0.0466 *
## chl          5.889e-02  2.727e-02   2.159   0.0308 *
## TOC         -2.590e-02  1.276e-02  -2.029   0.0424 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3068.5  on 2293  degrees of freedom
## Residual deviance: 3056.2  on 2290  degrees of freedom
## AIC: 3064.2
## 
## 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.300278e+00 -1.338911e-01
## TDS          1.433311e-07  1.951725e-05
## chl          5.520962e-03  1.124815e-01
## TOC         -5.097279e-02 -9.200715e-04

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.4888938 0.2724561 0.8746853
## TDS         1.0000098 1.0000001 1.0000195
## chl         1.0606592 1.0055362 1.1190516
## TOC         0.9744300 0.9503045 0.9990804

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
## 21  7.036752 227.4350 22305.57 10.333918 333.0735 554.8201 16.331693 45.38282
## 22  6.660212 168.2837 30944.36  5.858769 310.9309 523.6713 17.884235 77.04232
## 25  5.400302 140.7391 17266.59 10.056852 328.3582 472.8741 11.256381 56.93191
##         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
## 21 4.133423 0    0
## 22 3.749701 0    0
## 25 4.824786 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 599 383
## dudoan  979   3

Ma tran nham lan

table(kd$P, kd$pred)

##    
##       0   1
##   0 597   2
##   1 382   1

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.6089613

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.002610966

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.9966611

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.499636