Nhắc lại mô hình hồi quy đơn giản

    Ở nội dung này, chúng ta cùng nhắc lại các kiến thức cơ bản của một mô mình hồi quy tuyến tính đơn giản và cách đọc các kết quả từ mô hình hồi quy thông qua hàm lm(). Chúng ta sử dụng một bộ dữ liệu khác để thực hành lại mô hình hồi quy đơn giản, bộ dữ liệu thai kỳ (Gestation). Bộ dữ liệu trình bày mối tương quan giữa trọng lượng của thai nhi (theo tuần tuổi) và thời gian mang thai của sản phụ.

Francis Galton

    Trước khi đến với bài thực hành chúng ta cần nhắc lại lịch sử của hồi quy (regression). Thuật ngữ này bắt nguồn từ nhà khoa học Francis Galton (anh họ của Charles Darwin) vào thế kỷ 19 khi nhà khoa học này muốn tìm một phương pháp để miêu tả sự liên quan của chiều cao con cái với chiều cao của cha mẹ. Ông nhận thấy những người có cha mẹ có chiều cao lớn thì chiều cao cũng lớn và gần với giá trị trung bình, những người có chiều cao thấp thì chiều cao của cha mẹ cũng thấp nhưng gần với giá trị trung bình. Hiện tượng này được Galton gọi là “Regression toward the mediocrity”. Thuật ngữ regression được Galton lần đầu sử dụng và phổ biến cho tới ngày nay.

Dữ liệu về thai kỳ

    Để xem xét sự liên quan của tuần thai kỳ và trọng lượng. Chúng ta cùng xem xét tương quan giữa 2 biến này.

gestation <- c(34.7, 36, 29.3, 40.1, 35.7, 42.4, 40.3, 37.3, 40.9, 38.3, 38.5, 41.4, 39.7, 39.7, 41.1, 38, 38.7)
weight <- c(1895, 2030, 1040, 2835, 3090, 3827, 3260, 2690, 3285, 2920, 3430, 3657, 3685, 3345, 3260, 2680, 2005)
thaiky <- data.frame(gestation, weight)
corrplot(cor(thaiky), type = "upper", method = "number")

thaiky %>% 
  ggplot(aes(x = gestation, y = weight)) +
  geom_point(color = "red") + 
  geom_smooth(method = "lm") + 
  ggtitle("The correlation between Gestation and Weight") +
  xlab("Gestation") +
  ylab("Weight") +
  theme_economist()
## `geom_smooth()` using formula 'y ~ x'

Dạng phương trình hồi quy

Dành cho tổng thể

    Để biểu thị mô hình hồi quy tuyến tính cho tổng thể, chúng ta có thể viết như sau: \[y_{i} = \alpha+\beta x_{i}\]     Trong đó tham số \(\alpha\) là hệ số góc của phương trình hồi quy và \(\beta\) là hệ số góc (slope hay gradient) của phương trình hồi quy. Tuy nhiên phương trình này không thể nối tất cả các điểm \((x_{i}, y_{i})\) được và do đó chúng ta cần cộng thêm một giá trị bị lệch khỏi mô hình (giá trị nhiễu, giá trị phần dư của mô hình) \[y_{i} = \alpha+\beta x_{i}+\epsilon_{i}\]     Bỏ đi giá trị i nhỏ thì phương trình hồi quy cho tổng thể bây giờ là: \[y = \alpha+\beta x+\epsilon\]

Dành cho mẫu

    Phía trên là mô hình cho quần thể, đối với mẫu mà chúng ta xử lý được thì sử dụng phương trình sau: \[y=a+bx+e\]     Trong đó a và b là hai ước số của \(\beta\), phần dư \(e\) của mô hình là phần nhiễu không giải thích được. Tóm lại mô hình hồi quy được viết như sau:

Giá trị quan sát y = Giá trị dự đoán y + Phần dư(nhiễu) \[y = \hat{y}+e\]     Giá trị slope chỉ ra tương quan giữa y và x. Phản ánh tốc độ thay đổi của y theo sự thay đổi của x. Phương pháp tối thiểu hóa bình phương phần dư là phương pháp được sử dụng nhiều nhất để ước tính hệ số góc slope cho phương trình hồi quy. \[min\sum(y-\alpha-\beta x)^{2}\] \[\sum_{i}e_{i}=0\]     Bỏ qua phần nội dung tính toán. Người ta chứng minh được rằng, hệ số góc: \[b=\frac{\sum(x_{i}-\overline x)(y_{i}-\overline y)}{\sum(x_{i}-\overline x)^2}\]

    Hệ số chặn: \[a=\overline y -bx\]     Phương sai của phần dư: \[s^2=\frac{\sum(y_{i}-\hat y)^2}{n-2}\]

Ước tính tham số bằng R:

    Trong R, sử dụng hàm lm để ước tính tham số bằng phương pháp bình phương tối thiểu. Cú pháp chung của hàm này là: \[lm(y \sim x,data)\]     Trong đó y là biến phụ thuộc, x là biến độc lập. Đối với bài toán về thời gian thai kỳ, chúng ta có: \[Weight=\alpha+\beta *gestation+\epsilon\]

mymodel <- lm(weight ~ gestation, data = thaiky)
summary(mymodel)
## 
## Call:
## lm(formula = weight ~ gestation, data = thaiky)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -942.79 -175.44    2.14  160.61  751.76 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4915.45    1245.79  -3.946  0.00129 ** 
## gestation     203.18      32.37   6.276 1.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 408.8 on 15 degrees of freedom
## Multiple R-squared:  0.7242, Adjusted R-squared:  0.7058 
## F-statistic: 39.39 on 1 and 15 DF,  p-value: 1.485e-05
plot(weight ~ gestation, pch = 16, col = "blue", data = thaiky)
abline(lm(weight ~ gestation), col = "red")

    Phương trình hồi quy chúng ta có là: \(Weight=-4915.5+203.2*gestation\). Chúng ta để ý thấy rằng khi thay số thai kỳ là 40 và 41 vào thì độ lệch của weight là : 203 gram. Đây chính là giá trị slope. Do đó mô hình này có ý nghĩa như sau: Khi thai kỳ tăng lên 1 tuần thì trọng lượng trẻ sơ sinh tăng trung bình 203 gram

Một số giả định của hồi quy tuyến tính

    Trong mô hình hồi quy tuyến tính cũng như trong các phương pháp phân tích tương quan, luôn đặt ra một số giả định như sau:

  1. Mối liên quan giữa hai biến \(x\) và \(y\) là phải tuyến tính. Tuyến tính hiểu là tham số tuyến tính, ví dụ \(bx\) chứ không phải là \(b^{x}\).

2.Phần dư (nhiễu) của mô hình hồi quy phải tuân theo luật phân phối chuẩn với giá trị trung bình xấp xỉ 0 và phương sai \(s^{2}\). Viết ngắn gọn \(e\sim N(0, s^{2})\).

  1. Phương sai của phần dư không thay đổi theo giá trị của biến y (biến được dự đoán).

  2. Các giá trị của \(y_{i}\) và \(y_{i-1}\) là độc lập với nhau. Hiện tượng này khi quan sát dữ liệu các giá trị của lần thứ i và i-1 có liên quan đến nhau, gây ra một lỗi của mô hình gọi là hiện tượng tự tương quan (autocorrelation).

Kiểm tra các giả định trên bằng R

Để có thể kiểm tra các giả định trên của mô hình hồi quy tuyến tính. Sử dụng hàm autoplot() trong thư viện ggfortify:

autoplot(mymodel, prob.colour = "red")

  1. Ở biểu đồ thứ 1 (bên trái dòng 1) trình bày mối liên hệ giữa giá trị dự đoán (Fitted values) và giá trị phần dư (residuals). Các giá trị phần dư đa số đều xoay quanh giá trị 0, như vậy giả định về giá trị trung bình của phần dư là 0 thỏa mãn

  2. Ở biểu đồ thứ 2 (bên phải dòng 1) trình bày về sự phân phối của các giá trị phần dư chuẩn hóa, đa số các điểm của phần dư chuẩn hóa đều nằm trên đường chéo lý thuyết, chứng tỏ giá trị phần dư tuân theo luật phân phối chuẩn thỏa mãn.

  3. Ở biểu đồ thứ 3 (bên trái dòng 2) trình bày mối quan hệ giữa giá trị căn bậc hai của phần dư chuẩn hóa và giá trị dự đoán, chúng ta thấy các giá trị này không có liên quan đến nhau, do đó giả thuyết phương sai của y không liên quan đến giá trị x thỏa mãn

  4. Biểu đồ thứ 4 (bên phải dòng 2) trình bày mối liên hệ về các giá trị ngoại lai ảnh hưởng tới mô hình. Đa số các giá trị đều nằm trong phạm vi cho phép [-2;2], rất ít giá trị nằm ngoài khoảng này.

    Ở đây chúng ta quan sát trong bảng tóm tắt thống kê có một số các nội dung. Đó là 15 bậc tự do, vậy bậc tự do là gì?

Khái niệm bậc tự do

    Bậc tự do (Degrees Of Freedom) là con số có số lượng bằng n-1 khi sử dụng mẫu có quy mô n để ước lượng số bình quân và độ chênh lệch tiêu chuẩn của một tổng thể. Nó đề cập đến số lượng các giá trị độc lập tối đa của một hệ, là các giá trị có thể thay đổi tự do trong mẫu dữ liệu. Chúng ta có công thức tính bậc tự do như sau: \[D_{f}=N-1\]. Trong đó \(D_{f}\) là bậc tự do; \(N\) là kích thước của mẫu. Bậc tự do thường được sử dụng trong nhiều phương thức kiểm định giả thuyết trong thống kê, chẳng hạn như kiểm định Chi bình phương (đôi khi đọc là Khi bình phương). Khi sử dụng kiểm định Chi bình phương và kiểm tra tính hợp lệ của giả thuyết không, việc tính toán bậc tự do là rất quan trọng.

    Có hai dạng kiểm định Chi bình phương là:

  • Kiểm định tính độc lập đặt câu hỏi về các mối quan hệ của các biến, chẳng hạn như “Có mối quan hệ nào giữa giới tính thí sinh và điểm toán kì thi THPT 2019 không?”.

  • Kiểm định mức độ phù hợp với câu hỏi có dạng như “Nếu một đồng xu được tung 100 lần, liệu mặt ngửa sẽ xuất hiện 50 lần và mặt sấp sẽ xuất hiện 50 lần?”.

    Đối với các kiểm định này, bậc tự do được sử dụng để xác định liệu có giả thuyết không nào đó có thể bị bác bỏ không dựa trên tổng số biến và các mẫu trong kiểm định.

    Ví dụ, xem xét việc sinh viên lựa chọn khóa học, với 30 hoặc 40 sinh viên cỡ mẫu có thể không đủ lớn để tạo ra kết quả đáng tin cậy. Kết quả sử dụng cỡ mẫu 400 hoặc 500 sinh viên sẽ cho kết quả chính xác hơn.

Ví dụ về bậc tự do:

    Hãy xem xét một mẫu dữ liệu đơn giản bao gồm năm số nguyên dương. Các giá trị có thể là bất kì số nào sao cho không có mối quan hệ nào giữa chúng. Về mặt lí thuyết, mẫu dữ liệu này sẽ có năm bậc tự do.

    Giả sử biết bốn số trong mẫu là {3, 8, 5, 4} và giá trị trung bình mẫu dữ liệu là 6.

    Điều này có nghĩa là số thứ năm phải bằng 10 mà không thể là số nào khác, nó được gọi là không thể tự do thay đổi. Vì vậy, bậc tự do cho mẫu dữ liệu này là 4.

Mô hình hồi quy tuyến tính với biến phân nhóm

Biến dự đoán là biến nhị phân

Ở trường hợp này. Để bắt đầu, chúng ta sử dụng bộ dữ liệu nghiên cứu về lương (Salary) của các giáo sư được khảo sát ở Mỹ. Nghiên cứu được thực hiện trên 63 giáo sư, mục tiêu chính là tìm ra các yếu tố ảnh hưởng tới lương của giáo sư mỗi năm (USD). Các biến thành phần bao gồm:

  • TimeSincePhD: thời gian (năm) kể từ khi nhận bằng tiến sỹ
  • NPubs: số bài báo khoa học đã công bố
  • Sex: giới tính(nam = 0, nữ = 1)
  • Citations: Số lần trích dẫn các bài báo đã công bố
  • Salary: lương (USD)

Câu hỏi: Lương giáo sư có liên quan gì đến giới tính hay không, nói cách khác lương có bị ảnh hưởng bởi giới tính hay không? Phương pháp đầu tiên chúng ta nghĩ tới chính là kiểm định t-test dành cho 2 biến phân loại. Trước tiên chúng ta sẽ mô phỏng dữ liệu dành cho nam/nữ và sử dụng kiểm định t-test để kiểm tra kết quả.

datasalary <- read.csv("D:/Short Course Using R/Lesson 4/Hoi quy tuyen tinh/Income and PhDs.csv", header = TRUE, stringsAsFactors = FALSE)
# Coding
datasalary$Sex <- as.factor(datasalary$Sex)
ggplot(datasalary, aes(x = Sex, y = Salary, fill = Sex, col = Sex)) +
  geom_boxplot(alpha = 0.5)

Nhìn vào biểu đồ hộp box-plot chúng ta thấy dường như lương của nữ giáo sư thấp hơn so với lương của nam giáo sư. Nhưng chúng ta không rõ rằng sự khác biệt này có ý nghĩa thống kê hay không, hay đơn giản đây chỉ là một kết quả ngẫu nhiên. Sử dụng kiểm định t-test:

t.test(Salary ~ Sex, data = datasalary)
## 
##  Welch Two Sample t-test
## 
## data:  Salary by Sex
## t = 1.6548, df = 59.755, p-value = 0.1032
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -814.0454 8607.9480
## sample estimates:
## mean in group 0 mean in group 1 
##        56509.91        52612.96

Dựa vào kết quả thống kê trên, chúng ta thấy giá trị p-value = 0.1032 lớn hơn mức ý nghĩa alpha = 0.05, khoảng tin cậy 95% cũng dao động từ -814.0454 đến 8607.9489 tức là lương của nam giáo sư có thể caoo hơn nữ giáo sư 8607 nhưng cũng có thể là ít hơn 814.04. Như vậy sự dao động này trải dài từ âm tới dương, do đó kết quả trên cho thấy chưa có ý nghĩa thống kê về lương và giới tính của giáo sư.

Phân tích bằng phương pháp hồi quy tuyến tính

Giả sử chúng ta có y là biến Salary và x là biến giới tính. Mô hình hồi quy phát biểu: \[y_{i}=\alpha +\beta x_{i}+\epsilon_{i}\] Chúng ta chỉ có thêm một chú ý, đó là khi biến \(x_{i}=0\) thì sẽ là nam giới, còn khi \(x_{i}=1\) thì là nữ giới. Ở đây tham số \(\beta\) phản ánh mức độ khác biệt về lương giữa nam và nữ.

model.salary <- lm(Salary ~ Sex, data = datasalary)
summary(model.salary)
## 
## Call:
## lm(formula = Salary ~ Sex, data = datasalary)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -18571  -5730    -19   4859  26993 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    56510       1620  34.878   <2e-16 ***
## Sex1           -3897       2455  -1.587    0.118    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9585 on 60 degrees of freedom
## Multiple R-squared:  0.0403, Adjusted R-squared:  0.0243 
## F-statistic: 2.519 on 1 and 60 DF,  p-value: 0.1177

Giá trị Multiple R-squared ở đây bằng 0.0403 tức là giới tính chỉ giải thích được khoảng 4% giá trị phương sai của Lương.

Biến dự đoán là biến phân loại

Trên đây chúng ta mới chỉ xem xét khi biến dự đoán là biến nhị phân chỉ có hai giá trị, bây giờ chúng ta suy rộng ra biến dự đoán là biến phân loại, tức là nhiều hơn 2 giá trị. Cùng xem qua một bộ dữ liệu kinh điển liên quan đến các yếu tố ảnh hưởng tới cân nặng của trẻ sơ sinh. Đây là nghiên cứu trên 189 bà mẹ và trẻ sơ sinh, được thực hiện bởi trung tâm Y Khoa Baystate (Mỹ) vào năm 1986. Biến phụ thuộc ở đây là cân nặng của trẻ sơ sinh, tính bằng gram, ngoài ra các biến khác trong bộ dữ liệu bao gồm:

  1. age: Tuổi của bà mẹ

  2. lwt: Cân nặng của bà mẹ

  3. race: Sắc tộc của bà mẹ (1 = da trắng; 2 = da đen; 3 = loại khác)

  4. smoke: Bà mẹ có hút thuốc hay không (1 là có; 0 là không hút)

  5. ptl: Số lần sinh (premature labors)

  6. ht: Có tiền sử cao huyết áp hay không (1 là Có; 0 là không bị)

  7. ui: Tử cung có bị khó chịu hay không (1 là Có; 0 là không bị)

  8. ftv: Số lần đến khám bác sỹ

  9. bwt: Trọng lượng của trẻ sơ sinh khi mới sinh

Dữ liệu này có sẵn trong thư viện MASS, có tên là birthwt.

data("birthwt")
birthwt %>% 
  head() %>% 
  kable()
low age lwt race smoke ptl ht ui ftv bwt
85 0 19 182 2 0 0 0 1 0 2523
86 0 33 155 3 0 0 0 0 3 2551
87 0 20 105 1 1 0 0 0 1 2557
88 0 21 108 1 1 0 0 1 2 2594
89 0 18 107 1 1 0 0 1 0 2600
91 0 21 124 3 0 0 0 0 0 2622

Ở nội dung này, chúng ta quan tâm đến việc sắc tộc ảnh hưởng tới cân nặng của trẻ sơ sinh như thế nào? Dĩ nhiên, ở đây chúng ta thấy có nhiều nhóm khác nhau, phương pháp được nghĩ tới đầu tiên chính là kiểm định phương sai hay kiểm định ANOVA. Cùng trực quan dữ liệu:

birthwt$race <- as.factor(birthwt$race)
ggplot(birthwt, aes(x = reorder(race, bwt), y = bwt, fill = race, col = race)) + 
  geom_boxplot(col = "black", alpha = 0.6, notch = TRUE) + coord_flip()
## notch went outside hinges. Try setting notch=FALSE.

Biểu đồ trên cho thấy, nhóm các bà mẹ da trắng sinh con có trọng lượng lớn hơn so với nhóm các bà mẹ da đen và các nhóm khác. Chúng ta sử dụng mô hình hồi quy để biết sự khác biệt đó có ý nghĩa thống kê hay không?

model.race <- lm(bwt ~ race, data = birthwt)
summary(model.race)
## 
## Call:
## lm(formula = bwt ~ race, data = birthwt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2096.28  -502.72   -12.72   526.28  1887.28 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3102.72      72.92  42.548  < 2e-16 ***
## race2        -383.03     157.96  -2.425  0.01627 *  
## race3        -297.44     113.74  -2.615  0.00965 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 714.5 on 186 degrees of freedom
## Multiple R-squared:  0.05017,    Adjusted R-squared:  0.03996 
## F-statistic: 4.913 on 2 and 186 DF,  p-value: 0.008336

Giá trị p-value ở đây cho thấy, giá trị cân nặng của trẻ sơ sinh các bà mẹ da đen và nhóm khác đều thấp hơn so với nhóm bà mẹ da trắng. Ở đây chúng ta dùng giá trị tham chiếu là 1 (tức là giá trị của bà mẹ da trắng). Các giá trị thu thập được này đều có ý nghĩa thống kê, với trị số p-value = 0.016 và 0.00965. Tuy nhiên trong kết quả này thì không cung cấp sự khác biệt giữa nhóm da đen và nhóm khác (trừ nhóm da trắng). So sánh nhiều nhóm sau khi đã có kết quả chung trong thống kê gọi là phân tích hậu định “post-hoc comparison”.

Mức ý nghĩa của kiểm định giả thuyết

Không phải tất cả các kết quả kiểm định giả thuyết đều như nhau. Một thử nghiệm giả thuyết hoặc thử nghiệm có ý nghĩa thống kê thường có một mức ý nghĩa gắn liền với nó. Mức ý nghĩa này là một con số thường được biểu thị bằng chữ cái Hy Lạp alpha \((\alpha)\). Thông thường người ta chọn mức ý nghĩa alpha liên quan tới mức độ dương tính giả trong kiểm định giả thuyết thống kê, đa số các con số được chọn trong đoạn [0;1] và con số 0.05 hay 5% hay được chọn nhiều nhất, tức là chúng ta chấp nhận một giả thuyết được kiểm định thì sai số dương tính giả là 0.05 (đối với 1 giả thuyết.), xác suất chúng ta đúng là \(1-\alpha\). Khi kiểm định hai giả thuyết, xác suất chúng ta đúng cả hai là \((1-\alpha)^2\), và xác suất dương tính giả là \(1-(1-\alpha)^2\). Tương tự như vậy, khi chúng ta có \(k\) giả thuyết thống kê thì xác suất dương tính giả là \(1-(1-\alpha)^k\). Theo tập dữ liệu trên chúng ta có 3 nhóm, do đó sẽ có kiểm định cặp như sau (1-2; 1-3; 2-3). Và khi đó sai số dương tính giả có xác suất là 14%, không còn là 5% nữa. Điều này đưa tới cho chúng ta một câu hỏi là cần phải hiệu chỉnh trị số \(p-value\) như thế nào?

Sai số dương tính giả

Khi triển khai một kỹ thuật xét nghiệm mới, nhất là với các xét nghiệm định tính và bán định lượng, người quản lý Phòng xét nghiệm thường quan tâm đến các chỉ số sau: Độ nhạy và độ đặc hiệu của Phương pháp xét nghiệm. Khi đề cập đến Độ nhạy và Độ đặc hiệu có hai khái niệm là Dương tính giả và âm tính giả.

Độ nhạy của một xét nghiệm là tỷ lệ những trường hợp thực sự có bệnh và có kết quả xét nghiệm dương tính trong toàn bộ các trường hợp có bệnh. Công thức để tính độ nhạy như sau:

Độ nhạy = Số dương tính thật/(số dương tính thật + số âm tính giả)

Độ đặc hiệu của một xét nghiệm là tỷ lệ những trường hợp thực sự không có bệnh và có kết quả xét nghiệm âm tính trong toàn bộ các trường hợp không bị bệnh. Độ đặc hiệu được tính theo công thức sau:

Độ đặc hiệu = Số trường hợp âm tính thật/ (số trường hợp âm tính thật + số trường hợp dương tính giả)

Khi đề cập đến Độ nhạy và Độ đặc hiệu có hai khái niệm là Dương tính giả và âm tính giả. Vậy Dương tính giả và Âm tính giả là gì?

  • Dương tính giả là kết quả xét nghiệm Dương tính trong khi đó người đi kiểm tra thực sự không có bệnh.

  • Âm tính giả là kết quả xét nghiệm Âm tính trong khi đó người đi kiểm tra thực sự là có bệnh.

Như vậy, âm tính giả hay dương tính giả là tình trạng kết quả xét nghiệm thu được không chính xác với tình trạng bệnh.

Confusion matrix

Hiệu chỉnh p-value bằng Tukey

Chúng ta gọi lại mô hình hồi quy tuyến tính và hiệu chỉnh

comparison <- glht(model.race, mcp(race = "Tukey"))
summary(comparison)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = bwt ~ race, data = birthwt)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)  
## 2 - 1 == 0  -383.03     157.96  -2.425   0.0416 *
## 3 - 1 == 0  -297.44     113.74  -2.615   0.0254 *
## 3 - 2 == 0    85.59     165.09   0.518   0.8603  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(comparison)

Nếu khoảng tin cậy lệch hẳn ra khỏi đường tham chiếu 0 thì kết quả có ý nghĩa thống kê, còn nếu cắt đường tham chiếu thì kết quả không có ý nghĩa thống kê.