Mixed negative binomial

.Giới thiệu
.Mô tả dữ liệu
.Hồi quy poisson và negative binomial
.Mixed negative binomial với trường phái Bayes

.Giới thiệu

Trong thực tế, chúng ta hay gặp phải dữ liệu ở dạng số đếm (count data) như đếm số hoa ,quả, hạt trên cây trong 1 mùa vụ; số sản phẩm bán ra trong 1 ngày; số ngày nằm viện., etc. Đặc điểm cần lưu ý với count data là phân bố dữ liệu thường hay bị lệch phải, phương sai tăng khi giá trị trung bình tăng. Như vậy các phương pháp xử lý số liệu có giả định phân phối chuẩn (t-test, ANOVA, OLS regression, etc) khi áp dụng với dữ liệu count data sẽ không còn phù hợp nữa (trừ khi dữ liệu đã được chuyển đổi). Khi gặp dữ liệu dạng này chúng ta có thể có 2 lựa chọn: 1) chuyển đổi dự liệu về phân phối chuẩn; hoặc 2) thay đổi phương pháp xử lý (phi tham số, hồi quy tuyến tính hỗn hợp tổng quát- Generalized linear [mixed] model).

Ở bài hướng dẫn này chúng ta sẽ làm quen với một phương pháp khá phổ biến trong họ gia đình GLMM là Mixed negative binomial regression (tạm dịch; hồi quy nhị thức âm hỗn hợp) dùng để thay thế hồi quy OLS trong trường hợp biến phụ thuộc có dạng số đếm. Chúng tôi sẽ không đi sâu vào các công thức toán mà sẽ chú trọng giúp cho các bạn biết khi nào nên sử dụng mô hình negative binomial (NB), thực hành trên R, Stan theo trường phái Bayes và cách diễn giải kết quả.

Khi các bạn muốn xem xét ảnh hưởng của các biến độc lập \(\ x_{1}, x_{2},..x_{n}\) lên biến phụ thuộc dạng số đếm \(y\) thì chúng ta thường hay nghĩ ngay đến mô hình hồi quy Poisson. Mô hình Poisson có một giả định (assumption) rất quan trọng là: giá trị mean và variance bằng nhau. Tuy nhiên dữ liệu thực tế mà chúng ta thu thập ít khi thõa mãn giả định này. Khi variance > mean chúng ta gọi đó là over-dispersion (tạm dịch: quá phân tán), ngược lại, variance < mean thì gọi là under-dispersion .Khi dữ liệu gặp vấn đề over-dispersion thì đó chính là lúc chúng ta dùng mô hình negative binomial. Vậy làm cách nào để kiểm tra overdispersion ? Để trả lời câu hỏi đó, chúng ta hãy cùng thực hành với bộ dữ liệu “epil”.

Bộ dữ liệu này được thu thập từ 59 bệnh nhân bị động kinh. Các bệnh nhân được theo dõi số lần bị động kinh trong vòng 8 tuần (base). Sau đó bệnh nhân được phân chia ngẫu nhiên vào nhóm điều trị (treatment- dùng thuốc progabide ) và nhóm giả dược (placebo); Số lần bị động kinh được ghi nhận trong 2 tuần ở 8 tuần tiếp theo. Như vậy mỗi bệnh nhân sẽ được ghi nhận ở 4 thời điểm (period) và dạng dữ liệu này được gọi là dữ liệu tái đo lường (repeated measurement data). Ở bài thực hành này chúng ta sẽ sử dụng một mô hình hồi quy NB với random intercept để xem xét ảnh hưởng của biến lbase, điều trị (trt) và biến visit (được chuyển đổi từ biến period) với số lần bị động kinh trong thời gian theo dõi.

dat=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/MASS/epil.csv")
names(dat)[2] <- "sei"
dat$visit <- (2*dat$period-5)/10
head(dat, 5) %>% knitr:: kable()

X	sei	trt	base	age	V4	subject	period	lbase	lage	visit
1	5	placebo	11	31	0	1	1	-0.7563538	0.1142037	-0.3
2	3	placebo	11	31	0	1	2	-0.7563538	0.1142037	-0.1
3	3	placebo	11	31	0	1	3	-0.7563538	0.1142037	0.1
4	3	placebo	11	31	1	1	4	-0.7563538	0.1142037	0.3
5	3	placebo	11	30	0	2	1	-0.7563538	0.0814139	-0.3

.Mô tả dữ liệu

Dữ liệu có 236 quan sát, trong đó có 23 quan sát có giá trị zero. Mục đích của bài thực hành này là mô hình NB vì vậy chúng ta sẽ tạm thời bỏ qua vấn đề về zero.
Xem xét sự dao động của biến \(y\) ở các lần ghi nhận; và giữa progabide vs. giả dược.

box <-  dat %>% ggplot(., aes(x = trt, y = sei)) +
        geom_boxplot(aes(fill = as.factor(period)), alpha = 0.7,outlier.colour = "blue") +
        theme_bw() + labs(x = "treatment", y = "Count") +
        scale_fill_discrete(name = "Period")+
        theme(axis.text = element_text(size = 12, color = "black"),
              axis.title = element_text(size = 15, color = "black"))
        
box

Từ biểu đồ boxplot chúng ta thấy rằng một số bệnh nhân có số lần bị động kinh rất cao và hiện tượng này xuất hiện ở cả 2 nhóm dùng thuốc progabide và placebo.

Phân bố của lbase ở cả 2 nhóm không tuân theo phân bố chuẩn và không đồng dạng. Bệnh nhân được phân vào 2 nhóm (trt) một cách ngẫu nhiên nên không có sự khác biệt về lbase.

Kiểm tra phân bố của lbase.

tes.base <- dat %>% split(.$trt) %>% map(~ fBasics::shapiroTest(.$lbase))
tes.base

## $placebo
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9452
##   P VALUE:
##     0.0001702 
## 
## Description:
##  Mon Nov 06 00:45:25 2017 by user: Tien Tai
## 
## 
## $progabide
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9705
##   P VALUE:
##     0.00813 
## 
## Description:
##  Mon Nov 06 00:45:25 2017 by user: Tien Tai

wilcox.test(dat$lbase ~ dat$trt)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  dat$lbase by dat$trt
## W = 6432, p-value = 0.3284
## alternative hypothesis: true location shift is not equal to 0

dens.base <- dat %>% ggplot(., aes(x =lbase, color = trt)) + geom_line(stat = "density", size = 0.9) +
        labs(x = "lbase", y = "Density") + theme_bw() + scale_color_discrete(name = "Treatment") + 
        theme(axis.text = element_text(size = 12, color = "black"),
        axis.title = element_text(size = 15, color = "black")) +
        theme(panel.grid.minor = element_blank())
dens.base

.Hồi quy poisson và negative binomial

Hồi quy poisson là mô hình được sử dụng phổ biến cho count data. Số sự kiện (số hoa quả, số sản phẩm, etc) phải được đếm trong một đơn vị thời gian, không gian xác định (1 vụ mùa, 1 ngày, một ô đo đếm, etc) và các sự kiện đó phải độc lập với nhau.

Một biến ngẫu nhiên \(y\) với giá trị kì vọng \(\mu\) tuân theo phân bố poisson sẽ có hàm probability mass function (hàm khối xác suất) như sau:

\[P\left ( Y = y \right ) = \frac{e^{-\mu }\mu^y}{y!} \] Có thể nhận ra rằng phân phối poisson chỉ có 1 tham số là mean \(\mu\). Trong hồi quy poisson, mean \(\mu\) được mô tả bằng một hàm số mũ của các biến độc lập tuyến tính (linear predictor). Link function được sử dụng trong hồi quy poisson là log link. \[\mu = e^\eta\] \[\eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...+ \beta_n x_n\]

Đầu tiên chúng ta sẽ dùng mô hình poisson không có random effect để kiểm tra độ phân tán của dữ liệu (dispersion).

glm.pois <-  glm(sei ~ lbase + trt + visit , family = "poisson", data = dat) # fitting mô hình
summary(glm.pois)

## 
## Call:
## glm(formula = sei ~ lbase + trt + visit, family = "poisson", 
##     data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.4185  -1.4334  -0.2444   0.7166  10.9926  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.76538    0.03823  46.174  < 2e-16 ***
## lbase         1.17646    0.03086  38.121  < 2e-16 ***
## trtprogabide -0.10325    0.04532  -2.278  0.02272 *  
## visit        -0.29598    0.10148  -2.917  0.00354 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2517.83  on 235  degrees of freedom
## Residual deviance:  973.81  on 232  degrees of freedom
## AIC: 1747.7
## 
## Number of Fisher Scoring iterations: 5

library(COUNT)

## Warning: package 'COUNT' was built under R version 3.3.3

## Loading required package: msme

## Warning: package 'msme' was built under R version 3.3.3

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: lattice

## Loading required package: sandwich

P__disp(glm.pois)

## pearson.chi2   dispersion 
##  1120.845606     4.831231

Hàm *P__disp* cung cấp cho chúng ta Pearson \(\chi ^{2}\) statistic và Pearson dispersion statistic. Có thể tính thủ công Pearson dispersion statistic bằng cách lấy tổng bình phương của phần dư Pearson chia cho bậc tự do của mô hình (N quan sát - số biến trong mô hình bao gồm intercept).

poi.df <- df.residual(glm.pois) # 232
Pear.disp.sta <- sum(resid(glm.pois, type = "pearson")^2)/poi.df
Pear.disp.sta

## [1] 4.831231

Giá trị Pearson dispersion statistic > 1, có nghĩa là dữ liệu gặp vấn đề overdispersion. Xin có một lưu ý với các bạn là khi dữ liệu bị overdispersion, một giải pháp thay thế cho mô hình poisson là mô hình quasi-poisson.

glm.quasi <-  glm(sei ~ lbase + trt + visit , family = "quasipoisson", data = dat)
summary(glm.quasi) %>% coef()

##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)   1.7653831 0.08403771 21.007035 1.362195e-55
## lbase         1.1764626 0.06783359 17.343363 9.064533e-44
## trtprogabide -0.1032456 0.09961567 -1.036439 3.010764e-01
## visit        -0.2959814 0.22304668 -1.326993 1.858156e-01

Hệ số hồi quy của mô hình glm.poi và glm.quasi là giống nhau, standard error (SE) của các hệ số hồi quy ở mô hình glm.pois lớn hơn so với SE ở mô hình glm.pois.

Chúng ta có thể tính toán sự sai khác SE giữa 2 mô hình dựa vào Pearson dispersion statistic. Khi nhân SE của các hệ số hồi quy ở mô hình poisson với căn bậc hai của Pearson dispersion statistic chúng ta sẽ có được SE tương ứng trong mô hình quasi-poisson.

se.quasi <- summary(glm.quasi) %>% coef(. )%>% apply(., 1, "[", 2)  
se.poi <- summary(glm.pois) %>% coef(.) %>% apply(., 1, "[", 2)  
se.trans <- se.poi * sqrt(Pear.disp.sta) 
se.both = data.frame(se.poi, se.trans,  se.quasi) %>% round(., 5)
se.both

##               se.poi se.trans se.quasi
## (Intercept)  0.03823  0.08404  0.08404
## lbase        0.03086  0.06783  0.06783
## trtprogabide 0.04532  0.09962  0.09962
## visit        0.10148  0.22305  0.22305

Mixed negative binomial.

Mô hình NB có một số dạng riêng biệt, và ở phần này chúng ta sẽ tập trung dạng phổ biến của nó: quadratic negative binomial. Mô hình NB này có mean và variance như sau: \[\ mean = \mu\] \[variance = \mu + \alpha \mu^2\] Từ quadratic nhằm chỉ số mũ của mean \(\mu\) trong công thức tính variance. Mô hình negative binomial còn được gọi mô hình poisson-gamma distribution mixture model. Với biến \(y\) tuân theo phân phối NB thì chúng ta cần lưu ý những điểm sau:

Giả định \(y\) tuân theo phân phối poisson \(y \sim Poisson (\mu )\), với mean = variance = \(\mu\).
Mean \(\mu\) sẽ tuân theo phân phối gamma với mean = \(\mu\) và variance = \(\frac{\mu^2}{v}\) Khi NB model được mô hình hóa dưới dạng kết hợp của gamma-poisson thì variance của mô hình NB sẽ là bằng tổng variance của phân phối poisson và phân phối gamma. \[variance = \mu + \frac{\mu^2}{v}\] Thay \(\alpha\) = \(\frac{1}{v}\) thì chúng ta sẽ có được variance của mô hình NB \(variance = \mu + \alpha \mu^2\)

Ở đây \(v\) được gọi là shape parameter trong phân phối gamma và nó cũng là indirect dispersion parameter trong phân phối NB; \(\alpha\) là direct dispersion parameter. Như vậy đối với mô hình NB, có 2 tham số cần được ước tính là mean \(\mu\) và dispersion parameter (\(v\) hoặc \(\alpha\)) . Lưu ý răng khi sử dụng phần mềm R, thì \(v\) được gọi bằng tên thay thế là \(\theta\) ở trong phần kết quả của mô hình. Một số phần mềm khác như Stata, SAS không tính toán \(\theta\) (hay \(v\)) mà sẽ tính \(\alpha\) và gọi đó là dispersion parameter trong phần kết quả. Khi \(\theta \mapsto \ 0\), thì \(\alpha\rightarrow \infty\), và điều này cho thấy dữ liệu bị phân tán nghiêm trọng.

Mô hình mixed NB với random effect theo trường phái frequentist có thể được fit với package gamlss hoặc glmmADMB. Để đơn giản về mặt cú pháp, chúng ta sẽ dùng package glmmADMB. Với bộ dữ liệu epil, chúng ta sẽ xem xét ảnh hưởng của trt, lbase và visit lên số lần bị động kinh với random intercept cho từng bệnh nhân. Mô hình mixed NB với random intercept chúng ta cần ước tính có dạng: \[\mu_{ij} = e^\left ( \beta_{0} + \beta_{trt} x_{trt} + \beta_{lbase} x_{lbase} + \beta_{visit} x_{visit} + sub_i + \varepsilon_{ij} \right )\] \(sub_i\) thể hiện giá trị random intercept cho bệnh nhân thứ \(i^{th}\); \(e^{\varepsilon_{ij}}\) tuân theo phân bố gamma, có mean = 1 và variance là giá trị direct dispersion parameter \(\alpha\).

library(glmmADMB)

## 
## Attaching package: 'glmmADMB'

## The following object is masked from 'package:MASS':
## 
##     stepAIC

## The following object is masked from 'package:stats':
## 
##     step

dat$subject <-  as.factor(dat$subject)
nb.admb <-  glmmadmb(sei ~ trt+lbase+visit + (1|subject), zeroInflation = FALSE, family = "nbinom" , data = dat )
summary(nb.admb)

## 
## Call:
## glmmadmb(formula = sei ~ trt + lbase + visit + (1 | subject), 
##     data = dat, family = "nbinom", zeroInflation = FALSE)
## 
## AIC: 1264.7 
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.818      0.108   16.83   <2e-16 ***
## trtprogabide   -0.334      0.151   -2.21    0.027 *  
## lbase           1.011      0.101   10.03   <2e-16 ***
## visit          -0.273      0.167   -1.64    0.101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of observations: total=236, subject=59 
## Random effect variance(s):

## Warning in .local(x, sigma, ...): 'sigma' and 'rdig' arguments are present
## for compatibility only: ignored

## Group=subject
##             Variance StdDev
## (Intercept)   0.2394 0.4893
## 
## Negative binomial dispersion parameter: 7.4574 (std. err.: 1.7549)
## 
## Log-likelihood: -626.353

Tính incidence rate ratio (IRR) và khoảng tin cậy 95%.

IRR <-  cbind(est.IRR =  coef(nb.admb)  %>% exp(.),
              CI95 = confint(nb.admb) %>% exp(.)) %>% round(., 3)
IRR

##              est.IRR 2.5 % 97.5 %
## (Intercept)    6.159 4.984  7.612
## trtprogabide   0.716 0.533  0.963
## lbase          2.748 2.255  3.348
## visit          0.761 0.549  1.055

Trt và lbase ảnh hưởng có ý nghĩa thống kê lên số lần bị động kinh của các bệnh nhân trong thời gian theo dõi. Số lần bị động kinh ở bệnh nhân được điều trị bằng thuốc progabide giảm khoảng 28% so với các bệnh nhân dùng giả dược. Mô hình cho ta biết dispersion parameter, ở đây là \(\theta\) bằng 7.457. Variance của random intercept không quá lớn, điều này cho thấy có ít sự khác biệt giữa các bệnh nhân.

library(viridis)
new.df <-  data.frame(lbase = rep(seq(min(dat$lbase),max(dat$lbase), length.out =  50), times = 2), trt = rep(c(0,1), each = 50),visit = rep(mean(dat$visit), 100))
pred.count <-  predict(nb.admb, newdata = new.df, type = "link", se.fit = TRUE ) %>% map_df(~exp(.x))
new.df2 <- cbind(new.df, pred.count)
#
g.count <- ggplot(new.df2, aes(x = lbase, y = fit, color = as.factor(trt))) + geom_line(size = 1) +
  geom_ribbon( aes(ymin  = fit - 1.96*se.fit, ymax = fit + 1.96*se.fit, fill = as.factor(trt)), alpha = 0.6, color = "white") +
  labs(x = "lbase", y ="Predicted count")  +
  scale_fill_viridis(name = "Treatment", labels =c("Placebo","Progabide"),discrete=T,option="D", direction = 1) +
  scale_color_viridis(discrete= T,option="D", direction = 1) +
  guides(color =FALSE) + theme_bw() +
  theme(axis.title = element_text(size = 15, color = "black"), axis.text = element_text(size = 14, color = "black")) +
  theme(panel.grid.major = element_blank())
         
g.count