ngày 3 hồi quy tuyến tính
GS Nguyễn Văn Tuấn
1/7/2020
Package cần làm quen/sử dụng: BMA, visreg, GGally, table1, caret, relaimpo
Dữ liệu thực hành: “obesity data.csv”
Mục tiêu: đánh giá mối liên quan giữa pcfat với các biến gender, age, height, weight, bmi
Task 1: Đọc dữ liệu vào R, và gọi đối tượng là ob
ob = read.csv("~/Desktop/R GS Tuấn 2020/TDTU Datasets for 2020 Workshop/Obesity data.csv", 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.2Chọn biến liên quan
dat = ob[, c("gender", "age", "height", "weight", "bmi",
"pcfat")]Task 2: Khai thách dữ liệu (mô tả)
library(GGally); library(ggplot2); library(visreg); library(table1); library(BMA)## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-## Loading required package: survival## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.4-9)ggpairs(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`.
# Task 2: Mô hình hồi qui tuyến tính đơn giản Khác biệt gữa nam và nữ
m = lm(pcfat ~ gender, data=ob)
summary(m)##
## Call:
## lm(formula = pcfat ~ gender, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.0724 -3.2724 0.1484 3.6276 14.8439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.6724 0.1826 189.9 <2e-16 ***
## genderM -10.5163 0.3381 -31.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.362 on 1215 degrees of freedom
## Multiple R-squared: 0.4432, Adjusted R-squared: 0.4428
## F-statistic: 967.3 on 1 and 1215 DF, p-value: < 2.2e-16visreg(m)
# Khác biệt gữa nam và nữ, hiệu chỉnh cho BMI
m = lm(pcfat ~ gender + bmi, data=ob)
summary(m)##
## Call:
## lm(formula = pcfat ~ gender + bmi, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.4709 -2.4780 0.1773 2.6903 15.1761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.01035 0.85880 10.49 <2e-16 ***
## genderM -11.06631 0.25599 -43.23 <2e-16 ***
## bmi 1.15303 0.03809 30.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.049 on 1214 degrees of freedom
## Multiple R-squared: 0.6828, Adjusted R-squared: 0.6822
## F-statistic: 1306 on 2 and 1214 DF, p-value: < 2.2e-16Tương quan với bmi
m1 = lm(pcfat ~ bmi, data=ob)
m2 = lm(pcfat ~ bmi + I(bmi^2), data=ob)
m3 = lm(pcfat ~ bmi + I(bmi^2) + I(bmi^3), data=ob)
anova(m1, m2, m3)## Analysis of Variance Table
##
## Model 1: pcfat ~ bmi
## Model 2: pcfat ~ bmi + I(bmi^2)
## Model 3: pcfat ~ bmi + I(bmi^2) + I(bmi^3)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1215 50541
## 2 1214 49921 1 620.14 15.1175 0.0001065 ***
## 3 1213 49758 1 162.30 3.9565 0.0469163 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Chia thành 3 cửa số và hiển thị mô hình m1, m2, m3
par(mfrow=c(1,3))
visreg(m1)
visreg(m2)
visreg(m3)
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.