NGÀY 3 (23/10/2025)
DUC CHIEN NGO
2025-10-22
# Ngày 3 (23/10/2025)
# Việc 1 Phân tích phương sai
# 1.1 Tạo dữ liệu
# Nhập dữ liệu 4 nhóm cân nặng
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)
wt = c(A, B, C, D)
group = c(rep("A", 7), rep("B", 8), rep("C",6),rep("D",9))
data = data.frame(wt, group)
dim(data)## [1] 30 2# 1.2 Mô tả cân nặng giữa 4 nhóm
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~ wt | 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) |
|
|---|---|---|---|---|---|
| wt | |||||
| 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 sự khác biệt về cân nặng giữa 4 nhóm
av=aov(wt ~ group)
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 ' ' 1# 1.4 Thực hiện phân tích hậu định (post-hoc analysis)
TukeyHSD(av)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = wt ~ group)
##
## $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# Việc 2 Phân tích tương quan
# 2.1 Đọc dữ liệu "Demo data.csv" vào R và gọi dữ liệu là "df"
file.choose()## [1] "D:\\TAP HUAN R DAU\\Demo.csv"df = read.csv("D:\\TAP HUAN R DAU\\Demo.csv")
# 2.2 Mô tả đặc điểm cân nặng (weight) và chiều cao (height)
library(table1)
table1(~height+weight, data=df)| Overall (N=1217) |
|
|---|---|
| height | |
| Mean (SD) | 157 (7.98) |
| Median [Min, Max] | 155 [136, 185] |
| weight | |
| Mean (SD) | 55.1 (9.40) |
| Median [Min, Max] | 54.0 [34.0, 95.0] |
# 2.3 Vẽ biểu đồ tán xạ
library(ggplot2)
ggplot(data=df, aes(x=height, y=weight)) + geom_point() + labs(x="Chiều cao", y="Cân nặng")
# 2.4 Tiến hành phân tích tương quan định lượng
shapiro.test(df$height)##
## Shapiro-Wilk normality test
##
## data: df$height
## W = 0.9827, p-value = 7.179e-11shapiro.test(df$weight)##
## Shapiro-Wilk normality test
##
## data: df$weight
## W = 0.974, p-value = 5.471e-14cor.test(df$height, df$weight, method="spearman")## Warning in cor.test.default(df$height, df$weight, method = "spearman"): Cannot
## compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: df$height and df$weight
## S = 132686871, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.55832# 2.5 Tiến hành phân tích tương quan định lượng
shapiro.test(df$height)##
## Shapiro-Wilk normality test
##
## data: df$height
## W = 0.9827, p-value = 7.179e-11shapiro.test(df$pcfat)##
## Shapiro-Wilk normality test
##
## data: df$pcfat
## W = 0.98284, p-value = 8.193e-11cor.test(df$height, df$pcfat, method="spearman")## Warning in cor.test.default(df$height, df$pcfat, method = "spearman"): Cannot
## compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: df$height and df$pcfat
## S = 441808157, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4706643# Việc 3 Hồi quy tuyến tính
# 3.1 Đọc dữ liệu
library(gapminder)
data(gapminder)
vn = subset(gapminder, country == "Vietnam")
plot(vn$lifeExp ~ vn$year, pch=16, col="blue", xlab="Year",
ylab="Life Expextancy")
# 3.2 Thực hiện phân tích hồi qui tuyến tính
m = lm(lifeExp ~ year, data=vn)
summary(m)##
## 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) -1.272e+03 4.349e+01 -29.25 5.10e-11 ***
## year 6.716e-01 2.197e-02 30.57 3.29e-11 ***
## ---
## 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: 3.289e-11# Hoặc làm qua lessR
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':
##
## labelfit = reg(lifeExp ~ year, data=vn)
fit## >>> Suggestion
## # Create an R markdown file for interpretative output with Rmd = "file_name"
## reg(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 qui tuyến tính
library(ggfortify)
autoplot(m)
# 3.4 Viết phương trình đánh giá
# Phương trình: Life Exp = -1272 + 0.67*year
# 3.5 Sử dụng ChatGPT
# Nhập dữ liệu năm và tuổi thọ
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)
# Tạo khung dữ liệu
data <- data.frame(year, lifeExp)
# Xem dữ liệu
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.2# Vẽ đồ thị thể hiện xu hướng gia tăng tuổi thọ
plot(data$year, data$lifeExp,
type = "b", # vẽ đường nối giữa các điểm
col = "blue", # màu đường
pch = 19, # kiểu điểm tròn đặc
xlab = "Năm",
ylab = "Tuổi thọ trung bình (năm)",
main = "Gia tăng tuổi thọ của người Việt Nam (1952 - 2007)")
# Thêm đường xu hướng tuyến tính
model <- lm(lifeExp ~ year, data = data)
abline(model, col = "red", lwd = 2)
# Hiển thị kết quả mô hình hồi quy
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.00000000003527# Kiểm định sau mô hình
par(mfrow = c(2,2))
plot(model)
shapiro.test(residuals(model))##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.95642, p-value = 0.7317library(car)## Loading required package: carData##
## Attaching package: 'car'## The following objects are masked from 'package:lessR':
##
## bc, recode, spncvTest(model)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.1795134, Df = 1, p = 0.67179durbinWatsonTest(model)## lag Autocorrelation D-W Statistic p-value
## 1 0.386577 0.9453653 0.024
## Alternative hypothesis: rho != 0# Việc 4 Hồi quy tuyến tính đa biến
# 4.1 Nhập liệu
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)
head(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 2# 4.2 Đánh giá mối liên quan giữa Y và X1
model1 = lm(Y ~ X1, data = df)
summary(model1)##
## 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.0004737# 4.3 Đánh giá mối liên quan giữa Y và X2
model2 = lm(Y ~ X2, data = df)
summary(model2)##
## 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.01436# 4.4 Bạn muốn đánh giá mối liên quan độc lập giữa X2 và Y
model3 = lm(Y ~ X1 + X2, data = df)
summary(model3)##
## 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