資料來源:https://www.kaggle.com/datasets/fatemehmohammadinia/heart-attack-dataset-tarik-a-rashid

目的:資料給予我們心率、收縮壓、舒張壓、血糖、CK同工梅、肌鈣蛋白，我們去分析這些變數是否與心臟病有相關性。

資料筆數:1319筆

反應變數:Result(結果)

解釋變數:

Heart rate(心率)

Systolic blood pressure(收縮壓)

Diastolic blood pressure（舒張壓)

Blood sugar(血糖)

CK-MB(CK同工梅)

Troponin(肌鈣蛋白)

讀入檔案

heart <- read.csv("C:/Users/User/Desktop/Medicaldataset.csv")
head(heart)
##   Age Gender Heart.rate Systolic.blood.pressure Diastolic.blood.pressure
## 1  64      1         66                     160                       83
## 2  21      1         94                      98                       46
## 3  55      1         64                     160                       77
## 4  64      1         70                     120                       55
## 5  55      1         64                     112                       65
## 6  58      0         61                     112                       58
##   Blood.sugar CK.MB Troponin   Result
## 1         160  1.80    0.012 negative
## 2         296  6.75    1.060 positive
## 3         270  1.99    0.003 negative
## 4         270 13.87    0.122 positive
## 5         300  1.08    0.003 negative
## 6          87  1.83    0.004 negative

我們沒有用到年齡、性別，因此欄位選取時沒有選擇Age、Gender欄位

畫出各變數的矩陣散佈圖

data=heart[,c(2:9)]
data <- na.omit(data)
head(data)
##   Gender Heart.rate Systolic.blood.pressure Diastolic.blood.pressure
## 1      1         66                     160                       83
## 2      1         94                      98                       46
## 3      1         64                     160                       77
## 4      1         70                     120                       55
## 5      1         64                     112                       65
## 6      0         61                     112                       58
##   Blood.sugar CK.MB Troponin   Result
## 1         160  1.80    0.012 negative
## 2         296  6.75    1.060 positive
## 3         270  1.99    0.003 negative
## 4         270 13.87    0.122 positive
## 5         300  1.08    0.003 negative
## 6          87  1.83    0.004 negative
plot(data)

為了簡化程式碼，所以將反應變數以Y表示，解釋變數以X1~X6表示

X1=data$Heart.rate
X2=data$Systolic.blood.pressure
X3=data$Diastolic.blood.pressure
X4=data$Blood.sugar
X5=data$CK.MB
X6=data$Troponin
Y<-ifelse(data$Result == "positive", 1, 2)

配適最複雜模型m1，以vif檢查共線性

library(car)
## 載入需要的套件：carData
m1=lm(Y~X1+X2+X3+X4+X5+X6)
summary(m1)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6019 -0.4250 -0.2906  0.5439  1.5255 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.352e+00  7.857e-02  17.213   <2e-16 ***
## X1          -6.330e-05  2.484e-04  -0.255   0.7989    
## X2           4.775e-04  6.031e-04   0.792   0.4287    
## X3           5.412e-05  1.129e-03   0.048   0.9618    
## X4           3.094e-04  1.702e-04   1.817   0.0694 .  
## X5          -2.347e-03  2.750e-04  -8.535   <2e-16 ***
## X6          -9.915e-02  1.104e-02  -8.984   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4618 on 1312 degrees of freedom
## Multiple R-squared:  0.1047, Adjusted R-squared:  0.1006 
## F-statistic: 25.57 on 6 and 1312 DF,  p-value: < 2.2e-16
vif(m1)
##       X1       X2       X3       X4       X5       X6 
## 1.016593 1.533638 1.552041 1.005361 1.002978 1.003170

因為在vif檢查時，沒有vif≧10的值，因此選擇不刪除

運用向前選取法、向後選取法和逐步選取法，來選取最佳適配模型

none=lm(Y~1)
step(none,scope =list(upper=m1,lower=none),direction='forward' ) #向前選取法
## Start:  AIC=-1897.06
## Y ~ 1
## 
##        Df Sum of Sq    RSS     AIC
## + X6    1   16.4458 296.13 -1966.3
## + X5    1   14.8168 297.76 -1959.1
## <none>              312.58 -1897.1
## + X4    1    0.3416 312.24 -1896.5
## + X2    1    0.1356 312.44 -1895.6
## + X3    1    0.0292 312.55 -1895.2
## + X1    1    0.0150 312.56 -1895.1
## 
## Step:  AIC=-1966.35
## Y ~ X6
## 
##        Df Sum of Sq    RSS     AIC
## + X5    1   15.3247 280.81 -2034.4
## + X4    1    0.4490 295.68 -1966.4
## <none>              296.13 -1966.3
## + X2    1    0.2982 295.83 -1965.7
## + X3    1    0.1204 296.01 -1964.9
## + X1    1    0.0059 296.13 -1964.4
## 
## Step:  AIC=-2034.44
## Y ~ X6 + X5
## 
##        Df Sum of Sq    RSS     AIC
## + X4    1   0.72503 280.08 -2035.8
## <none>              280.81 -2034.4
## + X2    1   0.23483 280.57 -2033.5
## + X3    1   0.06657 280.74 -2032.8
## + X1    1   0.01619 280.79 -2032.5
## 
## Step:  AIC=-2035.85
## Y ~ X6 + X5 + X4
## 
##        Df Sum of Sq    RSS     AIC
## <none>              280.08 -2035.8
## + X2    1  0.218181 279.86 -2034.9
## + X3    1  0.078316 280.00 -2034.2
## + X1    1  0.012290 280.07 -2033.9
## 
## Call:
## lm(formula = Y ~ X6 + X5 + X4)
## 
## Coefficients:
## (Intercept)           X6           X5           X4  
##   1.4114709   -0.0986887   -0.0023512    0.0003134
step(m1,scope =list(upper=m1,lower=none),direction='backward' )  #向後選取法
## Start:  AIC=-2030.94
## Y ~ X1 + X2 + X3 + X4 + X5 + X6
## 
##        Df Sum of Sq    RSS     AIC
## - X3    1    0.0005 279.85 -2032.9
## - X1    1    0.0138 279.86 -2032.9
## - X2    1    0.1337 279.98 -2032.3
## <none>              279.85 -2030.9
## - X4    1    0.7046 280.56 -2029.6
## - X5    1   15.5386 295.39 -1961.7
## - X6    1   17.2173 297.07 -1954.2
## 
## Step:  AIC=-2032.94
## Y ~ X1 + X2 + X4 + X5 + X6
## 
##        Df Sum of Sq    RSS     AIC
## - X1    1    0.0134 279.86 -2034.9
## - X2    1    0.2193 280.07 -2033.9
## <none>              279.85 -2032.9
## - X4    1    0.7043 280.56 -2031.6
## - X5    1   15.5436 295.39 -1963.6
## - X6    1   17.2213 297.07 -1956.2
## 
## Step:  AIC=-2034.88
## Y ~ X2 + X4 + X5 + X6
## 
##        Df Sum of Sq    RSS     AIC
## - X2    1    0.2182 280.08 -2035.8
## <none>              279.86 -2034.9
## - X4    1    0.7084 280.57 -2033.5
## - X5    1   15.5350 295.40 -1965.6
## - X6    1   17.2338 297.10 -1958.1
## 
## Step:  AIC=-2035.85
## Y ~ X4 + X5 + X6
## 
##        Df Sum of Sq    RSS     AIC
## <none>              280.08 -2035.8
## - X4    1     0.725 280.81 -2034.4
## - X5    1    15.601 295.68 -1966.4
## - X6    1    17.099 297.18 -1959.7
## 
## Call:
## lm(formula = Y ~ X4 + X5 + X6)
## 
## Coefficients:
## (Intercept)           X4           X5           X6  
##   1.4114709    0.0003134   -0.0023512   -0.0986887
step(none,scope =list(upper=m1,lower=none),direction='both' )    #逐步選取法
## Start:  AIC=-1897.06
## Y ~ 1
## 
##        Df Sum of Sq    RSS     AIC
## + X6    1   16.4458 296.13 -1966.3
## + X5    1   14.8168 297.76 -1959.1
## <none>              312.58 -1897.1
## + X4    1    0.3416 312.24 -1896.5
## + X2    1    0.1356 312.44 -1895.6
## + X3    1    0.0292 312.55 -1895.2
## + X1    1    0.0150 312.56 -1895.1
## 
## Step:  AIC=-1966.35
## Y ~ X6
## 
##        Df Sum of Sq    RSS     AIC
## + X5    1   15.3247 280.81 -2034.4
## + X4    1    0.4490 295.68 -1966.4
## <none>              296.13 -1966.3
## + X2    1    0.2982 295.83 -1965.7
## + X3    1    0.1204 296.01 -1964.9
## + X1    1    0.0059 296.13 -1964.4
## - X6    1   16.4458 312.58 -1897.1
## 
## Step:  AIC=-2034.44
## Y ~ X6 + X5
## 
##        Df Sum of Sq    RSS     AIC
## + X4    1    0.7250 280.08 -2035.8
## <none>              280.81 -2034.4
## + X2    1    0.2348 280.57 -2033.5
## + X3    1    0.0666 280.74 -2032.8
## + X1    1    0.0162 280.79 -2032.5
## - X5    1   15.3247 296.13 -1966.3
## - X6    1   16.9537 297.76 -1959.1
## 
## Step:  AIC=-2035.85
## Y ~ X6 + X5 + X4
## 
##        Df Sum of Sq    RSS     AIC
## <none>              280.08 -2035.8
## + X2    1    0.2182 279.86 -2034.9
## - X4    1    0.7250 280.81 -2034.4
## + X3    1    0.0783 280.00 -2034.2
## + X1    1    0.0123 280.07 -2033.9
## - X5    1   15.6007 295.68 -1966.4
## - X6    1   17.0990 297.18 -1959.7
## 
## Call:
## lm(formula = Y ~ X6 + X5 + X4)
## 
## Coefficients:
## (Intercept)           X6           X5           X4  
##   1.4114709   -0.0986887   -0.0023512    0.0003134

★發現最佳適配模型是Y~X4+X5+X6 重新配m2

m2=lm(Y~X4+X5+X6)
summary(m2)
## 
## Call:
## lm(formula = Y ~ X4 + X5 + X6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5687 -0.4260 -0.2916  0.5500  1.5245 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.4114709  0.0283211  49.838   <2e-16 ***
## X4           0.0003134  0.0001699   1.845   0.0653 .  
## X5          -0.0023512  0.0002747  -8.558   <2e-16 ***
## X6          -0.0986887  0.0110144  -8.960   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4615 on 1315 degrees of freedom
## Multiple R-squared:  0.104,  Adjusted R-squared:  0.1019 
## F-statistic: 50.86 on 3 and 1315 DF,  p-value: < 2.2e-16

以下是殘差分析: 畫關係圖

par(mfrow=c(2,2))
plot(m2)

1.檢查是否具均齊性

\(H_0\):變異數具均齊性

\(H_1\):變異數不具均齊性

ncvTest(m2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 19.23273, Df = 1, p = 1.1571e-05
library(lmtest)
## 載入需要的套件：zoo
## 
## 載入套件：'zoo'
## 下列物件被遮斷自 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(m2)
## 
##  studentized Breusch-Pagan test
## 
## data:  m2
## BP = 195.12, df = 3, p-value < 2.2e-16

►ncvTest 檢定結果之p值<0.05，推翻虛無假設，表示變異數不具均齊性

►bptest 檢定結果之p值<0.05，推翻虛無假設，表示變異數不具均齊性

2.檢定殘差是否具常態性

\(H_0\):殘差具常態性

\(H_1\):殘差不具常態性

e=m2$residuals #殘差
qqnorm(e)
qqline(e)

shapiro.test(e)
## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.73941, p-value < 2.2e-16

►因為shapiro.test結果之p值<0.05，所以推翻虛無假設，表示殘差不具常態性

3.檢定殘差是否具一階自我相關

\(H_0\):殘差無一階自我相關

\(H_1\):殘差有一階自我相關

durbinWatsonTest(m2) #雙尾檢定
##  lag Autocorrelation D-W Statistic p-value
##    1      0.00107487      1.996704   0.936
##  Alternative hypothesis: rho != 0
dwtest(m2)           #單尾檢定
## 
##  Durbin-Watson test
## 
## data:  m2
## DW = 1.9967, p-value = 0.4763
## alternative hypothesis: true autocorrelation is greater than 0

►由於兩種檢定結果之p值>0.05，沒有推翻虛無假設，表示殘差不具一階自我相關

★因為不具均齊性、常態性，所以做Box-Cox轉換

做Box-Cox轉換

trans=powerTransform(m2)
summary(trans)
## bcPower Transformation to Normality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1   -2.0403          -2      -2.3234      -1.7572
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df       pval
## LR test, lambda = (0) 209.3967  1 < 2.22e-16
## 
## Likelihood ratio test that no transformation is needed
##                           LRT df       pval
## LR test, lambda = (1) 470.634  1 < 2.22e-16

\(\lambda\)=-2.0402917，以\(\lambda\)=-2.04取代來做newy

newy=bcPower(Y,-2.04)

以newy對單一變數來看是否顯著

summary(lm(newy~X4))
## 
## Call:
## lm(formula = newy ~ X4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1746 -0.1420 -0.1390  0.2279  0.2347 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.315e-01  1.094e-02   12.02   <2e-16 ***
## X4          7.972e-05  6.641e-05    1.20     0.23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1806 on 1317 degrees of freedom
## Multiple R-squared:  0.001093,   Adjusted R-squared:  0.0003345 
## F-statistic: 1.441 on 1 and 1317 DF,  p-value: 0.2302
summary(lm(newy~X5))
## 
## Call:
## lm(formula = newy ~ X5)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1558 -0.1537 -0.1396  0.2164  0.2208 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1561369  0.0051146  30.527  < 2e-16 ***
## X5          -0.0008491  0.0001049  -8.095 1.29e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1764 on 1317 degrees of freedom
## Multiple R-squared:  0.0474, Adjusted R-squared:  0.04668 
## F-statistic: 65.53 on 1 and 1317 DF,  p-value: 1.288e-15
summary(lm(newy~X6))
## 
## Call:
## lm(formula = newy ~ X6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1560 -0.1552 -0.1251  0.2150  0.5738 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.156123   0.005075  30.761   <2e-16 ***
## X6          -0.035894   0.004197  -8.552   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1759 on 1317 degrees of freedom
## Multiple R-squared:  0.05261,    Adjusted R-squared:  0.05189 
## F-statistic: 73.14 on 1 and 1317 DF,  p-value: < 2.2e-16

►發現變數X4顯著性不大，所以不用它

以newy做反應變數，X5、X6做解釋變數配適m3

m3=lm(newy~X5+X6)
summary(m3)
## 
## Call:
## lm(formula = newy ~ X5 + X6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1688 -0.1617 -0.1088  0.2032  0.5664 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1695149  0.0051905  32.659   <2e-16 ***
## X5          -0.0008636  0.0001019  -8.475   <2e-16 ***
## X6          -0.0364488  0.0040891  -8.914   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1714 on 1316 degrees of freedom
## Multiple R-squared:  0.1016, Adjusted R-squared:  0.1003 
## F-statistic: 74.45 on 2 and 1316 DF,  p-value: < 2.2e-16

畫關係圖

par(mfrow=c(2,2))
plot(m3)

1.檢查是否具均齊性

\(H_0\):變異數具均齊性

\(H_1\):變異數不具均齊性

library(car)
ncvTest(m3)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 19.39824, Df = 1, p = 1.061e-05
library(lmtest)
bptest(m3)
## 
##  studentized Breusch-Pagan test
## 
## data:  m3
## BP = 201.54, df = 2, p-value < 2.2e-16

►ncvTest 檢定結果之p值<0.05，推翻虛無假設，表示變異數不具均齊性

►bptest 檢定結果之p值<0.05，推翻虛無假設，表示變異數不具均齊性

2.檢定殘差是否具常態性

\(H_0\):殘差具常態性

\(H_1\):殘差不具常態性

e1=m3$residuals
qqnorm(e1)
qqline(e1)

shapiro.test(e1)
## 
##  Shapiro-Wilk normality test
## 
## data:  e1
## W = 0.717, p-value < 2.2e-16

►因為shapiro.test結果之p值<0.05，所以推翻虛無假設，表示殘差不具常態性

3.檢定殘差是否具一階自我相關

\(H_0\):殘差無一階自我相關

\(H_1\):殘差有一階自我相關

durbinWatsonTest(m3)
##  lag Autocorrelation D-W Statistic p-value
##    1     0.000214414      1.998404   0.954
##  Alternative hypothesis: rho != 0
dwtest(m3)
## 
##  Durbin-Watson test
## 
## data:  m3
## DW = 1.9984, p-value = 0.4882
## alternative hypothesis: true autocorrelation is greater than 0

►由於兩種檢定結果之p值均>0.05，沒有推翻虛無假設，表示殘差不具一階自我相關

因為m3模型X2變數不顯著，殘差不具均齊和常態性以及不具一階自我相關性

所以我們使用m3做離群值、影響點檢測

p=2    #解釋變數個數
n=nrow(data)
ei = m3$resid
hii = hatvalues(m3)
ri = rstandard(m3)
ti = rstudent(m3)
Di = cooks.distance(m3)
DFFITSi = dffits(m3)
DFBETASi = dfbetas(m3)
case1 = which(hii>2*(p+1)/n); names(case1)=NULL; case1
##  [1]    8   13   29   30   31   83   90   98  102  104  114  145  182  186  188
## [16]  206  240  279  334  394  422  428  430  434  435  436  446  449  452  467
## [31]  473  476  531  608  618  652  667  678  683  697  707  733  742  754  784
## [46]  795  797  807  809  810  864  921  925  931  936  973  989  998 1004 1029
## [61] 1048 1049 1095 1132 1145 1154 1210 1216 1225 1250 1253 1260 1281 1311 1317
case2 = which(abs(ri)>2); names(case2)=NULL; case2
## [1] 30
case3 = which(abs(ti)>qt(1-0.05/n, n-p-1)); names(case3)=NULL; case3
## integer(0)
case4 = which(Di>qf(0.5, p+1, n-(p+1))); names(case4)=NULL; case4
## integer(0)
case5 = which(abs(DFFITSi)>1); names(case5)=NULL; case5
## integer(0)
case6 = which(abs(DFBETASi)>1, arr.ind=T); names(case6)=NULL; case6
##      row col

我們看到case1和case2有重複的點(第30筆資料)，因此我們先刪除第30筆資料，重新配適m4

m4=lm(newy[-30]~X5[-30]+X6[-30])
summary(m4)
## 
## Call:
## lm(formula = newy[-30] ~ X5[-30] + X6[-30])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1694 -0.1623 -0.1066  0.2026  0.2596 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1702182  0.0051738  32.900   <2e-16 ***
## X5[-30]     -0.0008618  0.0001015  -8.491   <2e-16 ***
## X6[-30]     -0.0397313  0.0041848  -9.494   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1707 on 1315 degrees of freedom
## Multiple R-squared:  0.1084, Adjusted R-squared:  0.1071 
## F-statistic: 79.97 on 2 and 1315 DF,  p-value: < 2.2e-16
ncvTest(m4)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 25.32278, Df = 1, p = 4.8495e-07
bptest(m4)
## 
##  studentized Breusch-Pagan test
## 
## data:  m4
## BP = 329.58, df = 2, p-value < 2.2e-16
qqnorm(m4$residuals)
qqline(m4$residuals)

shapiro.test(m4$resid)
## 
##  Shapiro-Wilk normality test
## 
## data:  m4$resid
## W = 0.7178, p-value < 2.2e-16
durbinWatsonTest(m4) 
##  lag Autocorrelation D-W Statistic p-value
##    1    -0.003817371      2.006479   0.918
##  Alternative hypothesis: rho != 0
dwtest(m4)
## 
##  Durbin-Watson test
## 
## data:  m4
## DW = 2.0065, p-value = 0.5471
## alternative hypothesis: true autocorrelation is greater than 0
summary(m4)
## 
## Call:
## lm(formula = newy[-30] ~ X5[-30] + X6[-30])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1694 -0.1623 -0.1066  0.2026  0.2596 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.1702182  0.0051738  32.900   <2e-16 ***
## X5[-30]     -0.0008618  0.0001015  -8.491   <2e-16 ***
## X6[-30]     -0.0397313  0.0041848  -9.494   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1707 on 1315 degrees of freedom
## Multiple R-squared:  0.1084, Adjusted R-squared:  0.1071 
## F-statistic: 79.97 on 2 and 1315 DF,  p-value: < 2.2e-16

因為殘差分析出來的結果一樣不具均齊性、常態性且不具一階自我相關，但是殘差圖的結果m4比較好，所以我們選用m4當最終模型

\(\hat{heart}\) = 0.1702182 - 0.0008618\(^*\)\(CK.MB\) - 0.0397313\(^*\)\(Troponin\)

組員貢獻度:康凱綸 25%、薛育庭 25%、楊子謙 15%、劉亨峙 25 % 王瑋豪 10%