#Việc 1. Phân tích phương sai ##1.1.Nhập dữ liệu vào R và đặt tên dữ liệu là data

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 = factor(rep(c("A", "B", "C", "D"), times = c(length(A), length(B), length(C), length(D))))
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     D

##1.2 Mô tả cân nặng giữa 4 nhóm.

aggregate(weight ~ group, data = data,function(x) {c(Min = min(x),Max = max(x),Mean = mean(x),Median = median(x),SD = sd(x))})
##   group weight.Min weight.Max weight.Mean weight.Median weight.SD
## 1     A   4.000000  11.000000    7.428571      8.000000  2.370453
## 2     B   7.000000  24.000000   14.625000     13.000000  5.950690
## 3     C  11.000000  28.000000   21.500000     22.500000  6.090977
## 4     D   9.000000  26.000000   14.333333     13.000000  5.147815
print(aggregate(weight ~ group, data = data, function(x) {
  c(Min = min(x), Max = max(x), Mean = mean(x), Median = median(x), SD = sd(x))
}))
##   group weight.Min weight.Max weight.Mean weight.Median weight.SD
## 1     A   4.000000  11.000000    7.428571      8.000000  2.370453
## 2     B   7.000000  24.000000   14.625000     13.000000  5.950690
## 3     C  11.000000  28.000000   21.500000     22.500000  6.090977
## 4     D   9.000000  26.000000   14.333333     13.000000  5.147815
summary_stats <- aggregate(weight ~ group, data = data, function(x) {
  c(Min = min(x),
    Max = max(x),
    Mean = mean(x),
    Median = median(x),
    SD = sd(x))
})
print(summary_stats)
##   group weight.Min weight.Max weight.Mean weight.Median weight.SD
## 1     A   4.000000  11.000000    7.428571      8.000000  2.370453
## 2     B   7.000000  24.000000   14.625000     13.000000  5.950690
## 3     C  11.000000  28.000000   21.500000     22.500000  6.090977
## 4     D   9.000000  26.000000   14.333333     13.000000  5.147815
boxplot(weight ~ group, data = data,main = "So sánh cân nặng giữa 4 nhóm",xlab = "Nhóm",ylab = "Cân nặng",col = c("lightblue", "lightgreen", "lightpink", "lightyellow"))

##1.3 Phân tích sự khác biệt về cân nặng giữa 4 nhóm. Diễn giải kết quả.

av = aov(weight ~ 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

Kết luận: p < 0.05 => có sự khác biệt giữa các nhóm

##1.4 Thực hiện phân tích hậu định (post-hoc analysis) để xác định cụ thể nhóm có khác biệt về cân nặng. Diễn giải kết quả.

 TukeyHSD(av)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ 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
plot(TukeyHSD(av), las=2 , col="red")

Kết luận: Giá trị p ở kiểm định Turkey cho thấy sự khác biệt giữa B-A, C-A, D-C. Biểu đồ cũng cho thấy khoảng tin cậy về lênh lệch giữa các nhóm B-A, C-A, D-C gần như không chứa 0 => chênh lệch khác 0 => có khác biệt.

#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”

df=read.csv("F:\\VIEC_KHOA\\SU_KIEN\\2025_2026\\Taphuanphantichdulieu\\Datamophong\\Demo.csv")

##2.2 Mô tả đặc điểm cân nặng (weight) và chiều cao (height).

library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
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 (weight) và chiều cao (height). Nhận xét kết quả.

library(ggplot2)
ggplot(data=df, aes(x=height, y=weight)) + geom_point() + labs(x="Chiều cao", 
y="Cân nặng")

Kết luận: Chiều cao và cân nặng dường như có mối quan hệ, chiều cao tăng thì cân nặng tăng. ##2.4 Tiến hành phân tích tương quan định lượng mối liên quan giữa cân nặng (weight) và chiều cao (height). Nhận xét kết quả.

cor.test(df$height, df$weight)
## 
##  Pearson's product-moment correlation
## 
## data:  df$height and df$weight
## 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.5976667

Kết luận: p < 0.05 => có mối quan hệ với hệ số tương quan r = 0.5977. r lớn hơn 0 nên tương quan thuận chiều. ##2.5 Tiến hành phân tích tương quan định lượng mối liên quan giữa chiều cao (height) và tỉ trọng mỡ (pcfat). Nhận xét kết quả.

ggplot(data=df, aes(x=height, y=pcfat)) + geom_point() + labs(x="Chiều cao", 
y="Tỷ trọng mỡ")

cor.test(df$height, df$pcfat)
## 
##  Pearson's product-moment correlation
## 
## data:  df$height and df$pcfat
## t = -19.063, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5219407 -0.4353664
## sample estimates:
##        cor 
## -0.4798206

Kết luận: Kết luận: p < 0.05 => có mối quan hệ với hệ số tương quan r =-0.4797. r nhỏ hơn 0 nên tương quan ngược chiều Người cao thì tỷ lệ mỡ thấp.

#Việc 3. Hồi qui tuyến tính ##3.1 Đọc dữ liệu “gapminder” vào R từ gói lệnh gapminder. Tạo bộ dữ liệu “vn” gồm dữ liệu của Việt Nam. gapminder chứa dữ liệu về các quốc gia theo thời gian (GDP, tuổi thọ, dân số).

library(gapminder)
data(gapminder)
 vn = subset(gapminder, country == "Vietnam")
 head(vn)
## # A tibble: 6 × 6
##   country continent  year lifeExp      pop gdpPercap
##   <fct>   <fct>     <int>   <dbl>    <int>     <dbl>
## 1 Vietnam Asia       1952    40.4 26246839      605.
## 2 Vietnam Asia       1957    42.9 28998543      676.
## 3 Vietnam Asia       1962    45.4 33796140      772.
## 4 Vietnam Asia       1967    47.8 39463910      637.
## 5 Vietnam Asia       1972    50.3 44655014      700.
## 6 Vietnam Asia       1977    55.8 50533506      714.
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':
## 
##     label
plot(lifeExp~year,data=vn)

ggplot(data = vn, aes(x = year, y = lifeExp)) + geom_point(color = "blue") + labs(x = "Năm", y = "Tuổi thọ", title = "Tuổi thọ Việt Nam theo thời gian")

##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? cách 1: dùng hàm lm hàm summary dùng để hiển thị bảng tóm tắt chi tiết của mô hình hồi quy m

 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) -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.00000000003289

Kết luận: Trong giai đoạn 1952-2007, mỗi năm tuổi thọ tăng 0.672 Cách 2: dùng lessR

library(lessR)
 fit = 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 
## ----------------------------------

Vẽ đồ thị hồi quy

library(ggplot2)
ggplot(vn, aes(x = year, y = lifeExp)) +
 geom_point(size = 3, color = "steelblue") +
 geom_smooth(method = "lm", se = TRUE, color = "darkred") +
 theme_minimal(base_size = 14) +
 labs(title = "Linear Regression: lifeExp ~ year",
 x = "Năm",
 y = "Tuổi thọ (lifeExp)")
## `geom_smooth()` using formula = 'y ~ x'

##3.3 Kiểm tra các giả định của mô hình hồi qui tuyến tính

m = lm(lifeExp ~ year, data=vn)
library(ggfortify)
autoplot(m)

cách 2: dùng plot(m) par(mfrow = c(2, 2)) dùng để chia các hình vào các cột và hàng (ở đây là 2 cột, 2 hàng)

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

##3.4 Viết phương trình đánh giá gia tăng tuổi thọ của người Việt Nam trong khoảng 1952-2007. Nhận xét kết quả. lifeExp=−1272+0.6716×year

#Việc 4. Hồi qui tuyến tính đa biến ##4.1 Nhập dữ liệu vào R và đặt tên tập dữ liệu là df

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 ethylene oxide (Y) và hoạt tính bạc của xúc tác (X1).

m1 = lm(Y ~ X1, data=df)
summary(m1)
## 
## 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
plot(Y~X1,data=df)

ggplot(df, aes(x = X1, y = Y)) +
 geom_point(size = 3, color = "steelblue") +
 geom_smooth(method = "lm", se = TRUE, color = "darkred") +
 theme_minimal(base_size = 14) +
 labs(title = "Linear Regression: Y ~ X1",
 x = "hoạt tính bạc của xúc tác",
 y = "ethylene oxide")
## `geom_smooth()` using formula = 'y ~ x'

m2 = lm(Y ~ X2, data=df)
summary(m2)
## 
## 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
plot(Y~X2,data=df)

ggplot(df, aes(x = X2, y = Y)) +
 geom_point(size = 3, color = "steelblue") +
 geom_smooth(method = "lm", se = TRUE, color = "darkred") +
 theme_minimal(base_size = 14) +
 labs(title = "Linear Regression: Y ~ X2",
 x = "thời gian lưu",
 y = "ethylene oxide")
## `geom_smooth()` using formula = 'y ~ x'

##4.4 Đánh giá mối liên quan độc lập giữa thời gian lưu (X2) và ethylene oxide (Y). Nhận xét kết quả. Xem X2 có ảnh hưởng đến Y hay không, độc lập với X1 (không xét X1 trong mô hình này).