ngày 4
Thu Trang Nguyen
2024-09-11
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: # Việc 1: Đọc dữ liệu vào R
ob = read.csv ("Desktop/obesity data.csv")Việc 2: Đánh giá mối liên quan giữa tuổi và mật độ xương toàn thân
Biểu đồ histogram đánh giá phân bố tuổi và mật độ xương toàn thân
library (ggplot2)
library (gridExtra)
p = ggplot(data = ob, aes(x = age))
p1 = p + geom_histogram(fill = "blue", col = "white") + labs(x = "Tuổi (năm)", y = "Số người", title = "Phân bố tuổi")
p = ggplot(data = ob, aes(x = wbbmd))
p2 = p + geom_histogram(fill = "blue", col = "white") + labs(x = "Mật độ xương toàn thân (g/cm2)", y = "Số người", title = "Phân bố MĐX toàn thân")
grid.arrange(p1, p2, ncol = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Biểu đồ tán xạ đánh giá mối liên quan giữa tuổi và MĐX toàn thân
p = ggplot(data = ob, aes(x = age, y = wbbmd))
p + geom_point() + geom_smooth(method = "lm", formula = y~ x)
# Việc 2: Phân tích tương quan đánh giá mối liên quan giữa tuổi và MĐX
toàn thân
cor.test(ob$age, ob$wbbmd, method= "pearson")##
## Pearson's product-moment correlation
##
## data: ob$age and ob$wbbmd
## t = -17.154, df = 1215, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4856972 -0.3951677
## sample estimates:
## cor
## -0.4415556Việc 3: Mô hình tuyến tính đánh giá mối liên quan giữa tuổi và MĐX toàn thân
Thực hiện mô hình hồi quy
m.1 = lm(wbbmd ~ age, data = ob)
summary(m.1)##
## Call:
## lm(formula = wbbmd ~ age, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32749 -0.07268 -0.00533 0.06793 0.33178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1450766 0.0084638 135.29 <2e-16 ***
## age -0.0028914 0.0001686 -17.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1015 on 1215 degrees of freedom
## Multiple R-squared: 0.195, Adjusted R-squared: 0.1943
## F-statistic: 294.3 on 1 and 1215 DF, p-value: < 2.2e-16Đánh giá giả định của mô hình
par(mfrow = c(2,2))
plot(m.1)
library (ggfortify)
autoplot (m.1)
#Việc 4: Đánh gía mối liên quan giữa tuổi và giới tính với MĐX toàn thân ## Biểu đồ tán xạ đánh giá mối liên quan giữa tuổi với MĐX toàn thân theo giới tính
p = ggplot(data = ob, aes(x = age, y = wbbmd, fill = gender, col = gender))
p1 = p + geom_point() + geom_smooth() + labs(x = "Tuổi (năm)", y = "Mật độ xương toàn thân (g/cm2)") + ggtitle("Liên quan giữa tuổi và MĐX theo giới tính")
p1## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Nhận xét về ảnh hưởng của giới tính lên mối liên quan giữa tuổi với
MĐX toàn thân ## Thực hiện mô hình hồi quy tuyến tính
m.2 = lm(wbbmd ~ age + gender, data = ob)
summary(m.2)##
## Call:
## lm(formula = wbbmd ~ age + gender, data = ob)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36272 -0.06658 -0.00411 0.06549 0.34473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.118288 0.008636 129.485 <2e-16 ***
## age -0.002691 0.000164 -16.408 <2e-16 ***
## genderM 0.059417 0.006230 9.537 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09798 on 1214 degrees of freedom
## Multiple R-squared: 0.2511, Adjusted R-squared: 0.2498
## F-statistic: 203.5 on 2 and 1214 DF, p-value: < 2.2e-16Nhận xét: Trong cùng giới tính, khi tăng một tuổi thì giảm MĐX toàn thân 0.0027 g/m2, có ý nghĩa thống kê
Kiểm tra giả định mô hình
par(mfrow = c(2,2))
plot(m.2)
#
Việc 5: Xây dựng mô hình dự báo MĐX
library (readxl)
fx = as.data.frame(read_excel("Desktop/Osteo-data.xlsx"))
dim(fx)## [1] 2216 17summary(fx)## id sex age weight
## Min. : 1.0 Length:2216 Min. :57.00 Min. : 34.00
## 1st Qu.: 554.8 Class :character 1st Qu.:65.00 1st Qu.: 60.00
## Median :1108.5 Mode :character Median :70.00 Median : 69.00
## Mean :1108.5 Mean :70.89 Mean : 70.14
## 3rd Qu.:1662.2 3rd Qu.:76.00 3rd Qu.: 79.00
## Max. :2216.0 Max. :96.00 Max. :133.00
## NA's :53
## height prior_fx fnbmd smoking
## Min. :136.0 Min. :0.0000 Min. :0.2800 Min. :0.0000
## 1st Qu.:158.0 1st Qu.:0.0000 1st Qu.:0.7300 1st Qu.:0.0000
## Median :164.0 Median :0.0000 Median :0.8200 Median :0.0000
## Mean :164.9 Mean :0.1611 Mean :0.8287 Mean :0.4176
## 3rd Qu.:171.0 3rd Qu.:0.0000 3rd Qu.:0.9300 3rd Qu.:1.0000
## Max. :196.0 Max. :1.0000 Max. :1.5100 Max. :1.0000
## NA's :54 NA's :89 NA's :1
## parkinson rheum hypertension diabetes
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.000
## Median :0.00000 Median :0.00000 Median :1.0000 Median :0.000
## Mean :0.06498 Mean :0.03881 Mean :0.5063 Mean :0.111
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.000
## Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :1.000
##
## copd cancer cvd falls_n
## Min. :0.000 Min. :0.00000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.000 Median :0.00000 Median :0.0000 Median :0.0000
## Mean :0.111 Mean :0.08529 Mean :0.3872 Mean :0.2843
## 3rd Qu.:0.000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.0000
## Max. :1.000 Max. :1.00000 Max. :1.0000 Max. :2.0000
##
## fx
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.2595
## 3rd Qu.:1.0000
## Max. :1.0000
## Mô tả đặc điểm của dữ liệu theo giới tính
library (table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1 (~ age + weight + height + fnbmd + as.factor(prior_fx) + as.factor(falls_n) + as.factor(smoking) + as.factor(parkinson) + as.factor(rheum) + as.factor(hypertension) + as.factor(diabetes) + as.factor(copd) + as.factor(cancer) + as.factor(cvd) + as.factor(fx) | sex, data = fx)| Female (N=1358) |
Male (N=858) |
Overall (N=2216) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 71.2 (7.59) | 70.4 (6.44) | 70.9 (7.17) |
| Median [Min, Max] | 70.0 [57.0, 96.0] | 69.0 [59.0, 92.0] | 70.0 [57.0, 96.0] |
| weight | |||
| Mean (SD) | 64.9 (12.5) | 78.2 (12.7) | 70.1 (14.2) |
| Median [Min, Max] | 64.0 [34.0, 115] | 78.0 [45.0, 133] | 69.0 [34.0, 133] |
| Missing | 42 (3.1%) | 11 (1.3%) | 53 (2.4%) |
| height | |||
| Mean (SD) | 160 (6.37) | 173 (6.86) | 165 (9.35) |
| Median [Min, Max] | 160 [136, 181] | 173 [151, 196] | 164 [136, 196] |
| Missing | 41 (3.0%) | 13 (1.5%) | 54 (2.4%) |
| fnbmd | |||
| Mean (SD) | 0.777 (0.132) | 0.909 (0.153) | 0.829 (0.155) |
| Median [Min, Max] | 0.770 [0.280, 1.31] | 0.900 [0.340, 1.51] | 0.820 [0.280, 1.51] |
| Missing | 57 (4.2%) | 32 (3.7%) | 89 (4.0%) |
| as.factor(prior_fx) | |||
| 0 | 1125 (82.8%) | 734 (85.5%) | 1859 (83.9%) |
| 1 | 233 (17.2%) | 124 (14.5%) | 357 (16.1%) |
| as.factor(falls_n) | |||
| 0 | 1063 (78.3%) | 671 (78.2%) | 1734 (78.2%) |
| 1 | 206 (15.2%) | 128 (14.9%) | 334 (15.1%) |
| 2 | 89 (6.6%) | 59 (6.9%) | 148 (6.7%) |
| as.factor(smoking) | |||
| 0 | 962 (70.8%) | 328 (38.2%) | 1290 (58.2%) |
| 1 | 395 (29.1%) | 530 (61.8%) | 925 (41.7%) |
| Missing | 1 (0.1%) | 0 (0%) | 1 (0.0%) |
| as.factor(parkinson) | |||
| 0 | 1268 (93.4%) | 804 (93.7%) | 2072 (93.5%) |
| 1 | 90 (6.6%) | 54 (6.3%) | 144 (6.5%) |
| as.factor(rheum) | |||
| 0 | 1306 (96.2%) | 824 (96.0%) | 2130 (96.1%) |
| 1 | 52 (3.8%) | 34 (4.0%) | 86 (3.9%) |
| as.factor(hypertension) | |||
| 0 | 695 (51.2%) | 399 (46.5%) | 1094 (49.4%) |
| 1 | 663 (48.8%) | 459 (53.5%) | 1122 (50.6%) |
| as.factor(diabetes) | |||
| 0 | 1213 (89.3%) | 757 (88.2%) | 1970 (88.9%) |
| 1 | 145 (10.7%) | 101 (11.8%) | 246 (11.1%) |
| as.factor(copd) | |||
| 0 | 1211 (89.2%) | 759 (88.5%) | 1970 (88.9%) |
| 1 | 147 (10.8%) | 99 (11.5%) | 246 (11.1%) |
| as.factor(cancer) | |||
| 0 | 1235 (90.9%) | 792 (92.3%) | 2027 (91.5%) |
| 1 | 123 (9.1%) | 66 (7.7%) | 189 (8.5%) |
| as.factor(cvd) | |||
| 0 | 843 (62.1%) | 515 (60.0%) | 1358 (61.3%) |
| 1 | 515 (37.9%) | 343 (40.0%) | 858 (38.7%) |
| as.factor(fx) | |||
| 0 | 932 (68.6%) | 709 (82.6%) | 1641 (74.1%) |
| 1 | 426 (31.4%) | 149 (17.4%) | 575 (25.9%) |
Xây dựng mô hình dự báo MĐX tại cổ xương đùi
Phương pháp stepwise
fx.2 = na.omit(fx)
m.step = lm(fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking + parkinson + rheum + hypertension + diabetes + copd + cancer + cvd, data = fx.2)
step = step(m.step)## Start: AIC=-9100.52
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + copd + cancer +
## cvd
##
## Df Sum of Sq RSS AIC
## - copd 1 0.0000 28.641 -9102.5
## - cvd 1 0.0029 28.643 -9102.3
## - hypertension 1 0.0037 28.644 -9102.2
## - cancer 1 0.0056 28.646 -9102.1
## - parkinson 1 0.0074 28.648 -9102.0
## - diabetes 1 0.0086 28.649 -9101.9
## - falls_n 1 0.0134 28.654 -9101.5
## <none> 28.641 -9100.5
## - rheum 1 0.0433 28.684 -9099.3
## - height 1 0.1607 28.801 -9090.7
## - prior_fx 1 0.3350 28.976 -9077.9
## - smoking 1 0.3807 29.021 -9074.5
## - sex 1 0.8234 29.464 -9042.4
## - age 1 2.1037 30.744 -8952.2
## - weight 1 4.9616 33.602 -8763.7
##
## Step: AIC=-9102.52
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + cancer + cvd
##
## Df Sum of Sq RSS AIC
## - cvd 1 0.0029 28.643 -9104.3
## - hypertension 1 0.0037 28.644 -9104.2
## - cancer 1 0.0056 28.646 -9104.1
## - parkinson 1 0.0074 28.648 -9104.0
## - diabetes 1 0.0086 28.649 -9103.9
## - falls_n 1 0.0134 28.654 -9103.5
## <none> 28.641 -9102.5
## - rheum 1 0.0433 28.684 -9101.3
## - height 1 0.1607 28.801 -9092.7
## - prior_fx 1 0.3351 28.976 -9079.9
## - smoking 1 0.3810 29.022 -9076.5
## - sex 1 0.8236 29.464 -9044.4
## - age 1 2.1053 30.746 -8954.1
## - weight 1 4.9628 33.603 -8765.6
##
## Step: AIC=-9104.3
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + hypertension + diabetes + cancer
##
## Df Sum of Sq RSS AIC
## - hypertension 1 0.0045 28.648 -9106.0
## - cancer 1 0.0057 28.649 -9105.9
## - parkinson 1 0.0074 28.651 -9105.8
## - diabetes 1 0.0076 28.651 -9105.7
## - falls_n 1 0.0134 28.657 -9105.3
## <none> 28.643 -9104.3
## - rheum 1 0.0426 28.686 -9103.2
## - height 1 0.1604 28.804 -9094.5
## - prior_fx 1 0.3346 28.978 -9081.7
## - smoking 1 0.3826 29.026 -9078.2
## - sex 1 0.8236 29.467 -9046.2
## - age 1 2.1048 30.748 -8955.9
## - weight 1 4.9620 33.605 -8767.4
##
## Step: AIC=-9105.97
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + diabetes + cancer
##
## Df Sum of Sq RSS AIC
## - cancer 1 0.0057 28.654 -9107.5
## - diabetes 1 0.0066 28.655 -9107.5
## - parkinson 1 0.0077 28.656 -9107.4
## - falls_n 1 0.0130 28.661 -9107.0
## <none> 28.648 -9106.0
## - rheum 1 0.0428 28.691 -9104.8
## - height 1 0.1613 28.809 -9096.1
## - prior_fx 1 0.3334 28.981 -9083.4
## - smoking 1 0.3807 29.029 -9080.0
## - sex 1 0.8201 29.468 -9048.1
## - age 1 2.1035 30.751 -8957.7
## - weight 1 4.9596 33.608 -8769.3
##
## Step: AIC=-9107.55
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum + diabetes
##
## Df Sum of Sq RSS AIC
## - diabetes 1 0.0069 28.661 -9109.0
## - parkinson 1 0.0076 28.661 -9109.0
## - falls_n 1 0.0135 28.667 -9108.5
## <none> 28.654 -9107.5
## - rheum 1 0.0442 28.698 -9106.3
## - height 1 0.1628 28.817 -9097.5
## - prior_fx 1 0.3351 28.989 -9084.9
## - smoking 1 0.3854 29.039 -9081.2
## - sex 1 0.8219 29.476 -9049.6
## - age 1 2.0984 30.752 -8959.6
## - weight 1 4.9596 33.613 -8771.0
##
## Step: AIC=-9109.04
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## parkinson + rheum
##
## Df Sum of Sq RSS AIC
## - parkinson 1 0.0080 28.669 -9110.4
## - falls_n 1 0.0135 28.674 -9110.0
## <none> 28.661 -9109.0
## - rheum 1 0.0450 28.706 -9107.7
## - height 1 0.1614 28.822 -9099.1
## - prior_fx 1 0.3366 28.997 -9086.3
## - smoking 1 0.3844 29.045 -9082.8
## - sex 1 0.8262 29.487 -9050.8
## - age 1 2.1055 30.766 -8960.7
## - weight 1 4.9622 33.623 -8772.4
##
## Step: AIC=-9110.44
## fnbmd ~ age + sex + weight + height + prior_fx + falls_n + smoking +
## rheum
##
## Df Sum of Sq RSS AIC
## - falls_n 1 0.0132 28.682 -9111.5
## <none> 28.669 -9110.4
## - rheum 1 0.0454 28.714 -9109.1
## - height 1 0.1614 28.830 -9100.5
## - prior_fx 1 0.3345 29.003 -9087.8
## - smoking 1 0.3848 29.053 -9084.2
## - sex 1 0.8261 29.495 -9052.2
## - age 1 2.1123 30.781 -8961.7
## - weight 1 4.9672 33.636 -8773.5
##
## Step: AIC=-9111.47
## fnbmd ~ age + sex + weight + height + prior_fx + smoking + rheum
##
## Df Sum of Sq RSS AIC
## <none> 28.682 -9111.5
## - rheum 1 0.0467 28.729 -9110.0
## - height 1 0.1639 28.846 -9101.4
## - prior_fx 1 0.3383 29.020 -9088.6
## - smoking 1 0.3900 29.072 -9084.8
## - sex 1 0.8244 29.506 -9053.4
## - age 1 2.1100 30.792 -8962.9
## - weight 1 4.9608 33.643 -8775.1Phương pháp BMA
library (BMA)## 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.7-6)xvars = fx.2[, c("age", "sex", "weight", "height", "prior_fx", "falls_n", "smoking", "parkinson", "rheum", "hypertension", "diabetes", "copd", "cancer", "cvd")]
m.bma = bicreg(xvars, fx.2$fnbmd, strict = FALSE, OR = 20)
summary(m.bma)##
## Call:
## bicreg(x = xvars, y = fx.2$fnbmd, strict = FALSE, OR = 20)
##
##
## 3 models were selected
## Best 3 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100.0 0.626343 0.0977859 6.071e-01 6.098e-01 8.466e-01
## age 100.0 -0.004784 0.0003880 -4.765e-03 -4.778e-03 -4.995e-03
## sexMale 100.0 0.061554 0.0088666 6.021e-02 6.039e-02 7.690e-02
## weight 100.0 0.004329 0.0002379 4.306e-03 4.302e-03 4.603e-03
## height 92.1 0.001397 0.0005836 1.519e-03 1.503e-03 .
## prior_fx 100.0 -0.035389 0.0070901 -3.532e-02 -3.535e-02 -3.611e-02
## falls_n 0.0 0.000000 0.0000000 . . .
## smoking 100.0 -0.029091 0.0054107 -2.910e-02 -2.898e-02 -2.918e-02
## parkinson 0.0 0.000000 0.0000000 . . .
## rheum 10.0 0.002423 0.0083564 . 2.421e-02 .
## hypertension 0.0 0.000000 0.0000000 . . .
## diabetes 0.0 0.000000 0.0000000 . . .
## copd 0.0 0.000000 0.0000000 . . .
## cancer 0.0 0.000000 0.0000000 . . .
## cvd 0.0 0.000000 0.0000000 . . .
##
## nVar 6 7 5
## r2 0.434 0.435 0.431
## BIC -1.162e+03 -1.157e+03 -1.157e+03
## post prob 0.821 0.100 0.079Nhận xét: Nên chọn mô hình thứ nhất, hơn là mô hình thứ 3 (đưa ra lựa chọn vì mô hình 3 ít biến nhất) vì: post prob 82% so với 8% (lựa chọn mô hình chú trọng nhất post prob - xác xuất hậu định), trong khi biến height không khó thu thập (có thể cân nhắc lựa chọn mô hình 3 nếu biến không lựa chọn rất khó thu thập)
imageplot.bma(m.bma)
## Kiểm định giả định của mô hình
m.bmd = lm(fnbmd ~ age + sex + weight + height + prior_fx + smoking, data = fx.2)
summary(m.bmd)##
## Call:
## lm(formula = fnbmd ~ age + sex + weight + height + prior_fx +
## smoking, data = fx.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37425 -0.07568 -0.00663 0.07184 0.58750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6070710 0.0758296 8.006 1.94e-15 ***
## age -0.0047646 0.0003834 -12.428 < 2e-16 ***
## sexMale 0.0602131 0.0077533 7.766 1.25e-14 ***
## weight 0.0043063 0.0002252 19.126 < 2e-16 ***
## height 0.0015189 0.0004328 3.509 0.000459 ***
## prior_fx -0.0353237 0.0070861 -4.985 6.70e-07 ***
## smoking -0.0290965 0.0054097 -5.379 8.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1166 on 2114 degrees of freedom
## Multiple R-squared: 0.4341, Adjusted R-squared: 0.4325
## F-statistic: 270.3 on 6 and 2114 DF, p-value: < 2.2e-16par(mfrow = c(2,2))
plot(m.bmd)
## Viết công thức mô hình tối ưu
summary (m.bmd)##
## Call:
## lm(formula = fnbmd ~ age + sex + weight + height + prior_fx +
## smoking, data = fx.2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37425 -0.07568 -0.00663 0.07184 0.58750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6070710 0.0758296 8.006 1.94e-15 ***
## age -0.0047646 0.0003834 -12.428 < 2e-16 ***
## sexMale 0.0602131 0.0077533 7.766 1.25e-14 ***
## weight 0.0043063 0.0002252 19.126 < 2e-16 ***
## height 0.0015189 0.0004328 3.509 0.000459 ***
## prior_fx -0.0353237 0.0070861 -4.985 6.70e-07 ***
## smoking -0.0290965 0.0054097 -5.379 8.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1166 on 2114 degrees of freedom
## Multiple R-squared: 0.4341, Adjusted R-squared: 0.4325
## F-statistic: 270.3 on 6 and 2114 DF, p-value: < 2.2e-16