Giới thiệu package broom.mixed

1 Giới thiệu

Mixed model là một dạng mô hình hồi quy cho phép khảo sát đồng thời hiệu ứng chính (fixed effect) và hiệu ứng ngẫu nhiên (random effect) với nhiều ứng dụng trong nghiên cứu y học. Phương pháp này được hỗ trợ rất sâu và rộng trong ngôn ngữ R,đây là một lợi thế. Tuy nhiên, đây cũng chính là trở ngại quan trọng - vì có quá nhiều R packages để dựng mixed model, mỗi loại dùng cú pháp khác nhau, và cho ra kết quả với nội dung và cấu trúc khác nhau (về class, về các thành phần tham số, hình thức trình bày, suy diễn thống kê…), đa phần rối rắm khó diễn giải và bắt buộc dùng những phương pháp hậu kiểm, phân tích chuyên biệt và tương thích với mỗi package.

Ngày 21/10 vừa qua, package broom.mixed được công bố trên CRAN, và Nhi sẽ giới thiệu với các bạn về nó trong bài này. Đây là một công cụ tiện lợi, cho phép chúng ta diễn giải và khai thác nội dung của những mô hình Mixed model dựng bằng bất cứ packages nào theo cùng một cách thức. Như vậy, công dụng của broom.mixed tương tự như package broom mà Nhi từng giới thiệu, nhưng chuyên biệt cho mixed model.

Những tính năng tiện dụng của broom.mixed sẽ lần lượt được minh họa dưới đây

2 Phạm vi ứng dụng

Package broom.mixed tương thích với tất cả những định dạng mixed model từng xuất hiện trong ngôn ngữ R, bao gồm phái frequentist với họ mô hình lmer/glmer (package lme4), lme/gls/nlme (package nlme), glmmTMB glmmadmb (package glmmTMB, glmmADMB), và phái Bayes: brmsfit (package brms), stanreg (package rstanarm), MCMCglmm, thậm chí nó đọc được cả loại mô hình GAMLSS (package gamlss) vốn là kiểu định dạng mixed model khó chịu nhất.

Công dụng của package broom.mixed liên quan đến 3 phương pháp (function) :

1. glance: tóm tắt phẩm chất của mixed model, nó cung cấp thông tin khái quát về sự phù hợp của model với dữ liệu, với những tiêu chí như AIC, LOOIC, R2.

2. augment: phân tích sâu về tính chính xác của mô hình, cụ thể là sai số tồn lưu (residual error), giá trị dự báo (predicted values) trên chính dữ liệu hiện hành hay dữ liệu gốc, và các chỉ số khác (influence…).Dữ liệu xuất ra còn có thể dùng để trình bày biểu đồ để so sánh giá trị dự báo và quan sát thực tế.

3. Cuối cùng, phương pháp tidy; cho phép trình bày nội dung mô hình nhằm phục vụ cho suy diễn thống kê (estimate, s.e, CI (cho frequentist) hay HDI (phái Bayes), p_values Wald test…), đặc biệt người dùng có thể tùy chỉnh ngưỡng trên/dưới của CI, tính CI bằng bootstrap,chọn Mean hay median cho estimate… Những thông số này được chuyển thành tibble, nên có thể được dùng để vẽ biểu đồ dễ dàng bằng ggplot2

package	object	tidy	glance	augment	effects.fixed	effects.ran_vals	effects.ran_pars	effects.ran_coefs	confint.â…Ww.aldâ..	confint.â..profileâ..	confint.â..bootâ..	component.zi	component.disp	covstruct
lme4	glmer	y	y	y	y	y	y	y	NA	y	NA	NA	NA	?
lme4	lmer	y	y	y	y	y	y	y	NA	y	NA	NA	NA	?
nlme	lme	y	y	y	y	y	y	y	NA	n	NA	NA	?	?
nlme	gls	y	y	y	y	NA	NA	NA	NA	n	NA	NA	?	?
nlme	nlme	y	y	y	y	y	n	y	NA	n	NA	NA	?	?
glmmTMB	glmmTMB	y	y	y	y	y	y	n	NA		NA	y	y	?
glmmADMB	glmmadmb	y	y	y	y	y	y	n	NA		NA	y	?	?
brms	brmsfit	y	y	y	y	y	y	n	NA		NA	y	?	?
rstanarm	stanreg	y	y	y	y	y	y	n	NA		NA	NA	NA	?
MCMCglmm	MCMCglmm	y	y	y	y	y	y	n	NA		NA	?	?	?
TMB	TMB	y	n	n	n	n	n	n	NA		NA	NA	NA	?
INLA		n	n	n	n	n	n	n	NA		NA	?	?	?

3 Thí dụ minh họa

library(tidyverse)
library(broom.mixed)

Đầu tiên, Nhi sẽ thử broom.mixed trên 2 package thông dụng nhất cho mixed model, đó là nlme và lme4. Với package nlme, ta thử dựng một mô hình mixed model với random slope cho dataset “sleep study” quen thuộc;

Như các bạn thấy: nội dung nguyên thủy của mô hình nlme bao gồm 3 phần; hiệu ứng random được khảo sát dưới dạng within group residuals, còn kết quả hiệu ứng chính gần giống như mô hình glm. Ngoài ra còn có các tiêu chí phẩm chất mô hình như AIC, BIC và LogLikelihood. Kết quả này có quá nhiều chi tiết thừa thải, không có ích và thậm chí khó hiểu đối với nhiều bạn.

data("sleepstudy", package="lme4")

nlme_model <- nlme::lme(Reaction ~ Days, random =~ Days|Subject, data=sleepstudy)

summary(nlme_model)

## Linear mixed-effects model fit by REML
##  Data: sleepstudy 
##        AIC      BIC    logLik
##   1755.628 1774.719 -871.8141
## 
## Random effects:
##  Formula: ~Days | Subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 24.740241 (Intr)
## Days         5.922103 0.066 
## Residual    25.591843       
## 
## Fixed effects: Reaction ~ Days 
##                 Value Std.Error  DF  t-value p-value
## (Intercept) 251.40510  6.824516 161 36.83853       0
## Days         10.46729  1.545783 161  6.77151       0
##  Correlation: 
##      (Intr)
## Days -0.138
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.95355735 -0.46339976  0.02311783  0.46339621  5.17925089 
## 
## Number of Observations: 180
## Number of Groups: 18

Với hàm tidy, ta có thể tóm tắt những nội dung chính của mô hình như sau:

tidy(nlme_model)

## # A tibble: 6 x 8
##   effect  group   term        estimate std.error    df statistic   p.value
##   <chr>   <chr>   <chr>          <dbl>     <dbl> <dbl>     <dbl>     <dbl>
## 1 fixed   fixed   (Intercept) 251.          6.82   161     36.8   2.38e-80
## 2 fixed   fixed   Days         10.5         1.55   161      6.77  2.26e-10
## 3 ran_pa~ Subject sd_(Interc~  24.7        NA       NA     NA    NA       
## 4 ran_pa~ Subject cor_Days.(~   0.0656     NA       NA     NA    NA       
## 5 ran_pa~ Subject sd_Days       5.92       NA       NA     NA    NA       
## 6 ran_pa~ Residu~ sd_Observa~  25.6        NA       NA     NA    NA

Kết quả trình bày trong 1 bảng gọn gàng: các tham số của mô hình được chia thành nhóm: effect là nhóm hiệu ứng với Fixed = hiệu ứng chính, ran_pars là tham số cho hiệu ứng ngẫu nhiên, các hiệu ứng này lại được phân bố cho quần thể (fixed), Subject (cá thể) và residual: sai số không giải thích được, tương ứng với các terms trong mô hình bao gồm Intercept, biến độc lập, sd của random effects, correlation và rsd theo cá thể. Với fixed effect, kết quả gồm Mean, SE, df, kiểm định Wald và p.values.

Hình thức trình bày này dễ hiểu, trong sáng hơn nhiều so với nội dung gốc.

Bạn có thể thêm bớt arguments để lựa chọn: Chỉ xuất kết quả cho hiệu ứng chính, và thêm khoảng tin cậy 90% (hay một ngưỡng cao hơn, nếu bạn thích)

tidy(nlme_model, effects = "fixed", conf.int = T)

## # A tibble: 2 x 8
##   term      estimate std.error    df statistic  p.value conf.low conf.high
##   <chr>        <dbl>     <dbl> <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Interce~    251.       6.82   161     36.8  2.38e-80   238.       265. 
## 2 Days          10.5      1.55   161      6.77 2.26e-10     7.41      13.5

Không chỉ khảo sát hiệu ứng ngẫu nhiên dưới hình thức parameters, hàm tidy còn cho phép ước lượng giá trị random effect cho từng cá thể:

tidy(nlme_model, 
     effects = "ran_vals")%>%head()

## # A tibble: 6 x 3
##   level term        estimate
##   <chr> <chr>          <dbl>
## 1 308   (Intercept)     2.26
## 2 309   (Intercept)   -40.4 
## 3 310   (Intercept)   -39.0 
## 4 330   (Intercept)    23.7 
## 5 331   (Intercept)    22.3 
## 6 332   (Intercept)     9.04

Hàm glance sẽ xuất riêng giá trị LogLikelihood, AIC, BIC và sigma

glance(nlme_model)

## # A tibble: 1 x 4
##   sigma logLik   AIC   BIC
##   <dbl>  <dbl> <dbl> <dbl>
## 1  25.6  -872. 1756. 1775.

Hàm augment cho phép xuất giá trị dự báo cho chính dữ liệu gốc, hoặc áp dụng cho một dataset mớivà xuất kết quả dưới dạng dataframe: tương tự như hàm predict

augment(nlme_model, newdata=sample_frac(sleepstudy,0.2))%>%head()

## # A tibble: 6 x 5
##   .rownames Reaction  Days Subject .fitted
##   <chr>        <dbl> <dbl> <fct>     <dbl>
## 1 125           272.     4 351        286.
## 2 155           268.     4 370        287.
## 3 1             250.     0 308        254.
## 4 177           334.     6 372        334.
## 5 106           270.     5 349        284.
## 6 59            330.     8 332        342.

Dữ liệu này cho phép vẽ những biểu đồ như sau để kiểm tra trực quan hoạt động của mô hình:

augment(nlme_model, newdata=sample_frac(sleepstudy,0.2))%>%
  gather(Reaction,.fitted,key="Type",value="Score")%>%
  ggplot()+
  geom_density(aes(x=Score,fill=Type),alpha=0.5)+
  theme_bw()

augment(nlme_model, newdata=sample_frac(sleepstudy,0.2))%>%
  gather(Reaction,.fitted,key="Type",value="Score")%>%
  ggplot()+
  geom_smooth(aes(x=Days,y=Score,fill=Type,col=Type),alpha=0.2,method="glm")+
  theme_bw()

Bây giờ chúng ta qua package lme4, và thử dựng một logistic mixed model cho dataset “cbpp”

Như ta thấy, nội dung nguyên thủy của mô hình lme4 được trình bày hoàn toàn khác so với mô hình nlme, các tham số fixed effect được khảo sát bằng SE và test Z, trong khi random effect được trình bày dưới dạng variance và SD của intercept. Ngoài ra còn có thêm một correlation matrix cho fixed effects

data("cbpp", package="lme4")

lme4_logistic <- lme4::glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)

summary(lme4_logistic)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(incidence, size - incidence) ~ period + (1 | herd)
##    Data: cbpp
## 
##      AIC      BIC   logLik deviance df.resid 
##    194.1    204.2    -92.0    184.1       51 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3816 -0.7889 -0.2026  0.5142  2.8791 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  herd   (Intercept) 0.4123   0.6421  
## Number of obs: 56, groups:  herd, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3983     0.2312  -6.048 1.47e-09 ***
## period2      -0.9919     0.3032  -3.272 0.001068 ** 
## period3      -1.1282     0.3228  -3.495 0.000474 ***
## period4      -1.5797     0.4220  -3.743 0.000182 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) perid2 perid3
## period2 -0.363              
## period3 -0.340  0.280       
## period4 -0.260  0.213  0.198

Ta thử dùng hàm tidy, với 2 yêu cầu phụ:

1. Áp dụng hàm exponential để tính Odds-ratios.

2. Xuất 97.5% CI bằng phương pháp bootstrap .

Quy trình bootstrap mất một ít thời gian, sau đây là kết quả:

Như ta thấy, hàm tidy đã bỏ hết những chi tiết thừa, quy đồng nội dung mô hình về cùng một hình thức trình bày, cấu trúc bảng tương tự như trong mô hình nlme, các terms được phân thành loại effect và group; với giá trị đã được hoán chuyển thành Odds-ratio, Z-values giữ nguyên,estimate tương ứng vớitrung bình của Odds-ratio, SE và 97.5% CI của Odds-ratio được xác định bằng bootstrap. Kết quả này sẵn sàng để đưa vào báo cáo.

tidy(lme4_logistic,
     exponentiate=TRUE, 
     conf.int = T, 
     conf.level = 0.975, 
     conf.method ="boot")

## # A tibble: 5 x 9
##   effect group term  estimate std.error statistic  p.value conf.low
##   <chr>  <chr> <chr>    <dbl>     <dbl>     <dbl>    <dbl>    <dbl>
## 1 fixed  <NA>  (Int~    0.247    0.0571     -6.05  1.47e-9   0.139 
## 2 fixed  <NA>  peri~    0.371    0.112      -3.27  1.07e-3   0.189 
## 3 fixed  <NA>  peri~    0.324    0.104      -3.49  4.74e-4   0.138 
## 4 fixed  <NA>  peri~    0.206    0.0870     -3.74  1.82e-4   0.0677
## 5 ran_p~ herd  sd__~    0.642   NA          NA    NA         0.167 
## # ... with 1 more variable: conf.high <dbl>

Ta thử hàm glance để khảo sát phẩm chất mô hình:

glance(lme4_logistic)

## # A tibble: 1 x 6
##   sigma logLik   AIC   BIC deviance df.residual
##   <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1     1  -92.0  194.  204.     73.5          51

Và hàm augment để ước tính giá trị dự báo trên chính dữ liệu gốc:

augment(lme4_logistic,data=cbpp)%>%head()

##   herd incidence size period    .fitted     .resid    .fixed        .mu
## 1    1         2   14      1 -0.8087134 -1.4377078 -1.398343 0.30816472
## 2    1         3   12      2 -1.8006384  0.9884720 -2.390268 0.14177337
## 3    1         4    9      3 -1.9369296  2.3566883 -2.526559 0.12598556
## 4    1         0    5      4 -2.3884588 -0.9370227 -2.978088 0.08405701
## 5    2         3   22      1 -1.6974362 -0.2431732 -1.398343 0.15480041
## 6    2         1   18      2 -2.6893612 -0.1428294 -2.390268 0.06360406
##   .offset  .sqrtXwt  .sqrtrwt .weights     .wtres       .eta
## 1       0 1.7276542  8.103473       14 -1.3395656 -0.8087134
## 2       0 1.2083394  9.930984       12  1.0747970 -1.8006384
## 3       0 0.9954992  9.040690        9  2.8790881 -1.9369296
## 4       0 0.6204492  8.058678        5 -0.6773884 -2.3884588
## 5       0 1.6965905 12.967183       22 -0.2390730 -1.6974362
## 6       0 1.0354006 17.384575       18 -0.1399198 -2.6893612

Bây giờ chúng ta sẽ thử áp dụng broom.mixed cho mô hình Bayes, cụ thể là 2 package brms và rstanarm. Để tiết kiệm thời gian, Nhi sẽ không thực sự dùng brms nhưng tải dữ liệu một mô hình brms có sẵn trong package broom.mixed:

load(system.file("extdata", "brms_example.rda", package="broom.mixed"))

Đây là một mô hình mixed model với crossed random effect trên dataset mtcars

fit <- brms_crossedRE

fit$fit

## Inference for Stan model: 7fd04982032e6915d9435dfb9d716706.
## 2 chains, each with iter=100; warmup=50; thin=1; 
## post-warmup draws per chain=50, total post-warmup draws=100.
## 
##                           mean se_mean   sd    2.5%     25%    50%    75%
## b_Intercept              34.41    0.72 3.87   26.16   32.23  34.25  36.78
## b_wt                     -4.25    0.30 1.26   -6.90   -4.95  -4.25  -3.48
## sd_cyl__Intercept         4.88    0.30 2.98    0.76    2.47   4.38   6.62
## sd_gear__Intercept        4.70    0.78 3.23    0.27    2.43   4.13   6.47
## sd_gear__wt               1.72    0.46 1.45    0.05    0.71   1.29   2.36
## cor_gear__Intercept__wt  -0.38    0.07 0.57   -0.99   -0.83  -0.59  -0.07
## sigma                     2.53    0.07 0.40    1.77    2.25   2.47   2.77
## r_cyl[4,Intercept]        1.98    0.68 2.36   -2.30    0.33   2.15   3.58
## r_cyl[6,Intercept]       -0.92    0.64 2.27   -5.30   -2.21  -0.96   0.33
## r_cyl[8,Intercept]       -2.21    0.47 2.23   -6.95   -3.76  -1.90  -0.51
## r_gear[3,Intercept]      -3.68    0.80 4.43  -12.39   -7.55  -2.55  -0.19
## r_gear[4,Intercept]       0.99    0.34 2.52   -3.68   -0.66   0.85   2.44
## r_gear[5,Intercept]      -0.59    0.63 2.69   -6.62   -1.75  -0.11   1.00
## r_gear[3,wt]              0.96    0.31 1.27   -0.62    0.01   0.56   1.71
## r_gear[4,wt]             -0.30    0.18 1.01   -2.45   -0.98  -0.14   0.24
## r_gear[5,wt]             -0.29    0.26 1.18   -2.45   -0.92  -0.24   0.11
## lp__                    -98.19    1.27 3.85 -105.47 -101.12 -97.83 -95.01
##                          97.5% n_eff Rhat
## b_Intercept              41.78    29 1.12
## b_wt                     -1.94    17 1.12
## sd_cyl__Intercept        11.69   100 1.01
## sd_gear__Intercept       12.94    17 1.22
## sd_gear__wt               4.84    10 1.36
## cor_gear__Intercept__wt   0.85    64 1.00
## sigma                     3.27    32 1.11
## r_cyl[4,Intercept]        7.41    12 1.15
## r_cyl[6,Intercept]        2.87    13 1.18
## r_cyl[8,Intercept]        1.24    23 1.09
## r_gear[3,Intercept]       2.51    31 1.21
## r_gear[4,Intercept]       5.78    55 1.00
## r_gear[5,Intercept]       3.32    18 1.06
## r_gear[3,wt]              4.34    16 1.25
## r_gear[4,wt]              1.64    33 1.04
## r_gear[5,wt]              2.31    20 1.10
## lp__                    -92.24     9 1.57
## 
## Samples were drawn using NUTS(diag_e) at Tue Oct 16 08:42:29 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Khác với lme4 hay nlme, brms thực chất áp dụng suy diễn Bayes để tạo ra những chuỗi MCMC chứa phân phối hậu nghiệm cho mỗi tham số (terms, parameters) trong mô hình, và kết quả cũng được trình bày theo cách khác. Tuy nhiên khi áp dụng hàm tidy, kết quả mô hình Bayes cũng được tóm tắt gọn lại theo cùng một cách thức:

Kết quả này có thể được diễn giải tương tự như bất cứ mô hình nào, chỉ khác một chút về thuật ngữ: khoảng tin cậy trong thống kê Bayes được gọi là Highest posterior density interval(HPDI hay HDI).

brms_sum=tidy(fit,
     robust=TRUE, 
     conf.method="HPDinterval", 
     conf.level = 0.975)

brms_sum

## # A tibble: 7 x 8
##   effect  component group  term      estimate std.error conf.low conf.high
##   <chr>   <chr>     <chr>  <chr>        <dbl>     <dbl>    <dbl>     <dbl>
## 1 fixed   cond      <NA>   (Interce~   34.2       3.37   2.45e+1    42.1  
## 2 fixed   cond      <NA>   wt          -4.25      1.11  -7.15e+0    -1.37 
## 3 ran_pa~ cond      cyl    sd__(Int~    4.38      2.98   4.37e-1    12.2  
## 4 ran_pa~ cond      gear   sd__(Int~    4.13      2.90   1.96e-1    13.4  
## 5 ran_pa~ cond      gear   sd__wt       1.29      1.07   3.43e-3     7.18 
## 6 ran_pa~ cond      gear   cor__(In~   -0.595     0.435 -9.97e-1     0.895
## 7 ran_pa~ cond      Resid~ sd__Obse~    2.47      0.387  1.73e+0     3.47

Lợi thế của dữ liệu dạng bảng này đó là ta có thể đưa nó vào pipe để vẽ một biểu đồ tương tự như hình bên dưới:

brms_sum%>%ggplot(aes(estimate,term,xmin=conf.low,xmax=conf.high,col=estimate))+
  geom_errorbarh(height=0)+
  geom_vline(xintercept=0,lty=2,col="black")+
  geom_point(size=3)+
  scale_color_gradientn(colours = pals::coolwarm(10))+
  theme_bw()

Mô hình mixed model Bayes cho phép ước lượng được cả phân phối hậu nghiệm đến mức độ cá thể:

brms_ran_val=tidy(fit, effects = "ran_vals", 
     robust = T,
     conf.method="HPDinterval", 
     conf.level = 0.975) 

brms_ran_val

## # A tibble: 9 x 9
##   effect component group level term  estimate std.error conf.low conf.high
##   <chr>  <chr>     <chr> <chr> <chr>    <dbl>     <dbl>    <dbl>     <dbl>
## 1 ran_v~ cond      cyl   4     (Int~    2.15      2.43     -2.80      8.35
## 2 ran_v~ cond      cyl   6     (Int~   -0.963     1.95     -5.54      4.39
## 3 ran_v~ cond      cyl   8     (Int~   -1.90      2.39     -7.45      1.40
## 4 ran_v~ cond      gear  3     (Int~   -2.55      4.31    -15.2       3.34
## 5 ran_v~ cond      gear  4     (Int~    0.846     2.28     -4.88      7.16
## 6 ran_v~ cond      gear  5     (Int~   -0.110     2.10     -7.89      4.40
## 7 ran_v~ cond      gear  3     wt       0.565     1.01     -1.21      4.61
## 8 ran_v~ cond      gear  4     wt      -0.137     0.733    -2.58      3.06
## 9 ran_v~ cond      gear  5     wt      -0.243     0.661    -2.55      3.15

brms_ran_val%>%
  ggplot(aes(estimate,level,xmin=conf.low,xmax=conf.high,col=group))+
  geom_errorbarh(height=0)+
  geom_vline(xintercept=0,lty=2,col="black")+
    facet_grid(term~group,scale="free_x")+
  geom_point(size=3)+
  theme_bw()

Bây giờ ta thử một mô hình Bayes mixed model khác với package rstanarm:

library(rstanarm)

stan_fit <- stan_glmer(mpg ~ wt + (1|cyl) + (1+wt|gear), data = mtcars,
                  iter = 300, chains = 2)

## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 300 [  0%]  (Warmup)
## Iteration:  30 / 300 [ 10%]  (Warmup)
## Iteration:  60 / 300 [ 20%]  (Warmup)
## Iteration:  90 / 300 [ 30%]  (Warmup)
## Iteration: 120 / 300 [ 40%]  (Warmup)
## Iteration: 150 / 300 [ 50%]  (Warmup)
## Iteration: 151 / 300 [ 50%]  (Sampling)
## Iteration: 180 / 300 [ 60%]  (Sampling)
## Iteration: 210 / 300 [ 70%]  (Sampling)
## Iteration: 240 / 300 [ 80%]  (Sampling)
## Iteration: 270 / 300 [ 90%]  (Sampling)
## Iteration: 300 / 300 [100%]  (Sampling)
## 
##  Elapsed Time: 0.357 seconds (Warm-up)
##                0.109 seconds (Sampling)
##                0.466 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 300 [  0%]  (Warmup)
## Iteration:  30 / 300 [ 10%]  (Warmup)
## Iteration:  60 / 300 [ 20%]  (Warmup)
## Iteration:  90 / 300 [ 30%]  (Warmup)
## Iteration: 120 / 300 [ 40%]  (Warmup)
## Iteration: 150 / 300 [ 50%]  (Warmup)
## Iteration: 151 / 300 [ 50%]  (Sampling)
## Iteration: 180 / 300 [ 60%]  (Sampling)
## Iteration: 210 / 300 [ 70%]  (Sampling)
## Iteration: 240 / 300 [ 80%]  (Sampling)
## Iteration: 270 / 300 [ 90%]  (Sampling)
## Iteration: 300 / 300 [100%]  (Sampling)
## 
##  Elapsed Time: 0.39 seconds (Warm-up)
##                0.263 seconds (Sampling)
##                0.653 seconds (Total)

summary(stan_fit)

## 
## Model Info:
## 
##  function:     stan_glmer
##  family:       gaussian [identity]
##  formula:      mpg ~ wt + (1 | cyl) + (1 + wt | gear)
##  algorithm:    sampling
##  priors:       see help('prior_summary')
##  sample:       300 (posterior sample size)
##  observations: 32
##  groups:       cyl (3), gear (3)
## 
## Estimates:
##                                       mean   sd     2.5%   25%    50% 
## (Intercept)                           32.8    3.5   25.7   30.7   32.9
## wt                                    -4.0    1.3   -7.2   -4.6   -3.8
## b[(Intercept) cyl:4]                   2.6    2.3   -1.6    1.1    2.4
## b[(Intercept) cyl:6]                  -0.5    2.1   -4.9   -1.7   -0.3
## b[(Intercept) cyl:8]                  -1.7    2.3   -7.2   -2.8   -1.5
## b[(Intercept) gear:3]                 -1.1    2.4   -7.6   -0.9   -0.1
## b[wt gear:3]                           0.4    1.1   -1.2    0.0    0.1
## b[(Intercept) gear:4]                  0.6    1.5   -1.1   -0.1    0.1
## b[wt gear:4]                           0.0    0.9   -2.1   -0.2    0.0
## b[(Intercept) gear:5]                  0.2    1.6   -1.9   -0.2    0.0
## b[wt gear:5]                          -0.3    1.0   -2.6   -0.6   -0.1
## sigma                                  2.6    0.3    2.1    2.4    2.6
## Sigma[cyl:(Intercept),(Intercept)]    12.1   18.8    0.1    2.7    5.8
## Sigma[gear:(Intercept),(Intercept)]    3.2    7.4    0.0    0.1    0.3
## Sigma[gear:wt,(Intercept)]            -0.4    1.6   -5.0   -0.5    0.0
## Sigma[gear:wt,wt]                      1.5    3.9    0.0    0.0    0.2
## mean_PPD                              20.1    0.7   18.9   19.7   20.1
## log-posterior                       -108.4    3.5 -115.3 -110.7 -108.1
##                                       75%    97.5%
## (Intercept)                           34.7   39.7 
## wt                                    -3.3   -1.8 
## b[(Intercept) cyl:4]                   3.9    7.4 
## b[(Intercept) cyl:6]                   0.8    3.8 
## b[(Intercept) cyl:8]                  -0.4    2.1 
## b[(Intercept) gear:3]                  0.1    1.2 
## b[wt gear:3]                           0.6    3.4 
## b[(Intercept) gear:4]                  0.9    4.8 
## b[wt gear:4]                           0.2    1.4 
## b[(Intercept) gear:5]                  0.4    3.6 
## b[wt gear:5]                           0.0    1.5 
## sigma                                  2.8    3.3 
## Sigma[cyl:(Intercept),(Intercept)]    13.3   60.0 
## Sigma[gear:(Intercept),(Intercept)]    2.5   24.0 
## Sigma[gear:wt,(Intercept)]             0.0    1.2 
## Sigma[gear:wt,wt]                      1.5    9.2 
## mean_PPD                              20.7   21.4 
## log-posterior                       -105.8 -103.0 
## 
## Diagnostics:
##                                     mcse Rhat n_eff
## (Intercept)                         0.3  1.0  162  
## wt                                  0.2  1.0   36  
## b[(Intercept) cyl:4]                0.2  1.0   86  
## b[(Intercept) cyl:6]                0.2  1.0   89  
## b[(Intercept) cyl:8]                0.2  1.0  103  
## b[(Intercept) gear:3]               0.4  1.1   30  
## b[wt gear:3]                        0.2  1.1   28  
## b[(Intercept) gear:4]               0.2  1.0   62  
## b[wt gear:4]                        0.1  1.0   67  
## b[(Intercept) gear:5]               0.2  1.0   95  
## b[wt gear:5]                        0.1  1.0   62  
## sigma                               0.0  1.0  300  
## Sigma[cyl:(Intercept),(Intercept)]  2.2  1.0   72  
## Sigma[gear:(Intercept),(Intercept)] 1.3  1.0   33  
## Sigma[gear:wt,(Intercept)]          0.1  1.0  166  
## Sigma[gear:wt,wt]                   0.5  1.0   57  
## mean_PPD                            0.0  1.0  238  
## log-posterior                       0.3  1.0  129  
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Nội dung quá dài dòng này sẽ được tóm gọn bằng hàm tidy, theo cùng một cách như ta đã biết

tidy(stan_fit, conf.int = TRUE, conf.method = "HPDinterval", prob = 0.95)

## # A tibble: 2 x 5
##   term        estimate std.error conf.low conf.high
##   <chr>          <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)    32.9      3.00     26.6      38.1 
## 2 wt             -3.85     0.984    -6.10     -2.35

tidy(stan_fit, effects = "ran_pars")

## # A tibble: 5 x 3
##   term                    group    estimate
##   <chr>                   <chr>       <dbl>
## 1 sd_(Intercept).cyl      cyl         3.47 
## 2 sd_(Intercept).gear     gear        1.79 
## 3 sd_wt.gear              gear        1.23 
## 4 cor_(Intercept).wt.gear gear       -0.193
## 5 sd_Observation.Residual Residual    2.63

stan_fit%>%tidy(effects = "ran_vals",
                conf.int = TRUE, 
                conf.method = "HPDinterval")%>%
  group_by(group,level) %>%
  mutate(keep=any(estimate > 0))%>%
  ggplot(aes(estimate,level,
             xmin=conf.low,xmax=conf.high,
             col=factor(keep)))+
  geom_errorbarh(height=0)+
  geom_vline(xintercept=0,lty=2)+
  geom_point(size=3)+
    facet_grid(term~group,scale="free_x")+
  scale_colour_manual(values=c("blue","red"),guide=FALSE)+
  theme_bw()

Cuối cùng, ta thử dùng hàm tidy cho mô hình dựng bằng gamlss. Như đã nói, gamlss là một package khá kì dị, kết quả mà nó xuất ra không giống với bất cứ gói nào khác, và rất khó khai thác.Trong thí dụ này, ta dựng một mô hình BoxCox t với dataset abdomen.Ghi chú: đây không phải là mixed model.

data(abdom, package="gamlss.data") 

library(gamlss)

## Loading required package: splines

## Loading required package: gamlss.data

## Loading required package: gamlss.dist

## Loading required package: MASS

## 
## Attaching package: 'MASS'

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

## Loading required package: nlme

## 
## Attaching package: 'nlme'

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

## Loading required package: parallel

##  **********   GAMLSS Version 5.1-0  **********

## For more on GAMLSS look at http://www.gamlss.org/

## Type gamlssNews() to see new features/changes/bug fixes.

gamlss_mod<- gamlss(y~pb(x),sigma.fo=~pb(x),family=BCT,data=abdom)

## GAMLSS-RS iteration 1: Global Deviance = 4771.925 
## GAMLSS-RS iteration 2: Global Deviance = 4771.039 
## GAMLSS-RS iteration 3: Global Deviance = 4770.999 
## GAMLSS-RS iteration 4: Global Deviance = 4770.994 
## GAMLSS-RS iteration 5: Global Deviance = 4770.993

summary(gamlss_mod)

## ******************************************************************
## Family:  c("BCT", "Box-Cox t") 
## 
## Call:  gamlss(formula = y ~ pb(x), sigma.formula = ~pb(x),  
##     family = BCT, data = abdom) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -72.19470    1.33349  -54.14   <2e-16 ***
## pb(x)        11.05556    0.05785  191.12   <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.65086    0.10861 -24.406  < 2e-16 ***
## pb(x)       -0.01004    0.00380  -2.641  0.00849 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.1076     0.6301  -0.171    0.865
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.4903     0.4247   5.864 7.47e-09 ***
## ---
## 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:  610 
## Degrees of Freedom for the fit:  11.76162
##       Residual Deg. of Freedom:  598.2384 
##                       at cycle:  5 
##  
## Global Deviance:     4770.993 
##             AIC:     4794.516 
##             SBC:     4846.426 
## ******************************************************************

tidy(gamlss_mod)

##   parameter        term     estimate   std.error   statistic       p.value
## 1        mu (Intercept) -72.19470143 1.327904569 -54.3673869 7.314947e-235
## 2        mu       pb(x)  11.05556002 0.057722851 191.5283099  0.000000e+00
## 3     sigma (Intercept)  -2.65086130 0.108135709 -24.5142083  9.885987e-93
## 4     sigma       pb(x)  -0.01003582 0.003788057  -2.6493332  8.275390e-03
## 5        nu (Intercept)  -0.10755887 0.557700059  -0.1928615  8.471317e-01
## 6       tau (Intercept)   2.49027048 0.299096364   8.3259804  5.526415e-16

4 Nhận xét

Package broom.mixed là một công cụ rất hữu ích, nó giúp giải quyết được các vấn đề trở ngại khi thực hành mixed model bao gồm:

1. Sử dụng cùng một danh pháp thuật ngữ cho những thành phần trong nội dung mô hình, dù mô hình được dựng bằng package nào. Các terms được phân chia theo thứ bậc: effect (fixed hay random), nếu là random thì là parameters hay values ? kế đến là terms (intercept, các factors, predictors trong mô hình, sd hay cor_intercept), giúp dễ dàng phân biệt ý nghĩa của các thành phần này.

2. Giản lược những chi tiết thừa thải trong nội dung mô hình, chỉ giữ lại những thông tin thiết yếu nhất để suy diễn thống kê

3. Tách riêng 3 phần thông tin về phẩm chất mô hình, nội dung mô hình và giá trị dự báo

4. Quy đồng tất cả kết quả từ 8-9 package khác nhau về cùng một bảng kết quả có cấu trúc nhất quán.

5. Kết quả được xuất ra dưới dạng tibble, cho phép thực hiện các hàm phân tích sâu, sử dụng pipeline và vẽ biểu đồ bằng ggplot.

Bài thực hành đến đây là hết. Tạm biệt các bạn và hẹn gặp lại…