Phân tích hồi quy và tương quan

1 Phân tích tương quan

1.1 Ví dụ:

1.2 Hệ số tương quan tuyến tính (Correlation coefficient)

Hệ số tương quan đo lường mức độ quan hệ tuyến tính giữa hai biến, không phân biệt biến này phụ thuộc vào biến kia

1.3 Ðặc tính của hệ số tương quan

• Hệ số tương quan không có đơn vị
• Hệ số tương quan (r) nằm trong khoảng [-1,1]

• r > 0: tương quan dương
• r < 0: tương quan âm
• r = 0: không tương quan

1.4 Các phương pháp tính tương quan

### Pearson’s r - Đánh giá mức độ tương quan tuyến tính của 2 biến định lượng

\[ r_{xy} = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt[]{\sum_{i=1}^{n}(x_i - \bar{x})^2 \sum_{i=1}^{n}(y_i - \bar{y})^2}} \]

• Ví dụ: Sử dụng số liệu dân số và thu nhập của 50 bang của Mỹ

states <- state.x77[,c(1:2)] 
head(states)

##            Population Income
## Alabama          3615   3624
## Alaska            365   6315
## Arizona          2212   4530
## Arkansas         2110   3378
## California      21198   5114
## Colorado         2541   4884

• Công thức tính hệ số tương quan Pearson trên R: cor(df, method = “pearson”)

M <- cor(states)
M

##            Population    Income
## Population  1.0000000 0.2082276
## Income      0.2082276 1.0000000

1.4.1 Tương quan hạng (spearman)

• Đánh giá mức độ tương quan của 2 hạng của 2 biến (rank-ordered variables), sử dụng khi phân phối của tổng thể được giả sử không phải là phân phối chuẩn hoặc trong trường hợp có các giá trị quan sát bất thường (lớn quá hoặc nhỏ quá)

\[ 1-\frac{6\sum_i^n d_{i}^2}{n(n^2-1)}\]

di: hiệu hạng của 2 biến

\[ d_{i} = rgX_{i} - rgY_{i}\]

• Công thức tính trên R: cor(df, method = “spearman”)

• Minh họa bằng bảng tính:

states <- state.x77[,c(1:2)] 
head(states)

##            Population Income
## Alabama          3615   3624
## Alaska            365   6315
## Arizona          2212   4530
## Arkansas         2110   3378
## California      21198   5114
## Colorado         2541   4884

n <- dim(states)[1]

states %>% 
  as.data.frame() %>% 
  mutate(rgX = rank(Population, ties.method= "first"),
         rgY = rank(Income, ties.method= "first"),
         d = rgX - rgY,
         d2 = d^2) %>% 
  mutate(Spearman_cor = 1 - 6 * sum(d2)/(n * (n^2 -1))) %>% 
  datatable()

M <- cor(states[,c("Population", "Income")], method = "spearman")
M

##            Population    Income
## Population  1.0000000 0.1246098
## Income      0.1246098 1.0000000

1.4.2 Tương quan hạng (kendall)

• Đánh giá mức độ tương quan của 2 hạng của 2 biến (rank-ordered variables), hệ số này được sử dụng tương tự như spearman, thông thường hệ số này nhỏ hơn spearman

• Hệ số kendall ít dùng hơn so với 2 hệ số tương quan trên

• Công thức tính trên R: cor(df, method = “kendall”)

M <- cor(states, method = "kendall")
M

##            Population     Income
## Population 1.00000000 0.08408163
## Income     0.08408163 1.00000000

1.5 Kiểm định xem 2 biến có tương quan với nhau không

• Cũng như phương pháp tính, kiểm định cũng có 3 method: Pearson, Spearman, Kendall

• Giả thiết kiểm định:

  • H0: Không có tương quan (hệ số tương quan = 0)
  • Ha: Có tương quan

• Ví dụ: Sử dụng hàm cor.test của gói stats

states <- state.x77
states[,c(3,5)] %>% cor()

##            Illiteracy    Murder
## Illiteracy  1.0000000 0.7029752
## Murder      0.7029752 1.0000000

cor.test(states[,3], states[,5])

## 
##  Pearson's product-moment correlation
## 
## data:  states[, 3] and states[, 5]
## t = 6.8479, df = 48, p-value = 1.258e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5279280 0.8207295
## sample estimates:
##       cor 
## 0.7029752

Hàm này chỉ dùng được với 2 biến. Trong trường hợp muốn thực hiện với nhiều biến có thể sử dụng hàm corr.test của package psych

• Ví dụ sử dụng hàm corr.test

corr.test(states[,1:5], use = "complete", method = "pearson" )

## Call:corr.test(x = states[, 1:5], use = "complete", method = "pearson")
## Correlation matrix 
##            Population Income Illiteracy Life Exp Murder
## Population       1.00   0.21       0.11    -0.07   0.34
## Income           0.21   1.00      -0.44     0.34  -0.23
## Illiteracy       0.11  -0.44       1.00    -0.59   0.70
## Life Exp        -0.07   0.34      -0.59     1.00  -0.78
## Murder           0.34  -0.23       0.70    -0.78   1.00
## Sample Size 
## [1] 50
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##            Population Income Illiteracy Life Exp Murder
## Population       0.00   0.44       0.91     0.91   0.09
## Income           0.15   0.00       0.01     0.09   0.43
## Illiteracy       0.46   0.00       0.00     0.00   0.00
## Life Exp         0.64   0.02       0.00     0.00   0.00
## Murder           0.01   0.11       0.00     0.00   0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Hàm này cho ra kết quả của cả hệ số tương quan và xác suất kiểm định

Xác suất > mức ý nghĩa alpha (= 0.05) có thể kết luận hệ số tương quan = 0 với mức ý nghĩa alpha

1.6 Một số đồ thị trong phân tích tương quan

• corrplot (package: corrplot)

states <- state.x77
M <- cor(states)
corrplot(M, method = "circle")

• corrplot.mixed (package: corrplot)

corrplot.mixed(M)

• chart.Correlation (package: PerformanceAnalytics)

data(managers)
chart.Correlation(managers[,1:8], histogram=TRUE, pch="+")

Hàm này tương như: psych::pairs.panels, GGally::ggpairs, corrr::correlate

2 Phân tích hồi quy (regression analysis)

Sử dụng để dự báo 1 biến phụ thuộc (dependent variable, response variable) dựa vào 1 hay nhiều biến độc lập (independent variables, predictor variables, explanatory variables)

Ví dụ 1: Phân tích khi tăng 1 nhân viên thì lợi nhuận ngân hàng tăng hay giảm bao nhiêu tiền. Biến phụ thuộc là lợi nhuận ngân hàng, biến độc lập là số lượng nhân viên (simple linear)

Ví dụ 2: Tính toán xem khi tăng 1 cây ATM thì lợi nhuận ngân hàng tăng hay giảm bao nhiêu phần trăm. Biến phụ thuôc là log(lợi nhuận), biến độc lập có thể là số lượng máy ATM, số lượng máy ATM bình phương (Polynomial)

Ví dụ 3: Dự báo khả năng phát sinh nợ xấu của khách hàng. Biến phụ thuộc là khả năng phát sinh nợ xấu, biến độc lập ví dụ: tuổi, giới tính, trình độ học vấn … (Multiple linear)

Ví dụ 4: Dự báo giá cổ phiếu của ngân hàng tại các thời điểm trong tương lai. Biến phụ thuộc là giá cổ phiếu, biến độc lập có thể có là trễ của biến giá, hoặc 1 số yếu tố như GPD, lạm phát … (Time-series)

Có rất nhiều loại hồi quy như: Simple linear, Polynomial, Multiple linear, Multilevel, Multivariate, Logistic, Poisson, Cox proportional hazards, Time-series, Nonlinear, Nonparametric

Trong đó, 4 loại phổ biến hay được sử dụng là

type_model <- data.frame(Type = c("Simple linear", 
                    "Multiple linear",
                    "Logistic",
                    "Time-series"),
           `Typical use` = c("Predicting a quantitative response variable from a quantitative explanatory variable",
"Predicting a quantitative response variable from two or more explanatory variables",
"Predicting a categorical response variable from one or more explanatory variables",
"Modeling time-series data with correlated errors"))

type_model %>% datatable()

2.1 Những hàm thường dùng đối với hồi quy tuyến tính

data.frame(Function = c("summary()", "coefficients()", "confint()", "fitted()", "residuals()", "anova()", 
"vcov()", "AIC()", "plot()", "predict()"), 
Action =
c("Displays detailed results for the fitted model", "Lists the model parameters (intercept and slopes) for the fitted model", "Provides confidence intervals for the model parameters (95% by default)",
  "Lists the predicted values in a fitted model", "Lists the residual values in a fitted model", "Generates an ANOVA table for a fitted model, or an ANOVA table comparing two or more fitted models", "Lists the covariance matrix for model parameters", "Prints Akaike’s Information Criterion", "Generates diagnostic plots for evaluating the fit of a model", "Uses a fitted model to predict response values for a new dataset"))  %>%  datatable()

2.2 Hồi quy tuyến tính đơn biến (Simple linear)

Là mô hình hồi quy với 1 biến phụ thuộc, 1 biến độc lập. Hàm hồi quy có dạng Y = a + bX

Y là biến phụ thuộc, X là biến độc lập, a là hệ số chặn (intercept), b là hệ số góc (coefficient)

Ví dụ: Với bộ dữ liệu women, có 2 biến là chiều cao, cân nặng của phụ nữ. Muốn dự báo cân nặng của phụ nữ dựa vào chiều cao của họ ta xây dựng mô hình hồi quy đơn

data(women)

head(women)

##   height weight
## 1     58    115
## 2     59    117
## 3     60    120
## 4     61    123
## 5     62    126
## 6     63    129

women %>% 
  ggplot(aes(x = height, y = weight))+
  geom_point()+
  labs(x = "height (in inches)",
       y = "weight (in pounds)")

fit <- lm(weight ~ height, data=women)

Sử dụng phương pháp bình phương nhỏ nhất OLS (Ordinary least squares)

summary(fit)

## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

Như vậy phương trình hồi quy có dạng

\[ \hat{weight} = -87.51 + 3.45 * height \]

Dự báo (fitted)

fitted(fit)

##        1        2        3        4        5        6        7        8 
## 112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 
##        9       10       11       12       13       14       15 
## 140.1833 143.6333 147.0833 150.5333 153.9833 157.4333 160.8833

women$weight_pre <- fitted(fit)

Phần dư (sai lệch)

residuals(fit)

##           1           2           3           4           5           6 
##  2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 
##           7           8           9          10          11          12 
## -1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333 
##          13          14          15 
##  0.01666667  1.56666667  3.11666667

women$res <- residuals(fit)

women %>% 
  ggplot()+
  geom_point(aes(x = height, y = weight), color = "blue")+
  geom_point(aes(x = height, y = weight_pre), color = "red")+
  labs(x = "height (in inches)",
       y = "weight (in pounds)")

2.3 Hồi quy tuyến tính đa biến

Mô hình hồi quy dạng: \[\hat{Y} = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \beta_{3}X_{3} + ...\]

Ví dụ lấy bộ số liệu giết người ở các bang của Mỹ state.x77

states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])

states %>% head()

##            Murder Population Illiteracy Income Frost
## Alabama      15.1       3615        2.1   3624    20
## Alaska       11.3        365        1.5   6315   152
## Arizona       7.8       2212        1.8   4530    15
## Arkansas     10.1       2110        1.9   3378    65
## California   10.3      21198        1.1   5114    20
## Colorado      6.8       2541        0.7   4884   166

• Phân tích tương quan

chart.Correlation(states, histogram=TRUE, pch="+")

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
summary(fit)

## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
##     data = states)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7960 -1.6495 -0.0811  1.4815  7.6210 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
## Population  2.237e-04  9.052e-05   2.471   0.0173 *  
## Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
## Income      6.442e-05  6.837e-04   0.094   0.9253    
## Frost       5.813e-04  1.005e-02   0.058   0.9541    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared:  0.567,  Adjusted R-squared:  0.5285 
## F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

states$pre <- fitted(fit)
states$res <- residuals(fit)

Giá trị p-value Income và Frost > 0.05. Vì vậy, với mức ý nghĩa 0.05 có thể kết luận hệ số của 2 biến này bằng 0

• Khoảng tin cậy của hệ số hồi quy

confint(fit)

##                     2.5 %       97.5 %
## (Intercept) -6.552191e+00 9.0213182149
## Population   4.136397e-05 0.0004059867
## Illiteracy   2.381799e+00 5.9038743192
## Income      -1.312611e-03 0.0014414600
## Frost       -1.966781e-02 0.0208304170

Thay vì việc nhìn vào pvalue ở mô hình hồi quy, khoảng tin cậy của hệ số hồi quy cho phép đánh giá xem hệ số hồi quy có khác 0 hay không.

Khoảng tin cậy của Hệ số của Intercept, Income, Frost có cận dưới < 0, trong khi cận trên > 0 –> nhìn vào đây có thể biết được 3 hệ số này = 0 với mức ý nghĩa 0.05%

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

par(mfrow=c(1,1))

2.4 Kiểm định và lựa chọn mô hình

2.4.1 Dựa vào đồ thị phần dư

• Biểu đồ 1: Vẽ tương quan giữa phần dư và kết quả dự báo, giá trị phần dư càng ở quanh mức 0 kết quả dự báo càng tốt. Như vậy, tỷ lệ giết người đang được dự báo quá cao so với thực tế tại 2 bang Rhode Island và Masschusetts, dự báo quá thấp tại bang Nevada

• Biểu đồ 2: Kiểm tra xem phần dư có phân phối chuẩn N(0,1) hay không, kết quả cho thấy phần dư tại 3 bang Nevada và Rhode Island và Alaska đang ko theo quy luật phân phối chuẩn

• Biểu đồ 3: Đánh giá phương sai của phần dư có đồng nhất hay không

• Biều đồ 4: Cho phép phát hiện ra các outliers trong phần dư

2.4.2 Dựa vào các kiểm định

2.4.3 Giả thiết 1: Sai số ngẫu nhiên có phân phối chuẩn

2.4.4 Giả thiết 2: Kỳ vọng của sai số ngẫu nhiên tại mỗi giá trị bằng 0

res <- residuals(fit)

t.test(res, mu = 0)

## 
##  One Sample t-test
## 
## data:  res
## t = 6.5292e-16, df = 49, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6903919  0.6903919
## sample estimates:
##    mean of x 
## 2.243128e-16

p-value = 1: mô hình thỏa mãn giả thiết 2

2.4.5 Giả thiết 3: Phương sai của sai số ngẫu nhiên không đổi

H0: phương sai sai số không đổi Ha: Phương sai sai số thay đổi

ncvTest(fit)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.746514, Df = 1, p = 0.18632

p-value = 0.18, có thể kết luận phương sai sai số không đổi với mức ý nghĩa < 0.18

2.4.6 Giả thiết 4: Giữa các biến độc lập không có mối quan hệ đa cộng tuyến hoàn hảo

vif(fit)

## Population Illiteracy     Income      Frost 
##   1.245282   2.165848   1.345822   2.082547

vif < 10 cho thấy các biến độc lập không có đa cộng tuyến

2.5 Lựa chọn mô hình

Sau khi kiểm định các giả thiết của mô hình OLS, thử nghiệm trường hợp bỏ 2 biến Income, Frost ra khỏi mô hình

fit <- lm(Murder ~ Population + Illiteracy, data = states)
summary(fit)

## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = states)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7652 -1.6561 -0.0898  1.4570  7.6758 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.652e+00  8.101e-01   2.039  0.04713 *  
## Population  2.242e-04  7.984e-05   2.808  0.00724 ** 
## Illiteracy  4.081e+00  5.848e-01   6.978 8.83e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.481 on 47 degrees of freedom
## Multiple R-squared:  0.5668, Adjusted R-squared:  0.5484 
## F-statistic: 30.75 on 2 and 47 DF,  p-value: 2.893e-09

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

par(mfrow=c(1,1))