Hàm map của package purrr

1 Hàm map của package purrr

Vì R là một ngôn ngữ lập trình, sớm muộn thì người dùng cũng sẽ thay đổi cách viết code theo hướng “programing”, bao gồm function, loop và algorithm. Điều thú vị đó là trong 2 năm gần đây bản thân ngữ pháp R cũng có nhiều sự thay đổi, một trong số đó là package purrr của Lionel Henry và Hadley Wickham đã đưa ra thêm nhiều giải pháp tối ưu hóa hiệu quả của hàm và vòng lặp. Bài thực hành ngắn sau đây sẽ so sánh tốc độ của hàm map( ) của purrr với cách làm cổ điển với hàm apply( ) và vòng lặp “for loop”.

Package purr sau đó được đưa vào tidyverse để liên hợp với những công cụ khác như dplyr, broom, ggplot2… để tạo ra quy trình khép kín khi thao tác trên dữ liệu.

Hàm map( ) trong purrr cho phép áp dụng một quy trình lặp lại cho từng cột trong một cấu trúc dữ liệu.

Thí dụ với dữ liệu hypothyroid của viện Garvan (1988) như sau

library(tidyverse)

df=read.csv("https://www.openml.org/data/get_csv/53534/hypothyroid.arff",na.strings = "?")

df=df%>%dplyr::select(referral.source,
                      TSH,T3,TT4,T4U,FTI,binaryClass)%>%na.omit()

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

	referral.source	TSH	T3	TT4	T4U	FTI	binaryClass
1	SVHC	1.30	2.5	125	1.14	109	P
5	SVI	0.72	1.2	61	0.87	70	P
8	SVI	2.20	0.6	80	0.70	115	P
9	SVI	0.60	2.2	123	0.93	132	P
10	SVI	2.40	1.6	83	0.89	93	P
11	SVI	1.10	2.2	115	0.95	121	P

Ta có thể dùng hàm map để mô tả phân vị của 5 biến :output của hàm map là 1 list.

df%>%select(2:6)%>%
  map(~quantile(.,c(0.05,0.25,0.5,0.75,0.975),na.rm=T))

## $TSH
##     5%    25%    50%    75%  97.5% 
##  0.025  0.500  1.300  2.600 31.225 
## 
## $T3
##    5%   25%   50%   75% 97.5% 
##   0.8   1.5   2.0   2.3   3.9 
## 
## $TT4
##      5%     25%     50%     75%   97.5% 
##  61.000  88.000 103.000 124.000 194.225 
## 
## $T4U
##    5%   25%   50%   75% 97.5% 
## 0.740 0.870 0.975 1.090 1.520 
## 
## $FTI
##    5%   25%   50%   75% 97.5% 
##    67    93   107   124   188

Nếu bạn còn nhớ, các hàm summarise của dplyr cũng cho phép làm việc này, tuy nhiên chúng kém hiệu quả hơn vì chỉ cho phép dùng những function có output là 1 con số duy nhất. Hàm map ngược lại, cho phép làm nhiều function phức tạp với output có định dạng khác nhau, thí dụ double, character, list…

Ta sẽ thử làm 1 mô hình hồi quy logistic cho từng phân nhóm dữ liệu khác nhau chia theo bệnh viện, bằng 2 hàm map nối tiếp

df$Class=as.numeric(df$binaryClass)-1

formula=as.formula(Class~TSH+T4U+T3+TT4)

df%>%select(-binaryClass)%>%
  split(.$referral.source)%>%
  map(~glm(formula,
           data=.,
           family = "binomial"))%>%
  map(summary)

## $other
## 
## Call:
## glm(formula = formula, family = "binomial", data = .)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4698   0.1875   0.2660   0.3318   5.9140  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.957366   0.757980   3.902 9.55e-05 ***
## TSH         -0.124030   0.018072  -6.863 6.73e-12 ***
## T4U         -2.293738   0.923966  -2.482 0.013047 *  
## T3           0.110441   0.254549   0.434 0.664384    
## TT4          0.023534   0.006134   3.837 0.000125 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 790.64  on 1296  degrees of freedom
## Residual deviance: 459.70  on 1292  degrees of freedom
## AIC: 469.7
## 
## Number of Fisher Scoring iterations: 7
## 
## 
## $STMW
## 
## Call:
## glm(formula = formula, family = "binomial", data = .)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
##  -8.49    0.00    0.00    0.00    0.00  
## 
## Coefficients:
##               Estimate Std. Error   z value Pr(>|z|)    
## (Intercept) -1.610e+14  3.033e+07  -5308199   <2e-16 ***
## TSH         -1.417e+13  1.456e+05 -97279630   <2e-16 ***
## T4U          1.273e+15  2.545e+07  49999472   <2e-16 ***
## T3          -3.914e+14  9.599e+06 -40776417   <2e-16 ***
## TT4          7.544e+12  1.691e+05  44605834   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance:  39.703  on 99  degrees of freedom
## Residual deviance: 288.349  on 95  degrees of freedom
## AIC: 298.35
## 
## Number of Fisher Scoring iterations: 13
## 
## 
## $SVHC
## 
## Call:
## glm(formula = formula, family = "binomial", data = .)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.98147   0.01502   0.03047   0.05995   2.39122  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  7.20246    2.92375   2.463   0.0138 *  
## TSH         -0.89964    0.20366  -4.417 9.99e-06 ***
## T4U         -7.35965    4.73959  -1.553   0.1205    
## T3          -0.42942    0.67302  -0.638   0.5234    
## TT4          0.09341    0.04804   1.945   0.0518 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 98.831  on 366  degrees of freedom
## Residual deviance: 25.053  on 362  degrees of freedom
## AIC: 35.053
## 
## Number of Fisher Scoring iterations: 9
## 
## 
## $SVHD
## 
## Call:
## glm(formula = formula, family = "binomial", data = .)
## 
## Deviance Residuals: 
##        Min          1Q      Median          3Q         Max  
## -2.636e-05   2.100e-08   2.100e-08   2.100e-08   3.264e-05  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  2.632e+02  3.695e+05   0.001    0.999
## TSH         -2.082e+01  2.721e+04  -0.001    0.999
## T4U         -1.011e+02  3.952e+05   0.000    1.000
## T3          -1.939e+01  1.044e+05   0.000    1.000
## TT4          1.218e-01  3.406e+03   0.000    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2.0652e+01  on 35  degrees of freedom
## Residual deviance: 2.1416e-09  on 31  degrees of freedom
## AIC: 10
## 
## Number of Fisher Scoring iterations: 25
## 
## 
## $SVI
## 
## Call:
## glm(formula = formula, family = "binomial", data = .)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3195   0.0720   0.1052   0.1576   5.0401  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.94260    1.23472   3.193  0.00141 ** 
## TSH         -0.65188    0.07295  -8.936  < 2e-16 ***
## T4U         -1.60158    1.65787  -0.966  0.33402    
## T3           0.42126    0.44834   0.940  0.34743    
## TT4          0.02735    0.01117   2.449  0.01432 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 577.51  on 951  degrees of freedom
## Residual deviance: 163.06  on 947  degrees of freedom
## AIC: 173.06
## 
## Number of Fisher Scoring iterations: 8

2 Hàm map có thực sự tốt hơn vòng lặp và hàm apply ?

Trước kia, những quy trình lặp lại trong R thường được xử lý bằng 1 vòng lặp, hoặc họ hàm apply. Sau đây chúng ta sẽ làm 1 thí nghiệm benchmark để so sánh tốc độ của 3 phương pháp này,

Vấn đề giả định là ta muốn tính trung bình của mỗi biến trong một dataframe lớn với 100000 hàng và số biến tăng dần từ 2 đến 1024 biến.

Đầu tiên ta tạo ra 3 hàm cho 3 phương pháp

# Apply

colmean_apply=function(odf){
  out<-apply(odf,2,mean)%>%as.vector()
  return(out)
}

# For loop

colmean_loop=function(odf){
  out=vector(length=1000)
  for(i in seq_along(odf)){
    out[i]<-mean(odf[[i]])
  }
  return(out)
}

# Map

colmean_map=function(odf){
  out<-map_dbl(odf,mean)%>%as.vector()
  return(out)
}

Và tạo hàm benchmark với nội dung: tạo ra 1 dataframe mô phỏng phân phối uniform với nrow =số hàng và ncol = số cột (biến), sau đó lần lượt áp dụng 3 hàm như trên và xuất kết quả thời gian thi hành

run_benchmark <- function(ncol,nrow) {
  odf<-matrix(ncol=ncol,nrow=nrow)%>%as_tibble()
  set.seed(123)
  for(i in seq_along(odf)){
    odf[,i]<-runif(n=nrow,min=1,max=10)
  }
  
  res <- list(Apply=system.time(colmean_apply(odf)),
              loop=system.time(colmean_loop(odf)),
              map=system.time(colmean_map(odf))
              )
  
  res <- lapply(res, `[[`, "elapsed")
  
  res <- c(ncol= ncol, res)
  return(res)
}

all_times <- lapply(1:10, function(n) {
  run_benchmark(ncol=2^n,nrow=100000)
})

times <- lapply(all_times, as.data.frame)
times <- do.call(rbind, times)

knitr::kable(times)

ncol	Apply	loop	map
2	0.00	0.02	0.00
4	0.00	0.00	0.00
8	0.02	0.00	0.00
16	0.03	0.00	0.01
32	0.14	0.00	0.02
64	0.19	0.00	0.02
128	0.38	0.02	0.03
256	0.59	0.06	0.06
512	1.31	0.17	0.14
1024	553.20	55.39	2.32

library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

colnames(times)<-(c("n_var","Apply_func","For_Loop","purrr_Map"))

times%>%gather(Apply_func:purrr_Map,key="Method",value="Time")->long


long%>%ggplot(aes(x=n_var,y=Time,col=Method,group=Method))+
  geom_line(size=1)+
  geom_point(size=3)+
 scale_y_continuous(minor_breaks = NULL)+
  scale_x_continuous(trans = log2_trans(),
                     breaks = trans_breaks("log2", function(x) 2^x),
                     labels = trans_format("log2", math_format(2^.x)),
                     minor_breaks = NULL)+
  theme_bw()

Kết quả cho thấy bắt đầu từ 200 biến trở lên, hàm map vẫn chạy rất nhanh so với vòng lặp for và hàm apply. Hàm apply bắt đầu đuối từ 30 biến trở lên.