Giới thiệu package broom

1 Giới thiệu

Thân chào các bạn, hôm nay Nhi muốn giới thiệu với các bạn một R package rất tiện lợi, đó là broom. Công dụng của broom là trích xuất nội dung của một output object - kết quả của một mô hình thống kê và tóm tắt thông tin này vào một dataframe. Tên gọi broom chính là để diễn tả việc dọn dẹp từ một mô hình rối rắm thành một dataframe gọn gàng, sạch đẹp, sẵn sàng để đưa vào báo cáo, hoặc làm nguyên liệu để vẽ những biểu đồ.

Package broom do tác giả David Robinson tạo ra từ năm 2014, và trong 3 năm nay nhóm tác giả của broom đã âm thầm mở rộng khả năng của package này. Ở thời điểm hiện nay, broom rất mạnh, nó tương thích với hầu hết những package Hồi quy trong R, kể cả những mô hình phức tạp như gamlss, nlmer, lme4, mgcv, brms…

broom có 4 hàm quan trọng:

1. glance: hàm này trích xuất thông tin về phẩm chất mô hình, thí dụ goodness of fit, các trị số AIC, BIC, Rsq, F test…

2. tidy: hàm này trích xuất nội dung mô hình, thí dụ các tham số hồi quy, S.E, khoảng tin cậy, kiểm định t và p_value

3. augment: hàm này cho phép kiểm định mô hình trên dataset hiện thời hoặc mới, và cung cấp thông tin về độ chính xác của mô hình, thí dụ Residual error, …

4. bootstrap: broom có một hàm bootstrap của riêng nó, khi kết hợp với hàm tidy, sẽ cho phép bạn áp dụng bootstrap cho bất kì mô hình nào. Quy trình này cần sự hỗ trợ của dplyr

Tuy đã hiện diện từ rất lâu, nhưng broom chỉ thực sự gây chú ý từ năm 2016, khi Hadley Wickham hợp nhất broom với những bạn bè của nó là các package như dplyr, purr, tidyr, và ggplot2 để tạo thành một quy trình phân tích dữ liệu khép kín, trong đó broom vẫn giữ nhiệm vụ dọn dẹp và đóng gói dataframe nội dung mô hình, nhưng nó không làm việc một mình mà nhận được sự trợ thủ của dplyr cũng như kết quả mà broom xuất ra được chuyền qua cho ggplot một cách nhanh chóng nhờ các nhân tử “pipe”.

broom rất tiện lợi, nhưng có rất ít tutorial về nó ngoại trừ 3 thí dụ minh họa của chính tác giả và một vài đoạn code trong 2nd Ed. của quyển ggplot2 của Wickham. Do đó, Nhi viết bài hướng dẫn sau với thông điệp khuyến khích các bạn sử dụng broom mỗi khi bạn làm việc với mô hình hồi quy.

Trong bài, Nhi sẽ minh họa khả năng của broom với 4 thí dụ

1. Trích xuất thông tin từ một mô hình Logistic dựng bằng Gamlss, và bootstrap mô hình logistic này

2. Trích xuất thông tin từ một mô hình hồi quy Polynomial dựng bằng glm()

3. Bootstraping một mô hình GLM với family Student t chứa 3 tham số, dựng bằng Gamlss

4. Bootstraping một mô hình Cox-PH

2 Mô hình Logistic với gamlss

Ta dựng một mô hình logistic cho dataset Biopsy, nhưng sử dụng gamlss chứ không dùng hàm glm().

library(tidyverse)

library(gamlss)
nC<-detectCores()

library(broom)

library(ggridges)

df=read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/MASS/biopsy.csv")%>%as_tibble()%>%.[,c(3:12)]%>%na.omit()

names(df)=c("clumpthickness",
            "SizeUniformity",
            "ShapeUniformity",
            "Margin_adhesion",
            "EpiCellSize",
            "Barenuclei",
            "BlandChromatin",
            "NormalNucleoli",
            "Mitoses",
            "Class"
)

df2=df%>%mutate(.,Class=as.integer(.$Class)-1L)

glslogit=gamlss(Class~clumpthickness+
                     ShapeUniformity+
                     Margin_adhesion+
                     Barenuclei+
                     BlandChromatin,
                   data=df2,
                   family=BI(),
                   trace=F,
                parallel="multicore",
                ncpus = nC)


summary(glslogit)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = Class ~ clumpthickness + ShapeUniformity +  
##     Margin_adhesion + Barenuclei + BlandChromatin,  
##     family = BI(), data = df2, trace = F, parallel = "multicore",  
##     ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -9.74114    1.04989  -9.278  < 2e-16 ***
## clumpthickness   0.62576    0.13373   4.679 3.48e-06 ***
## ShapeUniformity  0.48994    0.15379   3.186 0.001510 ** 
## Margin_adhesion  0.33918    0.11221   3.023 0.002599 ** 
## Barenuclei       0.37330    0.09381   3.979 7.65e-05 ***
## BlandChromatin   0.55731    0.16341   3.411 0.000687 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  683 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  677 
##                       at cycle:  2 
##  
## Global Deviance:     112.566 
##             AIC:     124.566 
##             SBC:     151.725 
## ******************************************************************

Ta lần lượt thử hàm tidy trên mô hình này: Lưu ý: kết quả của hàm tidy luôn là 1 dataframe, do đây là mô hình logistic nên hàm augment không cho kết quả chính xác, bạn phải dùng mutate để tính probability từ predicted bằng hàm plogis.

tidy(glslogit)%>%knitr::kable()

parameter	term	estimate	std.error	statistic	p.value
mu	(Intercept)	-9.7411399	1.0498889	-9.278258	0.0000000
mu	clumpthickness	0.6257619	0.1337317	4.679233	0.0000035
mu	ShapeUniformity	0.4899369	0.1537879	3.185797	0.0015099
mu	Margin_adhesion	0.3391770	0.1122074	3.022767	0.0025993
mu	Barenuclei	0.3732985	0.0938082	3.979382	0.0000766
mu	BlandChromatin	0.5573061	0.1634060	3.410562	0.0006868

Với kết quả hàm tidy mà broom xuất ra, ta có thể tính vẽ biểu đồ Odds-ratio dễ dàng:

logitsum=glslogit%>%tidy()%>%mutate(OR=exp(estimate),
                                    UL=exp(estimate+1.96*std.error),
                                    LL=exp(estimate-1.96*std.error))

logitsum%>%dplyr::select(term,OR,UL,LL)%>%ggplot(aes(x=term,
                                                     y=OR,
                                                     ymin=LL,
                                                     ymax=UL,
                                                     col=term))+
  geom_hline(yintercept = c(1,2),linetype=2,col=c("blue","red"))+
  geom_pointrange(show.legend = F,size=1)+
  coord_flip()+
  theme_bw()

Trong R có một package khác là sjPlot cho phép làm những việc như thế này chỉ với 1 hàm plot_model, tuy nhiên sjPlot hoàn toàn bất lực với mô hình gamlss, và khi làm thủ công mọi chuyện, bạn sẽ kiểm soát tốt hơn tất cả mọi thứ, từ việc tính OR, 95%CI cho đến tùy chỉnh màu sắc cho ggplot; và làm thủ công không khó lắm khi broom đã giúp bạn làm phần khó nhất là chuyển model thành dataframe. Có dataframe là mọi việc đều khả thi.

Tiếp theo, ta sẽ thử bootstrap mô hình logistic trên đây và một lần nữa, xác định lại khoảng tin cậy cho OR dựa vào 100 mô hình bootstrap:

Để nối hàm bootstrap của broom vào hàm tidy, có 2 cách, Nhi giới thiệu cách dễ trước, đó là dùng hàm do của dplyr. Lợi thế của hàm bootstrap này đó là:

1. Bạn không cần phải dùng package boot, không cần viết hàm gì cả
2. Bạn có thể bootstrap mọi mô hình bạn muốn, từ lme4 cho đến survival…
3. Kết quả bootstrap được hàm tidy đóng gói luôn, rất gọn gàng và sẵn sàng để vẽ biểu đồ

Kết quả bootstrap như sau:

df2$Class=as.numeric(df2$Class)

bootlogit<- df2 %>%
  bootstrap(100) %>%
  do(tidy(gamlss(.$Class~.$clumpthickness+
                   .$ShapeUniformity+
                   .$Margin_adhesion+
                   .$Barenuclei+
                   .$BlandChromatin,
                 family=BI(),
                 trace=F,
                 parallel="multicore",
                 ncpus = nC)))

bootlogit<-mutate(bootlogit,OR=exp(estimate))

bootlogit%>%dplyr::select(term,OR)%>%
  ggplot(aes(y=term,
             x=OR))+
  geom_density_ridges(aes(fill=term),alpha=0.5,scale=1.5)+
  geom_vline(xintercept = c(1,2),linetype=2,col=c("blue","red"))+
  theme_bw()

muDF<-bootlogit%>%
  dplyr::select(replicate,term,estimate)%>%
  spread(key=term,value=estimate)

names(muDF)=c("replicate",
              "(Intercept)",
              "Barenuclei",
              "BlandChromatin",
              "clumpthickness",
              "Margin_adhesion",
              "ShapeUniformity")

predboot<-df%>%mutate(pred=NA,Iter=NA)

for(i in 1:nrow(muDF)){
  beta0=muDF$`(Intercept)`[i]
  beta1=muDF$Barenuclei[i]
  beta2=muDF$BlandChromatin[i]
  beta3=muDF$clumpthickness[i]
  beta4=muDF$Margin_adhesion[i]
  beta5=muDF$ShapeUniformity[i]
  tempdf<-df%>%mutate(pred=plogis(beta0+
                           ShapeUniformity*beta5+
                           Margin_adhesion*beta4+
                           Barenuclei*beta1+
                           BlandChromatin*beta2+
                           clumpthickness*beta3),
                         Iter=i)
  predboot=rbind(predboot,tempdf)%>%na.omit()
}


binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se=T,show.legend = F,...)
}

longdf=predboot%>%dplyr::select(Iter,
                         ShapeUniformity,
                         Margin_adhesion,
                         Barenuclei,
                         BlandChromatin,
                         clumpthickness,
                         pred,
                         Class)%>%
  gather(ShapeUniformity:clumpthickness,key="Feature",value="Score")

  longdf%>%ggplot(aes(x=Score,
                      y=pred,
                      group=as.factor(Iter),
                      fill=as.factor(Iter),
                      col=as.factor(Iter)))+
  binomial_smooth(alpha=0.1)+
  theme_bw()+
  facet_wrap(~Feature,scales="free",ncol=3)+
    scale_fill_brewer("PuRd")+
    scale_color_brewer("PuRd")

3 Mô hình Polynomial bậc 5 dùng GLM

Thí dụ tiếp theo, Nhi dựng một mô hình polynomial bậc 5 khảo sát tương quan giữa log(chiều cao) và tuổi của các bệnh nhi trong dataset Leukemia.

Leukemia%>%ggplot(aes(x=age,y=log(height)))+
  geom_point(col="blue",alpha=0.3)+
  geom_smooth(method="glm",formula = y ~ poly(x,5),col="red",fill="red",alpha=0.5)+
  theme_bw()

Mô hình đươc dựng đơn giản bằng hàm glm, với family = gaussian và link function = log:

lmod=glm(Leukemia,formula=height~poly(age,5),family=gaussian(link="log"))

summary(lmod)

## 
## Call:
## glm(formula = height ~ poly(age, 5), family = gaussian(link = "log"), 
##     data = Leukemia)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -21.9027   -3.4454   -0.0024    3.5268   24.6226  
## 
## Coefficients:
##                Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)    4.815358   0.001015 4743.948  < 2e-16 ***
## poly(age, 5)1  6.326672   0.044497  142.181  < 2e-16 ***
## poly(age, 5)2 -1.250836   0.044197  -28.301  < 2e-16 ***
## poly(age, 5)3 -0.134181   0.044419   -3.021  0.00255 ** 
## poly(age, 5)4 -0.124521   0.044173   -2.819  0.00487 ** 
## poly(age, 5)5  0.095655   0.042151    2.269  0.02335 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 29.90352)
## 
##     Null deviance: 699301  on 1987  degrees of freedom
## Residual deviance:  59269  on 1982  degrees of freedom
## AIC: 12405
## 
## Number of Fisher Scoring iterations: 4

Ta dùng 2 hàm glance và tidy của broom cho mô hình trên, như bạn thấy, ta có 2 dataframe rất gọn đẹp:

glance(lmod)%>%knitr::kable()

null.deviance	df.null	logLik	AIC	BIC	deviance	df.residual
699301.1	1987	-6195.434	12404.87	12444.03	59268.78	1982

tidy(lmod)%>%knitr::kable()

term	estimate	std.error	statistic	p.value
(Intercept)	4.8153583	0.0010151	4743.948410	0.0000000
poly(age, 5)1	6.3266725	0.0444972	142.181380	0.0000000
poly(age, 5)2	-1.2508362	0.0441974	-28.301153	0.0000000
poly(age, 5)3	-0.1341815	0.0444194	-3.020788	0.0025534
poly(age, 5)4	-0.1245206	0.0441734	-2.818906	0.0048666
poly(age, 5)5	0.0956553	0.0421512	2.269337	0.0233549

Khi có dataframe, bạn sẽ nghĩ ngay đến ggplot2, đúng như vậy, ta có thể vẽ rất nhiều biểu đồ từ dataframe mà broom cung cấp:

lmod%>%tidy()%>%
  ggplot(aes(x=term,
             y=estimate,
             ymin=estimate-1.96*std.error,
             ymax=estimate+1.96*std.error,
             col=term))+
  geom_hline(yintercept = 0,linetype=2,col="blue")+
  geom_pointrange(show.legend = F,size=1)+
  coord_flip()+
  theme_bw()

lmod%>%tidy()%>%
  ggplot(aes(x=term,
             y=p.value,
             col=term))+
  geom_hline(yintercept = 0.005,linetype=2,col="blue")+
  geom_point(show.legend = F,size=5)+
  coord_flip()+
  theme_bw()

Hàm augment cho phép tạo ra 1 dataframe chồng lên dataframe gốc, chứa các trị số như fitted, se.fit, residual error, cooksd…

diagDF<-augment(lmod)

diagDF%>%head()%>%knitr::kable()

## Warning in `[<-.data.frame`(`*tmp*`, , isn, value = structure(list(height =
## structure(c("105.5", : provided 13 variables to replace 9 variables

height	poly.age..5.	.fitted	.se.fit	.resid	.hat	.sigma	.cooksd	.std.resid
105.5	-0.0230743292	0.0110197426	0.0063424413	-0.0196888610	0.0209255870	4.659193	0.0022789	-0.0508689
106.5	-0.0200972669	0.0053346290	0.0119236175	-0.0207242879	0.0146627775	4.683920	0.0021049	-1.6933626
108.0	-0.0159293797	-0.0016066597	0.0167457510	-0.0179317966	0.0036937488	4.716927	0.0018470	-3.8241205
111.5	-0.0129523174	-0.0058376726	0.0182884993	-0.0138407745	-0.0039692970	4.739605	0.0016984	-2.8890149
113.5	-0.0087844302	-0.0107432200	0.0181540203	-0.0064736150	-0.0125556528	4.770389	0.0015880	-4.4651563
115.0	-0.0052119555	-0.0140028095	0.0162320005	0.0003761595	-0.0169252462	4.796055	0.0015701	-6.0320482

diagDF%>%ggplot(aes(x=.fitted,y=.resid,col=.resid))+
  geom_point(alpha=0.2,show.legend = F)+
  geom_smooth(alpha=0.2,show.legend = F,fill="blue",col="blue4")+
  theme_bw()+
  scale_color_gradient2(low="red",high="red",mid="blue",midpoint = 0)

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

Trong một số trường hợp, bạn muốn dựng mô hình riêng cho một số phân nhóm trong dữ liệu, ta có thể làm điều này dễ dàng khi kết hợp dplyr và broom: Thí dụ ta dựng cùng mô hình polynomial nói trên riêng cho 3 phân nhóm điều trị 1,2,3:

Kết quả vẫn là 1 dataframe:

models <- Leukemia %>%
  group_by(treatment) %>%
  do(mod = glm(data=.,
               formula=height~poly(age,5),
               family=gaussian(link="log")))

model_sum <- models %>% tidy(mod)

model_sum%>%knitr::kable()

treatment	term	estimate	std.error	statistic	p.value
1	(Intercept)	4.8229908	0.0011808	4084.6389927	0.0000000
1	poly(age, 5)1	5.1997563	0.0428299	121.4046879	0.0000000
1	poly(age, 5)2	-1.0255819	0.0424972	-24.1329397	0.0000000
1	poly(age, 5)3	-0.0275165	0.0427571	-0.6435539	0.5199745
1	poly(age, 5)4	-0.1476890	0.0424684	-3.4776209	0.0005220
1	poly(age, 5)5	0.0447927	0.0405215	1.1054077	0.2691808
2	(Intercept)	4.8055117	0.0025728	1867.8239016	0.0000000
2	poly(age, 5)1	2.6602309	0.0454282	58.5590702	0.0000000
2	poly(age, 5)2	-0.5174254	0.0451399	-11.4627170	0.0000000
2	poly(age, 5)3	-0.1365923	0.0457617	-2.9848593	0.0030564
2	poly(age, 5)4	0.0289002	0.0457650	0.6314922	0.5281710
2	poly(age, 5)5	0.1094441	0.0431043	2.5390551	0.0115899
3	(Intercept)	4.7922652	0.0024658	1943.4838759	0.0000000
3	poly(age, 5)1	2.4228498	0.0424623	57.0588894	0.0000000
3	poly(age, 5)2	-0.4721817	0.0420656	-11.2248949	0.0000000
3	poly(age, 5)3	-0.1269940	0.0419478	-3.0274272	0.0026755
3	poly(age, 5)4	-0.0435665	0.0419031	-1.0396963	0.2992991
3	poly(age, 5)5	0.0731085	0.0406323	1.7992699	0.0729579

Ta vẽ được 3 biểu đồ riêng biệt từ dataframe này:

model_sum%>%
  ggplot(aes(x=term,
             y=estimate,
             ymin=estimate-1.96*std.error,
             ymax=estimate+1.96*std.error,
             col=term))+
  geom_hline(yintercept = 0,linetype=2,col="blue")+
  geom_pointrange(show.legend = F)+
  coord_flip()+
  theme_bw()+facet_wrap(~treatment,ncol=1,scales="free")

4 Mô hình gamlss nhiều tham số

Ta trở lại package gamlss, lần này Nhi dựng một mô hình phức tạp hơn cho một dataset đơn giản. Bài toán ở đây là ước tính chiều dài catherter thông tim dựa vào chiều cao của bệnh nhân.

Mục tiêu của minh họa này là để cho thấy broom làm việc hiệu quả ra sao cho một mô hình gamlss với phân phối student t, với 3 mô hình con cho 3 tham số Mu(link=identity), sigma (link=log) và Nu (link=log).

data(heart, package="robustbase")

heart%>%ggplot(aes(x=height,y=clength))+
  geom_point(size=5,col="blue",alpha=0.3)+
  geom_smooth(method="lm",fill="blue",alpha=0.2)+
  theme_bw()

Trước tiên, ta dựng mô hình và dùng broom để trích xuất kết quả:

gmod=gamlss(data=heart,family=TF(),
            formula=clength ~ height,
            trace=FALSE,
            parallel="multicore",
            ncpus = nC            
            )


summary(gmod)

## ******************************************************************
## Family:  c("TF", "t Family") 
## 
## Call:  gamlss(formula = clength ~ height, family = TF(), data = heart,  
##     trace = FALSE, parallel = "multicore", ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.47898    3.73706   3.072   0.0118 *  
## height       0.61171    0.08909   6.866 4.37e-05 ***
## ---
## 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)   1.2605     0.2041   6.175 6.96e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  log 
## Nu Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 2.895e+01  2.887e-06 10028569   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  12 
## Degrees of Freedom for the fit:  4
##       Residual Deg. of Freedom:  8 
##                       at cycle:  8 
##  
## Global Deviance:     64.30629 
##             AIC:     72.30629 
##             SBC:     74.24591 
## ******************************************************************

tidy(gmod)%>%knitr::kable()

parameter	term	estimate	std.error	statistic	p.value
mu	(Intercept)	11.4789756	3.7370645	3.071656e+00	0.0118074
mu	height	0.6117124	0.0890944	6.865889e+00	0.0000437
sigma	(Intercept)	1.2604900	0.2041241	6.175115e+00	0.0000696
nu	(Intercept)	28.9499863	0.0000029	1.002857e+07	0.0000000

Tiếp theo, Nhi sẽ dùng broom để bootstrap 100 phiên bản mô hình trên đây:

bootgmod <- heart %>%
  bootstrap(100) %>%
  do(tidy(gamlss(family=TF(),
                 formula=.$clength~.$height,
                 trace=FALSE,
                 parallel="multicore",
                 ncpus = nC)
  ))

bootgmod%>%head()%>%knitr::kable()

replicate	parameter	term	estimate	std.error	statistic	p.value
1	mu	(Intercept)	18.2500034	0.0000108	1.693156e+06	0.0000000
1	mu	.$height	0.4999926	0.0000002	2.270082e+06	0.0000000
1	sigma	(Intercept)	-7.3949481	1.2906166	-5.729779e+00	0.0001322
1	nu	(Intercept)	-1.9103709	0.3075115	-6.212356e+00	0.0000660
2	mu	(Intercept)	8.0263124	4.1159903	1.950032e+00	0.0797441
2	mu	.$height	0.7048566	0.1012699	6.960180e+00	0.0000390

Có một cách khác để dùng hàm bootstrap của broom, cách này phức tạp hơn 1 chút vì phải thông qua 2 hàm nest và unnest của package tidy,

bootgmod2 <- heart %>%
  modelr::bootstrap(100) %>%
  mutate(mod = map(strap, ~ gamlss(data=.,family=TF(),
                                   formula=clength ~ height,
                                   trace=FALSE,
                                   parallel="multicore",
                                   ncpus = nC)))%>%
  unnest(map(mod, tidy))

bootgmod2%>%head()%>%knitr::kable()

.id	parameter	term	estimate	std.error	statistic	p.value
001	mu	(Intercept)	9.4097442	4.8973504	1.9213949	0.0836159
001	mu	height	0.6398555	0.1075773	5.9478701	0.0001416
001	sigma	(Intercept)	1.3662393	0.2041264	6.6931055	0.0000341
001	nu	(Intercept)	11.9264567	66.1976703	0.1801643	0.8602998
002	mu	(Intercept)	19.9369256	2.2153694	8.9993684	0.0000041
002	mu	height	0.4697114	0.0533093	8.8110663	0.0000050

Cấu trúc 2 dataframe là tương đương, chỉ khác ở biến .id cho cách làm thứ 2 thay vì biến replicate với hàm do, nhưng đổi lại tên của các predictor của cách 2 rõ ràng hơn.

Tiếp theo, giả sử Nhi đang dùng cách thứ 2, và muốn vẽ biểu đồ 100 mô hình tiên lượng cho Mu. Để làm việc này, trước hết ta đưa dataframe bootgmod2 vào hàm filter để lọc riêng giá trị parameter=mu, sau đó dùng hàm spread để tách Intercept và Height ra thành 2 biến (cột). Sau đó ta dùng 1 vòng lặp để tính Mu từ mô hình tuyến tính với Intercept và Height. Kết quả vòng lặp là 1 dataframe predboot.

Cuối cùng ta vẽ được biểu đồ như sau:

muDF<-bootgmod2%>%filter(parameter=="mu")%>%
  dplyr::select(.id:estimate)%>%
  spread(key=term,value=estimate)

predboot<-heart%>%mutate(pred=NA,Iter=NA)

for(i in 1:nrow(muDF)){
  intercept=muDF$`(Intercept)`[i]
  beta=muDF$height[i]
  tempdf<-heart%>%mutate(pred=heart$height*beta+intercept,
                         Iter=i)
  predboot=rbind(predboot,tempdf)%>%na.omit()
}

predboot%>%
  ggplot()+
  geom_point(aes(x=height,y=clength),alpha=0.1,size=4,col="blue4")+
  geom_point(aes(x=height,y=pred,group=Iter,col=Iter),alpha=0.02,size=3,show.legend = F)+
  geom_line(aes(x=height,y=pred,group=Iter,col=Iter),alpha=0.1,show.legend = F)+
  geom_smooth(method="lm",se=F,aes(x=height,y=pred),col="blue4")+
  scale_color_distiller(palette="PuRd")+
  theme_bw()

5 Mô hình Cox-PH

Trong thí dụ cuối cùng này, Nhi sẽ dùng broom cho một mô hình Cox-Ph trên dataset lung.

Hàm tidy cho ra một dataframe coxdf với cấu trúc hoàn toàn khác với mô hình glm. Từ dataframe này ta có thể vẽ được biểu đồ cho mô hình Cox:

library(survival)

data(lung)

coxmod <- survfit(coxph(Surv(time, status) ~ age + sex, data=lung))

glance(coxmod)%>%knitr::kable()

records	n.max	n.start	events	rmean	rmean.std.error	median	conf.low	conf.high
228	228	228	165	380.9381	20.27475	320	285	363

coxdf=tidy(coxmod)

coxdf%>%head()%>%knitr::kable()

time	n.risk	n.event	n.censor	estimate	std.error	conf.high	conf.low
5	228	1	0	0.9958207	0.0041893	1.0000000	0.9876776
11	227	3	0	0.9832688	0.0084466	0.9996823	0.9671247
12	224	1	0	0.9790670	0.0094750	0.9974188	0.9610529
13	223	2	0	0.9706383	0.0112867	0.9923495	0.9494021
15	221	1	0	0.9664127	0.0121065	0.9896182	0.9437513
26	220	1	0	0.9621802	0.0128837	0.9867862	0.9381877

coxmod%>%tidy()%>%
  ggplot(aes(time, estimate)) +
  geom_line(col="red4") +
  geom_ribbon(aes(ymin = conf.low, 
                  ymax = conf.high),
              fill="red",
              alpha = 0.4)+
  theme_bw()

Ta có thể bootstrap mô hình Cox này dễ dàng bằng broom:

bootcox<- lung %>%
  bootstrap(200) %>%
  do(tidy(survfit(coxph(Surv(time, status) ~ age + sex, .))))

bootcox%>%head%>%knitr::kable()

replicate	time	n.risk	n.event	n.censor	estimate	std.error	conf.high	conf.low
1	5	228	1	0	0.9957728	0.0042371	1.0000000	0.9875375
1	11	227	4	0	0.9788329	0.0095779	0.9973815	0.9606293
1	12	223	1	0	0.9745817	0.0105243	0.9948934	0.9546847
1	13	222	1	0	0.9703219	0.0114035	0.9922532	0.9488754
1	30	221	1	0	0.9660584	0.0122284	0.9894918	0.9431799
1	31	220	1	0	0.9617874	0.0130107	0.9866289	0.9375713

Đây là biểu đồ của 200 mô hình Cox từ bootstrap:

bootcox%>%
  ggplot(aes(time, estimate,group=replicate,col=replicate)) +
  geom_line(alpha=0.1,show.legend = F) +
  scale_color_distiller(palette="RdBu")+
  theme_bw()

Hoặc giản dị hơn: Chỉ trình bày trung vị và 97.5% CI của giá trị dự báo từ 200 mô hình bootstrap;

bootcoxsum <- bootcox %>%
  group_by(time) %>%
  summarize(median = median(estimate),
            LL= quantile(estimate, .025),
            UL = quantile(estimate, .975))

bootcoxsum%>%ggplot()+
  geom_ribbon(aes(x=time,
                  ymin = LL, 
                  ymax = UL),
              fill="red",
              alpha =0.1)+
  geom_line(aes(time, median),col="red4",size=1)+
  theme_bw()

6 Tổng kết

Bài thực hành đến đây là chấm dứt. Tuy Nhi chỉ mới khảo sát package broom với 3 loại mô hình: GLM, Gamlss và Survival; tuy nhiên như đã giới thiệu, broom còn tương thích với rất nhiều package hồi quy khác, bao gồm nlme, lme4, gam, mgcv, rstanarm và brms… và với mỗi package nó đều làm được chức năng đóng gói mô hình bằng hàm tidy và tái chọn mẫu bằng hàm bootstrap.

Như ta thấy, khi tải broom đồng thời với tidyverse chúng ta đã kết nối những công dụng tiện lợi của dplyr, broom và ggplot2 trong một quy trình khép kín và rất tiện lợi.Nhi khuyến khích các bạn sử dụng quy trình này khi phân tích dữ liệu bằng mô hình hồi quy, bất kể bạn đang dùng algorithm nào.

Tạm biệt các bạn và hẹn gặp lại.