Sử dụng package Gamlss

Bài 6: Thành phần đa thức

1 Giới thiệu

Một thế giới phi tuyến tính

Trong hình dưới đây, tuy mọi thứ có vẻ lộn xộn khó hiểu, nhưng thật ra Nhi chỉ muốn chuyển đến các bạn một thông điệp đơn giản : thế giới mà chúng ta đang sống là một thế giới phi tuyến tính. Dòng chảy của một con sông, quỹ đạo chuyển động của một vật thể, của ánh sáng, âm nhạc, quá trình tăng trưởng của em bé và sự thoái hóa của người già, nhịp sinh học, giá chứng khoán … đều biểu hiện thành những đường cong. Do đó, ta có thể hình dung rằng mình không thể chỉ dùng duy nhất mô hình hồi quy tuyến tính để mô tả về thế giới này.

Hồi quy đa thức bậc cao trong gamlss

Chúng ta không quên mục tiêu của hành trình, đó là khám phá package gamlss. Chặng thứ 6 hôm nay ta sẽ tìm hiểu về thành phần đa thức trong mô hình GAMLSS. Hồi quy đa thức là chức năng thú vị nhất mà package gamlss cung cấp, vì những chiêu thức trong đó cho phép ta mô hình hóa được những quy luật phi tuyến tính.

Có 4 thành phần đa thức Nhi sẽ trình bày, theo thứ tự về độ tinh tế của chúng. Đó là hàm đa thức cơ bản (basis polynomial), phương pháp trực giao (orthogonal polynomial), fractional polynomial, phương pháp phân đoạn (Piecewise polynomial) bao gồm truncated polynomial và B-spline piecewise polynomial. Phương pháp đi sau mạnh, tinh tế hơn trước, và hoàn thiện cho phương pháp đi trước.

2 Mô hình tăng trưởng BMI ở trẻ em

Ta sẽ bắt đầu bằng một bài toán cụ thể ; mục tiêu là dựng một mô hình mô tả tăng trưởng của chỉ số BMI ở trẻ em từ 0-2 tuổi. Dataset nằm ngay trong package gamlss và có nguồn gốc từ Hà Lan.

library(tidyverse)
library(gamlss)

data(dbbmi)

data=subset(dbbmi,age<2)

Đồ thị cho thấy giữa BMI (biến kết quả) và Tuổi (predictor) có một quan hệ phi tuyến tính.

ggplot(data,aes(x=age*12,y=bmi))+
  geom_point(alpha=0.05,color="red",size=2)+
  theme_bw()+
  geom_smooth(fill="#f74c4c",color="red4",alpha=0.5,size=1)+
  scale_x_continuous("Age (month)")+
  scale_y_continuous("BMI")

Thói quen thường gặp là ta thử áp dụng một số hàm chuyển dạng cho X và/hoặc Y. Nhưng bằng một phép thử đơn giản ta sẽ thấy ngay là không có cách chuyển dạng nào cho phép đưa quan hệ BMI và tuổi về tuyến tính.

page=c("age^(-0.5)", "log(age)", "age^(0.5)", "age")
pbmi=c("bmi^(-0.5)", "log(bmi+1)", "bmi^(0.5)")

op <- par(mfrow = c(3, 4), mar = par("mar") + c(0, 1, 0, 0),
          pch = "+", cex = 0.45, cex.lab = 2, cex.axis = 1.6)

for (i in 1:3) {
yy <- with(data,eval(parse(text = pbmi[i])))
for (j in 1:4) {
xx <- with(data, eval(parse(text = page[j])))

plot(yy ~ xx, xlab = page[j], ylab = pbmi[i],col="blue")
}
}

par(op)

3 Hàm đa thức cơ bản

Mô hình hồi quy đa thức bậc cao là một giải pháp cho phép tăng khả năng (capacity) thích ứng với dữ liệu của mô hình tuyến tính. Đa thức bậc P có dạng :

\[ Y ~ f(x)=\beta _{0}+\beta _{1}x+\beta _{2}x^2+...\beta _{p}x^p \]

Thành phần đa thức bậc p có thể được đưa vào matrix của model bằng nhiều cách, đơn giản nhất là sử dụng hàm đa thức cơ bản I(x^p)

Thí dụ một mô hình bậc 4 sử dụng hàm đa thức cơ bản

mbp=gamlss(bmi~age+                                  I(age^2)+I(age^3)+I(age^4),data=data)%>%summary()

## GAMLSS-RS iteration 1: Global Deviance = 6309.12 
## GAMLSS-RS iteration 2: Global Deviance = 6309.12 
## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = bmi ~ age + I(age^2) + I(age^3) +  
##     I(age^4), data = data) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.4744     0.1498  89.948  < 2e-16 ***
## age          14.8251     1.1681  12.691  < 2e-16 ***
## I(age^2)    -18.7890     2.4567  -7.648 3.26e-14 ***
## I(age^3)      9.7747     1.9068   5.126 3.26e-07 ***
## I(age^4)     -1.8792     0.4896  -3.838 0.000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27615    0.01639   16.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1861 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  1855 
##                       at cycle:  2 
##  
## Global Deviance:     6309.12 
##             AIC:     6321.12 
##             SBC:     6354.293 
## ******************************************************************

Trong thí dụ sau đây, Nhi dùng hàm I(age^p) với p từ 2 đến 6 để tạo ra một mô hình đa thức bậc 6. Mỗi thành phần được cộng thêm là 1 hàm đa thức cơ bản (basis polynomial).

Hình sau đây mô tả 6 hàm đa thức cơ bản cho biến số Age, như ta thấy, những hàm đa thức cơ bản không hữu hiệu, vì kết quả của chúng không thể được kiểm soát, kết quả này có thể rất cao đến mức phi lý (BMI không thể cao đến 40 hay 60 ?) cho những giá trị Age cao nhất và ngược lại, vô cùng thấp cho những giá trị Age thấp nhất trên thang đo.

polybasis<- with(data,
          model.matrix(formula(~age+
                                      I(age^2)+I(age^3)+I(age^4)+I(age^5)+I(age^6)
                               )
                       )
          )[,-1]%>%as.data.frame()

names(polybasis)=c("1","2","3","4","5","6")

polybasis%>%mutate(.,Age=.$`1`)%>%
  gather(`1`:`6`,key="Degree",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Degree))+
  geom_smooth()+
  theme_bw()+
  ggtitle("Basis polynomial functions")

Do đó, thay vì ta sử dụng trực tiếp hàm đa thức cơ bản, ta dùng phương pháp trực giao (orthogonal) qua hàm poly(x,p) với p là bậc đa thức, thí dụ poly(x,3) sẽ tạo ra đa thức bậc 3

4 Đa thức trực giao (Orthogonal polynomial)

Cơ chế của phương pháp orthogonal đó là nó sẽ tạo ra một hàm tuyến tính với các biến số là kết quả của hàm đa thức cơ bản, và áp dụng trọng số khác nhau tùy theo tham số của hàm tuyến tính. Vi lý do này, mặc dù đồ thị của hàm orthogonal polynomial là một đường cong, thực chất mô hình đa thức trực giao chính là 1 mô hình tuyến tính và mỗi thành phần đa thức bậc p bên trong mô hình có thể được xem như một thành phần tuyến tính cộng thêm.

Hình vẽ sau đây trình bày đặc tính của 6 thành phần đa thức được đưa vào mô hình bậc 6 qua hàm poly().

poly=with(data,
                 model.matrix(formula(~age+
                                        poly(age,6)
                 )
                 )
)[,-1]%>%as.data.frame()

names(poly)=c("Age","1","2","3","4","5","6")

poly%>%gather(`1`:`6`,key="Poly",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Poly))+
  geom_smooth()+
  theme_bw()

poly%>%gather(`1`:`6`,key="Poly",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Poly))+
  geom_smooth()+
  facet_wrap(~Poly,ncol=2)

Bây giờ Ta sẽ sử dụng gamlss để thăm dò 5 mô hình đa thức bậc 2-6, dùng hàm poly (trực giao).

Kết quả BIC cho thấy mô hình bậc 4 là phù hợp nhất.

mp2=gamlss(bmi~poly(age,2),data=data)

## GAMLSS-RS iteration 1: Global Deviance = 6425.805 
## GAMLSS-RS iteration 2: Global Deviance = 6425.805

mp3=gamlss(bmi~poly(age,3),data=data)

## GAMLSS-RS iteration 1: Global Deviance = 6323.795 
## GAMLSS-RS iteration 2: Global Deviance = 6323.795

mp4=gamlss(bmi~poly(age,4),data=data)

## GAMLSS-RS iteration 1: Global Deviance = 6309.12 
## GAMLSS-RS iteration 2: Global Deviance = 6309.12

mp5=gamlss(bmi~poly(age,5),data=data)

## GAMLSS-RS iteration 1: Global Deviance = 6306.534 
## GAMLSS-RS iteration 2: Global Deviance = 6306.534

mp6=gamlss(bmi~poly(age,6),data=data)

## GAMLSS-RS iteration 1: Global Deviance = 6306.444 
## GAMLSS-RS iteration 2: Global Deviance = 6306.444

GAIC(mp2,mp3,mp4,mp5,mp6,k=log(nrow(data)))

##     df      AIC
## mp4  6 6354.293
## mp5  7 6359.236
## mp3  5 6361.440
## mp6  8 6366.675
## mp2  4 6455.920

Ta có thể xem nội dung bên trong mô hình.Các bạn lưu ý rằng nội dung mô hình với hàm poly hoàn toàn khác với mô hình với hàm I() cơ bản:

summary(mp4)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = bmi ~ poly(age, 4), data = data) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    16.57191    0.03055 542.394  < 2e-16 ***
## poly(age, 4)1  23.15257    1.31805  17.566  < 2e-16 ***
## poly(age, 4)2 -30.18798    1.31805 -22.904  < 2e-16 ***
## poly(age, 4)3  13.55008    1.31805  10.280  < 2e-16 ***
## poly(age, 4)4  -5.05919    1.31805  -3.838 0.000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27615    0.01639   16.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1861 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  1855 
##                       at cycle:  2 
##  
## Global Deviance:     6309.12 
##             AIC:     6321.12 
##             SBC:     6354.293 
## ******************************************************************

Nếu ta vẽ đồ thị của tất cả các mô hình từ bậc 1 (tuyến tính) đến bậc 4, ta có thể thấy: trong mô hình đa thức trực giao, chỉ có 2 thành phần thực sự quan trọng (quyết định giá trị kết quả), đó là hằng số (đường thẳng nằm ngang) và thành phần tuyến tính của hàm đa thức cơ bản (basis), các hàm basis bậc cao hơn chỉ đóng vai trò hiệu chỉnh (uốn cong) với hiệu ứng rất nhỏ.

Ta cũng có thể thấy rằng mô hình đa thức trực giao bậc 4 rất phù hợp với dữ liệu:

b=coef(mp4)
P<-model.matrix(with(data,formula(~poly(age,4))))

Fit=(rep(b,nrow(P)) *t(P))%>%t(.)%>%as.data.frame()
names(Fit)=c("Intercept","1","2","3","4")
Fit=Fit%>%mutate(.,Fitted=fitted(mp4),Age=data$age,Truth=data$bmi)

Fit%>%gather(Intercept:Fitted,key="Poly",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Poly))+
  geom_smooth()+
  theme_bw()

## Warning: attributes are not identical across measure variables; they will
## be dropped

## `geom_smooth()` using method = 'gam'

Fit%>%gather(Intercept,Fitted,Truth,key="Poly",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Poly))+
  geom_smooth(fill="gold",aes(linetype=Poly))+
  theme_bw()+scale_color_manual(values=c("red","#9727f9","black"))

## Warning: attributes are not identical across measure variables; they will
## be dropped

## `geom_smooth()` using method = 'gam'

5 Fractional polynomial

Khái niệm fractional polynomial do Royston và Altman giới thiệu năm 1994. Mô hình đa thức phân số chỉ có 3 cấp:

\[ Y ~ f(x)=\beta _{0}+\beta _{1}x^p1+\beta _{2}x^p2+\beta _{3}x^p3 \]

Trong đó mỗi tham số fractional polynomial x^p có thể nhận bất cứ giá trị nào trong một tập quy ước bao gồm (-2, -1, -0.5, 0, 0.5, 1,2 và 3). Nếu có 2 tham số và chúng giống nhau, tham số thứ 2 sẽ được nhân cho log(x) \[ \beta _{1j}x^pj ; \beta _{2j}x^pj.log(x) \]

Tương tự nếu có 3 tham số giống nhau, mô hình sẽ gồm:

\[\beta _{1j}x^pj ; \beta _{2j}x^pj.log(x) \ và \ \beta _{3j}x^pj.log(x)^2\]

Gamlss cho phép dựng mô hình fractional polynomial với hàm fp(x,npoly=1,2 hoặc 3)

Nhi thăm dò 3 mô hình fractional polynomial. Kết quả cho thấy mô hình fractional bậc 2 là phù hợp nhất, thậm chí nó tốt hơn mô hình orthogonal bậc 4

m1f=gamlss(bmi~fp(age,1), sigma.fo=~1, data=data, trace=FALSE)
m2f=gamlss(bmi~fp(age,2), sigma.fo=~1, data=data, trace=FALSE)
m3f=gamlss(bmi~fp(age,3), sigma.fo=~1, data=data, trace=FALSE)

GAIC(m1f,m2f,m3f,k=log(nrow(data)))

##     df      AIC
## m2f  6 6354.218
## m3f  8 6364.829
## m1f  4 6488.368

summary(m2f)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = bmi ~ fp(age, 2), sigma.formula = ~1,  
##     data = data, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.57191    0.03055   542.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27613    0.01639   16.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  1861 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  1855 
##                       at cycle:  2 
##  
## Global Deviance:     6309.045 
##             AIC:     6321.045 
##             SBC:     6354.218 
## ******************************************************************

getSmo(m2f)

## 
## Call:
## lm(formula = y ~ x.fp, weights = w)
## 
## Coefficients:
## (Intercept)        x.fp1        x.fp2  
##      -8.507        9.333       -5.134

GAIC(mp4,m1f,m2f,m3f,k=log(nrow(data)))

##     df      AIC
## m2f  6 6354.218
## mp4  6 6354.293
## m3f  8 6364.829
## m1f  4 6488.368

data%>%mutate(.,Truth=.$bmi,fpol2=fitted(m2f),fpol3=fitted(m3f),poly4=fitted(mp4),fpol1=fitted(m1f))%>%
  gather(.,Truth:fpol1,key="Model",value="BMI")%>%
  ggplot()+
  geom_point(aes(x=age,y=bmi),color="grey",alpha=0.1)+
  geom_smooth(aes(x=age,y=BMI,color=Model),se=F)+
  theme_bw()

## Warning: attributes are not identical across measure variables; they will
## be dropped

## `geom_smooth()` using method = 'gam'

Hình sau đây trình bày các thành phần đa thức fractional cho biến số Age

td=data%>%mutate(.,`0`=log(.$age),
                    `1`=(.$age),
                 `-2`=(.$age)^(-2),
                 `-1`=(.$age)^(-1),
                 `-0.5`=(.$age)^(-0.5),
                 `0.5`=(.$age)^0.5,
                 `2`=(.$age)^2,
                 `3`=(.$age)^3,
                 )
                        
td%>%gather(`0`:`3`,key="Poly",value="Value")%>%
  ggplot()+
  geom_smooth(aes(x=age,y=Value,color=Poly),se=F)+
  theme_bw()+scale_y_continuous(limits=c(0,24))

6 Đa thức phân đoạn cơ bản

Ý tưởng mô hình đa thức phân đoạn khởi nguồn từ việc quan sát hiện tượng thay đổi khuynh hướng đột ngột, rõ rệt giá trị biến kết quả Y tại 1 số vị trí nhất định của X. Hiện tượng này tạo ra sự đứt gãy trên đồ thị, và vị trí X được gọi là nút (knots) hay điểm gãy (breakpoints).

Giả sử ta chia thang đo X thành 2 đoạn trước và sau 1 knot, và dựng 2 mô hình tuyến tính trên 2 đoạn bị chặn này, 2 đoạn thẳng sẽ gắn kết với nhau tại vị trí knot tạo ra 1 đồ thị hình chữ V

Bây giờ nếu ta sử dụng nhiều knot hơn và chia thang đo X thành nhiều đoạn hơn, và làm tương tự, ta sẽ thấy những phân đoạn đồ thị của các mô hình tuyến tính sẽ nối với nhau và tạo ra ảo giác về một đồ thị của mô hình phi tuyến tính.

td=data%>%mutate(.,t1=ifelse(.$age>0.0 &.$age<=0.25,1,0),
                 t2=ifelse(.$age>0.25 & .$age<=0.5,1,0),
                 t3=ifelse(.$age>0.5 & .$age<=0.75,1,0),
                 t4=ifelse(.$age>0.75,1,0),
                 )

td%>%gather(t1:t4,key="breakpoint",value="val")%>%
  mutate(.,response=.$age*val)%>%
  mutate(.,response=ifelse(.$response>0,.$response,NA))%>%
  na.omit(.)%>%
  ggplot()+
  geom_smooth(method="lm",se=F,aes(x=age,y=bmi*val,color=breakpoint))+
  geom_smooth(se=F,aes(x=age,y=bmi),color="black",linetype=2)+
  theme_bw()+
  geom_vline(xintercept=c(0.25,0.5,0.75),linetype=2)+
  scale_color_manual(values=c("orange","red","blue","purple"))

## `geom_smooth()` using method = 'gam'

Bây giờ thay vì dựng mô hình tuyến tính, ta dựng mô hình đa thức bậc cao, thí dụ p=2, ta sẽ có những đường cong nối với nhau. Mô hình có dạng:

\[ Y ~ f(x)=\beta _{00} +\beta _{01}x+ (\beta _{10} +\beta _{11}(x-b)).H(x-b) \]

Trong đó b là giá trị của X tại điểm gãy (knot, breakpoint)

H là hàm Heaviside với 2 giá trị : H(t) =1 nếu t>=0 và H(t)=0 nếu t<0

Hàm đa thức phân đoạn bậc D với K knots được viết tổng quát như sau :

\[ f(x)=\sum_{j=0}^{D}\beta _{0j}x^j+\sum_{k=1}^{K}\sum_{j=0}^{D}\beta _{kj}(x-b_{k})^jH(x-b_{k}) \]

td=data%>%mutate(.,t1=ifelse(.$age>0.0 &.$age<=2,1,0),
                 t2=ifelse(.$age>0.25 & .$age<=2,1,0),
                 t3=ifelse(.$age>0.5 & .$age<=2,1,0),
                 t4=ifelse(.$age>0.75&.$age<=2,1,0),
)

td=td%>%gather(t1:t4,key="breakpoint",value="val")%>%
  mutate(.,`0`=val,
         `1`=age*val,
         `2`=(age+age^2)*val,
         `3`=(age+age^2+age^3)*val,
         `4`=(age+age^2+age^3+age^4)*val
         )

td%>%gather(`0`:`4`,key="Poly",value="Value")%>%
  mutate(.,Value=ifelse(.$Value>0,.$Value,NA))%>%
  na.omit(.)%>%
  ggplot()+
  geom_smooth(aes(x=age,y=Value,color=breakpoint),se=F)+
  theme_bw()+
  facet_wrap(~Poly,ncol=2,scales="free")

## `geom_smooth()` using method = 'gam'

Mô hình đa thức phân đoạn cơ bản (hồi quy đa thức trên các khoảng chặn) linh động hơn so với 2 phương pháp trực giao hoặc phân số, tuy nhiên nó cũng chịu chung nhược điểm với hàm đa thức cơ bản đó là tính bất ổn định. Do đó, Nhi sẽ bước qua phương pháp B-splines như sau đây :

7 B-Spline piece wise polynomial

Việc áp dụng B-splines cho phương pháp đa thức phân đoạn trên khoảng chặn cũng có ý nghĩa tương tự như khi tác động của đa thức trực giao so với hàm đa thức cơ bản. B-spline là chữ viết tắt của basis spline.

B-spline là một thành phần đa thức nhưng đồng thời cũng là 1 thành phần bù trừ (hiệu chỉnh) trong mô hình gamlss. B-spline được xác định bởi các hàm mang tính cục bộ có miền xác định trong khoảng 2+D knots. Thí dụ hàm spline bậc 3 (D=3) có các hàm cơ bản được xác định trong giới hạn 2+3=5 knots.

Tùy theo bậc của đa thức phân đoạn (giá trị D), ta có :

Hằng số cục bộ (D=0)

2 hàm tuyến tính cục bộ cho x (D=1)

3 hàm bậc 2 cục bộ (D=2)

4 hàm bậc 3 cục bộ (D=3)

D+1 hàm đa thức cục bộ (D cao hơn hoặc = 4)

Hình sau đây trình bày đồ thị của 4 hàm B-spline cục bộ cho D=4

library(splines)

Pmx=bs(data$age, 
     df=4, 
     intercept=T)%>%
  as.data.frame()

Pmx%>%mutate(.,Age=data$age)%>%gather(`1`:`4`,key="Element",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Element))+
  geom_smooth()+
  theme_bw()

## `geom_smooth()` using method = 'gam'

Ta có thể thấy rằng vị trí các knots cho hàm B-spline không bắt buộc phải cách đều nhau, cho thấy phương pháp này cục kì linh động trong việc bám sát dữ liệu (đây là 1 điều tốt, trừ khi dữ liệu có outliers…). Số knots sẽ xác định số lượng hàm Bspline cơ bản.

m=gamlss(bmi~bs(data$age, 
                df=4, 
                intercept=T)-1,
         data=data,
         trace=F)

P=model.matrix(with(data,
                      formula(~bs(data$age, 
                                  df=4,
                                  intercept=T)-1)
                    ))
b=coef(m)
Fit=t(rep(b,nrow(data))*t(P))%>%as.data.frame()
names(Fit)=c("1","2","3","4")

Fit=Fit%>%mutate(.,Fitted=fitted(m),BMI=data$bmi,Age=data$age)

knots<-c(1, attr(Pmx, "knots"),nrow(Pmx))

Fit%>%gather(`1`:`4`,BMI,Fitted,key="Element",value="Value")%>%
  ggplot(aes(x=Age,y=Value,color=Element))+
  geom_smooth()+
  theme_bw()+scale_y_continuous(limits=c(0,24))

## Warning: attributes are not identical across measure variables; they will
## be dropped

## `geom_smooth()` using method = 'gam'

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_smooth).

Mỗi hàm Bspline cục bộ là 1 cột trong matrix tạm gọi là B, ta có thể dùng matrix B này như một tập biến trong mô hình hồi quy tuyến tính có dạng

\[ \hat{y} = \hat{\beta }B \]

Tham số hồi quy “hat beta” có thể được diễn giải như một mối quan hệ tuyến tính « linh động » giữa biến kết quả Y và biến độc lập X tại các khoảng cục bộ khác nhau.

Khi dùng hàm summary() cho 1 mô hình đa thức phân đoạn Bspline, thí dụ bậc 3, ta sẽ có kết quả của hệ số hồi quy beta cho Intercept, và 3 hàm Bspline cục bộ của X.

Sau đây Nhi sẽ thăm dò 5 mô hình bậc 1,3,4,5,6 sử dụng B-spline. Kết quả BIC cho thấy mô hình bậc 4 là thích hợp nhất. Đây cũng là mô hình tối ưu nhất từ đầu đến giờ.

mpw1=gamlss(bmi ~ bs(age), data = data, trace = FALSE)
mpw3=gamlss(bmi ~ bs(age,df=3), data = data, trace = FALSE)
mpw4=gamlss(bmi ~ bs(age,df=4), data = data, trace = FALSE)
mpw5=gamlss(bmi ~ bs(age,df=5), data = data, trace = FALSE)
mpw6=gamlss(bmi ~ bs(age,df=6), data = data, trace = FALSE)

GAIC(mpw1,mpw3,mpw4,mpw5,mpw6,m2f,k=log(nrow(data)))

##      df      AIC
## mpw4  6 6353.208
## m2f   6 6354.218
## mpw5  7 6356.679
## mpw1  5 6361.440
## mpw3  5 6361.440
## mpw6  8 6363.885

data%>%mutate(.,Truth=.$bmi,
              fpol2=fitted(m2f),
              pwPol4=fitted(mpw4),
              poly4=fitted(mp4)
              )%>%
  gather(.,Truth:poly4,key="Model",value="BMI")%>%
  ggplot()+
  geom_point(aes(x=age,y=bmi),color="grey",alpha=0.2)+
  geom_smooth(aes(x=age,y=BMI,color=Model),se=F)+
  theme_bw()

summary(mpw4)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = bmi ~ bs(age, df = 4), data = data,  
##     trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       13.8457     0.1277  108.39   <2e-16 ***
## bs(age, df = 4)1   3.3480     0.2660   12.59   <2e-16 ***
## bs(age, df = 4)2   4.1413     0.2370   17.48   <2e-16 ***
## bs(age, df = 4)3   2.8756     0.3160    9.10   <2e-16 ***
## bs(age, df = 4)4   2.4744     0.2294   10.79   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27586    0.01639   16.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1861 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  1855 
##                       at cycle:  2 
##  
## Global Deviance:     6308.034 
##             AIC:     6320.034 
##             SBC:     6353.208 
## ******************************************************************

Ta có thể áp dụng B-spline cho bất cứ họ phân phối nào. Trong thí dụ dưới đây, sau khi thăm dò họ phân phối của BMI ta thấy T family cho ra kết quả tốt hơn so với Gaussian.

scan=fitDist(bmi,data=data,type="realAll",k=log(nrow(data)),parallel="multicore",ncpus = 4)

scan$fits

##        TF        LO        GT        PE       ST3       ST4       ST2 
##  7086.504  7091.743  7092.074  7092.708  7092.807  7093.724  7093.839 
##       ST5       ST1      BCTo       JSU      EGB2    exGAUS      SEP1 
##  7093.924  7093.983  7093.986  7094.157  7094.464  7094.865  7097.386 
##      SEP2    SHASHo      SEP3     BCPEo      SEP4       SN1     BCCGo 
##  7098.103  7099.077  7099.344  7100.182  7100.235  7100.532  7101.386 
##       SN2        GG        GA       GIG     LOGNO        IG    IGAMMA 
##  7101.451  7102.686  7105.440  7112.969  7121.784  7122.838  7145.882 
##      WEI3        RG        GU       EXP   PARETO2 
##  7327.650  7422.185  7574.382 14179.821 14187.352

scan

## 
## Family:  c("TF", "t Family") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = y, family = DIST[i], parallel = "multicore",  
##     ncpus = 4, data = sys.parent()) 
## 
## Mu Coefficients:
## [1]  16.57
## Sigma Coefficients:
## [1]  0.415
## Nu Coefficients:
## [1]  2.765
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   1858 
## Global Deviance:     7063.92 
##             AIC:     7069.92 
##             SBC:     7086.5

gamlss.family(TF)

## 
## GAMLSS Family: TF t Family 
## Link function for mu   : identity 
## Link function for sigma: log 
## Link function for nu   : log

library(viridis)

set.seed(1234)

p1=rTF(2000, mu = 16.57, sigma =exp(0.415), nu=exp(2.765))%>%
  as_tibble()%>%
  ggplot(aes(x=value))+
  geom_histogram(aes(y=..density..,fill=..density..,color=..density..), alpha=0.7,show.legend = F)+
  geom_density(alpha=.1,color="black",fill="gold",linetype=2)+
  theme_bw()+
  scale_fill_viridis(discrete=F,option="C",direction = -1)+
  scale_color_viridis(discrete=F,option="C",direction = -1)+
  ggtitle("Simulated TF (mu=16.57,s=exp(0.415),nu=exp(2.765))")

p2=rNO(2000, mu = 16.57, sigma =1.62)%>%
  as_tibble()%>%
  ggplot(aes(x=value))+
  geom_histogram(aes(y=..density..,fill=..density..,color=..density..), alpha=0.7,show.legend = F)+
  geom_density(alpha=.1,color="black",fill="gold",linetype=2)+
  theme_bw()+
  scale_fill_viridis(discrete=F,option="C",direction = -1)+
  scale_color_viridis(discrete=F,option="C",direction = -1)+
  ggtitle("Simulated Gaussian(mu=16.57,s=1.62)")

p3=data%>%ggplot(aes(x=bmi))+
  geom_histogram(aes(y=..density..,fill=..density..,color=..density..), alpha=0.7,show.legend = F)+
  geom_density(alpha=.1,color="black",fill="gold",linetype=2)+
  theme_bw()+
  scale_fill_viridis(discrete=F,option="C",direction = -1)+
  scale_color_viridis(discrete=F,option="C",direction = -1)+
  ggtitle("True observations")

gridExtra::grid.arrange(p1,p2,p3,ncol=1)

Mô hình B-spline bậc 4 với phân phối T family, lưu ý là phân phối TF có đến 3 tham số gồm Mu, Sigma và Nu

mpw4tf=gamlss(bmi ~ bs(age,df=4), data = data, trace = FALSE,family=TF)

summary(mpw4tf)

## ******************************************************************
## Family:  c("TF", "t Family") 
## 
## Call:  gamlss(formula = bmi ~ bs(age, df = 4), family = TF,  
##     data = data, trace = FALSE) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       13.8018     0.1239 111.410   <2e-16 ***
## bs(age, df = 4)1   3.3729     0.2584  13.055   <2e-16 ***
## bs(age, df = 4)2   4.1989     0.2324  18.069   <2e-16 ***
## bs(age, df = 4)3   2.8320     0.3077   9.203   <2e-16 ***
## bs(age, df = 4)4   2.5281     0.2222  11.377   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.16839    0.02473   6.808 1.33e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  log 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.3785     0.2006   11.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1861 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  1854 
##                       at cycle:  5 
##  
## Global Deviance:     6260.921 
##             AIC:     6274.921 
##             SBC:     6313.623 
## ******************************************************************

GAIC(mpw4,mpw4tf)

##        df      AIC
## mpw4tf  7 6274.921
## mpw4    6 6320.034

library(viridisLite)
library(gamlss.util)

## Warning: package 'gamlss.util' was built under R version 3.4.1

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

plotSimpleGamlss(y=bmi,x=age,formula=bmi ~ bs(age,df=4),family=TF(),data=data,cols=viridis(option="A",n=100),x.val=c(0.25,0.5,0.75,1,1.25,1.5,1.75,2),N=1000)

## Warning in predict.gamlss(object, newdata = newdata, what = "mu", type = type, : There is a discrepancy  between the original and the re-fit 
##  used to achieve 'safe' predictions 
## 

8 Kết luận:

Các thành phần đa thức bậc cao là giải pháp hữu hiệu để cải thiện khả năng của mô hình. Ta có ít nhất là 3 lựa chọn là orthogonal, fractional và B-spline piece wise, thứ tự tăng dần theo độ mạnh và tinh tế.

Tuy nhiên cần cẩn trọng vì mô hình đa thức bậc cao dễ bị overfitting, và khó diễn giải nên không phù hợp cho nghiên cứu diễn dịch và dataset đa biến hoặc có xuất hiện outlier.

Xin cảm ơn và hẹn gặp lại.