Sử dụng package Gamlss

Bài 5: Thành phần tuyến tính

1 Giới thiệu về thành phần tuyến tính

Đây là chặng thứ 5 trên hành trình khám phá packae gamlss. Để nhắc lại, mô hình gamlss có 2 đặc tính quan trọng : thứ nhất là chuyên biệt cho tham số của họ phân phối xác định, thứ hai là có cấu trúc nhiều thành phần, bao gồm 3 phần chính là : tham số, bù trừ và ngẫu nhiên. Phần tham số lại chia ra thành 2 loại bao gồm tuyến tính và phi tuyến tính.

Bài hôm nay Nhi sẽ giới thiệu sơ lược về thành phần tham số tuyến tính (các yếu tố tuyến tính này có thể gia giảm nên còn được gọi là phần cộng thêm : linear parametric additive term). Phần phi tuyến sẽ được giới thiệu trong 1 bài khác.

Bạn sẽ thấy trong hình sau đây 3 món ăn mà Nhi thích: Maki sushi, Kem hoa hồng của tiệm Armorino và Phở Việt Nam. Có thể các bạn không tin, nhưng chúng ta có thể tìm thấy những ý tưởng về thống kê và mô hình hồi quy từ những món ăn và trong lúc nấu ăn. Thí dụ có lần trong khi làm món sushi, Nhi liên tưởng tới data transformation và piecewise regression (có một số giáo trình thống kê có tựa đề là Cookbook).

Cả 3 món ăn trong hình đều có một đặc điểm chung, đó là công thức chế biến (công thức mô hình) đòi hỏi bạn kết hợp các thành phần khác nhau (các biến số trong dữ liệu đầu vào) bằng phép cộng (như món sushi) hoặc kết hợp giữa nguyên liệu và gia vị (tương tác) để tạo ra kết quả cuối cùng (response).

Đôi khi cùng 1 dữ liệu đầu vào nhưng có khả năng tạo ra rất nhiều mô hình khác nhau bằng cách gia giảm, thêm bớt các biến số bên trong mô hình. Tại cửa hàng Armorino, bạn đứng trước một tủ kính với hàng chục loại kem, và có thể chọn 2,3, thậm chí 5 mùi khác nhau. Chị bán kem sẽ dán từng lớp kem với nhau trên cornet để sau đó bạn có thể vui vẻ dạo phố với cành hoa hồng trên tay và chụp selfie với nó. Bằng cách thức tương tự, mọi mô hình hồi quy đã được tạo ra trong R, dù bạn sử dụng package gamlss hay glm.

Công thức thành phần tuyến tính và toán tử

Các thành phần tuyến tính được đưa vào mô hình thông qua công thức (formula) và liên kết với nhau bằng 2 toán tử : cộng thêm (+) và tương tác (*).

Hai toán tử này có nguồn gốc lâu đời (đề xuất bởi Wilkinson và Rogers năm 1973 và đưa vào ngôn ngữ lập trình thống kê S (là tiền thân của R) bởi Chambers và Hastie năm 1992). Do được quy ước sẵn trong R nên hầu hết (không phải tất cả) các package hồi quy sử dụng formula đều dùng + và * với ý nghĩa như nhau.

Các bạn từng làm ANOVA hay hồi quy với hàm glm() hoặc lm() đều quen thuộc với công thức và 2 toán tử +, *, nhưng ở đây Nhi sẽ nhắc lại về cú pháp của công thức hồi quy cho các bạn mới làm quen với R :

Toán tử +

Toán tử + có công dụng kết hợp 2 tham số hồi quy (predictors) với nhau bằng phép cộng : Thí dụ : Y ~ X1 + X2 + … Xn

Nếu X là 1 biến liên tục, toán tử + sẽ gán vector X vào matrix của mô hình hồi quy (cộng thêm 1 cột duy nhất). Với 2 predictors X1+X2, một mặt phẳng hồi quy sẽ được tạo ra.

Nếu ta đưa vào 1 yếu tố (factor) F với nhiều level, F sẽ bị cắt ra thành nhiều biến giả (dummy variable) được mã hóa bằng giá trị 0 và 1, mỗi cột dummy variable tương ứng với 1 level.

Khi mô hình có chứa intercept, và ta muốn khảo sát một yếu tố F có k levels thì chỉ có (k-1) level được mã hóa theo dạng dummy variable, level đầu tiên của F sẽ bị giản lược và không được trình bày trong mô hình sau cùng.

Ta có thể tổ chức lại nội dung của F bằng hàm contrasts( ).

library(tidyverse)

F1=gl(2,1,25)
levels(F1)=c("Group1", "Group2")

F2=gl(3,2,25)
levels(F2)=c("CondA", "CondB","CondC")

X1=rnorm(25,10,2)

X2=rnorm(25,5,1.5)

head(model.matrix(~X1),10)%>%knitr::kable()

(Intercept)	X1
1	9.594486
1	8.797346
1	7.458943
1	12.099089
1	7.917185
1	9.890659
1	7.157442
1	9.588284
1	10.800068
1	9.328071

head(model.matrix(~X1+X2),10)%>%knitr::kable()

(Intercept)	X1	X2
1	9.594486	4.366550
1	8.797346	4.403565
1	7.458943	6.300203
1	12.099089	6.416712
1	7.917185	6.041099
1	9.890659	4.058703
1	7.157442	7.891579
1	9.588284	5.119930
1	10.800068	7.052355
1	9.328071	2.636837

head(model.matrix(~F2),10)%>%knitr::kable()

(Intercept)	F2CondB	F2CondC
1	0	0
1	0	0
1	1	0
1	1	0
1	0	1
1	0	1
1	0	0
1	0	0
1	1	0
1	1	0

head(model.matrix(~F2-1),10)%>%knitr::kable()

F2CondA	F2CondB	F2CondC
1	0	0
1	0	0
0	1	0
0	1	0
0	0	1
0	0	1
1	0	0
1	0	0
0	1	0
0	1	0

Toán tử tương tác

Toán tử * có công dụng biểu diễn một sự tương tác giữa 2 hay nhiều predictors, chúng có thể là biến liên tục hay yếu tố , thí dụ :

Y ~ X1*X2

Y ~ X*F

Y ~ F1*F2

Lưu ý rằng khi viết X1*X2, thực ra ta muốn xét cả 3 hiệu ứng : 2 hiệu ứng riêng phần của X1 và X2, và hiệu ứng tương tác giữa chúng, do đó :

Y ~ X1*X2 tương đương : Y ~ X1 +X2 + X1:X2

Tương tác giữa 2 factors F1 (k1 level) và F2 (k2 level) sẽ tạo ra trong matrix (k1-1)*(k2-1)

và (k1-1) + (k2-1) hay (k1*k2)-1 cột dummy variable cho hiệu ứng chính và hiệu ứng tương tác.

Thí dụ F1 có 2 level (a,b), F2 có 3 level là (1,2,3) thì sẽ có 5 dummy variable gồm F1b, F22, F23, F1B :F22, F1B :F23

head(model.matrix(~X1*F2),10)%>%knitr::kable()

(Intercept)	X1	F2CondB	F2CondC	X1:F2CondB	X1:F2CondC
1	9.594486	0	0	0.000000	0.000000
1	8.797346	0	0	0.000000	0.000000
1	7.458943	1	0	7.458943	0.000000
1	12.099089	1	0	12.099089	0.000000
1	7.917185	0	1	0.000000	7.917185
1	9.890659	0	1	0.000000	9.890659
1	7.157442	0	0	0.000000	0.000000
1	9.588284	0	0	0.000000	0.000000
1	10.800068	1	0	10.800068	0.000000
1	9.328071	1	0	9.328071	0.000000

head(model.matrix(~F1+F2),10)%>%knitr::kable()

(Intercept)	F1Group2	F2CondB	F2CondC
1	0	0	0
1	1	0	0
1	0	1	0
1	1	1	0
1	0	0	1
1	1	0	1
1	0	0	0
1	1	0	0
1	0	1	0
1	1	1	0

head(model.matrix(~F1*F2),10)%>%knitr::kable()

(Intercept)	F1Group2	F2CondB	F2CondC	F1Group2:F2CondB	F1Group2:F2CondC
1	0	0	0	0	0
1	1	0	0	0	0
1	0	1	0	0	0
1	1	1	0	1	0
1	0	0	1	0	0
1	1	0	1	0	1
1	0	0	0	0	0
1	1	0	0	0	0
1	0	1	0	0	0
1	1	1	0	1	0

head(model.matrix(~X1*X2*F1*F2),10)%>%knitr::kable()

(Intercept)	X1	X2	F1Group2	F2CondB	F2CondC	X1:X2	X1:F1Group2	X2:F1Group2	X1:F2CondB	X1:F2CondC	X2:F2CondB	X2:F2CondC	F1Group2:F2CondB	F1Group2:F2CondC	X1:X2:F1Group2	X1:X2:F2CondB	X1:X2:F2CondC	X1:F1Group2:F2CondB	X1:F1Group2:F2CondC	X2:F1Group2:F2CondB	X2:F1Group2:F2CondC	X1:X2:F1Group2:F2CondB	X1:X2:F1Group2:F2CondC
1	9.594486	4.366550	0	0	0	41.89480	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0	0	0.00000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	8.797346	4.403565	1	0	0	38.73968	8.797346	4.403565	0.000000	0.000000	0.000000	0.000000	0	0	38.73968	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	7.458943	6.300203	0	1	0	46.99286	0.000000	0.000000	7.458943	0.000000	6.300203	0.000000	0	0	0.00000	46.99286	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	12.099089	6.416712	1	1	0	77.63637	12.099089	6.416712	12.099089	0.000000	6.416712	0.000000	1	0	77.63637	77.63637	0.00000	12.099089	0.000000	6.416712	0.000000	77.63637	0.00000
1	7.917185	6.041099	0	0	1	47.82849	0.000000	0.000000	0.000000	7.917185	0.000000	6.041099	0	0	0.00000	0.00000	47.82849	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	9.890659	4.058703	1	0	1	40.14325	9.890659	4.058703	0.000000	9.890659	0.000000	4.058703	0	1	40.14325	0.00000	40.14325	0.000000	9.890659	0.000000	4.058703	0.00000	40.14325
1	7.157442	7.891579	0	0	0	56.48353	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0	0	0.00000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	9.588284	5.119930	1	0	0	49.09134	9.588284	5.119930	0.000000	0.000000	0.000000	0.000000	0	0	49.09134	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	10.800068	7.052355	0	1	0	76.16591	0.000000	0.000000	10.800068	0.000000	7.052355	0.000000	0	0	0.00000	76.16591	0.00000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000
1	9.328071	2.636837	1	1	0	24.59660	9.328071	2.636837	9.328071	0.000000	2.636837	0.000000	1	0	24.59660	24.59660	0.00000	9.328071	0.000000	2.636837	0.000000	24.59660	0.00000

Khi ta muốn xét tương tác giữa 2 biến liên tục X và 1 yếu tố F (k level), trong matrix sẽ gồm : 1 cột cho vector X, (k-1) cột cho các dummy hiệu ứng chính của (k-1) level trong F và (k-1) cột cho hiệu ứng tương tác giữa X và các dummy variable của (k-1) level trong F.

Liên hệ giữa mô hình tuyến tính và phân tích định lượng

Như chúng ta đã biết, có liên hệ chặt chẽ giữa mô hình tuyến tính và các phương pháp thống kê cổ điển, và những thiết kế nghiên cứu có thể được biểu diễn thông qua mô hình, thí dụ :

Y ~ 1 là một mô hình chỉ chứa intercept, nó tương đương với việc khảo sát giá trị trung bình của Y, hoặc với hàm ý chuẩn bị 1 mô hình cơ bản (baseline model) dùng vào việc đánh giá phẩm chất của những mô hình khác ; thí dụ trong hồi quy logistic, một model chỉ bắt đầu đáng tin nếu nó có thể đánh bại quy luật binnomial ngẫu nhiên.

Y ~ X với X là 1 biến liên tục : tương đương với hồi quy đơn biến hay khảo sát tương quan tuyến tính giữa Y và X

Y ~ F với F là 1 factor từ 2 level trở lên : tương đương với One way ANOVA (phân tích phương sai 1 yếu tố).

Y ~ X1 + X2 : tương đương với khảo sát quan hệ riêng phần giữa Y,X1 và X2

Y ~ X*F : tương đương với ANCOVA

Y ~ F1*F2…Fn : tương đương với n ways ANOVA

Thí dụ minh họa

Nhi lấy thí dụ dataset cigarette, đây là 1 bài toán hồi quy với câu hỏi là liệu việc đánh thuế lên thuốc lá có làm giảm doanh số của sản phẩm này hay không ? Tác giả đã ghi nhận giá thuốc lá trong vòng 30 năm từ 1963 tại 51 bang của Hoa kỳ, đồng thời ghi nhận các biến như dân số, dân số >16 tuổi, thu nhập bình quân đầu người, giá thuốc lá. Biến kết quả là mức tiêu thụ thuốc lá.

Trong khi thăm dò quan hệ giữa biến kết quả (mức tiêu thụ thuốc lá tính bằng gói trên đầu người) và những biến số khác trong dữ liệu, trên một không gian 2 chiều, chúng ta thực ra có thể cảm nhận được về những mô hình hồi quy đơn biến, ANOVA, mô hình có hiệu ứng tương tác…

cigarette=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/plm/Cigar.csv")

cigarette$state=factor(cigarette$state)
cigarette$year=factor(cigarette$year-62)

library(viridis)

cigarette%>%ggplot()+
  geom_boxplot(aes(x=as.integer(year),y=sales,fill=year,color=year),show.legend = F,alpha=0.2)+
  theme_bw(10)+
  scale_fill_viridis(discrete=T,option="C",direction = -1)+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  geom_smooth(aes(x=as.integer(year),y=sales),show.legend = F,color="darkviolet",fill="gold")+
  scale_x_continuous("30 years period")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  ggtitle("Model: Sale~ Year (as dummy variables)")

cigarette%>%ggplot()+
  geom_path(aes(x=as.integer(year),y=sales,color=reorder(state,sales)),show.legend = F)+
  theme_bw(10)+
  scale_colour_viridis(discrete=T,option="C",direction = -1)+
  scale_x_continuous("30 years period")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  ggtitle("Model: Sale~ Year*State (both as factors)")

cigarette%>%ggplot()+
  geom_jitter(aes(x=as.integer(year),y=sales,color=sales),size=1,alpha=0.5)+
  geom_smooth(aes(x=as.integer(year),y=sales),show.legend = F,color="purple",fill="gold")+
  theme_bw(10)+
  scale_x_continuous("30 years period")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=F,option="B",direction = -1)+
  ggtitle("Model: Sale~ Year (as discrete number)")

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales,color=sales),size=1,alpha=0.5)+
  geom_smooth(aes(x=price,y=sales),show.legend = F,color="purple",fill="gold")+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=F,option="B",direction = -1)+
  ggtitle("Model: Sale~ Price per pack")

cigarette%>%ggplot()+
  geom_point(aes(x=pimin,y=sales,color=sales),size=1,alpha=0.2)+
  geom_smooth(aes(x=pimin,y=sales),color="purple",fill="gold")+
  theme_bw(10)+
  scale_x_continuous("minimum price in adjoining states per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=F,option="B",direction = -1)+
  ggtitle("Model: Sale~ minimum price")

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales),size=1,alpha=0.2,show.legend = F,color="grey")+
  geom_smooth(aes(x=price,y=sales,color=year),show.legend = F,se=F,method="lm")+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  ggtitle("Model: Sale~ Price:Year")

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales),size=1,alpha=0.2,show.legend = F,color="grey")+
  geom_smooth(aes(x=price,y=sales,color=state),show.legend = F,se=F)+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  ggtitle("Model: Sale~ Price*State")

cigarette%>%ggplot()+
  geom_jitter(aes(x=log(pop16),y=sales,color=sales),size=1,alpha=0.5)+
  geom_smooth(aes(x=log(pop16),y=sales),show.legend = F,color="purple",fill="gold")+
  theme_bw(10)+
  scale_x_continuous("Population above 16 yrs old in Log scale")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=F,option="B",direction = -1)+
  ggtitle("Model: Sale~ log(Population>16 yrs)")

cigarette%>%ggplot()+
  geom_jitter(aes(x=ndi,y=sales,color=sales),size=1,alpha=0.5)+
  geom_smooth(aes(x=ndi,y=sales),show.legend = F,color="purple",fill="gold")+
  theme_bw(10)+
  scale_x_continuous("Per capita disposable income")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=F,option="B",direction = -1)+
  ggtitle("Model: Sale~ Per capita disposable income")

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales,color= factor(round(log(cigarette$ndi),0))),size=1,alpha=0.2,show.legend = F)+
  geom_smooth(aes(x=price,y=sales,color= factor(round(log(cigarette$ndi),0))),show.legend = F,se=F,method="lm")+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  ggtitle("Model: Sale~ Price * Per capita disposable income")

cigarette%>%ggplot()+
  geom_boxplot(aes(x=reorder(state,sales),y=sales,fill=reorder(state,sales),color=reorder(state,sales)),show.legend = F,alpha=0.5)+
  theme_bw(10)+
  scale_fill_viridis(discrete=T,option="B",direction = -1)+
  scale_color_viridis(discrete=T,option="B",direction = -1)+
  scale_x_discrete("51 States")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  ggtitle("One_way ANOVA: Sale~ State (as dummy variables)")

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales,color=year),size=1,alpha=0.2,show.legend = F)+
  geom_smooth(aes(x=price,y=sales,color=year),show.legend = F,se=F)+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  facet_wrap(~year,scales="free")+
  ggtitle("Model: Sale~ Price + random(1|Year)")+
  theme_bw(3)

cigarette%>%ggplot()+
  geom_jitter(aes(x=price,y=sales,color=state),size=1,alpha=0.2,show.legend = F)+
  geom_smooth(aes(x=price,y=sales,color=state),show.legend = F,se=F)+
  theme_bw(10)+
  scale_x_continuous("Price per pack")+
  scale_y_continuous("Cigarette sales (packs per capita)")+
  scale_color_viridis(discrete=T,option="C",direction = -1)+
  facet_wrap(~state,scales="free")+
  ggtitle("Model: Sale~ Price + Year + random(1|State)")+
  theme_bw(3)

p2=cigarette%>%ggplot(aes(x=sales))+
  geom_histogram(aes(y=..density..,fill=..density..,color=..density..),binwidth = 1, 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")

2 Thăm dò phân phối

Như trong bài 2, bước đầu tiên khi dựng mô hình trong gamlss là thăm dò về đặc tính phân phối của biến kết quả. Bằng hàm fitdist, Nhi nhanh chóng xác định được phân phối phù hợp nhất cho biến sales là họ phân phối Skewed T type 4. Họ phân phối này được mô tả bởi 4 tham số là Mu, Sigma, Nu và Tau (có thể hình dung là Mean, Sd, skewness và kurtosis). Sử dụng hàm rST4, Nhi mô phỏng một mẫu ngẫu nhiên với Mu, Sigma, Tau và Nu của mô hình chỉ chứa Intercept), kết quả mô phỏng khá phù hợp với quans át trên thực tế :

library(gamlss)

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

scan$fits

##      ST4    BCPEo     BCTo      ST1      ST5      JSU      ST2      ST3 
## 12800.57 12801.47 12812.24 12813.80 12819.56 12829.01 12829.37 12840.16 
##     SEP4   SHASHo       TF       GT     SEP3     SEP1       PE     SEP2 
## 12848.93 12849.70 12851.47 12857.84 12881.36 12888.94 12901.18 12909.46 
##   exGAUS   IGAMMA       RG    BCCGo       GG      GIG    LOGNO       LO 
## 12915.38 13009.24 13011.65 13012.41 13016.37 13016.47 13037.61 13043.48 
##       IG       GA      SN2      SN1     WEI3       GU      EXP  PARETO2 
## 13051.96 13118.36 13175.23 13413.97 13634.60 14491.46 16070.11 16077.34

gamlss.family(ST4)

## 
## GAMLSS Family: ST4 skew t type 4 
## Link function for mu   : identity 
## Link function for sigma: log 
## Link function for nu   : log 
## Link function for tau  : log

set.seed(123456)
p1=rST4(1000, mu = 119, sigma =exp(2.833), nu=exp(1.892), tau=exp(0.6194))%>%
  as_tibble()%>%
  ggplot(aes(x=value))+
  geom_histogram(aes(y=..density..,fill=..density..,color=..density..),binwidth = 1, 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 ST4 (mu=119,s=exp(2.83),nu=exp(1.89),tau=exp(0.62))")

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

3 Mô hình gamlss với các thành phần tuyến tính

Nhi sẽ lần lượt dựng 5 mô hình, các bạn có thể xem xét formula và nội dung của mô hình xuất ra để có thể hình dung về hoạt động của toán tử + và *, đồng thời diễn dịch các thành phần tuyến tính, tương tác trong mô hình.

Mô hình chỉ chứa hằng số

nC<-detectCores()

m0=gamlss(data=cigarette,
          formula=sales~1,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
          )

summary(m0)

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  gamlss(formula = sales ~ 1, family = ST4, data = cigarette,  
##     trace = F, parallel = "multicore", ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 119.0192     0.5679   209.6   <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)  2.83347    0.03661    77.4   <2e-16 ***
## ---
## 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)   1.8924     0.2484   7.619 4.72e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61952    0.07959   7.784 1.37e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  4
##       Residual Deg. of Freedom:  1376 
##                       at cycle:  10 
##  
## Global Deviance:     12771.65 
##             AIC:     12779.65 
##             SBC:     12800.57 
## ******************************************************************

Mô hình đơn biến

m1=gamlss(data=cigarette,
          formula=sales~price,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
)

summary(m1)

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  gamlss(formula = sales ~ price, family = ST4, data = cigarette,  
##     trace = F, parallel = "multicore", ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 130.61543    1.01783  128.33   <2e-16 ***
## price        -0.16791    0.01213  -13.85   <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)  2.74743    0.03589   76.56   <2e-16 ***
## ---
## 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)   1.7963     0.2181   8.238 4.04e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.57867    0.07716   7.499 1.15e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  5
##       Residual Deg. of Freedom:  1375 
##                       at cycle:  9 
##  
## Global Deviance:     12595.12 
##             AIC:     12605.12 
##             SBC:     12631.27 
## ******************************************************************

Mô hình 2 biến tương tác

m2=gamlss(data=cigarette,
          formula=sales~price*ndi,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
)

summary(m2)

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  gamlss(formula = sales ~ price * ndi, family = ST4,  
##     data = cigarette, trace = F, parallel = "multicore",  
##     ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.172e+02  1.792e+00  65.424  < 2e-16 ***
## price       -3.751e-01  4.989e-02  -7.518 9.97e-14 ***
## ndi          5.528e-03  3.274e-04  16.884  < 2e-16 ***
## price:ndi   -2.051e-05  2.444e-06  -8.395  < 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)  2.61914    0.03565   73.47   <2e-16 ***
## ---
## 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)   1.8638     0.2096   8.891   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.50719    0.07556   6.712 2.79e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  1373 
##                       at cycle:  9 
##  
## Global Deviance:     12309.88 
##             AIC:     12323.88 
##             SBC:     12360.48 
## ******************************************************************

ANOVA với Year là fixed factor

m3=gamlss(data=cigarette,
          formula=sales~year,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
)

summary(m3)

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  gamlss(formula = sales ~ year, family = ST4, data = cigarette,  
##     trace = F, parallel = "multicore", ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 119.6801     2.9393  40.717  < 2e-16 ***
## year2        -3.7095     4.1785  -0.888 0.374833    
## year3        -1.0348     4.1252  -0.251 0.801979    
## year4        -3.2866     4.2544  -0.773 0.439944    
## year5        -1.1830     4.1496  -0.285 0.775627    
## year6        -2.2424     4.0820  -0.549 0.582861    
## year7        -3.8627     3.9638  -0.975 0.329980    
## year8        -5.5466     3.8691  -1.434 0.151927    
## year9        -3.4533     3.8295  -0.902 0.367348    
## year10       -0.2822     3.9038  -0.072 0.942380    
## year11        0.4806     3.8979   0.123 0.901881    
## year12        4.1379     3.9596   1.045 0.296196    
## year13        5.1056     3.9123   1.305 0.192113    
## year14        7.7548     3.8941   1.991 0.046640 *  
## year15        8.9079     3.8886   2.291 0.022131 *  
## year16       10.1401     3.8375   2.642 0.008327 ** 
## year17        8.4816     3.8384   2.210 0.027295 *  
## year18       10.1201     3.8488   2.629 0.008651 ** 
## year19       10.1538     3.8050   2.669 0.007710 ** 
## year20        8.8099     3.8237   2.304 0.021375 *  
## year21        4.9056     3.7674   1.302 0.193091    
## year22        0.5878     3.7874   0.155 0.876692    
## year23       -0.8223     3.8064  -0.216 0.828998    
## year24       -3.1214     3.8586  -0.809 0.418676    
## year25       -6.9478     3.9006  -1.781 0.075105 .  
## year26      -10.3361     4.0035  -2.582 0.009934 ** 
## year27      -14.6591     3.9756  -3.687 0.000236 ***
## year28      -20.0512     3.9548  -5.070 4.53e-07 ***
## year29      -21.3457     4.0198  -5.310 1.28e-07 ***
## year30      -23.3993     4.0704  -5.749 1.11e-08 ***
## ---
## 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.65762    0.03569   74.46   <2e-16 ***
## ---
## 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)   1.8054     0.1968   9.173   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.48745    0.07488    6.51 1.06e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  33
##       Residual Deg. of Freedom:  1347 
##                       at cycle:  9 
##  
## Global Deviance:     12453.97 
##             AIC:     12519.97 
##             SBC:     12692.55 
## ******************************************************************

Mô hình tương tác 3 biến

m4=gamlss(data=cigarette,
          formula=sales~price*ndi*pop16,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
)

summary(m4)

## Warning in summary.gamlss(m4): summary: vcov has failed, option qr is used instead

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  
## gamlss(formula = sales ~ price * ndi * pop16, family = ST4, data = cigarette,  
##     trace = F, parallel = "multicore", ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.129e+02  2.203e+00  51.257  < 2e-16 ***
## price           -5.889e-01  5.821e-02 -10.117  < 2e-16 ***
## ndi              7.827e-03  3.939e-04  19.868  < 2e-16 ***
## pop16            1.326e-03  5.991e-04   2.214  0.02702 *  
## price:ndi       -1.957e-05  3.207e-06  -6.102 1.36e-09 ***
## price:pop16      5.100e-05  1.845e-05   2.764  0.00578 ** 
## ndi:pop16       -5.575e-07  1.004e-07  -5.552 3.39e-08 ***
## price:ndi:pop16 -1.606e-10  7.742e-10  -0.207  0.83569    
## ---
## 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.58474    0.02693   95.97   <2e-16 ***
## ---
## 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)   1.9334     0.1861   10.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4833     0.0734   6.585 6.44e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  11
##       Residual Deg. of Freedom:  1369 
##                       at cycle:  9 
##  
## Global Deviance:     12225.52 
##             AIC:     12247.52 
##             SBC:     12305.05 
## ******************************************************************

Mô hình có cả biến liên tục và yếu tố = ANCOVA

m5=gamlss(data=cigarette,
          formula=sales~price+year,
          family=ST4,
          trace=F,
          parallel="multicore",
          ncpus = detectCores()
)

summary(m5)

## ******************************************************************
## Family:  c("ST4", "skew t type 4") 
## 
## Call:  gamlss(formula = sales ~ price + year, family = ST4,  
##     data = cigarette, trace = F, parallel = "multicore",  
##     ncpus = detectCores()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 143.73434    3.30681  43.466  < 2e-16 ***
## price        -0.88847    0.06856 -12.959  < 2e-16 ***
## year2        -3.04465    3.90047  -0.781 0.435183    
## year3         0.20743    3.81303   0.054 0.956624    
## year4        -0.97454    4.03306  -0.242 0.809098    
## year5         2.16177    3.95856   0.546 0.585087    
## year6         2.75561    3.85874   0.714 0.475275    
## year7         2.37619    3.80034   0.625 0.531908    
## year8         3.03405    3.74560   0.810 0.418065    
## year9         7.20534    3.67231   1.962 0.049960 *  
## year10       11.90394    3.66587   3.247 0.001194 ** 
## year11       13.34395    3.70133   3.605 0.000323 ***
## year12       18.14435    3.76781   4.816 1.63e-06 ***
## year13       21.79917    3.82046   5.706 1.42e-08 ***
## year14       27.05878    3.85073   7.027 3.34e-12 ***
## year15       29.40588    3.87276   7.593 5.82e-14 ***
## year16       35.41725    3.99930   8.856  < 2e-16 ***
## year17       35.78210    4.06322   8.806  < 2e-16 ***
## year18       40.03731    4.17714   9.585  < 2e-16 ***
## year19       43.25537    4.30948  10.037  < 2e-16 ***
## year20       47.58785    4.57468  10.402  < 2e-16 ***
## year21       54.18075    5.11089  10.601  < 2e-16 ***
## year22       59.10690    5.61826  10.520  < 2e-16 ***
## year23       62.65829    5.92569  10.574  < 2e-16 ***
## year24       67.16567    6.38709  10.516  < 2e-16 ***
## year25       69.60638    6.82172  10.204  < 2e-16 ***
## year26       73.33025    7.29444  10.053  < 2e-16 ***
## year27       78.53767    7.94158   9.889  < 2e-16 ***
## year28       85.35976    8.73088   9.777  < 2e-16 ***
## year29       91.03478    9.24477   9.847  < 2e-16 ***
## year30      109.03117   10.85850  10.041  < 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)  2.53185    0.03788   66.85   <2e-16 ***
## ---
## 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)   1.4000     0.1501    9.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.43097    0.07401   5.823 7.21e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1380 
## Degrees of Freedom for the fit:  34
##       Residual Deg. of Freedom:  1346 
##                       at cycle:  9 
##  
## Global Deviance:     12272.55 
##             AIC:     12340.55 
##             SBC:     12518.36 
## ******************************************************************

4 Kết luận

Đây là 1 bài rất ngắn và nội dung của nó không có gì mới lạ đối với các bạn đã quen sử dụng R. Chúng ta sẽ dừng lại ở đây. Thực ra công dụng chính của Gamlss là cho mô hình tiên lượng chứ rất hiếm khi nó được sử dụng với mục đích diễn dịch, và output gamlss không tương thích với tất cả hàm phân tích ANOVA mà ta từng biết nên để làm được suy diễn thống kê, hoặc là ta phải tự viết hàm để tính F statistics và làm ANOVA, post-hoc test một cách thủ công, hoặc là ta phải suy diễn thống kê dựa vào marginal effect với bootstrap (cũng thủ công luôn). Thông điệp chính của bài đó là cú pháp cho phần tuyến tính trong gamlss hoàn toàn giống như của hàm glm() hoặc lm() và nhiều package hồi quy khác. Bài này chỉ là khởi động cho những thành phần phức tạp hơn mà Nhi sẽ đề cập trong lần tới, bao gồm đa thức bậc cao và smoothing splines.

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