D4_exercise
Ha Tan Duc
2022-09-18
Việc 1: Đọc dữ liệu Insurance dataset.xlsx vào R và gọi
đối tượng là ins
ins = readxl::read_excel("D:\\OneDrive\\Statistical courses\\Can Tho University of Medicine_Sep2022\\Data for practice\\Insurance dataset.xlsx")
names (ins)## [1] "age" "sex" "bmi" "children" "smoker" "region" "charge"head (ins)## # A tibble: 6 x 7
## age sex bmi children smoker region charge
## <dbl> <chr> <dbl> <dbl> <chr> <chr> <dbl>
## 1 19 female 27.9 0 yes southwest 16885.
## 2 18 male 33.8 1 no southeast 1726.
## 3 28 male 33 3 no southeast 4449.
## 4 33 male 22.7 0 no northwest 21984.
## 5 32 male 28.9 0 no northwest 3867.
## 6 31 female 25.7 0 no southeast 3757.- Vẽ biểu đồ histogram của biến
chargevà ghi lại một nhận xét
library (ggplot2)
ggplot(data= ins, aes(x= charge)) + geom_histogram(aes(y= ..density..), color="white", fill="blue") + geom_density(col="red")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- Vẽ biểu đồ boxplot với
y = chargevàx = smoker. Fit mô hình hồi quycharge = alpha + beta*smoker. Diễn giải ý nghĩa của mô hình này.
ggplot(data= ins, aes(x= smoker, y= charge, fill= smoker)) + geom_boxplot()
Fit mô hình hồi quy charge = alpha + beta*smoker. Diễn
giải ý nghĩa của mô hình
m1 = lm (charge ~ smoker, data = ins)
summary (m1)##
## Call:
## lm(formula = charge ~ smoker, data = ins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19221 -5042 -919 3705 31720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8434.3 229.0 36.83 <2e-16 ***
## smokeryes 23616.0 506.1 46.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7470 on 1336 degrees of freedom
## Multiple R-squared: 0.6198, Adjusted R-squared: 0.6195
## F-statistic: 2178 on 1 and 1336 DF, p-value: < 2.2e-16Việc 2: Phân tích các mô hình sau đây:
Mô hình 1: charge= α + β*sex
Mô hình 2: charge = α + β*sex + γ*smoker
So sánh tham số β giữa 2 mô hình. Diễn giải ý nghĩa của tham số β của mô hình 1 và mô hình 2.
summary (lm (charge ~ sex, data = ins))##
## Call:
## lm(formula = charge ~ sex, data = ins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12835 -8435 -3980 3476 51201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12569.6 470.1 26.740 <2e-16 ***
## sexmale 1387.2 661.3 2.098 0.0361 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12090 on 1336 degrees of freedom
## Multiple R-squared: 0.003282, Adjusted R-squared: 0.002536
## F-statistic: 4.4 on 1 and 1336 DF, p-value: 0.03613summary (lm (charge ~ sex + smoker, data= ins))##
## Call:
## lm(formula = charge ~ sex + smoker, data = ins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19193 -5074 -909 3739 31682
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8466.04 303.54 27.89 <2e-16 ***
## sexmale -65.38 409.81 -0.16 0.873
## smokeryes 23622.13 507.74 46.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7473 on 1335 degrees of freedom
## Multiple R-squared: 0.6198, Adjusted R-squared: 0.6192
## F-statistic: 1088 on 2 and 1335 DF, p-value: < 2.2e-16Việc 3: Phân tích mô hình tương tác giữa sex và smoker
Mô hình 3:
charge= α + β*sex + γ*smoker + δ*sex*smoker
Viết ra phương trình của mô hình 3, và diễn giải ý nghĩa của mô hình.
summary (lm (charge ~ sex + smoker + sex*smoker, data = ins))##
## Call:
## lm(formula = charge ~ sex + smoker + sex * smoker, data = ins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20213 -5166 -920 3655 33091
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8762.3 318.6 27.504 < 2e-16 ***
## sexmale -675.1 457.0 -1.477 0.13988
## smokeryes 21916.7 764.4 28.673 < 2e-16 ***
## sexmale:smokeryes 3038.1 1020.2 2.978 0.00295 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7451 on 1334 degrees of freedom
## Multiple R-squared: 0.6223, Adjusted R-squared: 0.6214
## F-statistic: 732.6 on 3 and 1334 DF, p-value: < 2.2e-16Việc 4:
Đọc dữ liệu SAT scores dataset.xlsx vào R và gọi đối
tượng là sat. Dữ liệu này là một nghiên cứu về mối liên
quan giữa điểm ‘grade point average’ (gpa) sau 3 semester ở đại học, và
điểm lúc theo học trung học. Có 224 học sinh tham gia nghiên cứu. Dữ
liệu bao gồm các biến số sau đây:
Obs mã số của học sinh gpa grade point
average (0=low, 4=high) hsm điểm trung học môn toán (0-10)
hss điểm trung học môn khoa học (0-10) hse
điểm trung học môn tiếng Anh (0-10) satm điểm SAT môn toán
satv điểm SAT về verbal knowledge sex giới
tính (1/2 = male/female)
sat = readxl::read_excel("D:\\OneDrive\\Statistical courses\\Can Tho University of Medicine_Sep2022\\Data for practice\\SAT scores dataset.xlsx")
names (sat)## [1] "Obs" "gpa" "hsm" "hss" "hse" "satm" "satv" "sex"head (sat)## # A tibble: 6 x 8
## Obs gpa hsm hss hse satm satv sex
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 3.32 10 10 10 670 600 1
## 2 2 2.26 6 8 5 700 640 1
## 3 3 2.35 8 6 8 640 530 1
## 4 4 2.08 9 10 7 670 600 1
## 5 5 3.38 8 9 8 540 580 1
## 6 6 3.29 10 8 8 760 630 1Hoán chuyển các biến hsm, hss, hse, satm và satv sao cho có giá trị trung bình là 0 và độ lệch chuẩn là 1 (z-score):
sat$hsm = scale(sat$hsm)
sat$hss = scale(sat$hss)
sat$hse = scale(sat$hse)
sat$satm = scale(sat$satm)
sat$satv = scale(sat$satv)
sat = data.frame (sat) #sau dùng hàm scale cần định nghĩa lại dataset thành dataframe- Vẽ biểu đồ tương quan giữa các biến
gpa,hsm,hss,hse,satmvàsatv. Dùng ma trận các biểu đồ tương quan.
library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphatemp = dplyr::select(sat, c(gpa, hsm, hss, hse, satm, satv))
pairs.panels (temp)
- Phân tích mô hình
gpa = α + β*hsm + ε. Chú ý hệ số hồi qui liên quan đếnhsmvà diễn giải ý nghĩa của hệ số.
summary (lm (gpa ~ hsm, data= sat))##
## Call:
## lm(formula = gpa ~ hsm, data = sat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.12610 -0.36370 0.08651 0.46351 1.63911
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.63522 0.04696 56.119 < 2e-16 ***
## hsm 0.34020 0.04706 7.229 7.77e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7028 on 222 degrees of freedom
## Multiple R-squared: 0.1905, Adjusted R-squared: 0.1869
## F-statistic: 52.25 on 1 and 222 DF, p-value: 7.774e-12- Phân tích mô hình
gpa = α + β*satm + ε. Chú ý hệ số hồi qui liên quan đếnsatmvà diễn giải ý nghĩa của hệ số.
summary (lm (gpa ~ satm, data= sat))##
## Call:
## lm(formula = gpa ~ satm, data = sat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.59405 -0.37876 0.08366 0.55984 1.39948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.63522 0.05051 52.170 < 2e-16 ***
## satm 0.19618 0.05063 3.875 0.00014 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.756 on 222 degrees of freedom
## Multiple R-squared: 0.06336, Adjusted R-squared: 0.05914
## F-statistic: 15.02 on 1 and 222 DF, p-value: 0.0001402- Phân tích mô hình
gpa = α + β*satm + γ*hsm + ε. Chú ý hệ số hồi qui liên quan đếnsatmvàhsmvà diễn giải tại sao kết quả này khác với kết quả ở b. và c.
summary (lm (gpa ~ satm + hsm, data= sat))##
## Call:
## lm(formula = gpa ~ satm + hsm, data = sat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1435 -0.3292 0.0793 0.4591 1.6169
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.63522 0.04696 56.118 < 2e-16 ***
## satm 0.05275 0.05281 0.999 0.319
## hsm 0.31628 0.05281 5.990 8.45e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7028 on 221 degrees of freedom
## Multiple R-squared: 0.1942, Adjusted R-squared: 0.1869
## F-statistic: 26.63 on 2 and 221 DF, p-value: 4.365e-11Việc 5: Tìm mô hình tối ưu
Biến phụ thuộc của dữ liệu này là gpa. Các biến tiên
lượng là hsm, hss, hse,
satm, satv và sex. Bạn hãy tìm
các biến liên quan và từ đó xây dựng mô hình tiên lượng GPA dựa vào các
biến liên quan.
library (BMA)## Loading required package: survival## Loading required package: leaps## Warning: package 'leaps' was built under R version 4.0.5## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.0.5##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Warning: package 'inline' was built under R version 4.0.5## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-1)xvars = sat [, c('hsm', 'hss', 'hse', 'satm', 'satv')]
r = bicreg (x= xvars, y = sat$gpa, maxCol= 31)
summary (r)##
## Call:
## bicreg(x = xvars, y = sat$gpa, maxCol = 31)
##
##
## 5 models were selected
## Best 5 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100.0 2.6352232 0.04690 2.63522 2.63522 2.63522
## hsm 100.0 0.3239816 0.05433 0.34020 0.29932 0.28777
## hss 13.8 0.0125371 0.03792 . . 0.09109
## hse 18.0 0.0164316 0.04154 . 0.09148 .
## satm 6.4 0.0033760 0.01858 . . .
## satv 4.0 0.0005964 0.01016 . . .
##
## nVar 1 2 2
## r2 0.191 0.202 0.200
## BIC -41.93648 -39.59809 -39.06568
## post prob 0.578 0.180 0.138
## model 4 model 5
## Intercept 2.63522 2.63522
## hsm 0.31628 0.33695
## hss . .
## hse . .
## satm 0.05275 .
## satv . 0.01473
##
## nVar 2 2
## r2 0.194 0.191
## BIC -37.53438 -36.61894
## post prob 0.064 0.040Việc 6: Đánh giá tầm quan trọng của các biến số
Dùng phương pháp “relative importance analysis”, các bạn hãy đánh giá tầm quan trọng của các biến tiên lượng GPA.
m = lm(gpa ~ hsm + hss + hse + satm + satv, data= sat)
relaimpo::calc.relimp(m, type="lmg", rela=T, rank=T)## Response variable: gpa
## Total response variance: 0.6074565
## Analysis based on 224 observations
##
## 5 Regressors:
## hsm hss hse satm satv
## Proportion of variance explained by model: 21.15%
## Metrics are normalized to sum to 100% (rela=TRUE).
##
## Relative importance metrics:
##
## lmg
## hsm 0.50102158
## hss 0.19336795
## hse 0.15632091
## satm 0.13173217
## satv 0.01755738
##
## Average coefficients for different model sizes:
##
## 1X 2Xs 3Xs 4Xs 5Xs
## hsm 0.34020499 0.31008030 0.28463658 0.26096075 0.23919178
## hss 0.25675243 0.18853263 0.13503334 0.09207860 0.06102693
## hse 0.22524611 0.15734863 0.11853408 0.09595220 0.08337422
## satm 0.19618484 0.14162783 0.11384806 0.09612545 0.08152776
## satv 0.08923328 0.01815515 -0.01558741 -0.03133515 -0.03777113