Chương 11 Hàm hồi quy phi tuyến

Cho đến giờ ta mặc định hàm hồi quy là tuyến tính, nghĩa là tham số nghiêng của hàm hồi quy là hằng số. Điều này hàm ý, ảnh hưởng lên \(Y\) của một đơn vị thay đổi của \(X\) không phụ thuộc vào giá trị của \(X\). Nếu ảnh hưởng này thật sự phụ thuộc vào giá trị của \(X\), ta cần phải sử dụng hàm hồi quy phi tuyến.

11.1 Chiến lược chung

Ta xem xét quan hệ giữa thu nhập sinh viên và điểm thi.

# prepare the data
data(CASchools)
CASchools$size <- CASchools$students/CASchools$teachers
CASchools$score <- (CASchools$read + CASchools$math) / 2
cor(CASchools$income, CASchools$score)
## [1] 0.7124308

Ta thấy có sự tương quan tuyến tính dương ở hai biến này: thu nhập trên trung bình thì điểm thi trên trung bình. Liệu hồi quy tuyến tính có ước lượng được quan hệ dữ liệu.

# fit a simple linear model
linear_model<- lm(score ~ income, data = CASchools)

# plot the observations
plot(CASchools$income, CASchools$score,
     col = "steelblue",
     pch = 20,
     xlab = "District Income (thousands of dollars)", 
     ylab = "Test Score",
     cex.main = 0.9,
     main = "Test Score vs. District Income and a Linear OLS Regression Function")

# add the regression line to the plot
abline(linear_model, 
       col = "red", 
       lwd = 2)

Ta thấy rằng khi thu nhập cao thì đường hồi quy đang ước lượng “lố” quan hệ đúng nhưng khi thu nhập trung bình thì đường hồi quy lại ước lượng “chưa đủ”.

Ta xem xét một mô hình bậc hai như sau.

\[\text{Điểm thi} = \beta_0 + \beta_1 \text{Thu nhập}_i + \beta_2 \text{Thu nhập}_i^2 + u_i\]

Khi đó \(\text{Thu nhập}_i^2\) được dùng như một biến giải thích khác cho điểm thi. Ta xem xét kết quả trong R như sau.

income 3.851***
(0.268)
I(income2) -0.042***
(0.005)
Constant 607.302***
(2.902)
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.

Mô hình này cho phép ta kiểm định quan hệ giữa thu nhập và điểm thi là tuyến tính hay phi tuyến bậc hai. Nói cách khác:

\[\begin{cases} H_0: \beta_2 = 0 \\ H_A: \beta_2 \ne 0 \end{cases}\]

Ta thấy rằng \(H_0\) bị bác bỏ ở bất kỳ mức ý nghĩa thông thường nào, do đó ta kết luận quan hệ giữa hai biến là phi tuyến. Điều này cũng trùng khớp với hình vẽ dưới đây.

# draw a scatterplot of the observations for income and test score
plot(CASchools$income, CASchools$score,
     col  = "steelblue",
     pch = 20,
     xlab = "District Income (thousands of dollars)",
     ylab = "Test Score",
     main = "Estimated Linear and Quadratic Regression Functions")

# add a linear function to the plot
abline(linear_model, col = "black", lwd = 2)

# add quatratic function to the plot
order_id <- order(CASchools$income)

lines(x = CASchools$income[order_id], 
      y = fitted(quadratic_model)[order_id],
      col = "red", 
      lwd = 2) 

11.2 Hàm phi tuyến đối với biến đơn nhất

11.2.1 Hàm đa thức

Từ ý tưởng hàm bậc hai, ta có thể mở rộng ra thành hàm đa thức tổng quát:

\[Y_i = \beta_0 + \beta_1 X_i + \beta_2 X_i^2 + \dots + \beta_r X_i^r + u_i\]

Trong R ta dùng hàm poly() để thể hiện số bậc mô hình. Chẳng hạn với mô hình bậc 3 ta code như sau.

# estimate a cubic model
cubic_model <- lm(score ~ poly(income, degree = 3, raw = TRUE), data = CASchools)

Ta có thể dùng kiểm định thống kê \(F\) để xác nhận mô hình tuyến tính hay phi tuyến đến bậc nào đó. Chẳng hạn, ta quan tâm đến bậc 3 mô hình.

# test the hypothesis of a linear model against quadratic or polynomial
# alternatives

# set up hypothesis matrix
R <- rbind(c(0, 0, 1, 0),
            c(0, 0, 0, 1))

# do the test
linearHypothesis(cubic_model,
                 hypothesis.matrix = R,
                 white.adj = "hc1")
## Linear hypothesis test
## 
## Hypothesis:
## poly(income, degree = 3, raw = TRUE)2 = 0
## poly(income, degree = 3, raw = TRUE)3 = 0
## 
## Model 1: restricted model
## Model 2: score ~ poly(income, degree = 3, raw = TRUE)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F    Pr(>F)    
## 1    418                        
## 2    416  2 37.691 9.043e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ở đây ta đang xem xét giả thuyết \(H_0: \beta_2 = \beta_3 = 0\) bằng cách lợi dụng phép toán ma trận:

\[\begin{aligned} \mathbf{R}\beta &= s \\ \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix} &= \begin{pmatrix} 0 \\ 0 \end{pmatrix} \\ \begin{pmatrix} \beta_2 \\ \beta_3 \end{pmatrix} &= \begin{pmatrix} 0 \\ 0 \end{pmatrix} \end{aligned}\]

Bởi vì hàm linearHypothesis() sử dụng các vector \(0\) nên việc dùng một ma trận \(\mathbf{R}\) sẽ rút ngắn được đoạn code. Ta thấy \(p\)-value khá nhỏ và do đó ta bác bỏ \(H_0\).

Trong thực tế, để xác định được bậc mô hình, ta nên kiểm định \(t\) nhiều lần từ một số bậc lớn nhất nào đó \(r\). Ta có thể thấy trong đoạn code dưới đây.

summary(cubic_model)
## 
## Call:
## lm(formula = score ~ poly(income, degree = 3, raw = TRUE), data = CASchools)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -44.28  -9.21   0.20   8.32  31.16 
## 
## Coefficients:
##                                         Estimate Std. Error t value
## (Intercept)                            6.001e+02  5.830e+00 102.937
## poly(income, degree = 3, raw = TRUE)1  5.019e+00  8.595e-01   5.839
## poly(income, degree = 3, raw = TRUE)2 -9.581e-02  3.736e-02  -2.564
## poly(income, degree = 3, raw = TRUE)3  6.855e-04  4.720e-04   1.452
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## poly(income, degree = 3, raw = TRUE)1 1.06e-08 ***
## poly(income, degree = 3, raw = TRUE)2   0.0107 *  
## poly(income, degree = 3, raw = TRUE)3   0.1471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.71 on 416 degrees of freedom
## Multiple R-squared:  0.5584, Adjusted R-squared:  0.5552 
## F-statistic: 175.4 on 3 and 416 DF,  p-value: < 2.2e-16
# test the hypothesis using robust standard errors
coeftest(cubic_model, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##                                          Estimate  Std. Error  t value
## (Intercept)                            6.0008e+02  5.1021e+00 117.6150
## poly(income, degree = 3, raw = TRUE)1  5.0187e+00  7.0735e-01   7.0950
## poly(income, degree = 3, raw = TRUE)2 -9.5805e-02  2.8954e-02  -3.3089
## poly(income, degree = 3, raw = TRUE)3  6.8549e-04  3.4706e-04   1.9751
##                                        Pr(>|t|)    
## (Intercept)                           < 2.2e-16 ***
## poly(income, degree = 3, raw = TRUE)1 5.606e-12 ***
## poly(income, degree = 3, raw = TRUE)2  0.001018 ** 
## poly(income, degree = 3, raw = TRUE)3  0.048918 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Khi sử dụng ước lượng vững, ta thấy kết quả có sự thay đổi ở mức ý nghĩa của hệ số biến bậc 3. Điều này có nghĩa ta bác bỏ giả thuyết \(H_0\) hàm hồi quy là bậc hai với \(H_A\) hàm hồi quy bậc ba tại mức ý nghĩa \(5\%\).

11.2.2 Giải thích hệ số hồi quy

Chẳng hạn mô hình hồi quy có dạng

\[\hat{\text{Điểm thi}} = 607.3 + 3.85 \times \text{Thu nhập} - 0.0423 \times \text{Thu nhập}^2\]

Như vậy khi tăng thu nhập từ \(10\) lên \(11\) thì điểm thi tăng \(2.96\) điểm, nhưng khi thu nhập tăng từ \(40\) lên \(41\) thì điểm thi chỉ tăng \(0.42\). Cho thấy độ nghiêng của hàm hồi quy dốc hơn ở mức thu nhập thấp và thoải hơn ở mức thu nhập cao.

# compute and assign the quadratic model
quadriatic_model <- lm(score ~ income + I(income^2), data = CASchools)

# set up data for prediction
new_data <- data.frame(income = c(10, 11))

# do the prediction
Y_hat <- predict(quadriatic_model, newdata = new_data)

# compute the difference
diff(Y_hat)
##        2 
## 2.962517

11.2.3 Hàm Logarithms

Một cách khác để quy cách hàm số phi tuyến là dùng hàm lograthims lên biến \(Y\) hoặc/và \(X\). Chuyển hàm logarithms thay đổi các biến sang phần trăm thay đổi. Có nhiều cách để xây dựng mô hình theo cách tiếp cận logarithm.

Trường hợp 1: logarithm \(X\)

Mô hình trở thành:

\[Y_i =\beta_0 + \beta_1 \times \ln(X_i) + u_i\]

# estimate a level-log model
LinearLog_model <- lm(score ~ log(income), data = CASchools)

# compute robust summary
coeftest(LinearLog_model, 
         vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 557.8323     3.8399 145.271 < 2.2e-16 ***
## log(income)  36.4197     1.3969  26.071 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta vẽ đường hồi quy như sau.

# draw a scatterplot
plot(score ~ income, 
     col = "steelblue",
     pch = 20,
     data = CASchools,
     main = "Linear-Log Regression Line")

# add the linear-log regression line
order_id  <- order(CASchools$income)

lines(CASchools$income[order_id],
      fitted(LinearLog_model)[order_id], 
      col = "red", 
      lwd = 2)

Ta có thể giải thích \(\hat{\beta}_1\) như sau. \(1\%\) tăng trưởng thu nhập thì điểm thi tăng \(0.01 \times 36.42 = 0.36\) điểm.

Trường hợp 2: logarithm \(Y\)

Mô hình trong trường hợp này như sau.

\[\ln(Y_i) = \beta_0 + \beta_1 \times X_i + u_i\]

# estimate a log-linear model 
LogLinear_model <- lm(log(score) ~ income, data = CASchools)

# obtain a robust coefficient summary
coeftest(LogLinear_model, 
         vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 6.43936234 0.00289382 2225.210 < 2.2e-16 ***
## income      0.00284407 0.00017509   16.244 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Khi thu nhập tăng 1 đơn vị thu nhập, điểm thi tăng \((100\times 0.00284)\% = 0.284\%\).

Khi dùng hàm chuyển dạng cho biến \(Y\) ta cần cẩn trọng. Mô hình hồi quy sẽ cho ta ước lượng của \(\ln(Y)\), một cách thông thường, ta có thể chuyển về \(\hat{Y}\) bằng hàm exp(). Khi đó, mô hình có dạng:

\[Y_i = e^{\beta_0 + \beta_1 X_i}e^{u_i}\]

Từ đó, ta có:

\[E(Y_i|X_i) = E(e^{\beta_0 + \beta_1 X_i}e^{u_i}|X_i) = e^{\beta_0 + \beta_1 X_i}E(e^{u_i}|X_i)\]

\(E(u_i|X_i)=0\) nên \(E(e^{u_i}|X_i) = 1\), nói cách khác, \(E(e^{u_i}|X_i) \ne 0\). Do đó ước lượng \(\hat{Y} = e^{\hat{\beta}_0 + \hat{\beta_1}X_i}\) sẽ bị thiên lệch do thiếu mất thành tố \(E(e^{u_i}|X_i)\). Hình vẽ sau đây sẽ thể hiện sự thiên lệch này.

# draw a scatterplot
order_id  <- order(CASchools$income)
par(mfrow = c(1,2))
plot(score ~ income, 
     col = "steelblue",
     pch = 20,
     data = CASchools,
     main = "Regression Line with Log Transformation")
abline(linear_model, col = "black", lwd = 2)
# add the linear-log regression line
lines(CASchools$income[order_id],
      exp(fitted(LogLinear_model)[order_id]), 
      col = "red", 
      lwd = 2)
# draw with different y-axis
plot(log(score) ~ income, 
     col = "steelblue", 
     pch = 20, 
     data = CASchools,
     main = "Log-Linear Regression Function")
lines(CASchools$income[order_id], 
      fitted(LogLinear_model)[order_id], 
      col = "red", 
      lwd = 2)

Một cách giải quyết đó là cố gắng ước lượng thành tố thiếu \(E(e^{u_i}|X_i)\). Tuy nhiên nếu \(u_i\) có hiện tượng heteroskedasticity thì việc ước lượng càng phức tạp hơn.

Một cách giải quyết khác đó là giữ nguyên \(\ln(Y)\). Trong tài chính, điều này có thể chấp nhận được khi phân tích diễn biến giá tài sản, nghĩa là phiên bản \(\ln(Y)\) có ý nghĩa kinh tế.

Trường hợp 2: logarithm \(X\)\(Y\)

Hồi quy log-log có dạng như sau

\[\ln(Y_i) = \beta_0 + \beta_1 \ln (X_i) + u_i\]

# estimate the log-log model
LogLog_model <- lm(log(score) ~ log(income), data = CASchools)

# print robust coefficient summary to the console
coeftest(LogLog_model, 
         vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 6.3363494  0.0059246 1069.501 < 2.2e-16 ***
## log(income) 0.0554190  0.0021446   25.841 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Như vậy, \(1\%\) thay đổi thu nhập dẫn đến \(\hat{\beta}_1 \%\) thay đổi điểm số.

# generate a scatterplot
plot(log(score) ~ income, 
     col = "steelblue", 
     pch = 20, 
     data = CASchools,
     main = "Log-Linear Regression Function")

# add the log-linear regression line
order_id  <- order(CASchools$income)

# add the log-log regression line
lines(sort(CASchools$income), 
      fitted(LogLog_model)[order(CASchools$income)], 
      col = "red", 
      lwd = 2)

Ta có thể mở rộng biến đổi logarithms bằng cách kết hợp với đa thức, ta có hàm polylog, chẳng hạn như sau.

\[\text{Điểm thi}_i = \beta_0 + \beta_1 \ln(\text{Thu nhập}_i) + \beta_2 \ln(\text{Thu nhập}_i)^2 + \beta_3 \ln(\text{Thu nhập}_i)^3 + u_i\]

# estimate the polylog model
polyLog_model <- lm(score ~ log(income) + I(log(income)^2) + I(log(income)^3), 
                    data = CASchools)

# print robust summary to the console
coeftest(polyLog_model, 
         vcov = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##                  Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      486.1341    79.3825  6.1239 2.115e-09 ***
## log(income)      113.3820    87.8837  1.2901    0.1977    
## I(log(income)^2) -26.9111    31.7457 -0.8477    0.3971    
## I(log(income)^3)   3.0632     3.7369  0.8197    0.4128    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sau tất cả, ta đánh giá các mô hình đã được xem xét thông qua \(\bar{R}^2\).

# compute the adj. R^2 for the nonlinear models
adj_R2 <-rbind("Quadratic" = summary(quadratic_model)$adj.r.squared,
               "Cubic" = summary(cubic_model)$adj.r.squared,
               "LinearLog" = summary(LinearLog_model)$adj.r.squared,
               "LogLinear" = summary(LogLinear_model)$adj.r.squared,
               "LogLog" = summary(LogLog_model)$adj.r.squared,
               "PolyLog" = summary(polyLog_model)$adj.r.squared)
Adjusted \(R^2\)
Quadratic 0.5540444
Cubic 0.5552279
LinearLog 0.5614605
LogLinear 0.4970106
LogLog 0.5567251
PolyLog 0.5599944

Ta thấy các hệ số \(R^2\) hiệu chỉnh gần xấp xỉ nhau. Ta có thể so sánh đồ thị của các mô hình này, chẳng hạn như sau.

# generate a scatterplot
plot(score ~ income, 
     data = CASchools,
     col = "steelblue", 
     pch = 20,
     main = "Linear-Log and Cubic Regression Functions")

# add the linear-log regression line
order_id  <- order(CASchools$income)

lines(CASchools$income[order_id],
      fitted(LinearLog_model)[order_id], 
      col = "darkgreen", 
      lwd = 2)

# add the cubic regression line
lines(x = CASchools$income[order_id], 
      y = fitted(cubic_model)[order_id],
      col = "darkred", 
      lwd = 2)

Cả hai đường đều gần như nhau. Tuy nhiên ta thấy mô hình linear-log được ưu tiên hơn vì ít tham số hơn, ta không cần những hàm bậc cao.

11.3 Sự tương tác giữa các biến độc lập

Có nhiều vấn đề thực tiễn liên quan đến ảnh hưởng lên biến \(Y\) của sự thay đổi một biến \(X_i\) lại phụ thuộc vào giá trị của một biến \(X_j\) (\(j \ne i\)) khác. Chẳng hạn, liệu các khu vực nhiều học sinh học tiếng Anh có được lợi gì về mặt điểm số không từ việc giảm quy mô lớp học. Thay vì câu hỏi hồi quy thông thường: điểm số sẽ bị tác động như thế nào nếu giảm quy mô lớp học. Để đánh giá vấn đề này, ta cần thêm vào mô hình một sự tương tác giữa các biến \(X\).

Ta xét ba trường hợp:

  • sự tương tác giữa hai biến nhị phân.
  • sự tương tác giữa một biến nhị phân và một biến liên tục.
  • sự tương tác giữa hai biến liên tục.

11.3.1 Sự tương tác giữa hai biến nhị phân

Mô hình có dạng

\[Y_i = \beta_0 + \beta_1 D_{1i} + \beta_2 D_{2i} + u_i\]

Giả sử:

\[\begin{aligned} Y_i &= \ln(\text{Thu nhập})_i \\ D_{1i} &= \begin{cases} 1 \text{ nếu người thứ i có bằng cao đẳng} \\ 0 \end{cases} \\ D_{2i} &= \begin{cases} 1 \text{ nếu người thứ i là nữ} \\ 0 \end{cases} \end{aligned}\]

Ta biết rằng \(\beta_1\) đo lường sự khác biệt trong logarithm thu nhập trung bình giữa hai nhóm người có bằng cao đẳng và người không, còn \(\beta_2\) đo lường sự khác biệt trong logarithm thu nhập trung bình giữa nam và nữ. Tuy nhiên mô hình này không cho phép ta đánh giá sự tác động của một người có đặc điểm của cả \(D_1\)\(D_2\). Do đó ta điều chỉnh mô hình đề xuất thành:

\[Y_i = \beta_0 + \beta_1 D_{1i} + \beta_2 D_{2i} + \beta_3 (D_{1i} \times D_{2i}) + u_i\]

\(D_{1i} \times D_{2i}\) được gọi là biến tương tác. Ta thấy:

\[\begin{aligned} &E(Y_i|D_{1i} = 0, D_{2i}= d_2) = \beta_0 + \beta_2 \times d_2 \\ &E(Y_i|D_{1i} = 1, D_{2i}= d_2) = \beta_0 +\beta_1+\beta_2\times d_2 +\beta_3 \times d_2 \\&E(Y_i|D_{1i} = 1, D_{2i}= d_2) -E(Y_i|D_{1i} = 0, D_{2i}= d_2) = \beta_1+\beta_3 \times d_2 \end{aligned}\]

\[\begin{aligned} &E(Y_i|D_{1i} = d_1, D_{2i}= 0) = \beta_0 + \beta_1 \times d_1 \\ &E(Y_i|D_{1i} = d_1, D_{2i}= 1) = \beta_0 +\beta_1 \times d_1+\beta_2 +\beta_3 \times d_1 \\&E(Y_i|D_{1i} = d_1, D_{2i}= 1) -E(Y_i|D_{1i} = d_1, D_{2i}= 0) = \beta_2+\beta_3 \times d_1 \end{aligned}\]

Theo đó ta sẽ giải thích được ý nghĩa hệ số \(\beta_3\).

Trong R, ta khảo sát sự tương tác giữa biến STRPctEL. Đặt:

\[\begin{aligned} &HiSTR = \begin{cases} 1 \text{ nếu STR } \ge 20 \\ 0 \end{cases} \\ &HiEL = \begin{cases} 1 \text{ nếu PctEL } \ge 10 \\ 0 \end{cases} \end{aligned}\]

# append HiSTR to CASchools
CASchools$HiSTR <- as.numeric(CASchools$size >= 20)

# append HiEL to CASchools
CASchools$HiEL <- as.numeric(CASchools$english >= 10)

Ta ước lượng mô hình.

# estimate the model with a binary interaction term
bi_model <- lm(score ~ HiSTR * HiEL, data = CASchools)

# print a robust summary of the coefficients
coeftest(bi_model, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 664.1433     1.3881 478.4589 < 2.2e-16 ***
## HiSTR        -1.9078     1.9322  -0.9874    0.3240    
## HiEL        -18.3155     2.3340  -7.8472 3.634e-14 ***
## HiSTR:HiEL   -3.2601     3.1189  -1.0453    0.2965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11.3.2 Sự tương tác giữa một biến liên tục và biến nhị phân

Đặt \(X_i\) thể hiện số năm kinh nghiệm của người thứ \(i\), là biến ngẫu nhiên liên tục. Ta chuyển biến \(D_{1i}\) thành biến \(D_i\). Mô hình mới đề xuất:

\[Y_i = \beta_0 + \beta_1 X_i + \beta_2 D_i + u_i\]

Ta thêm vào ảnh hưởng tương tác bằng cách thêm biến \((X_i \times D_i)\) vào mô hình. Điều này sẽ làm độ nghiêng đường hồi quy phụ thuộc vào biến \(D_i\). Có ba trường hợp:

  • cùng độ nghiêng nhưng khác intercept: \(Y_i = \beta_0 + \beta_1 X_i + \beta_2 D_i + u_i\)
  • cùng intercept nhưng khác độ nghiêng: \(Y_i = \beta_0 + \beta_1 X_i + \beta_2 (X_i\times D_i) + u_i\)
  • khác intercept khác độ nghiêng: \(Y_i = \beta_0 + \beta_1 X_i + \beta_2 D_i + \beta_3 (X_i\times D_i) u_i\)
# generate artificial data
set.seed(1)

X <- runif(200,0, 15)
D <- sample(0:1, 200, replace = T)
Y <- 450 +  150 * X + 500 * D + 50 * (X * D) + rnorm(200, sd = 300)

# divide plotting area accordingly
m <- rbind(c(1, 2), c(3, 0))
graphics::layout(m)

# estimate the models and plot the regression lines

# 1. (baseline model)
plot(X, log(Y),
     pch = 20,
     col = "steelblue",
     main = "Different Intercepts, Same Slope")

mod1_coef <- lm(log(Y) ~ X + D)$coefficients

abline(coef = c(mod1_coef[1], mod1_coef[2]), 
       col = "red",
       lwd = 1.5)

abline(coef = c(mod1_coef[1] + mod1_coef[3], mod1_coef[2]), 
       col = "green",
       lwd = 1.5)
       
# 2. (baseline model + interaction term)
plot(X, log(Y),
     pch = 20,
     col = "steelblue",
     main = "Different Intercepts, Different Slopes")

mod2_coef <- lm(log(Y) ~ X + D + X:D)$coefficients

abline(coef = c(mod2_coef[1], mod2_coef[2]), 
       col = "red",
       lwd = 1.5)

abline(coef = c(mod2_coef[1] + mod2_coef[3], mod2_coef[2] + mod2_coef[4]), 
       col = "green",
       lwd = 1.5)

# 3. (omission of D as regressor + interaction term)
plot(X, log(Y),
     pch = 20,
     col = "steelblue",
     main = "Same Intercept, Different Slopes")

mod3_coef <- lm(log(Y) ~ X + X:D)$coefficients

abline(coef = c(mod3_coef[1], mod3_coef[2]), 
       col = "red",
       lwd = 1.5)

abline(coef = c(mod3_coef[1], mod3_coef[2] + mod3_coef[3]), 
       col = "green",
       lwd = 1.5)

Ta áp dụng trong R bằng cách đánh giá tương tác giữa sizeHiEL theo dạng thứ ba.

# estimate the model
bci_model <- lm(score ~ size + HiEL + size * HiEL, data = CASchools)

# print robust summary of coefficients to the console
coeftest(bci_model, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 682.24584   11.86781 57.4871   <2e-16 ***
## size         -0.96846    0.58910 -1.6440   0.1009    
## HiEL          5.63914   19.51456  0.2890   0.7727    
## size:HiEL    -1.27661    0.96692 -1.3203   0.1875    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# identify observations with PctEL >= 10
id <- CASchools$english >= 10

# plot observations with HiEL = 0 as red dots
plot(CASchools$size[!id], CASchools$score[!id],
     xlim = c(0, 27),
     ylim = c(600, 720),
     pch = 20,
     col = "red",
     main = "",
     xlab = "Class Size",
     ylab = "Test Score")

# plot observations with HiEL = 1 as green dots
points(CASchools$size[id], CASchools$score[id],
     pch = 20,
     col = "green")

# read out estimated coefficients of bci_model
coefs <- bci_model$coefficients

# draw the estimated regression line for HiEL = 0
abline(coef = c(coefs[1], coefs[2]),
       col = "red",
       lwd = 1.5)

# draw the estimated regression line for HiEL = 1
abline(coef = c(coefs[1] + coefs[3], coefs[2] + coefs[4]),
       col = "green", 
       lwd = 1.5 )

# add a legend to the plot
legend("topright", 
       pch = c(20, 20), 
       col = c("red", "green"), 
       legend = c("HiEL = 0", "HiEL = 1"))

11.3.3 Sự tương tác giữa hai biến liên tục

Sự tương tác giữa hai biến liên tục \(X_1\)\(X_2\) được đánh giá thông qua \(X_1 \times X_2\). Khi đó, mô hình trở thành:

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 (X_{1i}\times X_{2i}) + u_i\]

Đánh giá đạo hàm theo \(X_1\)\(X_2\) ta sẽ có:

\[\Delta Y_i = (\beta_1 +\beta_3 X_2)\Delta X_1 +(\beta_2 + \beta_3 X_1)\Delta X_2 + \beta_3 \Delta X_1 \Delta X_2\]

Áp dụng trong R, ta đánh giá sizeenglish.

# estimate regression model including the interaction between 'PctEL' and 'size'
cci_model <- lm(score ~ size + english + english * size, data = CASchools) 

# print a summary to the console
coeftest(cci_model, vcov. = vcovHC, type = "HC1")
## 
## t test of coefficients:
## 
##                 Estimate  Std. Error t value Pr(>|t|)    
## (Intercept)  686.3385268  11.7593466 58.3654  < 2e-16 ***
## size          -1.1170184   0.5875136 -1.9013  0.05796 .  
## english       -0.6729119   0.3741231 -1.7986  0.07280 .  
## size:english   0.0011618   0.0185357  0.0627  0.95005    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11.3.4 Phân tích dữ liệu Cầu Economic Journals

Trong phần này, ta phân tích dữ liệu Journals trong package AER, gồm các quan sát của 180 tạp chí khoa học trong năm 2000. Ta đo lường giá tiền cho mỗi citation và tính toán tuổi của journal và số lượng ký tự cho mỗi journal.

data("Journals")
# define and rename variables
Journals$PricePerCitation <- Journals$price/Journals$citations
Journals$Age <- 2000 - Journals$foundingyear
Journals$Characters <- Journals$charpp * Journals$pages/10^6
Journals$Subscriptions <- Journals$subs

Vùng giá trị của PricePerCitation khá rộng và dàn trải.

# compute summary statistics for price per citation
summary(Journals$PricePerCitation)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.005223  0.464495  1.320513  2.548455  3.440171 24.459459

Ta tiến hành quy cách mô hình theo 4 cách, tất cả đều ở dạng Log-Log. Quy cách như vậy sẽ tiện lợi theo cách giải thích ý nghĩa kinh tế của các hệ số về độ co giãn. Nhằm loại bỏ khả năng omitted variable bias, ta quy cách các loại như sau.

  1. \(\ln (Subscriptions_i) = \beta_0 + \beta_1 \ln (PricePerCitation_i) + u_i\)
  2. \(\ln (Subscription_i) = \beta_0 +\beta_1 \ln(PricePerCitation_i) + \beta_4\ln(Age_i) + \beta_6 \ln(Characters_i) + u_i\)
  3. \(\ln (Subscription_i) = \beta_0 +\beta_1 \ln(PricePerCitation_i) + \beta_2 \ln(PricePerCitation_i)^2+\beta_3 \ln(PricePerCitation_i)^3 + \beta_4\ln(Age_i) + \beta_5[\ln(Age_i)\times \ln(PricePerCitation_i)] + \beta_6 \ln(Characters_i) + u_i\)
  4. \(\ln (Subscription_i) = \beta_0 +\beta_1 \ln(PricePerCitation_i) + \beta_4\ln(Age_i) +\beta_5[\ln(Age_i)\times \ln(PricePerCitation_i)]+ \beta_6 \ln(Characters_i) + u_i\)
# Estimate models (I) - (IV)
Journals_mod1 <- lm(log(Subscriptions) ~ log(PricePerCitation), 
                    data = Journals)

Journals_mod2 <- lm(log(Subscriptions) ~ log(PricePerCitation) 
                    + log(Age) + log(Characters), 
                    data = Journals)

Journals_mod3 <- lm(log(Subscriptions) ~ 
                    log(PricePerCitation) + I(log(PricePerCitation)^2) 
                    + I(log(PricePerCitation)^3) + log(Age) 
                    + log(Age):log(PricePerCitation) + log(Characters), 
                    data = Journals)

Journals_mod4 <- lm(log(Subscriptions) ~ 
                    log(PricePerCitation) + log(Age) 
                    + log(Age):log(PricePerCitation) + 
                    log(Characters), 
                    data = Journals)

Sử dụng summary() ta đạt được kết quả sau.

  1. \(\ln (\hat{Subscriptions}_i) = 4.77 - 0.53 \ln (PricePerCitation_i)\)
  2. \(\begin{aligned}\ln (\hat{Subscriptions}_i) = &3.21 - 0.41 \ln(PricePerCitation_i) + 0.42\ln(Age_i) \\&+ 0.21 \ln(Characters_i) \end{aligned}\)
  3. \(\begin{aligned}\ln (\hat{Subscriptions}_i) = &3.41 -0.96 \ln(PricePerCitation_i) + 0.02 \ln(PricePerCitation_i)^2\\&+0.004 \ln(PricePerCitation_i)^3 + 0.37\ln(Age_i) \\&+ 0.16[\ln(Age_i)\times \ln(PricePerCitation_i)] + 0.23 \ln(Characters_i)\end{aligned}\)
  4. \(\begin{aligned}\ln (\hat{Subscriptions}_i) = &3.43 -0.90 \ln(PricePerCitation_i) + 0.37\ln(Age_i) \\&+0.14[\ln(Age_i)\times \ln(PricePerCitation_i)]+ 0.23 \ln(Characters_i)\end{aligned}\)

Ta có thể dùng kiểm định \(F\) để kiểm tra quy cách polylog đối với biến \(\ln(PricePerCitation_i)\).

# F-Test for significance of cubic terms
linearHypothesis(Journals_mod3, 
                 c("I(log(PricePerCitation)^2)=0", "I(log(PricePerCitation)^3)=0"),
                 vcov. = vcovHC, type = "HC1")
## Linear hypothesis test
## 
## Hypothesis:
## I(log(PricePerCitation)^2) = 0
## I(log(PricePerCitation)^3) = 0
## 
## Model 1: restricted model
## Model 2: log(Subscriptions) ~ log(PricePerCitation) + I(log(PricePerCitation)^2) + 
##     I(log(PricePerCitation)^3) + log(Age) + log(Age):log(PricePerCitation) + 
##     log(Characters)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df      F Pr(>F)
## 1    175                 
## 2    173  2 0.1943 0.8236

Ta không thể bác bỏ \(H_0:\beta_3 = \beta_4 = 0\) đối với mô hình (III). Kết quả mô hình được thể hiện trong bảng dưới đây.

log(Subscriptions)
(I) (II) (III) (IV)
(1) (2) (3) (4)
log(PricePerCitation) -0.533*** -0.408*** -0.961*** -0.899***
(0.034) (0.044) (0.160) (0.145)
I(log(PricePerCitation)2) 0.017
(0.025)
I(log(PricePerCitation)3) 0.004
(0.006)
log(Age) 0.424*** 0.373*** 0.374***
(0.119) (0.118) (0.118)
log(Characters) 0.206** 0.235** 0.229**
(0.098) (0.098) (0.096)
log(PricePerCitation):log(Age) 0.156*** 0.141***
(0.052) (0.040)
Constant 4.766*** 3.207*** 3.408*** 3.434***
(0.055) (0.380) (0.374) (0.367)
Observations 180 180 180 180
R2 0.557 0.613 0.635 0.634
Adjusted R2 0.555 0.607 0.622 0.626
Residual Std. Error 0.750 (df = 178) 0.705 (df = 176) 0.691 (df = 173) 0.688 (df = 175)
F Statistic 224.037*** (df = 1; 178) 93.009*** (df = 3; 176) 50.149*** (df = 6; 173) 75.749*** (df = 4; 175)
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.

Hình hoạ các mô hình như sau.

# divide plotting area
m <- rbind(c(1, 2), c(3, 0))
graphics::layout(m)

# scatterplot
plot(Journals$PricePerCitation, 
     Journals$Subscriptions, 
     pch = 20, 
     col = "steelblue",
     ylab = "Subscriptions",
     xlab = "ln(Price per ciation)",
     main = "(a)")

# log-log scatterplot and estimated regression line (I)
plot(log(Journals$PricePerCitation), 
     log(Journals$Subscriptions), 
     pch = 20, 
     col = "steelblue",
     ylab = "ln(Subscriptions)",
     xlab = "ln(Price per ciation)",
     main = "(b)")

abline(Journals_mod1,
       lwd = 1.5)

# log-log scatterplot and regression lines (IV) for Age = 5 and Age = 80
plot(log(Journals$PricePerCitation), 
     log(Journals$Subscriptions), 
     pch = 20, 
     col = "steelblue",
     ylab = "ln(Subscriptions)",
     xlab = "ln(Price per ciation)",
     main = "(c)")

JM4C <-Journals_mod4$coefficients

# Age = 80
abline(coef = c(JM4C[1] + JM4C[3] * log(80), 
                JM4C[2] + JM4C[5] * log(80)),
       col = "darkred",
       lwd = 1.5)

# Age = 5
abline(coef = c(JM4C[1] + JM4C[3] * log(5), 
                JM4C[2] + JM4C[5] * log(5)),
       col = "darkgreen",
       lwd = 1.5)

Những kết luận có thể được rút ra:

  • Cầu của tạp chí co giãn nhiều đối với các tạp chí trẻ tuổi.
  • Việc không bác bỏ \(H_0\) của kiểm định \(F\) đối với mô hình (III) thống nhất với quan hệ tuyến tính giữa log(subscriptions)log(price).
  • Cầu cao hơn đối với tạp chí nhiều ký tự, với giá và tuổi không thay đổi.

Cầu tạp chí phi co giãn với giá: ta thấy mô hình (IV), kể cả tạp chí trẻ tuổi (\(Age=5\)) ta thấy ước lượng độ co giãn giá \(-0.899+0.374\times \ln(5)+ 0.141\times [\ln(1) \times \ln(5)] \approx -0.3\), nghĩa là một phần trăm tăng giá dẫn đến cầu giảm chỉ \(0.3\) phần trăm. Kết quả này không có gì ngạc nhiên vì nguồn tiêu thụ đầu ra của các tạp chí thường là các thư viện.

11.4 Bài tập

11.4.1 Bài 1: Hệ số tương quan và Phi tuyến 1

Xem xét mô hình đơn giản:

\[\hat{medv}_i = 34.554 - 0.95 \times lstat_i\] với \(medv\) là trung vị giá nhà và \(lstat\) là phần trăm hộ gia đình với tình trạng kinh tế thấp, trong bộ dữ liệu Boston.

  1. Tính hệ số tương quan giữa \(medv\)\(lstat\) và lưu vào biến corr.
  2. Đồ thị hoá medvlstat và thêm vào đường hồi quy mod. Nhận xét.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiJkYXRhKFwiQm9zdG9uXCIpXG5tb2QgPC0gbG0obWVkdiB+IGxzdGF0ICwgZGF0YSA9IEJvc3RvbilcbiMgY29tcHV0ZSB0aGUgY29ycmVsYXRpb24gYmV0d2VlbiBtZWR2IGFuZCBsc3RhdFxuXG5cbiMgcGxvdCBtZWR2IGFnYWluc3QgbHN0YXQgYW5kIGFkZCB0aGUgcmVncmVzc2lvbiBsaW5lIiwic29sdXRpb24iOiIjIGNvbXB1dGUgdGhlIGNvcnJlbGF0aW9uIGJldHdlZW4gbWVkdiBhbmQgbHN0YXRcbmNvcnIgPC0gY29yKEJvc3RvbiRtZWR2LCBCb3N0b24kbHN0YXQpXG5cbiMgcGxvdCBtZWR2IGFnYWluc3QgbHN0YXQgYW5kIGFkZCB0aGUgcmVncmVzc2lvbiBsaW5lXG5wbG90KG1lZHYgfiBsc3RhdCwgZGF0YSA9IEJvc3RvbilcbmFibGluZShyZWcgPSBtb2QsIGNvbCA9IFwicmVkXCIpIn0=

11.4.2 Bài 2: Hệ số tương quan và Phi tuyến 2

Ta xem xét mối quan hệ sau đây.

\[medv_i = \beta_0 + \beta_1 \times \log(lstat_i) + u_i\]

  1. Thực hiện hồi quy và lưu trữ vào biến log_mod.
  2. Mô tả điểm rải và thêm vào đường hồi quy. So sánh với kết quả bài trước.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbmR1Y3QgdGhlIHJlZ3Jlc3Npb24gYW5kIGFzc2lnbiBpdCB0byBtb2RfbG9nXG5cblxuIyBkcmF3IGEgc2NhdHRlcnBsb3QgYW5kIGFkZCB0aGUgcmVncmVzc2lvbiBsaW5lIiwic29sdXRpb24iOiIjIGNvbmR1Y3QgdGhlIHJlZ3Jlc3Npb24gYW5kIGFzc2lnbiBpdCB0byBtb2RfbG9nXG5tb2RfbG9nIDwtIGxtKG1lZHYgfiBsb2cobHN0YXQpLCBkYXRhID0gQm9zdG9uKVxuXG4jIGRyYXcgYSBzY2F0dGVycGxvdCBhbmQgYWRkIHRoZSByZWdyZXNzaW9uIGxpbmVcbnBsb3QobWVkdiB+IGxvZyhsc3RhdCksIGRhdGEgPSBCb3N0b24pXG5hYmxpbmUobW9kX2xvZywgY29sID0gXCJyZWRcIikifQ==

11.4.3 Bài 3: Bậc đa thức tối ưu

Ta thấy ở bài tập trước quy cách \(medv_i = \beta_0 + \beta_1 \times \log(lstat_i) + u_i\) là một sự lựa chọn hợp lý. Tuy nhiên, đa thức bậc cao đối với \(\log(lstat_i)\) có thể phù hợp hơn.

Giả sử bậc cao nhất xem xét là \(r=4\), sử dụng for() để chọn ra bậc tối ưu theo cách sau đây.

  1. Ước lượng mô hình, mod, bắt đầu từ \(r=4\).
  2. Lưu biến \(p\)-value vững của các tham số liên quan và so sánh với mức ý nghĩa \(0.05\).
  3. Nếu không thể bác bỏ mô hình, lặp lại bước (i) và (ii) đối với bậc thấp hơn.
  4. Dừng lại cho đến khi chọn ra được bậc tối ưu.

Tính \(R^2\) của mô hình được chọn và phân bổ vào R2.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGZpbmQgdGhlIG9wdGltYWwgcG9seW5vbWlhbCBvcmRlciBvZiB0aGUgcG9seWxvZyBtb2RlbFxuXG5cblxuXG5cblxuXG5cblxuIyBleHRyYWN0IHRoZSBSXjIgZnJvbSB0aGUgc2VsZWN0ZWQgbW9kZWwgYW5kIGFzc2lnbiBpdCB0byBSMiIsInNvbHV0aW9uIjoiIyBmaW5kIHRoZSBvcHRpbWFsIHBvbHlub21pYWwgb3JkZXIgb2YgdGhlIHBvbHlsb2cgbW9kZWxcbmZvcihpIGluIDQ6MSl7XG5tb2QgIDwtIGxtKG1lZHYgfiBwb2x5KGxvZyhsc3RhdCksIGksIHJhdyA9IFQpLCBkYXRhID0gQm9zdG9uKVxucHZhbCA8LSBjb2VmdGVzdChtb2QsIHZjb3YgPSB2Y292SEMpWyhpKzEpLCA0XVxuaWYocHZhbCA8IDAuMDUpe1xuICBwcmludChpKVxuICBicmVha1xuICB9XG59XG5cbiMgZXh0cmFjdCB0aGUgUl4yIGZyb20gdGhlIHNlbGVjdGVkIG1vZGVsIGFuZCBhc3NpZ24gaXQgdG8gUjJcblIyIDwtIHN1bW1hcnkobW9kKSRyLnNxdWFyZWQifQ==

11.4.4 Bài 4: Tương tác giữa các biến độc lập 1

Xem mô hình hồi quy

\[medv_i = \beta_0 + \beta_1 \times chas_i + \beta_2 \times old_i + \beta_3 \times (chas_i \times old_i) + u_i\]

trong đó \(chas_i\)\(old_i\) là các biến giả. Đối với biến đầu, mang giá trị \(1\) nếu sông Charles (một con sông lân cận Boston) đi ngang qua khu vự \(i\), biến sau mang giá trị \(1\) nếu \(age \ge 95\).

  1. Tạo lập biến giả old.
  2. Thực hiện hồi quy trong mod_bb.
  3. Tổng hợp hệ số hồi quy vững.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGdlbmVyYXRlIHRoZSBiaW5hcnkgdmFyaWFibGUgYG9sZGAgYW5kIGFwcGVuZCBpdCB0byB0aGUgZGF0YXNldFxuXG5cbiMgY29uZHVjdCB0aGUgcmVncmVzc2lvbiBhbmQgYXNzaWduIGl0IHRvIGBtb2RfYmJgXG5cblxuIyBwcmludCBhIHJvYnVzdCBzdW1tYXJ5IHRvIHRoZSBjb25zb2xlIiwic29sdXRpb24iOiIjIGdlbmVyYXRlIHRoZSBiaW5hcnkgdmFyaWFibGUgYG9sZGAgYW5kIGFwcGVuZCBpdCB0byB0aGUgZGF0YXNldFxuQm9zdG9uJG9sZCA8LSBhcy5udW1lcmljKEJvc3RvbiRhZ2UgPj0gOTUpXG5cbiMgY29uZHVjdCB0aGUgcmVncmVzc2lvbiBhbmQgYXNzaWduIGl0IHRvIGBtb2RfYmJgXG5tb2RfYmIgPC0gbG0obWVkdiB+IGNoYXMqb2xkLCBkYXRhID0gQm9zdG9uKVxuXG4jIHByaW50IGEgcm9idXN0IHN1bW1hcnkgdG8gdGhlIGNvbnNvbGVcbmNvZWZ0ZXN0KG1vZF9iYiwgdmNvdi4gPSB2Y292SEMpIn0=

11.4.5 Bài 5: Tương tác giữa các biến độc lập 2

Bây giờ xem xét mô hình hồi quy

\[medv_i = \beta_0 + \beta_1 \times indus_i + \beta_2 \times old_i + \beta_3 (indus_i \times old_i) + u_i\]

Dùng ?Boston để tham khảo ý nghĩa các biến trong mô hình. Biến old đã được append vào bộ dữ liệu Boston.

  1. Ước lượng mô hình hồi quy và lưu vào biến mod_bc.
  2. Lấy các hệ số hồi quy lưu vào biến params.
  3. Vẽ đường hồi quy medvindus cho hai trường hợp của old.
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIGNvbmR1Y3QgdGhlIHJlZ3Jlc3Npb24gYW5kIGFzc2lnbiBpdCB0byBtb2RfYmMuXG5cblxuIyBleHRyYWN0IHRoZSBlc3RpbWF0ZWQgbW9kZWwgY29lZmZpY2llbnRzIGFuZCBhc3NpZ24gdGhlbSB0byBwYXJhbXMuXG5cblxuIyBwbG90IG1lZHYgYWdhaW5zdCBpbmR1cyBhbmQgYWRkIHRoZSByZWdyZXNzaW9uIGxpbmVzLiIsInNvbHV0aW9uIjoiIyBjb25kdWN0IHRoZSByZWdyZXNzaW9uIGFuZCBhc3NpZ24gaXQgdG8gbW9kX2JjLlxubW9kX2JjIDwtIGxtKG1lZHYgfiBpbmR1cypvbGQsIGRhdGEgPSBCb3N0b24pXG5cbiMgZXh0cmFjdCB0aGUgZXN0aW1hdGVkIG1vZGVsIGNvZWZmaWNpZW50cyBhbmQgYXNzaWduIHRoZW0gdG8gcGFyYW1zLlxucGFyYW1zIDwtIGNvZWYobW9kX2JjKVxuXG4jIHBsb3QgbWVkdiBhZ2FpbnN0IGluZHVzIGFuZCBhZGQgdGhlIHJlZ3Jlc3Npb24gbGluZXMuXG5wbG90KG1lZHYgfiBpbmR1cywgZGF0YSA9IEJvc3RvbilcbmFibGluZShhID0gcGFyYW1zWzFdLCBiID0gcGFyYW1zWzJdLCBjb2wgPSBcInJlZFwiKVxuYWJsaW5lKGEgPSBwYXJhbXNbMV0gKyBwYXJhbXNbM10sIGIgPSBwYXJhbXNbMl0gKyBwYXJhbXNbNF0sIGNvbCA9IFwiZGFya2JsdWVcIikifQ==