BT ngay 3
Trieu Thi Vy Vy
2025-10-24
Ngày 3: Hồi qui tuyến tính
Việc 1. Phân tích phương sai
1.1. Đọc dữ liệu vào R
A = c(8, 9, 11, 4, 7, 8, 5)
B = c(7, 17, 10, 14, 12, 24, 11, 22)
C = c(28, 21, 26, 11, 24, 19)
D = c(26, 16, 13, 12, 9, 10, 11, 17, 15)
weight = c(A, B, C, D)
group = c(rep("A", 7), rep("B", 8), rep("C", 6), rep("D", 9))
data = data.frame(weight, group)
data## weight group
## 1 8 A
## 2 9 A
## 3 11 A
## 4 4 A
## 5 7 A
## 6 8 A
## 7 5 A
## 8 7 B
## 9 17 B
## 10 10 B
## 11 14 B
## 12 12 B
## 13 24 B
## 14 11 B
## 15 22 B
## 16 28 C
## 17 21 C
## 18 26 C
## 19 11 C
## 20 24 C
## 21 19 C
## 22 26 D
## 23 16 D
## 24 13 D
## 25 12 D
## 26 9 D
## 27 10 D
## 28 11 D
## 29 17 D
## 30 15 D1.2. Mô tả
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~ weight | group, data = data, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]"))| A (N=7) |
B (N=8) |
C (N=6) |
D (N=9) |
Overall (N=30) |
|
|---|---|---|---|---|---|
| weight | |||||
| Mean (SD) | 7.43 (2.37) | 14.6 (5.95) | 21.5 (6.09) | 14.3 (5.15) | 14.2 (6.75) |
| Median [Q1, Q3] | 8.00 [6.00, 8.50] | 13.0 [10.8, 18.3] | 22.5 [19.5, 25.5] | 13.0 [11.0, 16.0] | 12.0 [9.25, 18.5] |
1.3. Phân tích khác biệt về cân nặng giữa các nhóm
av = aov(weight ~ group, data = data)
summary(av)## Df Sum Sq Mean Sq F value Pr(>F)
## group 3 642.3 214.09 8.197 0.000528 ***
## Residuals 26 679.1 26.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.4. Phân tích hậu định
TukeyHSD(av)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = weight ~ group, data = data)
##
## $group
## diff lwr upr p adj
## B-A 7.1964286 -0.05969765 14.4525548 0.0525014
## C-A 14.0714286 6.27132726 21.8715299 0.0002134
## D-A 6.9047619 -0.16073856 13.9702624 0.0571911
## C-B 6.8750000 -0.69675602 14.4467560 0.0850381
## D-B -0.2916667 -7.10424368 6.5209103 0.9994049
## D-C -7.1666667 -14.55594392 0.2226106 0.0597131#Vẽ biểu đồ
plot(TukeyHSD(av), las = 2, col = "blue")
Diễn giải kết quả: Trong bốn nhóm A, B, C, D, chỉ nhóm C có trung bình khác biệt rõ rệt so với nhóm A. Các nhóm còn lại không có khác biệt đáng kể về giá trị trung bình.
2. Phân tích tương quan
2.1 Đọc dữ liệu Demo data.csv vào R
df = read.csv("F:\\ME\\R\\Demo.csv")
dim(df)## [1] 1217 7head(df)## X id age gender weight height pcfat
## 1 1 1 53 F 49 150 37.3
## 2 2 2 65 M 52 165 16.8
## 3 3 3 64 F 57 157 34.0
## 4 4 4 56 F 53 156 33.8
## 5 5 5 54 M 51 160 14.8
## 6 6 6 52 F 47 153 32.22.2. Mô tả cân nặng và chiều cao
library(table1)
table1(~ weight + height, data = df)| Overall (N=1217) |
|
|---|---|
| weight | |
| Mean (SD) | 55.1 (9.40) |
| Median [Min, Max] | 54.0 [34.0, 95.0] |
| height | |
| Mean (SD) | 157 (7.98) |
| Median [Min, Max] | 155 [136, 185] |
2.3. Vẽ biểu đồ tán xạ đánh giá mối liên quan giữa cân nặng và chiều cao
plot(height ~ weight, data = df)
library(ggplot2)
ggplot(data = df, aes(x = weight, y = height)) + geom_point() + geom_smooth()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
2.4. Đánh giá tương quan giữa cân nặng và chiều cao
cor.test(df$weight, df$height)##
## Pearson's product-moment correlation
##
## data: df$weight and df$height
## t = 25.984, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5602911 0.6326135
## sample estimates:
## cor
## 0.59766672.5. Đánh giá mối liên quan giữa chiều cao và tỉ trọng mỡ
library(lessR)##
## lessR 4.4.5 feedback: gerbing@pdx.edu
## --------------------------------------------------------------
## > d <- Read("") Read data file, many formats available, e.g., Excel
## d is default data frame, data= in analysis routines optional
##
## Many examples of reading, writing, and manipulating data,
## graphics, testing means and proportions, regression, factor analysis,
## customization, forecasting, and aggregation from pivot tables
## Enter: browseVignettes("lessR")
##
## View lessR updates, now including time series forecasting
## Enter: news(package="lessR")
##
## Interactive data analysis
## Enter: interact()##
## Attaching package: 'lessR'## The following object is masked from 'package:table1':
##
## labelPlot(height, pcfat, data = df, fit = "lm")
##
##
## >>> Suggestions or enter: style(suggest=FALSE)
## Plot(height, pcfat, enhance=TRUE) # many options
## Plot(height, pcfat, color="red") # exterior edge color of points
## Plot(height, pcfat, MD_cut=6) # Mahalanobis distance from center > 6 is an outlier
##
##
## >>> Pearson's product-moment correlation
##
## Number of paired values with neither missing, n = 1217
## Sample Correlation of height and pcfat: r = -0.480
##
## Hypothesis Test of 0 Correlation: t = -19.063, df = 1215, p-value = 0.000
## 95% Confidence Interval for Correlation: -0.522 to -0.435
##
##
## Line: b0 = 99.311633 b1 = -0.432014 Linear Model MSE = 39.747926 Rsq = 0.230
## Nhận xét kết quả: Kết quả phân tích cho thấy có mối tương quan âm có ý nghĩa thống kê giữa chiều cao và tỷ trọng mỡ cơ thể (r ≈ -0.58, p < 0.001).Điều này có nghĩa là khi chiều cao tăng, tỷ lệ mỡ cơ thể có xu hướng giảm. Mức tương quan ở mức trung bình – khá mạnh, chứng tỏ hai biến có mối quan hệ đáng kể về mặt sinh học và thống kê.
Việc 3. Hồi qui tuyến tính
3.1 Đọc dữ liệu gapminder
library(gapminder)
data(gapminder)
vn = subset(gapminder, country == "Vietnam")
library(lessR)
Plot(lifeExp, year, data = vn, xlab = "Life expectancy (years)", ylab = "Year")
##
## >>> Suggestions or enter: style(suggest=FALSE)
## Plot(lifeExp, year, enhance=TRUE) # many options
## Plot(lifeExp, year, color="red") # exterior edge color of points
## Plot(lifeExp, year, fit="lm", fit_se=c(.90,.99)) # fit line, stnd errors
## Plot(lifeExp, year, MD_cut=6) # Mahalanobis distance from center > 6 is an outlier
##
##
## >>> Pearson's product-moment correlation
##
## Number of paired values with neither missing, n = 12
## Sample Correlation of lifeExp and year: r = 0.995
##
## Hypothesis Test of 0 Correlation: t = 30.569, df = 10, p-value = 0.000
## 95% Confidence Interval for Correlation: 0.981 to 0.999
## 3.2. Thực hiện phân tích hồi qui tuyến tính để xác định xem mỗi năm, trong thời gian từ 1952-2007, tuổi thọ của người Việt Nam gia tăng như thế nào?
m.1 = lm(lifeExp ~ year, data = vn)
summary(m.1)##
## Call:
## lm(formula = lifeExp ~ year, data = vn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1884 -0.5840 0.1335 0.7396 1.7873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1271.98315 43.49240 -29.25 0.0000000000510 ***
## year 0.67162 0.02197 30.57 0.0000000000329 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared: 0.9894, Adjusted R-squared: 0.9884
## F-statistic: 934.5 on 1 and 10 DF, p-value: 0.00000000003289m.2 <- lm(lifeExp ~ year, data = vn)summary(m.2)##
## Call:
## lm(formula = lifeExp ~ year, data = vn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1884 -0.5840 0.1335 0.7396 1.7873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1271.98315 43.49240 -29.25 0.0000000000510 ***
## year 0.67162 0.02197 30.57 0.0000000000329 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.314 on 10 degrees of freedom
## Multiple R-squared: 0.9894, Adjusted R-squared: 0.9884
## F-statistic: 934.5 on 1 and 10 DF, p-value: 0.00000000003289m.2 = Regression(lifeExp ~ year, data = vn)
m.2## >>> Suggestion
## # Create an R markdown file for interpretative output with Rmd = "file_name"
## Regression(my_formula=lifeExp ~ year, data=vn, Rmd="eg")
##
##
## BACKGROUND
##
## Data Frame: vn
##
## Response Variable: lifeExp
## Predictor Variable: year
##
## Number of cases (rows) of data: 12
## Number of cases retained for analysis: 12
##
##
## BASIC ANALYSIS
##
## Estimate Std Err t-value p-value Lower 95% Upper 95%
## (Intercept) -1271.9832 43.4924 -29.246 0.000 -1368.8903 -1175.0760
## year 0.6716 0.0220 30.569 0.000 0.6227 0.7206
##
## Standard deviation of lifeExp: 12.17233
##
## Standard deviation of residuals: 1.31365 for df=10
## 95% range of residuals: 5.85399 = 2 * (2.228 * 1.31365)
##
## R-squared: 0.989 Adjusted R-squared: 0.988 PRESS R-squared: 0.984
##
## Null hypothesis of all 0 population slope coefficients:
## F-statistic: 934.455 df: 1 and 10 p-value: 0.000
##
## -- Analysis of Variance
##
## df Sum Sq Mean Sq F-value p-value
## Model 1 1612.5653 1612.5653 934.4554 0.000
## Residuals 10 17.2567 1.7257
## lifeExp 11 1629.8221 148.1656
##
##
## K-FOLD CROSS-VALIDATION
##
##
## RELATIONS AMONG THE VARIABLES
##
## lifeExp year
## lifeExp 1.00 0.99
## year 0.99 1.00
##
##
## RESIDUALS AND INFLUENCE
##
## -- Data, Fitted, Residual, Studentized Residual, Dffits, Cook's Distance
## [sorted by Cook's Distance]
## [n_res_rows = 12, out of 12 ]
## -------------------------------------------------------------
## year lifeExp fitted resid rstdnt dffits cooks
## 12 2007 74.2490 75.9489 -1.6999 -1.6742 -1.0827 0.4965
## 1 1952 40.4120 39.0101 1.4019 1.3167 0.8515 0.3377
## 5 1972 50.2540 52.4424 -2.1884 -2.0016 -0.6637 0.1694
## 9 1992 67.6620 65.8747 1.7873 1.5563 0.5937 0.1543
## 10 1997 70.6720 69.2328 1.4392 1.2327 0.5559 0.1469
## 4 1967 47.8380 49.0843 -1.2463 -1.0172 -0.3880 0.0750
## 2 1957 42.8870 42.3682 0.5188 0.4300 0.2316 0.0292
## 11 2002 73.0170 72.5908 0.4262 0.3520 0.1896 0.0197
## 3 1962 45.3630 45.7262 -0.3632 -0.2891 -0.1304 0.0094
## 7 1982 58.8160 59.1585 -0.3425 -0.2596 -0.0792 0.0035
## 8 1987 62.8200 62.5166 0.3034 0.2315 0.0768 0.0032
## 6 1977 55.7640 55.8005 -0.0365 -0.0275 -0.0084 0.0000
##
##
## PREDICTION ERROR
##
## -- Data, Predicted, Standard Error of Prediction, 95% Prediction Intervals
## [sorted by lower bound of prediction interval]
## ----------------------------------------------
##
## year lifeExp pred s_pred pi.lwr pi.upr width
## 1 1952 40.4120 39.0101 1.4948 35.6794 42.3408 6.6614
## 2 1957 42.8870 42.3682 1.4539 39.1286 45.6077 6.4790
## 3 1962 45.3630 45.7262 1.4203 42.5616 48.8909 6.3293
## 4 1967 47.8380 49.0843 1.3946 45.9770 52.1917 6.2147
## 5 1972 50.2540 52.4424 1.3772 49.3738 55.5109 6.1371
## 6 1977 55.7640 55.8005 1.3684 52.7515 58.8494 6.0979
## 7 1982 58.8160 59.1585 1.3684 56.1096 62.2075 6.0979
## 8 1987 62.8200 62.5166 1.3772 59.4481 65.5852 6.1371
## 9 1992 67.6620 65.8747 1.3946 62.7673 68.9820 6.2147
## 10 1997 70.6720 69.2328 1.4203 66.0681 72.3974 6.3293
## 11 2002 73.0170 72.5908 1.4539 69.3513 75.8304 6.4790
## 12 2007 74.2490 75.9489 1.4948 72.6182 79.2796 6.6614
##
## ----------------------------------
## Plot 1: Distribution of Residuals
## Plot 2: Residuals vs Fitted Values
## ----------------------------------3.3. Kiểm tra các giả định của mô hình hồi quy tuyến tính
library(ggfortify)
autoplot(m.1)
3.5. Sử dụng chatGPT viết code để đánh giá gia tăng tuổi thọ người Việt Nam 1952-2007
year <- c(1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, 2002, 2007)
lifeExp <- c(40.4, 42.9, 45.4, 47.8, 50.3, 55.8, 58.8, 62.8, 67.7, 70.7, 73.0, 74.2)
data <- data.frame(year, lifeExp)print(data)## year lifeExp
## 1 1952 40.4
## 2 1957 42.9
## 3 1962 45.4
## 4 1967 47.8
## 5 1972 50.3
## 6 1977 55.8
## 7 1982 58.8
## 8 1987 62.8
## 9 1992 67.7
## 10 1997 70.7
## 11 2002 73.0
## 12 2007 74.2summary(data$lifeExp)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 40.40 47.20 57.30 57.48 68.45 74.20plot(year, lifeExp,
type = "b", pch = 19, col = "blue",
main = "Gia tăng tuổi thọ theo thời gian (1952–2007)",
xlab = "Năm", ylab = "Tuổi thọ trung bình (năm)")
model <- lm(lifeExp ~ year, data = data)
summary(model)##
## Call:
## lm(formula = lifeExp ~ year, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1494 -0.5944 0.1387 0.7324 1.8268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1271.13492 43.77173 -29.04 0.0000000000547 ***
## year 0.67119 0.02211 30.35 0.0000000000353 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.322 on 10 degrees of freedom
## Multiple R-squared: 0.9893, Adjusted R-squared: 0.9882
## F-statistic: 921.4 on 1 and 10 DF, p-value: 0.00000000003527model <- lm(lifeExp ~ year, data = df)plot(year, lifeExp,
type = "b", pch = 19, col = "blue",
main = "Gia tăng tuổi thọ theo thời gian (1952–2007)",
xlab = "Năm", ylab = "Tuổi thọ trung bình (năm)")
abline(model, col = "red", lwd = 2)
Gia tăng tuổi thọ trung bình mỗi năm
coef(model)## (Intercept) year
## -1271.1349184 0.6711888increase_per_year <- coef(model)[2]
cat("Tuổi thọ tăng trung bình mỗi năm:", round(increase_per_year, 3), "năm/năm\n")## Tuổi thọ tăng trung bình mỗi năm: 0.671 năm/nămDự báo tuổi thọ năm 2012
predict(model, newdata = data.frame(year = 2012))## 1
## 79.29697Nhận xét kết quả: Kết quả hồi quy cho thấy tuổi thọ trung bình của người dân tăng ổn định và tuyến tính theo thời gian, với tốc độ tăng khoảng 0.67 năm mỗi năm trong giai đoạn 1952–2007.Mô hình có khả năng dự báo tốt, và dự đoán rằng đến năm 2012, tuổi thọ trung bình đạt khoảng 79,3 năm.Điều này phản ánh sự cải thiện đáng kể trong điều kiện sống, chăm sóc y tế và kinh tế – xã hội qua các thập kỷ.
Việc 4. Hồi quy tuyến tính đa biến
4.1. Nhập dữ liệu vào R
Y = c(12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6)
X1 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
X2 = c(7, 4, 4, 6, 4, 2, 1, 1, 1, 0)
df = data.frame(Y, X1, X2)
dim(df)## [1] 10 3head(df)## Y X1 X2
## 1 12.1 0 7
## 2 11.9 1 4
## 3 10.2 2 4
## 4 8.0 3 6
## 5 7.7 4 4
## 6 5.3 5 24.2 Đánh giá mối liên quan giữa ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1).
m.3 = lm(Y ~ X1, data = df)
summary(m.3)##
## Call:
## lm(formula = Y ~ X1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1606 -1.0735 0.1742 0.8621 2.0970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.8545 0.8283 14.312 0.000000554 ***
## X1 -0.8788 0.1552 -5.664 0.000474 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.409 on 8 degrees of freedom
## Multiple R-squared: 0.8004, Adjusted R-squared: 0.7755
## F-statistic: 32.08 on 1 and 8 DF, p-value: 0.00047374.3 Đánh giá mối liên quan giữa ethylene oxide (Y) và thời gian lưu (X2)
m.4 = lm(Y ~ X2, data = df)
summary(m.4)##
## Call:
## lm(formula = Y ~ X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.702 -1.533 -0.034 1.667 3.066
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0980 1.1222 4.543 0.00189 **
## X2 0.9340 0.2999 3.114 0.01436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.121 on 8 degrees of freedom
## Multiple R-squared: 0.548, Adjusted R-squared: 0.4915
## F-statistic: 9.698 on 1 and 8 DF, p-value: 0.014364.4. Đánh giá mối liên quan độc lập giữa thời gian lưu (X2) và ethylene oxide (Y)
m.5 = lm(Y ~ X1 + X2, data = df)
summary(m.5)##
## Call:
## lm(formula = Y ~ X1 + X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.46078 -0.33384 0.00026 0.81856 1.98476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7076 2.9785 4.938 0.00168 **
## X1 -1.2042 0.3614 -3.332 0.01255 *
## X2 -0.4629 0.4642 -0.997 0.35187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.41 on 7 degrees of freedom
## Multiple R-squared: 0.8252, Adjusted R-squared: 0.7753
## F-statistic: 16.53 on 2 and 7 DF, p-value: 0.002232