GAMLSS
Pham Dinh Khanh
May 18, 2018
1 GAMLSS model
GAMLSS model (Generalized Addictive Model for Location, Scale and Shape (GAMLSS)) là model hồi qui phân phối thuộc lớp các model học có giám sát. Model được xây dựng dựa trên giả định về phân phối của biến được dự báo sẽ phụ thuộc vào các tham số biến đổi theo biến giải thích sử dụng phương pháp nào trong các phương pháp linear, non-linear hoặc smooth function. Trong model GAMLSS giả định rằng phân phối của biến được dự báo (response variable hoặc target variable) sẽ có tham số phân phối đăc trưng mô phỏng được các hình dạng phân phối bất kì chẳng hạn như: Đuôi dày, đuôi mỏng, đỉnh nhọn, đỉnh dẹp, phân phối lệch trái hay lệch phải,…. và các phân phối này có thể được mô hình hóa theo các dạng hàm linear, non-linear, smoothing function mà tả sẽ trình bày bên dưới về phương trình hồi qui của các tham số theo biến giải thích (hoặc biến độc lập).
Trái với các model GLM và GAM khi dựa trên giả thuyết về họ các hàm phân phối dựa trên số mũ đã được thay thế bởi mội giả thuyết lỏng và linh hoạt hơn đó là họ phân phối tổng quát dựa trên hệ số skew, kurtosis liên tục hoặc rời rạc. Giả định về hàm mật độ xác xuất (density function) của GAM sẽ cho các quan sát độc lập \(y_{i}\) với \(i = \overline{1,n}\) có dạng: \[f(y_{i}| \mu_{i}, \sigma_{i}, \nu_{i}, \tau_{i})\] trong điều kiện phân phối có tham số đặc trưng là \((\mu_{i}, \sigma_{i}, \nu_{i}, \tau_{i})\). Trong đó \(\mu_{i}, \sigma_{i}\) là 2 tham số đặc tả lần lượt cho location (vị trí hoặc trung bình) và scale (độ rộng hoặc phương sai). Các tham số còn lại đặc tả về hình dạng như \(\nu_{i}\) là skewness (độ lệch) và \(\tau_{i}\) là kurtosis (độ nhọn). Với 4 tham số này model có thể mô tả được phân phối đặc trưng của hầu hết các hình dạng phân phối. Các tham số này được ước lượng theo biến độc lập như sau:
\[{\begin{aligned}g_{1}(\mu )=\eta _{1}=X_{1}\beta _{1}+\sum _{{j=1}}^{{J_{1}}}{h}_{{j1}}(x_{{j1}})\\g_{2}(\sigma )=\eta _{2}=X_{2}\beta _{2}+\sum _{{j=1}}^{{J_{2}}}{h}_{{j2}}(x_{{j2}})\\g_{3}(\nu )=\eta _{3}=X_{3}\beta _{3}+\sum _{{j=1}}^{{J_{3}}}{h}_{{j3}}(x_{{j3}})\\g_{4}(\tau )=\eta _{4}=X_{4}\beta _{4}+\sum _{{j=1}}^{{J_{4}}}{h}_{{j4}}(x_{{j4}})\end{aligned}}\] Trong đó :
\(\eta_{i}\) là các vector ước lượng của các tham số có độ dài bằng n. Mỗi một phần tử trong vector này đại diện cho giá trị của tham số phân phối ứng với một quan sát.
\(\beta_{k}^{T}=(\beta_{1k},\beta_{2k},\ldots ,\beta_{J'_{k}k})\) là vector các hệ số ước lượng trong phương trình hồi qui các tham số theo các biến giải thích. Vector này có độ dài là \(J'_{k}\).
\(X_{k}\) là matrix các biến giải thích có kích thước \(n X J'_{k}\). \(X_{k}\beta_{k}\) đại diện cho thành phần hồi qui có hệ thống.
\(h_{jk}(x_{jk})\) là hàm smoothing phi tham số của các biến giải thích \(x_{jk}\). Thành phần này đại diện cho thành phần hồi qui phi hệ thống. Thành phần này giúp các tham số phân phối giải thích được các biến động không theo qui luật.
Như vậy ta có thể thấy GAMLSS sẽ xây dựng các phương trình hồi qui độc lập nhau để dự báo các tham số phân phối gồm \(\mu, \sigma, \nu, \tau\) ứng với mỗi quan sát. Từ các ước lượng này sẽ tìm ra được qui luật phân phối của biến được dự báo và dự báo chúng.
2 Ứng dụng GAMLSS
Bộ dữ liệu được sử dụng trong ví dụ này là abdom, đây là bộ dữ liệu đo lường chu vi bất thường từ bào thai trong độ tuổi từ 12 đến 42 tuần. Trong đó:
y: Chu vi bào thai.
x: Độ tuổi của thai nhi.
Cú pháp model
gamlss(formula = formula(data), sigma.formula = ~1,
nu.formula = ~1, tau.formula = ~1, family = NO(),
data = sys.parent(), weights = NULL,
contrasts = NULL, method = RS(), start.from = NULL,
mu.start = NULL, sigma.start = NULL,
nu.start = NULL, tau.start = NULL,
mu.fix = FALSE, sigma.fix = FALSE, nu.fix = FALSE,
tau.fix = FALSE, control = gamlss.control(...),
i.control = glim.control(...), ...)Diễn giải arguments
formula: Là phương trình hồi qui. Phương trình hồi qui có thể sử dụng các phương pháp nonparametrics smoothing terms như penalised beta splines, smoothing splines, loess smoothing và random terms thông qua các hàm lần lượt là pb(), cs(), lo() và ra(). Chẳng hạn như ta có một phương trình hồi qui
y~cs(x,df=5) + x1+x2+x3thì có nghĩa rằng nhân tử được sử dụng để smoothing là đa thức bậc 5 của x theo phương pháp smoothing splines. Những phần tử làm mịn có thể được thêm vào để tạo ra một mặt phẳng fitted với các biến predictor với response.sigma.formula: Là phương trình sử dụng để fitting model cho tham số \(\sigma\). Chẳng hạn chúng ta có thể viết:
sigma.fo = ~cs(x,df=5).nu.formula: Phương trình dùng để fitting model đối với tham số nu. Chẳng hạn
nu.fo = ~x.tau.formula: Phương trình dùng để fitting model đối với tham số tau. Chẳng hạn
tau.fo = ~cs(x,df=12).family: Là một gamlss.family object, được sử dụng để xác định phân phối và hàm liên kết cửa rất nhiều các tham số khác nhau.
data: Dữ liệu sử dụng để training.
weights: Là vector trọng số sử dụng để đánh trọng số trong phân tích likelihood.
contrasts: Danh sách các so sánh được sử dụng cho các nhân tố xuất hiện như là biến trong phương trình hồi qui.
method: Algorithms được sử dụng để hồi qui model. Có thể là RS(), CG() hoặc mixed(). Trong đó RS() là phương pháp Rigby and Stasinopoulos. CG() là phương pháp Cole and Green. và mixed(2,10) sẽ sử dụng cả 2 phương pháp trên với RS lặp lại 2 lần và CG lập lại 10 lần.
start.from: Model GAMLSS sử dụng để fitted model.
mu.start, sigma.start, nu.start, tau.start: Giá trị khởi tạo cho các tham số \(\mu, \sigma, \nu, \tau\).
control: Thiết lập tham số kiểm soát các vòng lặp bên ngoài.
Ví dụ bên dưới ta sẽ sử dụng phương pháp non-parametric smoothing line là penalised beta splines thông qua hàm pb() để hồi qui phương trình. Công thức được dùng để fitting model đối với tham số \(\sigma\) cũng là penalised beta splines được sử dụng thông qua hàm sigma.formula = ~pb(x).
library(gamlss)
library(dplyr)
library(ggplot2)
library(hrbrthemes)
data(abdom)
str(abdom)## 'data.frame': 610 obs. of 2 variables:
## $ y: int 59 64 56 61 74 60 75 63 62 67 ...
## $ x: num 12.3 12.3 12.3 12.4 12.7 ...gam <- gamlss(y~pb(x), sigma.fo = ~pb(x), family = BCT,
data = abdom, method = mixed(1,20))## GAMLSS-RS iteration 1: Global Deviance = 4771.925
## GAMLSS-CG iteration 1: Global Deviance = 4771.013
## GAMLSS-CG iteration 2: Global Deviance = 4770.994
## GAMLSS-CG iteration 3: Global Deviance = 4770.994plot(gam)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.001237747
## variance = 1.001867
## coef. of skewness = -0.004513695
## coef. of kurtosis = 2.994576
## Filliben correlation coefficient = 0.999324
## ******************************************************************Ta có thể thấy phần dư của model có dạng phân phối chuẩn khi biểu đồ Normal Q-Q Plot có phân phối thực nghiệm và phân phối thực tế số dư của mẫu nằm rất sát với đường chéo chính. Đồ thị hàm mật độ xác xuất cũng cho thấy hình dạng chuông chuẩn.
theme_set(theme_minimal())
abdom %>% ggplot(aes(x,y)) +
geom_point(size = 2, color = 'grey') +
geom_line(aes(y = fitted(gam,"mu")), linetype = 2, size = 0.5, color = 'red') +
labs(
title = 'Circumference fetuses according to gestational age',
subtitle = "GAMLSS model \nDatasource: Kings College Hospital, London",
x = "Gestational Age",
y = "Circumference Fetuses"
)
So sánh với model hồi qui tuyến tính thông thường:
lm <- lm(y~x, data = abdom)
plot(lm)
Ta có thể thấy đối với model hồi qui tuyến tính thông thường các điểm outlier sẽ không có khả năng tự hiệu chỉnh các tham số phân phối trong hàm mật độ xác xuất để fit dữ liệu. Do đó đường Normal Q-Q plot sẽ xuất hiện các điểm outlier như tại các điểm số 523, 524, 592 và sai số ngẫu nhiên của model sẽ không có dạng phân phối chuẩn.
theme_set(theme_minimal())
abdom %>% ggplot(aes(x,y)) +
geom_point(size = 2, color = 'grey') +
geom_line(aes(y = fitted(gam,"mu")), linetype = 2, size = 0.5, color = 'red') +
geom_line(aes(y = fitted(lm)), linetype = 4, size = 0.5, color = 'blue') +
labs(
title = 'Circumference fetuses according to gestational age',
subtitle = "GAMLSS model \nDatasource: Kings College Hospital, London",
x = "Gestational Age",
y = "Circumference Fetuses"
)
So sánh sai số ngẫu nhiên:
#Viết hàm tính sai số ngẫu nhiên
f_error <- function(x, pred){
n = length(x)
mae = mean(abs(x-pred))/n
mape = mean(sum(abs(x-pred))/(n*mean(x)))
mse = mean((x-pred)^2)
rmse = sqrt(mse)
data.frame(n,mae, mape, mse, rmse)
}
f_error(abdom$y, fitted(gam,'mu'))f_error(abdom$y, fitted(lm))Kết quả cho thấy mape của gamlss model chỉ là 4.44% trong khi của linear model là 4.73%.
Kiểm tra tính tự tương quan và tương quan riêng phần của các sai số ngẫu nhiên model gam.
acfResid(gam)
acfResid(lm)
gamlss đã loại bỏ được yếu tố tự tương quan của sai số ngẫu nhiên trong khi linear model không loại bỏ được.
3 Lớp các hàm fit đa thức phân số
Bên cạnh các hàm smoothing như pb(), cs(), ra(), lo() chúng ta cũng có thể sử dụng lớp các hàm fit dạng đa thức phân số như bfp(), fp(), pp(). Các hàm này sẽ hỗ trợ hàm gamlss trong quá trình fitting model. Trong đó:
bfp tạo ra một đa thức phân số dựa trên bậc của các phần tử được xác định trước để fit model chẳng hạn bfp(x, powers = c(1,2)) thì phương trình sử dụng để fit model sẽ có dạng: \[a + b_{1} x + b_{2} x^{2}\].
fp sẽ không trực tiếp đưa ra phương trình hồi qui mà sẽ lựa chọn sự kết hợp của bao nhiêu bậc trong các tham số của vector degree = c(-2, -1, -0.5, 0, 0.5, 1, 2, 3) để lựa chọn ra model phù hợp nhất. Chẳng hạn ta lựa chọn fp(x,npoly = 2) thì model sẽ dựa trên sự kết hợp của các cặp 2 bậc đa thức trong vector degree và tìm ra kết hợp có kết quả fit nhất.
pp có thể sử dụng để fit đa thức phân số có dạng:
\[a + b_{1} x^{p_{1}} + b_{2} x^ {p_{2}}\] trong đó các giá trị \(p_{1}\) và \(p_{2}\) là tùy chọn.
Cú pháp:
bfp(x, powers = c(1,2), shift = NULL, scale = NULL)
fp(x, npoly = 2, shift = NULL, scale = NULL)
pp(x, start = list(), shift = NULL, scale = NULL)Giải thích Arguments:
x: Biến giải thích được sử dụng để backfiting model thông qua hàm bfp, fp hoặc pp.
powers: Bậc của x trong phương trình đa thức.
shift: Giá trị lần dịch chuyển biến x. Mặc định là 0 nếu x dương hoặc bằng giá trị dương nhỏ nhất của x - min(x).
scale: Một số dương sử dụng để scale biến x.
npoly: Một số dương biểu thị bao nhiêu phân số đa thức được sử dụng trong fitting model. Có thể nhận các giá trị 1,2 hoặc 3 với 2 là giá trị mặc định.
start: List các giá trị khởi đầu cho quá trình tìm bậc phù hợp để tối ưu hóa phi tuyến.
Bên dưới ta sẽ áp dụng phương trình gamlss đối với hàm đa thức phân số bằng cách sử dụng các hàm smoothing gồm bfp(), fp(), pp().
#fits đa thức với bậc là 1 và 0.5
gambfp <- gamlss(y~bfp(x,c(1,0.5)), data = abdom)## GAMLSS-RS iteration 1: Global Deviance = 4948.285
## GAMLSS-RS iteration 2: Global Deviance = 4948.285summary(gambfp)## ******************************************************************
## Family: c("NO", "Normal")
##
## Call: gamlss(formula = y ~ bfp(x, c(1, 0.5)), data = abdom)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -222.19 22.27 -9.977 < 2e-16 ***
## bfp(x, c(1, 0.5))1 38.02 8.71 4.365 1.50e-05 ***
## bfp(x, c(1, 0.5))2 211.75 28.13 7.527 1.89e-13 ***
## ---
## 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) 2.63703 0.02863 92.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 610
## Degrees of Freedom for the fit: 4
## Residual Deg. of Freedom: 606
## at cycle: 2
##
## Global Deviance: 4948.285
## AIC: 4956.285
## SBC: 4973.939
## ******************************************************************#tìm ra sự kết hợp đa thức phân số với 1 phần tử trong vector degree.
gamfp1 <- gamlss(y~fp(x,1), data = abdom) ## GAMLSS-RS iteration 1: Global Deviance = 4967.042
## GAMLSS-RS iteration 2: Global Deviance = 4967.042gamfp1$mu.coefSmo## [[1]]
##
## Call:
## lm(formula = y ~ x.fp, weights = w)
##
## Coefficients:
## (Intercept) x.fp
## -544.9 334.2#tìm ra sự kết hợp của đa thức phân số 2 phần tử trong vector degree.
gamfp2 <- gamlss(y~fp(x,2), data = abdom)## GAMLSS-RS iteration 1: Global Deviance = 4941.099
## GAMLSS-RS iteration 2: Global Deviance = 4941.099gamfp2$mu.coefSmo## [[1]]
##
## Call:
## lm(formula = y ~ x.fp, weights = w)
##
## Coefficients:
## (Intercept) x.fp1 x.fp2
## -314.2734 123.1273 -0.8205#tìm ra sự kết hợp của đa thức phân số 3 phần tử trong vector degree.
gamfp3 <- gamlss(y~fp(x,3), data = abdom) ## GAMLSS-RS iteration 1: Global Deviance = 4933.937
## GAMLSS-RS iteration 2: Global Deviance = 4933.937gamfp3$mu.coefSmo## [[1]]
##
## Call:
## lm(formula = y ~ x.fp, weights = w)
##
## Coefficients:
## (Intercept) x.fp1 x.fp2 x.fp3
## -121.645 -107.252 14.507 -7.601#sử dụng phương trình đa thức định nghĩa sẵn bậc các phần từ là c(1,3)
gampp <- gamlss(y~pp(x,c(1,3)), data = abdom)## GAMLSS-RS iteration 1: Global Deviance = 4924.698
## GAMLSS-RS iteration 2: Global Deviance = 4924.698gampp$mu.coefSmo## [[1]]
## p1 p2 .lin1 .lin2 .lin3
## 7.731353e-01 7.638203e+01 -3.632312e+02 1.690460e+02 -3.695033e-47Biểu diễn đồ thị các đường hồi qui:
abdom %>% ggplot(aes(x = x, y = y)) +
geom_point(size = 2, color = 'grey') +
geom_line(aes(y=fitted(gambfp,"mu")), color = 'red', size = 0.5) +
geom_line(aes(y=fitted(gamfp3,"mu")), color = 'blue', size = 0.5) +
geom_line(aes(y=fitted(gampp,"mu")), color = 'black', size = 0.5) +
labs(
title = 'Circumference fetuses according to gestational age',
subtitle = "GAMLSS model \nDatasource: Kings College Hospital, London",
x = "Gestational Age",
y = "Circumference Fetuses"
)
Ta thấy tùy thuộc vào lựa chọn tinh chính hàm smoothing mà kết quả chúng ta thu được sẽ là khác nhau. Tuy nhiên ta cũng có thể khẳng định rằng đồ thị của đường fitted value là linh hoạt hơn so với phương pháp linear regression khi đường hồi qui có thể uốn lượn theo sự thay đổi của dữ liệu thực tế chẳng hạn như với đường màu đen của phương pháp smoothing pp(x, c(1,3)) ở độ tuổi từ > 40 tuần.
Một ví dụ khác đối tìm mối liên hệ mức độ tiêu thụ năng lượng (mpg) của các loại xe và mã lực (horsepower) trong bộ dữ liệu Auto thuộc package ISLR.
theme_set(theme_minimal())
library(ISLR)
attach(Auto)
Auto %>% ggplot(aes(x = horsepower, y = mpg)) +
geom_point(size = 2, color = 'grey')
Ta có thể thấy mối liên hệ này có dạng đa thức phân số. Do đó nếu lựa chọn mối quan hệ tuyến tính trong tình huống này là không hợp lý. Chúng ta sẽ sử dụng phương trình hồi qui có dạng đa thức phân số.
Hồi qui model:
mod1 <- gamlss(mpg~bfp(horsepower, c(1,2)), data = Auto)## GAMLSS-RS iteration 1: Global Deviance = 2266.354
## GAMLSS-RS iteration 2: Global Deviance = 2266.354mod2 <- gamlss(mpg~fp(horsepower, 1), data = Auto)## GAMLSS-RS iteration 1: Global Deviance = 2279.813
## GAMLSS-RS iteration 2: Global Deviance = 2279.813mod3 <- gamlss(mpg~fp(horsepower, 2), data = Auto)## GAMLSS-RS iteration 1: Global Deviance = 2261.314
## GAMLSS-RS iteration 2: Global Deviance = 2261.314mod4 <- gamlss(mpg~fp(horsepower, 3), data = Auto)## GAMLSS-RS iteration 1: Global Deviance = 2250.981
## GAMLSS-RS iteration 2: Global Deviance = 2250.981Visualize các đường hồi qui:
Auto %>% ggplot(aes(x = horsepower, y = mpg)) +
geom_point(size = 2, color = 'grey') +
geom_line(aes(y=fitted(mod2,"mu")), color = 'red', size = 1) +
geom_line(aes(y=fitted(mod3,"mu")), color = 'blue', size = 1) +
geom_line(aes(y=fitted(mod4,"mu")), color = 'grey', size = 1)
Các đường hồi dạng đa thức phân số đã phản ánh chuẩn xác hơn so mối quan hệ giữa mã lực và mức độ tiêu thụ năng lượng.
4 Khai phá gamlss model
Chúng ta sẽ sử dụng hàm model.frame() để khai phái các giá trị đầu vào cho phương trình hồi qui của một model nào đó. Chẳng hạn:
model.frame(mod1) %>% head()Ta sẽ thấy ở từng bước dữ liệu được sử dụng trong việc ước lượng mpg là bao nhiêu cho các biến bfp(horsepower, c(1, 2)).1 và bfp(horsepower, c(1, 2)).2.
Hoặc chúng ta muốn xem thêm những biến dummy là input của phương trình hồi qui thì sử dụng lệnh: model.matrix()
model.matrix(mod1) %>% head()## (Intercept) bfp(horsepower, c(1, 2))1 bfp(horsepower, c(1, 2))2
## 1 1 1.30 1.6900
## 2 1 1.65 2.7225
## 3 1 1.50 2.2500
## 4 1 1.50 2.2500
## 5 1 1.40 1.9600
## 6 1 1.98 3.9204Muốn tìm hiểu các thành phần trong phương trình hồi qui sử dụng lệnh terms:
terms(mod1)## mpg ~ bfp(horsepower, c(1, 2))
## attr(,"variables")
## list(mpg, bfp(horsepower, c(1, 2)))
## attr(,"factors")
## bfp(horsepower, c(1, 2))
## mpg 0
## bfp(horsepower, c(1, 2)) 1
## attr(,"term.labels")
## [1] "bfp(horsepower, c(1, 2))"
## attr(,"specials")
## attr(,"specials")$cs
## NULL
##
## attr(,"specials")$scs
## NULL
##
## attr(,"specials")$ps
## NULL
##
## attr(,"specials")$pb
## NULL
##
## attr(,"specials")$cy
## NULL
##
## attr(,"specials")$tp
## NULL
##
## attr(,"specials")$pvc
## NULL
##
## attr(,"specials")$pbm
## NULL
##
## attr(,"specials")$pbj
## NULL
##
## attr(,"specials")$pbo
## NULL
##
## attr(,"specials")$pbz
## NULL
##
## attr(,"specials")$pbc
## NULL
##
## attr(,"specials")$pbts
## NULL
##
## attr(,"specials")$pbp
## NULL
##
## attr(,"specials")$pcat
## NULL
##
## attr(,"specials")$pbq
## NULL
##
## attr(,"specials")$gmrf
## NULL
##
## attr(,"specials")$mrfa
## NULL
##
## attr(,"specials")$mrf
## NULL
##
## attr(,"specials")$sap
## NULL
##
## attr(,"specials")$sap3
## NULL
##
## attr(,"specials")$krig
## NULL
##
## attr(,"specials")$lo
## NULL
##
## attr(,"specials")$random
## NULL
##
## attr(,"specials")$re
## NULL
##
## attr(,"specials")$re4
## NULL
##
## attr(,"specials")$fp
## NULL
##
## attr(,"specials")$pp
## NULL
##
## attr(,"specials")$nl
## NULL
##
## attr(,"specials")$ri
## NULL
##
## attr(,"specials")$boost
## NULL
##
## attr(,"specials")$fk
## NULL
##
## attr(,"specials")$own
## NULL
##
## attr(,"specials")$test
## NULL
##
## attr(,"specials")$test0
## NULL
##
## attr(,"specials")$test1
## NULL
##
## attr(,"specials")$arma
## NULL
##
## attr(,"specials")$rw
## NULL
##
## attr(,"specials")$ar
## NULL
##
## attr(,"specials")$seas
## NULL
##
## attr(,"specials")$srw
## NULL
##
## attr(,"specials")$sar
## NULL
##
## attr(,"specials")$la
## NULL
##
## attr(,"specials")$tr
## NULL
##
## attr(,"specials")$ga
## NULL
##
## attr(,"specials")$ba
## NULL
##
## attr(,"specials")$mm
## NULL
##
## attr(,"specials")$nn
## NULL
##
## attr(,"specials")$sv
## NULL
##
## attr(,"specials")$ma
## NULL
##
## attr(,"specials")$pr
## NULL
##
## attr(,"specials")$pc
## NULL
##
## attr(,"specials")$h2o
## NULL
##
## attr(,"specials")$pa
## NULL
##
## attr(,"specials")$gnet
## NULL
##
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(mpg, bfp(horsepower, c(1, 2)))
## attr(,"dataClasses")
## mpg bfp(horsepower, c(1, 2))
## "numeric" "nmatrix.2"Link tài liệu tham khảo: