Hồi quy phi tuyến tính là gì
|
xDuLieu ⮞Dữ liệu đa biến ⮞Tương quan & Hồi quy ⮞Hồi quy đa biến phi tuyến Show
Hồi quy đa biến phi tuyếnTrong thực tế, phần lớn tương quan giữa biến phụ thuộc và biến độc lập là phi tuyến. Mặc dù vậy, mô hình tuyến tính vẫn được áp dụng rộng rãi vì đơn giản, dễ xử lý; hơn nữa nó vẫn đáp ứng nhiều yêu cầu, hữu ích khi sử dụng. Tuy nhiên, trong một số trường hợp, mô hình này không đáp ứng được yêu cầu của thực tế. Khi ấy chúng ta phải sử dụng các mối tương quan phi tuyến. Trong phần này, chúng ta đề cập đến một số phương pháp giải quyết vấn đề này. Mở rộng mô hình tuyến tínhNguyên tắcMột trong những xu hướng phổ biến là mở rộng mô hình tuyến tính. Mô hình tuyến tính tổng quát (generalized linear model hay GLM) có dạng tương tự như mô hình tuyến tính, nghĩa là tương quan giữa biến phụ thuộc `Y` và các biến độc lập `X_i` vẫn được trình bày dưới dạng:
Tuy nhiên ở đây `X_i` có thể là :
Ghi chú : Mô hình tuyến tính tổng quát chỉ áp dụng khi có sự tuyến tính với các hệ số `b_i`, nghĩa là ta có thể áp dụng cho `Y=b_i e^X`, nhưng không thể áp dụng cho `Y=e^(b_i X)` vì khi ấy `b_i` không còn tuyến tính nữa so với `Y`. Thí dụĐể khảo sát ảnh hưởng của vận tốc tàu chở công ten nơ đến chi phí nhiên liệu sử dụng, người ta ghi nhận hai thông số trên cho một số tàu có trọng tải 8000 công ten nơ quy chuẩn (TEU: 20-foot equivalent unit). Kết quả ghi nhận như Bảng 1 sau (tập tin chi-phi-nhien-lieu.csv): Bảng 1 Vận tốc và chi phí nhiên liệu của tàu chở công ten nơ
Vận tốc : hải lý ; Chi phí nhiên liệu : tấn/ngày Mô hình tuyến tínhĐể xem xét mối tương quan giữa hai đại lượng trên ta dùng R. Sau khi nhập dữ liệu này vào R với tên bảng dữ liệu là nhl, ta khảo sát sơ bộ bảng này bằng đoạn lệnh: library(ggplot2) ggplot(nhl, aes(Van_Toc, CP_NL)) + geom_point(color = "red") + labs(x = "Vận tốc (hải lý)", y = "Chi phí nhiên liệu (tấn/ngày)")Ta thu được kết quả như trên Hình 1. Hình 1 Tương quan giữa vận tốc và chi phí nhiên liệu Qua Hình 1, ta thấy mối tương quan giữa hai đại lượng này không hoàn toàn tuyến tính. Để tìm hiểu kỹ hơn, ta vẽ biểu đồ phần dư khi dùng mô hình tuyến tính với đoạn lệnh: kqa <- lm(CP_NL ~ Van_Toc, data = nhl)pdu <- resid(kqa) pdu <- resid(kqa) cp <- nhl$CP_NL dfpdu <- data.frame(cp, pdu) ggplot(dfpdu, aes(cp, pdu)) + geom_point() + geom_abline(intercept = 0, slope = 0, linetype = "64", color = "blue", size = 1) + labs(x = "Chi phí nhiên liệu (tấn/ngày)", y = "Sai lệch (tấn/ngày)") ta thu được kết quả trên Hình 2. Hình 2 Biểu đồ phần dư cho mô hình tuyến tính Qua Hình 2, ta thấy phần dư không phân bố ngẫu nhiên mà theo một quy luật nào đó. Những điểm trên cho thấy sử dụng mô hình tuyến tính cho trường hợp này chưa hợp lý. Ta khảo sát một số mô hình phi tuyến như sau. Mô hình đa thứcTrong mô hình đa thức (polynomial), tương quan giữa biến độc lập `X` và biến phụ thuộc `Y` được biểu diễn bằng phương trình:
trong đó `m` là bậc của đa thức. Để xác định bậc của mối tương quan trong thí dụ trên, ta thử với đa thức bậc 5 trong đoạn lệnh sau: kqb <- lm(CP_NL ~ poly(Van_Toc, 5), data = nhl)summary(kqb) Kết quả thu được là : > summary(kqb) Call: lm(formula = CP_NL ~ poly(Van_Toc, 5), data = nhl) Residuals: Min 1Q Median 3Q Max -2.8268 -1.3719 0.1089 0.9478 5.2570 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 107.3250 0.2870 373.917 <2e-16 *** poly(Van_Toc, 5)1 182.0215 1.8153 100.269 <2e-16 *** poly(Van_Toc, 5)2 31.2746 1.8153 17.228 <2e-16 *** poly(Van_Toc, 5)3 0.2782 1.8153 0.153 0.879 poly(Van_Toc, 5)4 -0.1359 1.8153 -0.075 0.941 poly(Van_Toc, 5)5 0.8313 1.8153 0.458 0.650 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 1.815 on 34 degrees of freedom Multiple R-squared: 0.9967, Adjusted R-squared: 0.9962 F-statistic: 2070 on 5 and 34 DF, p-value: < 2.2e-16Kết quả này cho thấy đa thức chỉ có nghĩa đến bậc 2, các thừa số của bậc lớn hơn 2 không có nghĩa. Ta cũng lưu ý là kết quả này chưa cho ta viết được phương trình hồi quy. Để thực hiện việc này, ta sử dụng đoạn lệnh sau: kqc <- lm(CP_NL ~ Van_Toc + I(Van_Toc^2), data = nhl)summary(kqc) Kết quả nhận được là : > summary(kqc) Call: lm(formula = CP_NL ~ Van_Toc + I(Van_Toc^2), data = nhl) Residuals: Min 1Q Median 3Q Max -2.9960 -1.4686 0.0916 1.0673 4.8852 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 318.50904 30.19048 10.55 1.04e-12 *** Van_Toc -39.30734 3.17406 -12.38 9.97e-15 *** I(Van_Toc^2) 1.48630 0.08299 17.91 < 2e-16 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 1.746 on 37 degrees of freedom Multiple R-squared: 0.9967, Adjusted R-squared: 0.9965 F-statistic: 5593 on 2 and 37 DF, p-value: < 2.2e-16Dựa vào đây, ta có thể biểu diễn mối tương quan giữa chi phí nhiên liệu (CP_NL) và vận tốc (Van_Toc) bằng phương trình hồi quy sau: CP_NL = 318,50904 39,30734 Van_Toc + 1,48630 Van_Toc 2 Phương trình này có độ tương thích rất cao với dữ liệu thực tế với `R^2=0,9967` Để vẽ biểu đồ cho mô hình đa thức ta sử dụng đoạn lệnh sau: vt <- nhl$Van_Toc cpmh <- predict(kqc) dfmh <- data.frame(vt, cpmh) ggplot(data = nhl, aes(x=Van_Toc, y=CP_NL)) + geom_point(color = "red", size = 2) + geom_line(data = dfmh, aes(x = vt, y = cpmh), color = "blue", size = 1) + labs(x = "Vận tốc (hải lý)", y = "Chi phí nhiên liệu (tấn/ngày)"Trong đoạn lệnh trên, 2 dòng đầu tạo hai biến là vt (vận tốc) và cpmh (chi phí nhiên liệu xác định theo mô hình đa thức). Dòng thứ ba tạo bảng dữ liệu dfmh từ hai biến trên. Các dòng tiếp theo dùng để vẽ các điểm và đường cong hồi quy. Kết quả được thể hiện trên Hình 3. Hình 3 Tương quan giữa chi phí nhiên liệu và vận tốc theo mô hình đa thức bậc 2 Nếu ta thử nghiệm với đa thức bậc cao hơn, 3 chẳng hạn, bằng đoạn lệnh: kqd <- lm(CP_NL ~ Van_Toc + I(Van_Toc^2) + I(Van_Toc^3), data = nhl)summary(kqd) Ta thu được kết quả sau : > summary(kqd) Call: lm(formula = CP_NL ~ Van_Toc + I(Van_Toc^2) + I(Van_Toc^3), data = nhl) Residuals: Min 1Q Median 3Q Max -2.9512 -1.4639 0.0607 1.0957 4.8917 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 273.15587 290.13902 0.941 0.353 Van_Toc -32.17510 45.48673 -0.707 0.484 I(Van_Toc^2) 1.11477 2.36500 0.471 0.640 I(Van_Toc^3) 0.00641 0.04078 0.157 0.876 Residual standard error: 1.77 on 36 degrees of freedom Multiple R-squared: 0.9967, Adjusted R-squared: 0.9964 F-statistic: 3630 on 3 and 36 DF, p-value: < 2.2e-16Ta thấy tất cả các hệ số của phương trình hồi quy đều không có ý nghĩa thống kê. Mô hình hàm mũTrong mô hình hàm mũ, tương quan giữa biến độc lập `X` và biến phụ thuộc `Y` được biểu diễn bằng phương trình: `Y=a e^(bX) + epsilon`(3) Bỏ qua sai lệch `epsilon`, ta có thể viết lại là : `ln Y=ln a +bX+epsilon`(4) Ghi chú : Trong R, `log X` dùng để chỉ logarit tự nhiên, `log10 X` dùng để chỉ logarit thập phân. Vì thế khi dùng ký hiệu của R, phương trình (4) viết lại là `log Y=log a + bX + epsilon`(5) Vậy để xác định phương trình hồi quy theo mô hình hàm mũ của thí dụ, ta sử dụng đoạn lệnh sau: kqe <- lm(log(CP_NL) ~ Van_Toc, data = nhl)summary(kqe) Kết quả nhận được là : > summary(kqe) Call: lm(formula = log(CP_NL) ~ Van_Toc, data = nhl) Residuals: Min 1Q Median 3Q Max -0.03775 -0.01106 -0.00051 0.01316 0.03418 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.77457 0.03004 59.08 <2e-16 *** Van_Toc 0.15332 0.00160 95.85 <2e-16 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 0.01668 on 38 degrees of freedom Multiple R-squared: 0.9959, Adjusted R-squared: 0.9958 F-statistic: 9187 on 1 and 38 DF, p-value: < 2.2e-16Từ kết quả này, mô hình hồi quy có phương trình là: log Y = 1,77457 + 0,15332 Van_Toc + ε hay : Y = 5,897745 e0,15332 Van_Toc + ε Phương trình này có độ tương thích rất cao với dữ liệu thực tế với `R^2=0,9959`. Mô hình lũy thừaTrong mô hình lũy thừa, tương quan giữa biến độc lập `X` và biến phụ thuộc `Y` được biểu diễn bằng phương trình: `Y=aX^b +e`(6) Bỏ qua sai lệch `e`, ta có thể viết lại là : `ln Y=ln a + b ln X + e`(7) Khi sử dụng quy định của R thì (17) viết lại là: `log Y=log a + b log X + e `(8) Vậy để xác định phương trình hồi quy theo mô hình lũy thừa của thí dụ, ta sử dụng đoạn lệnh sau: kqf <- lm(log(CP_NL) ~ log(Van_Toc), data = nhl)summary(kqf) Và kết quả thu được là : > summary(kqf) Call: lm(formula = log(CP_NL) ~ log(Van_Toc), data = nhl) Residuals: Min 1Q Median 3Q Max -0.043797 -0.014458 -0.000699 0.011756 0.057970 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.83671 0.12913 -29.71 <2e-16 *** log(Van_Toc) 2.89887 0.04413 65.69 <2e-16 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 0.02428 on 38 degrees of freedom Multiple R-squared: 0.9913, Adjusted R-squared: 0.991 F-statistic: 4316 on 1 and 38 DF, p-value: < 2.2e-16Dựa vào kết quả này, mô hình lũy thừa của thí dụ được viết như sau: log Y = 3,83671 + 2,89887 log Van_Toc + e hay : Y = 0,021564 X2,89887 Van_Toc + e Phương trình này có độ tương thích rất cao với dữ liệu thực tế với `R^2=0,9913`. Hồi quy tuyến tính nhiều đoạnMột số đặc điểmMô hình hồi quy tuyến tính nhiều đoạn (piecewise hay multiphase) có một số đặc điểm sau:
Mô hình hồi quy tuyến tính nhiều đoạn có thể biểu diễn bằng công thức sau: `Y=b_0+b_1X+d_1(X-C_1)_+ + d_2(X-C_2)_+ +...+d_K(X-C_K)_+ + e`(9) hay :
trong đó `K` là số điểm nút, `C_i` là giá trị của biến độc lập `X` tại điểm nút `i`, còn "hàm" `(x-C_i)_+` được định nghĩa như sau:
Hàm `(x-C_i)_+` còn có thể biểu diễn dưới dạng khác là `max(0, x-C_i)` và cũng được gọi là hàm bản lề (hinge function). Thí dụTrong một công ty, người ta nhận thấy giữa thời gian thực tập (TG) và năng suất (NS) của công nhân có mối quan hệ như được trình bày ở Bảng 2. Bảng 2 Thời gian thực tập và năng suất của công nhân
TG : thời gian thực tập (giờ) ; NS : năng suất (sản phẩm/ngày) Mô hình hồi quy tuyến tínhTa sử dụng R để lập mô hình hồi quy tuyến tính cho dữ liệu trên bằng đoạn lệnh sau: TG <- c(3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30) NS <- c(150,165, 155,160, 145,170, 205,265, 340,405, 465,520, 585,590, 585,595, 590) ns <- data.frame(TG, NS) kq1 <- lm(NS ~ TG, data = ns) summary(kq1)Trong đoạn lệnh trên, hai dòng đầu tạo hai biến TG và NS cùng với các giá trị của chúng; dòng thứ ba tạo bảng dữ liệu từ hai biến trên; dòng thứ tư dùng lệnh lm để phân tích hồi quy tuyến tính với NS là biến phụ thuộc và TG là biến độc lập; kết quả phân tích này được lưu vào biến kq1; và ở dòng thứ năm, ta dùng lệnh summary để xem xét nội dung kết quả này. Ta có: > summary(kq1) Call: lm(formula = NS ~ TG, data = ns) Residuals: Min 1Q Median 3Q Max -81.877 -35.384 0.071 38.285 79.097 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49.475 22.466 2.202 0.0437 * TG 20.747 1.305 15.899 8.51e-11 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 46.57 on 15 degrees of freedom Multiple R-squared: 0.944, Adjusted R-squared: 0.9402 F-statistic: 252.8 on 1 and 15 DF, p-value: 8.511e-11Kết quả này cho thấy mô hình hồi quy tương thích tốt với số liệu thực tế với `R^2=0,944`. Tuy nhiên, nếu ta vẽ biểu đồ xy cho hai biến trên bằng đoạn lệnh: library(ggplot2) ggplot(ns, aes(x=TG, y=NS)) + geom_point(colour="red", size=I(3)) + labs(x="Thời gian thực tập (giờ)", y="Năng suất (sản phẩm/ngày)") + xlim(0, 30) + ylim(100,600)thì ta thu được Hình 4. Hình 4 Sự phụ thuộc của năng suất vào thời gian thực tập. Qua Hình 4, ta thấy năng suất phụ thuộc vào thời gian thực tập theo một số dạng khác nhau. Khi thời gian thực tập nhỏ (khoảng dưới 8 giờ) hay lớn (khoảng trên 22 giờ) thì năng suất gần như không bị ảnh hưởng vào thời gian thực tập. Nhưng khi thời gian thực tập trung bình (khoảng 8 đến 22 giờ) năng suất phụ thuộc đáng kể vào biến độc lập này. Ngoài ra, khi xem xét biểu đồ phần dư bằng đoạn lệnh: nspre <- predict(kq1) nsres <- resid(kq1) pdu <- data.frame(nspre, nsres) ggplot(pdu, aes(x=nspre, y=nsres)) + geom_point(colour="red", size=I(2)) + geom_abline(intercept = 0, slope = 0, colour="blue", linetype="55", size=1) + labs(x="Năng suất (sản phẩm/giờ)", y="Sai lệch (sản phẩm/giờ)")ta thu được Hình 5. Hình 5 Biểu đồ phần dư của năng suất theo mô hình tuyến tính Ta thấy phần dư không phân phối một cách ngẫu nhiên mà theo quy luật nhất định. Những điều trên cho thấy mô hình tuyến tính không phản ánh tốt tình hình thực tế. Quan sát Hình 4, ta thấy rằng trong ba vùng của biến độc lập, ta có thể sử dụng mô hình tuyến tính cho từng vùng. Hồi quy tuyến tính nhiều đoạnĐể hỗ trợ cho việc tìm mô hình tuyến tính trong nhiều đoạn, R có phụ kiện "segmented". Sau khi tải và cài đặt phụ kiện này, ta thực hiện đoạn lệnh sau: library(segmented) nsseg <- segmented(kq1, seg.Z = ~ TG, psi = c(6, 18)) summary(nsseg)Trong đoạn lệnh trên, ta dùng lệnh segmented để thực hiện hồi quy nhiều đoạn từ kết quả hồi quy tuyến tính thông thường kq1. Sự phân đoạn dựa vào đối số TG: seg.Z = ~ TG. Các điểm nút dự kiến ở vào khoảng 6 giờ và 22 giờ được đưa vào đối số psi. Hai giá trị này là điểm xuất phát cho quá trình dò tìm điểm nút của hàm segmented. Kết quả hồi quy được lưu vào biến nsseg và dùng hàm summary để xem nội dung. Kết quả xử lý cho ta: > summary(nsseg) ***Regression Model with Segmented Relationship(s)*** Call: segmented.lm(obj = kq1, seg.Z = ~TG, psi = c(6, 18)) Estimated Break-Point(s): Est. St.Err psi1.TG 8.588 0.293 psi2.TG 21.851 0.289 Meaningful coefficients of the linear terms: Estimate Std. Error t value Pr(>|t|) (Intercept) 150.429 10.247 14.681 1.43e-08 *** TG 1.286 1.779 0.723 0.485 U1.TG 30.714 1.989 15.440 NA U2.TG -31.250 1.475 -21.182 NA --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 7.443 on 11 degrees of freedom Multiple R-Squared: 0.999, Adjusted R-squared: 0.9985 Convergence attained in 3 iterations with relative change -1.678982e-15Kết quả này cho ta thấy :
Mô hình nhiều đoạn này được biểu diễn trên Hình 6 Hình 6 Mô hình tuyến tính nhiều đoạn Ghi chú : Trong mô hình này, ta thấy hệ số góc ở hai bên của điểm nút không bằng nhau, nghĩa là vi phân cấp 1 ở hai bên điểm này không bằng nhau. Như vậy mô hình này chỉ liên tục ở cấp 0 (còn được ký hiệu là `C^0`). Nếu ta xem xét phần dư của mô hình này bằng đoạn lệnh nspre2 <- predict(nsseg) nsres2 <- resid(nsseg) pdu2 <- data.frame(nspre2, nsres2) ggplot(pdu2, aes(x=nspre2, y=nsres2)) + geom_point(colour="red", size=I(2)) + geom_abline(intercept = 0, slope = 0, colour="blue", linetype="55", size=1) + labs(x="Năng suất (sản phẩm/giờ)", y="Sai lệch (sản phẩm/giờ)")ta thu được Hình 7 Hình 7 Biểu đồ phần dư của năng suất theo mô hình tuyến tính nhiều đoạn Khi so sánh với Hình 5, sự phân phối phần dư của mô hình mới có tính ngẫu nhiên cao hơn. So sánh hai mô hìnhKhi so sánh giá trị của phần dư trong hai mô hình, ta có: > resid(kq1) 1 2 3 4 5 6 38.28484345 32.53809920 1.79135495 -13.95538930 -49.70213356 -45.44887781 7 8 9 10 11 12 -51.94236631 -33.43585481 0.07065669 23.57716819 42.08367969 55.59019119 13 14 15 16 17 79.09670269 42.60321419 -3.89027431 -35.38376282 -81.87725132 > resid(nsseg) 1 2 3 4 5 6 7 -4.285714 9.428571 -1.857143 1.857143 -14.428571 9.285714 -1.666667 8 9 10 11 12 13 14 -5.666667 5.333333 6.333333 2.333333 -6.666667 -1.000000 2.500000 15 16 17 -4.000000 4.500000 -2.000000Hay bằng Hình 8 sau : Hình 8 Biểu đồ so sánh phần dư của năng suất theo hai mô hình Qua những so sánh trên, ta thấy mô hình nhiều đoạn có phần dư nhỏ hơn, phân phối có tính ngẫu nhiên cao hơn, nên tính tương thích và độ tin cậy cũng cao hơn. Hồi quy bằng đường cong khớpCũng như hồi quy tuyến tính nhiều đoạn, ta cũng chia vùng biến thiên của biến độc lập làm một số đoạn. Nhưng ở đây, trong mỗi đoạn, mối quan hệ giữa biến độc lập và biến phụ thuộc là phi tuyến và được biểu diễn bằng một đoạn cong. Toàn mô hình được biểu diễn bằng một đường cong khớp (spline) là tổng hợp các đoạn cong trên. Các mô hình thuộc nhóm này được ứng dụng rất rộng rãi trong lĩnh vực đồ họa. Để việc khảo sát và trình bày được đơn giản, trong phần sau chúng ta chỉ xem xét trường hợp có một biến độc lập. Đường cong khớp dạng đa thứcTrong mô hình đường cong khớp đa thức, phương trình hồi quy được trình bày dưới dạng:
trong đó `m` là bậc của đa thức, `K` là số điểm nút, `C_i` là giá trị của biến độc lập tại điểm nút, hàm `(x-C_i)_+^m` là hàm ngắt (truncated) được định nghĩa như sau:
Thí dụ mô hình đường cong khớp bậc 2 có 3 điểm nút có phương trình như sau: `Y=b_0+b_1X+b_2X^2+d_1(X-C_1)_+^2 + d_2(X-C_2)_+^2 + d_3(X-C_3)_+^2`(14) Trong các phương trình (12) và (14), `b_j` là các hệ số của đa thức ban đầu, `d_i` là các hệ số của hàm ngắt. Trong hồi quy bằng đường cong khớp, ta phải xác định các hệ số này. Trong các mô hình hồi quy này, dạng đa thức bậc 3 là dạng thường dùng nhất vì số hệ số không nhiều, vẫn đảm bảo sự liên tục đến cấp 2 (`C^2`). Ngoài ra, mô hình tuyến tính nhiều đoạn như ta vừa khảo sát ở trên là một trường hợp riêng của đường cong khớp đa thức với bậc 1. Hàm cơ sở Một hướng khác để xem xét mô hình đường cong khớp là xem mô hình là tổ hợp tuyến tính của một số hàm cơ sở `h_i(X)` và trình bày lại dưới dạng:
Theo cách trình bày này, mô hình có :
Lấy thí dụ đường cong khớp bậc 2 với 3 điểm nút, mô hình có:
Đường cong khớp BMô hình đường cong khớp đa thức có một khuyết điểm lớn là có sự khác biệt đáng kể (không ổn định) khi tính toán. Chỉ cần một khác biệt nhỏ trong số liệu, thí dụ thay đổi số chữ số khi làm tròn cũng có thể dẫn tới sự khác biệt đáng kể về mô hình. Điều này đặc biệt nghiêm trọng đổi với những vùng biên của biến độc lập (giá trị nhỏ nhất hay lớn nhất). Một trong những cách khắc phục sự không ổn định này là thay mô hình đa thức bằng mô hình B-spline. Mô hình B-spline cho `K` điểm nút có dạng như sau:
trong đó `B_i(X)` là các hàm cơ sở. Việc xác định các hàm này tương đối phức tạp, không được đề cập ở đây. Đường cong khớp tự nhiênĐể khắc phục sự không ổn định, ta có thể sử dụng mô hình đường cong khớp tự nhiên (natural spline). Đây là mô hình dạng đa thức với một số điều chỉnh sau:
Nhìn chung, mức độ ổn định của mô hình đường cong khớp tự nhiên cao hơn so với mô hình B-spline, vì vậy được sử dụng khá rộng rãi. Các mô hình đường cong khớp tự nhiên phổ biến nhất có bậc 3. Thí dụNhiệt độ trung bình của tỉnh T trong 730 ngày được ghi nhận trong tập tin nhiet-do.csv. Ta dùng R để phân tích xu hướng thay đổi nhiệt độ trung bình của tỉnh này trong 730 ngày. R có phụ kiện spline, hỗ trợ chúng ta thực hiện phân tích một số mô hình đường cong khớp. Sau khi nhập số liệu vào R thành bảng dữ liệu ndo với hai biến là Ngay và Nhiet_Do, ta trình bày xu hướng thay đổi nhiệt độ của tỉnh T bằng đoạn lệnh sau. library(splines) library(ggplot2) ggplot(ndo, aes(Ngay, Nhiet_Do)) + stat_smooth(method = "lm", formula = y ~ ns(x,10), color = "red") + geom_point(size = 0.8) + xlab("Ngày") + ylab("Nhiệt độ (C)")Theo đoạn lệnh trên, mô hình hồi quy bằng đương cong khớp tự nhiên được xây dựng bằng hàm ns của phụ kiện splines với 10 điểm nút. Kết quả được thể hiện trên Hình 9. Hình 9 Biểu đồ nhiệt độ tỉnh T |
