Trích xuất feature dữ liệu chuỗi

1 Giới thiệu

Chuỗi thời gian (time series), là một cấu trúc dữ liệu gồm nhiều giá trị liên tiếp theo trình tự thời gian. Một cách phổ quát hơn, Nhi thích sử dụng thuật ngữ “chuỗi dữ liệu” (sequence, sequential data) cho bất cứ dữ liệu nào có cấu trúc liên tục và theo trình tự, không bắt buộc phải là những con số và phải có yếu tố thời gian biểu kiến (thí dụ chuỗi sự kiện/trạng thái, văn bản, gen, protein, … ).

Trong các bài toán Machine learning, dữ liệu chuỗi có thể được hoán chuyển rất đa dạng ở đầu vào, theo 2 trường phái khác nhau: hoặc trực tiếp sử dụng dữ liệu thô (arrays/tensor từ 1D cho đến 5D) - cách làm này thường gặp trong deep learning (CNN, LSTM, …) hoặc gián tiếp bằng cách trích xuất các features từ dữ liệu gốc rồi dùng các features này như input data cho những algorithm như SVM, XGBoost hay Random Forest. Bài thực hành hôm nay sẽ bàn về trường phái thứ 2, trích xuất features, với nội dung cụ thể là giới thiệu về R package tsfeatures. Công cụ này cho phép thực hiện rất nhiều quy trình trích xuất features - vốn trước kia chỉ có thể làm bằng Matlab.

Bạn có thể cài đặt bản chính thức từ CRAN như sau:

# Cài đặt từ CRAN
install.packages('tsfeatures', dependencies = TRUE)

Bạn cũng có thể cài bản thử nghiệm từ Github như sau:

# Cần package "devtools"
devtools::install_github("robjhyndman/tsfeatures")

2 Một vài thí dụ về time series

Sau đây là một số dataset kinh điển giúp các bạn hình dung về time series,

library(tsfeatures)
library(tidyverse)

1. Data sunpost : timeseries có yếu tố thời gian biểu kiến (mốc thời gian là năm, tháng, ngày, giờ xác định)

annual_sunpost = read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv",
                          row.names = 1)%>%
  ts(data=.$value,
     start=min(.$time),
     end=max(.$time),
     frequency=1)

head(annual_sunpost,20)

## Time Series:
## Start = 1700 
## End = 1719 
## Frequency = 1 
##  [1]  5 11 16 23 36 58 29 20 10  8  3  0  0  2 11 27 47 63 60 39

plot(annual_sunpost,type="l",col="red")

2. Dữ liệu AirPassenger: là một thí dụ về time series có yếu tố chu kỳ

airpass = read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/datasets/AirPassengers.csv",
                   row.names = 1)%>%
  ts(data=.$value,
      class="ts",
      start=1949,end=1960,
      frequency=12)

airpass

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417

plot(airpass,type="l",col="blue")

3. Một đoạn tín hiệu điện tim:

Đây là thí dụ về timeseries không có mốc thời gian biểu kiến (có thể chỉ đơn giản là 1 vector số)

ecg = read.csv("ecg0606_1.csv",header = F)%>%as.ts(.$V1)

plot(ecg,type="l",col="black")

4. Dữ liệu chuỗi đa chiều:

8 băng Tín hiệu cảm biến gyro và gia tốc ghi nhận tư thế, chuyển động của cơ thể trong 5 trạng thái giấc ngủ khác nhau: có thể xem như một tensor 3D.

sleep_df = read.csv("origin_comp.csv")

Awk_df = sleep_df%>%filter(Stade == "Awk")%>%.[c(1:100),]
N1_df = sleep_df%>%filter(Stade == "N1")%>%.[c(1:100),]
N2_df = sleep_df%>%filter(Stade == "N2")%>%.[c(1:100),]
N3_df = sleep_df%>%filter(Stade == "N3")%>%.[c(1:100),]
Rem_df = sleep_df%>%filter(Stade == "REM")%>%.[c(1:100),]

df=rbind(Awk_df,N1_df,N2_df,N3_df,Rem_df)%>%select(-idx)

df$idx = rep(c(1:100),5)

df%>%gather(1:8,key="band",value="amplitude")%>%
  ggplot(aes(x=idx,y=amplitude,col=band))+
  geom_line(alpha=0.6,size=1)+
  theme_bw(5)+
  facet_wrap(~Stade, ncol=1, scales="free")

df%>%gather(1:8,key="band",value="amplitude")%>%
  ggplot(aes(x=idx,y=amplitude,col=band))+
  geom_line(alpha=0.6,size=0.8)+
  theme_bw(8)+
  facet_grid(band~Stade, scales="free_y")

Hình vẽ trên chỉ mới trình bày 1 phân đoạn ngắn (10 giây tín hiệu với tần số lấy mẫu 10 Hz),khi một chuỗi dài được chia đều thành nhiều phân đoạn, nó sẽ trở thành một tensor/array 4D, thậm chí 5D, bao gồm chiều thời gian. Giả định nếu muốn phân biệt 5 trạng thái giấc ngủ khác nhau (bài toán multiclass classification), ta khó lòng sử dụng trực tiếp dữ liệu gốc, ngay cả bằng mô tả trực quan.Đây là một không gian dữ liệu quá lớn, nhiều thông tin tiềm ẩn bên trong và khó phân tích. Kỹ thuật trích xuất features có thể rút gọn chiều dữ liệu xuống bằng cách sử dụng các chỉ số cho phép mô tả về đặc tính, khuynh hướng của chuỗi dữ liệu gốc, package tsfeatures cung cấp một số hàm tiện lợi để làm việc này.

3 Trích xuất feature cho 1 chuỗi bằng hàm chuyên dụng

về mặt kỹ thuật, mỗi chuỗi dữ liệu cơ bản là một vector (1D array) nên ta có thể áp dụng tất cả những hàm thông thường, bao gồm thống kê, tính toán số học, hoán chuyển, ..Nếu hàm này xuất ra kết quả là một con số, kết quả này có thể được xem như một feature,

Thí dụ đơn giản nhất, hàm median cho biết vị trí trung tâm của chuỗi:

median(annual_sunpost)

## [1] 39

ggplot(data=annual_sunpost)+
  geom_point(aes(y=annual_sunpost,x=c(1700:1988)),col="red",alpha=0.5)+
  geom_line(aes(y=annual_sunpost,x=c(1700:1988)),col="grey",alpha=0.5)+
  geom_hline(yintercept = median(annual_sunpost), color = "blue", linetype=2)+
  scale_x_continuous('Year')+
  theme_bw()

Package tsfeature cung cấp khoảng 30 hàm khác nhau tương ứng với nhiều thuộc tính chuyên biệt cho time series,

https://cran.r-project.org/web/packages/tsfeatures/vignettes/tsfeatures.html#stl_features

Thí dụ, hàm acf_features cung cấp hệ số tự tương quan và tổng bình phương của 10 autocorrelation coefficients đầu tiên

acf_features(ecg)%>%knitr::kable()

	x
x_acf1	0.9690683
x_acf10	2.9172409
diff1_acf1	0.9187728
diff1_acf10	2.5582842
diff2_acf1	0.7313610
diff2_acf10	1.6636852

Hàm stl_features trích xuất 8 thuộc tính về trend và seasonality của time series dựa vào phân tích STL decomposition.

stl_features(ecg)%>%knitr::kable()

	x
nperiods	0.0000000
seasonal_period	1.0000000
trend	0.0043295
spike	0.0000000
linearity	0.4771984
curvature	-0.3072355
e_acf1	0.9690477
e_acf10	2.9124859

Hàm heterogeneity trích xuất 4 thuộc tính nhằm khảo sát heterogeneity dựa vào tổng bình phương và R2 của mô hình GARCH và ARCH:

heterogeneity (airpass)%>%knitr::kable()

	x
arch_acf	0.2294789
garch_acf	0.2278266
arch_r2	0.1891338
garch_r2	0.1894199

4 Trích xuất features cho hàng loạt chuỗi

Hàm tsfeatures cho phép áp dụng một hay nhiều phương pháp trích xuất features khác nhau cho 1 list gồm nhiều chuỗi khác nhau. Nó làm đơn giản code mà không cần đến vòng lặp, hàm apply hay hàm map.

Trước hết ta đưa các chuỗi cần xử lý vào 1 list (có thể dùng hàm split trên 1 dataframe), thí dụ ta muốn xử lý 3 chuỗi anorm_raw ở 3 trạng thái Awk, N2 và REM :

ts_list = list(Awk_df$anorm_raw,
               N2_df$anorm_raw,
               Rem_df$anorm_raw)

str(ts_list)

## List of 3
##  $ : num [1:100] 1.02 1.02 1.02 1.02 1.02 ...
##  $ : num [1:100] 1.01 1.01 1.01 1.01 1.01 ...
##  $ : num [1:100] 1.01 1.02 1.02 1.01 1.01 ...

Có thể áp dụng 1 hay hàng loạt phương pháp, tên của phương pháp là tên của hàm có trong Environment, có thể là 1 hàm có sẵn trong package hay 1 hàm đặc biệt do bạn tạo ra:

Thí dụ quy trình sau trích xuất các features: entropy, stability, nonlinearity, lumpiness, crossing_points cho list 3 series:

tsfeatures(ts_list, features = c("entropy","stability","nonlinearity", 
                                 "lumpiness", "crossing_points"))%>%knitr::kable()

entropy	stability	nonlinearity	lumpiness	crossing_points
0.9866066	0.3323816	1.1801256	0.4820703	53
0.9516680	0.0166262	0.3401689	0.3215552	62
0.8107511	0.0340071	0.1445946	0.2417979	78

Kết quả xuất ra là 1 data frame, mỗi hàng tương ứng 1 series trong list, mỗi cột là 1 features. Nếu hàm xuất ra nhiều features, chúng sẽ được ghép nối tiếp theo nhiều cột

Như đã nói, ta có thể viết hàm tùy thích, thí dụ Nhi tạo ra 1 hàm có tên là boxcox_stl, nó sẽ làm 1 quy trình phức tạp gồm: chuẩn hóa (scaling) chuỗi thành thang đo 0:1, sau đó hoán chuyển BoxCox, cuối cùng áp dụng phương pháp STL để xuất ra 8 features:

boxcox_stl <- function(x,...) {
  minmax = (x- min(x)) /(max(x)-min(x))
  lambda <- forecast::BoxCox.lambda(minmax)
  y <- forecast::BoxCox(minmax, lambda)
  c(stl_features(y, s.window='periodic', robust=TRUE), lambda=lambda)
}

tsfeatures(ts_list, features = 'boxcox_stl')%>%knitr::kable()

nperiods	seasonal_period	trend	spike	linearity	curvature	e_acf1	e_acf10	lambda
0	1	0.2707876	1e-07	0.4579602	0.3788559	-0.0796860	0.0908502	1.077253
0	1	0.0000000	2e-07	0.0668875	-0.0161874	-0.3999478	0.3523617	1.018802
0	1	0.0122354	1e-07	0.0182492	0.1391286	-0.6922963	1.7365545	1.421700

5 Trích xuất features cho dữ liệu chuỗi đa chiều

Ta sẽ làm 1 quy trình phức tạp hơn, thí dụ ta có dataframe sleep_df gồm 8 băng tín hiệu, ghi trong 5 trạng thái giấc ngủ khác nhau, mỗi trạng thái lại được chia thành 10 phân đoạn khác nhau có độ dài 100 mẫu (10 giây, 10Hz), cấu trúc này tương ứng với 1 tensor/array 4D. Ta muốn trích xuất entropy của từng phân đoạn

sleep_df$idx = rep(c(1:50),each = 100)

sleep_df%>%sample_n(10)%>%knitr::kable()

anorm_raw	ax_raw	ay_raw	az_raw	gnorm_raw	gx_raw	gy_raw	gz_raw	idx	Stade
1.009713	-0.7873840	0.6317135	0.0219420	0.8806421	-0.7765406	0.2447864	0.3354096	28	N2
1.012446	-0.7616880	0.6274110	-0.2263485	0.3145053	-0.1477296	0.2104879	0.1810646	45	REM
1.014557	-0.7628377	0.6295167	-0.2260740	0.1365469	-0.0245079	0.0478859	0.0807093	43	REM
1.009033	-0.7902835	0.6268615	0.0253600	1.4076002	-0.8641931	0.4715389	0.9966136	24	N2
1.022717	0.0953163	0.6772870	-0.7603557	1.0109192	-0.4437156	-0.3700511	-0.3308760	10	Awk
1.014243	-0.7630310	0.6286010	-0.2265625	0.2718565	-0.1934616	-0.1610826	0.0190981	46	REM
1.013054	-0.7928260	0.6295577	0.0367433	0.8897557	0.7910407	-0.2252344	-0.0259980	26	N2
1.011761	-0.7847290	0.6381230	0.0256960	1.3894277	-0.9327906	0.6944809	0.7603336	25	N2
1.012252	-0.7820437	0.6424767	0.0168457	1.0432809	0.5661931	0.5255283	-0.6980000	34	N3
1.011412	-0.7983400	0.6207885	0.0150755	0.6220690	0.5096639	-0.0982016	-0.3315114	13	N1

sleep_df%>%group_by(idx,Stade)%>%tally()%>%knitr::kable()

idx	Stade	n
1	Awk	100
2	Awk	100
3	Awk	100
4	Awk	100
5	Awk	100
6	Awk	100
7	Awk	100
8	Awk	100
9	Awk	100
10	Awk	100
11	N1	100
12	N1	100
13	N1	100
14	N1	100
15	N1	100
16	N1	100
17	N1	100
18	N1	100
19	N1	100
20	N1	100
21	N2	100
22	N2	100
23	N2	100
24	N2	100
25	N2	100
26	N2	100
27	N2	100
28	N2	100
29	N2	100
30	N2	100
31	N3	100
32	N3	100
33	N3	100
34	N3	100
35	N3	100
36	N3	100
37	N3	100
38	N3	100
39	N3	100
40	N3	100
41	REM	100
42	REM	100
43	REM	100
44	REM	100
45	REM	100
46	REM	100
47	REM	100
48	REM	100
49	REM	100
50	REM	100

Để làm việc này, dự kiến ta phải dùng hàm map và 1 hàm cải biên để làm việc trên list.

Ta tạo ra hàm ent_func, nó sẽ chuyển dataframe nhỏ gồm 8 băng dữ liệu thành 1 list, sau đó áp dụng hàm ts_features (hay trực tiếp hàm entropy cũng được) để trích xuất kết quả entropy cho mỗi băng dữ liệu.

ent_func = function(.x){
  ts_list = as.list(.x[,-c(9,10)])
  out = tsfeatures(ts_list,
             features="entropy")%>%as.vector()
  return(out)
 }

Tiếp theo, ta dùng hàm map_df để chạy hàm ent_func cho từng phân đoạn từ 1 đến 50 (có 50 phân đoạn, mỗi nhóm label có 10 phân đoạn và mỗi phân đoạn dài 100 mẫu) và ghi kết quả trong dataframe out_df

sleep_df%>%split(.$idx)%>%
  map_df(.,~data_frame(ent=ent_func(.x)%>%.$entropy,
                       para = colnames(df[,-c(9,10)]),
                       stade=.x$Stade[[1]]))->out_df

out_df$idx = factor(c(1:50))

out_df%>%sample_n(10)%>%knitr::kable()

ent	para	stade	idx
0.9112369	gx_raw	N1	42
0.8880870	gy_raw	N2	49
0.9242151	gx_raw	N2	6
0.8515206	gz_raw	N2	18
0.9102416	az_raw	REM	48
0.6979388	ax_raw	N2	18
0.9169486	ax_raw	N1	22
0.8572580	gx_raw	N1	10
0.8757794	ax_raw	REM	44
0.8719307	gy_raw	N2	15

Lúc này ta có thể dùng entropy như 1 feature để phân biệt 5 label sleep stades:

out_df%>%ggplot()+
  geom_density(aes(x=ent,fill=stade),alpha=0.7)+
  theme_bw(8)+
  facet_wrap(~para,ncol=3,scales = "free")

out_df%>%ggplot()+
  geom_boxplot(aes(x=stade,y=ent,fill=stade),alpha=0.7)+
  theme_bw(8)+
  coord_flip()+
  facet_wrap(~para,ncol=2,scales = "free")

Thay vì xuất 1 features duy nhất là entropy, ta có thể xuất nhiều features hơn, như sau đây - Nhi dùng quy trình stl_features để xuất ra 8 features, sau đó chọn ra 4 là trend,spike,linearity và curvature. Lưu ý, hàm stl_func xuất kết quả là 1 dataframe thay vì vector như hàm ent_func

stl_func = function(.x){
  ts_list = as.list(.x[,-c(9,10)])
  out = tsfeatures(ts_list,
             features="stl_features")
  return(out)
 }

Kết hợp hàm split và map_df, ta có thể chạy áp dụng hàm stl_func vừa tạo cho hàng loạt phân đoạn, và lưu kết quả lại thành dataframe out_df2

sleep_df%>%split(.$idx)%>%
  map_df(.,~stl_func(.x)%>%
           select(trend,spike,linearity,curvature))->out_df2

out_df2%>%sample_n(10)%>%knitr::kable()

trend	spike	linearity	curvature
0.2821990	0.0000906	-1.0794391	1.8113425
0.0584019	0.0001293	1.1788674	-0.1435161
0.4926281	0.0001104	-0.3078516	-1.1142797
0.3645284	0.0000748	2.3945986	-1.9783793
0.8308025	0.0000071	-0.7858780	0.7821673
0.6655960	0.0000375	5.5228314	-3.0745624
0.7475569	0.0000328	1.9326001	3.3317111
0.0084289	0.0002730	-0.0255262	-0.6787094
0.3947238	0.0008709	0.0167956	0.8797239
0.4312240	0.0000725	0.4181652	0.5867820

Ta ghép 2 dataframe out_df và out_df2 lại thành 1, lúc này ta có 5 features trong 1 dataframe

stl_out_df = bind_cols(out_df2,out_df)

dataframe này có thể được dùng như input data cho bài toán classfication, hoặc unsupervised learning.

Ta có thể thăm dò trực quan đặc tính của dữ liệu features bằng heatmap:

library(viridis)

stl_out_df%>%
  dplyr::select(-c(para,stade,idx))%>%
  scale(.)%>%
  as_tibble()%>%
  mutate(Stade=out_df$stade,
         Id=factor(rep(c(1:5),80)),
         para=out_df$para)%>%
  gather(trend:ent,
         key="Feature",value="Value")%>%
  ggplot(aes(x=reorder(Id,Value),
             y=reorder(Feature,Value),
             fill=Value))+
  geom_tile(show.legend=T)+
  theme_bw(7)+
  theme(axis.text.y=element_blank(),
        axis.text.x = element_text(angle =45,hjust=1,vjust=1))+
  coord_flip()+
  scale_y_discrete("Features")+
  scale_x_discrete("Fragments",breaks=NULL)+
  facet_grid(Stade~para,shrink = T)+
  scale_fill_viridis(option="A",begin=1,end=0)

Hoặc bằng violin plot

stl_out_df%>%
  dplyr::select(-c(para,stade,idx))%>%
  scale(.)%>%
  as_tibble()%>%
  mutate(stade=out_df$stade,
         para=out_df$para)%>%
  gather(1:5,key="Feature",value="Value")%>%
  ggplot()+
  geom_violin(aes(x=stade,y=Value,fill=stade),alpha=0.5)+
  facet_grid(Feature~para)+
  theme_bw(7)+
  theme(axis.text.y=element_blank(),
        axis.text.x = element_text(angle =45,hjust=1,vjust=1))

6 Tổng kết

Bài thực hành đến đây là hết, Nhi hy vọng package tsfeature sẽ mang lại nhiều tiên ích cho các bạn nghiên cứu sinh khi làm việc với dữ liệu chuỗi. Chúc các bạn thực hành vui và hẹn gặp lại.