1 Giới thiệu

Ngoài giờ làm chính thức, nghề tay trái của Nhi là tư vấn về phân tích số liệu, nói một cách bình dân là đi làm thông kê « dạo ». Công việc này có thể rất nhẹ nhàng hoặc cực khổ, tùy vào tình huống, đôi khi những đòi hỏi kì quặc và khó khăn của một vài thân chủ lại là động lực để Nhi đi tìm tòi học hỏi kiến thức, kỹ năng mới.

Cách đây 3 ngày có môt thân chủ đưa ra yêu cầu rất thú vị : Ông ta nghĩ ra một phương pháp điều trị mới cho triệu chứng thiếu oxy trong khi ngủ do hẹp đường thở trên, bằng cách gắn 1 số điện cực kích thích một vài nhóm cơ hô hấp, bệnh nhân lập tức cải thiện được lưu lượng khí thở và phục hồi độ bão hòa Oxygen (SpO2). Ông ta muốn Nhi tư vấn nên dùng công cụ thống kê nào để chứng minh là kích thích điện thực sự đã làm cải thiện SpO2. Vấn đề là SpO2 được theo dõi liên tục trong khi ngủ như 1 time series data với tần suất lấy mẫu là 1 Hz. Đây là 1 bài toán chứng minh liên hệ nhân quả rất hiếm gặp trong nghiên cứu y học.

Tìm kiếm trên mạng, Nhi phát hiện ra bài toán tương tự có tồn tại trong những lĩnh vực hoàn toàn khác, đó là Marketing, quảng cáo, Kinh tế tài chính. Một giải pháp rất hiệu quả đã được giới thiệu bởi Brodersen KH (nhân viên của Google) từ năm 2015.

Brodersen KH, Gallusser F, Koehler J, Remy N, Scott SL. Inferring causal impact using Bayesian structural time-series models. Annals of Applied Statistics, 2015, Vol. 9, No. 1, 247-274. http://research.google.com/pubs/pub41854.html

Một cách đơn giản, tác giả áp dụng 1 mô hình cấu trúc Bayes trên một phần của Time series trước khi 1 can thiệp bắt đầu, mô hình này đưa ra dự báo về diễn tiến của Time serie nếu như không có can thiệp nào xảy ra, từ đó so sánh với kết quả trên thực tế sau can thiệp. Kết quả phân tích là 1 phân phối hậu nghiệm về độ lớn của hiệu ứng can thiệp lên diễn tiến của time series này.

Bằng cách này, người ta có thể khảo sát doanh số của 1 sản phẩm SAU khi bắt đầu 1 chiến dịch quảng cáo trên truyền thông?

2 Bài toán mô phỏng

Nhi quyết định thử nghiệm phương pháp này trên 1 bài toán mô phỏng như sau:

Giả sử ta muốn khảo sát hiệu quả điều trị của 1 loại thuốc điều trị hen phế quản trên 1 bệnh nhân. Kết quả là khả năng thông khí được đo bằng lưu lượng đỉnh (Peakflow). Chỉ số Peakflow được đo mỗi ngày vào 1 giờ nhất định trong 2 tháng trước và sau khi điều trị bằng loại thuốc này. Ta có 1 time series 120 điểm, điểm khởi đầu điều trị nằm ở ngày 61.

Số liệu được mô phỏng như sau:

library(tidyverse)
set.seed(123)

peakflow=c(round(rnorm(60,160,25)+arima.sim(model=list(ar=0.3),n=60),0),
           round(rnorm(60,350,40)+arima.sim(model=list(ar=0.3),n=60),0))

df=as.data.frame(peakflow)%>%mutate(.,SBP=round(rnorm(120,100,2)+arima.sim(model=list(ar=0.3),n=120),0),
                                    Treatment=c(rep("Before",59),rep("After",61)))

df%>%ggplot(aes(x=c(1:120),color=Treatment,group=1))+
              geom_rect(xmin=0,xmax=60,ymin=min(df$peakflow[1:60]),
                        ymax=max(df$peakflow[1:60]),fill="red",alpha=0.005,color=NA)+
              geom_rect(xmin=61,xmax=120,ymin=max(df$peakflow[1:60]),
                        ymax=max(df$peakflow[61:120]),fill="skyblue",alpha=0.01,color=NA)+
              geom_path(aes(y=peakflow),linetype=1,size=1)+
              scale_x_continuous("Following up time (Day)")+
              scale_y_continuous("Value",breaks=c(0,100,150,200,250,300,350,400,450,500))+
              geom_vline(xintercept=60,size=1,linetype=2)+
              theme_bw()+scale_color_manual(values=c("blue3","red3"))

Nhi tạo ra thêm 1 biến số nữa là Huyết áp tâm thu, với giả định biến này không có liên quan gì đến loại thuốc trị hen nói trên.

3 Package CausalImpact

Tác giả của method này đã tạo ra hẳn 1 package rất đơn giản, dễ sử dụng mà bạn có thể tải về từ CRAN: đó là package CausalImpact

https://google.github.io/CausalImpact/

Kết quả cho thấy loại thuốc trị hen này đã làm cải thiện chức năng thông khí của bệnh nhân : hiệu ứng tuyệt đối trung bình là 187 (180 – 194), hiệu ứng tương đối: tăng trung bình +116% (từ 111 tới 120%), xác suất hậu nghiệm của quan hệ nhân quả là 99.9% (p=0.001)

Tác giả còn cung cấp 1 văn bản diễn dịch rất cụ thể cho các bạn từ kết quả.

Và 1 hình vẽ rất đẹp cho thấy diễn tiến của biến kết quả trong 2 kịch bản : có và không có điều trị. Thật thú vị, điều này giống như bạn có thể du hành xuyên thời gian và nhìn thấy cả 2 hiện tượng khác nhau tại cùng 1 thời điểm ở tương lai cho cùng 1 bệnh nhân tùy theo quyết định của bạn trong quá khứ (đây là điều mà không có thiết kế thử nghiệm lâm sàng nào làm được, đơn giản là ta không thể quan sát cùng lúc hiệu ứng Placebo lẫn hiệu ứng của thuốc trên cùng 1 bệnh nhân một khi người này đã được phân nhóm).

Lưu ý là ngoài hiệu quả trung bình, còn có cả hiệu quả tích lũy (dùng trong trường hợp biến kết quả là 1 tần suất sự kiện hay số đếm).

library(CausalImpact)

model1=CausalImpact(df[,-3], c(1,60), c(61,120))

summary(model1)
## Posterior inference {CausalImpact}
## 
##                          Average        Cumulative    
## Actual                   349            20915         
## Prediction (s.d.)        162 (3.7)      9701 (219.6)  
## 95% CI                   [155, 169]     [9279, 10143] 
##                                                       
## Absolute effect (s.d.)   187 (3.7)      11214 (219.6) 
## 95% CI                   [180, 194]     [10772, 11636]
##                                                       
## Relative effect (s.d.)   116% (2.3%)    116% (2.3%)   
## 95% CI                   [111%, 120%]   [111%, 120%]  
## 
## Posterior tail-area probability p:   0.001
## Posterior prob. of a causal effect:  99.9%
## 
## For more details, type: summary(impact, "report")
plot(model1)

summary(model1, "report")
## Analysis report {CausalImpact}
## 
## 
## During the post-intervention period, the response variable had an average value of approx. 348.58. By contrast, in the absence of an intervention, we would have expected an average response of 161.68. The 95% interval of this counterfactual prediction is [154.64, 169.05]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 186.90 with a 95% interval of [179.53, 193.94]. For a discussion of the significance of this effect, see below.
## 
## Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of 20.91K. By contrast, had the intervention not taken place, we would have expected a sum of 9.70K. The 95% interval of this prediction is [9.28K, 10.14K].
## 
## The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +116%. The 95% interval of this percentage is [+111%, +120%].
## 
## This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (186.90) to the original goal of the underlying intervention.
## 
## The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.001). This means the causal effect can be considered statistically significant.

Nếu tò mò, bạn có thể thử dựng mô hình cho biến SBP, kết quả sẽ âm tính !

4 Kết luận:

Ngoài công việc chuyên môn, chúng ta nên quan sát và học hỏi thêm từ các chuyên ngành, lĩnh vực khác như Kinh tế học, Xã hội học, Tâm lý học. bên đó người ta có những phương pháp thống kê độc đáo mà nếu mình lĩnh hội được thì có thể mang áp dụng cho chuyên môn của mình.

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

LS0tDQp0aXRsZTogIkdp4bubaSB0aGnhu4d1IHBhY2thZ2UgQ2F1c2FsSW1wYWN0Ig0Kc3VidGl0bGU6ICJTdXkgZGnhu4VuIG5ow6JuIHF14bqjIGLhurFuZyBtw7QgaMOsbmggQmF5ZXNpYW4gc3RydWN0dXJhbCINCmF1dGhvcjogIkxlIE5nb2MgS2hhIE5oaSINCmRhdGU6ICIxMyBKdWx5IDIwMTciDQpvdXRwdXQ6DQogIGh0bWxfZG9jdW1lbnQ6IA0KICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUNCiAgICBudW1iZXJfc2VjdGlvbnM6IHllcw0KICAgIHRoZW1lOiAiZGVmYXVsdCINCiAgICB0b2M6IFRSVUUNCiAgICB0b2NfZmxvYXQ6IFRSVUUNCi0tLQ0KDQpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0NCmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkNCmBgYA0KDQohW10ocGtmbG93LnBuZykNCg0KIyBHaeG7m2kgdGhp4buHdQ0KDQpOZ2/DoGkgZ2nhu50gbMOgbSBjaMOtbmggdGjhu6ljLCBuZ2jhu4EgdGF5IHRyw6FpIGPhu6dhIE5oaSBsw6AgdMawIHbhuqVuIHbhu4EgcGjDom4gdMOtY2ggc+G7kSBsaeG7h3UsIG7Ds2kgbeG7mXQgY8OhY2ggYsOsbmggZMOibiBsw6AgxJFpIGzDoG0gdGjDtG5nIGvDqiDCqyBk4bqhbyDCuy4gQ8O0bmcgdmnhu4djIG7DoHkgY8OzIHRo4buDIHLhuqV0IG5o4bq5IG5ow6BuZyBob+G6t2MgY+G7sWMga2jhu5UsIHTDuXkgdsOgbyB0w6xuaCBodeG7kW5nLCDEkcO0aSBraGkgbmjhu69uZyDEkcOyaSBo4buPaSBrw6wgcXXhurdjIHbDoCBraMOzIGtoxINuIGPhu6dhIG3hu5l0IHbDoGkgdGjDom4gY2jhu6cgbOG6oWkgbMOgIMSR4buZbmcgbOG7sWMgxJHhu4MgTmhpIMSRaSB0w6xtIHTDsmkgaOG7jWMgaOG7j2kga2nhur9uIHRo4bupYywga+G7uSBuxINuZyBt4bubaS4gDQoNCkPDoWNoIMSRw6J5IDMgbmfDoHkgY8OzIG3DtHQgdGjDom4gY2jhu6cgxJHGsGEgcmEgecOqdSBj4bqndSBy4bqldCB0aMO6IHbhu4sgOiDDlG5nIHRhIG5naMSpIHJhIG3hu5l0IHBoxrDGoW5nIHBow6FwIMSRaeG7gXUgdHLhu4sgbeG7m2kgY2hvIHRyaeG7h3UgY2jhu6luZyB0aGnhur91IG94eSB0cm9uZyBraGkgbmfhu6cgZG8gaOG6uXAgxJHGsOG7nW5nIHRo4bufIHRyw6puLCBi4bqxbmcgY8OhY2ggZ+G6r24gMSBz4buRIMSRaeG7h24gY+G7sWMga8OtY2ggdGjDrWNoIG3hu5l0IHbDoGkgbmjDs20gY8ahIGjDtCBo4bqlcCwgYuG7h25oIG5ow6JuIGzhuq1wIHThu6ljIGPhuqNpIHRoaeG7h24gxJHGsOG7o2MgbMawdSBsxrDhu6NuZyBraMOtIHRo4bufIHbDoCBwaOG7pWMgaOG7k2kgxJHhu5kgYsOjbyBow7JhIE94eWdlbiAoU3BPMikuIMOUbmcgdGEgbXXhu5FuIE5oaSB0xrAgduG6pW4gbsOqbiBkw7luZyBjw7RuZyBj4bulIHRo4buRbmcga8OqIG7DoG8gxJHhu4MgY2jhu6luZyBtaW5oIGzDoCBrw61jaCB0aMOtY2ggxJFp4buHbiB0aOG7sWMgc+G7sSDEkcOjIGzDoG0gY+G6o2kgdGhp4buHbiBTcE8yLiBW4bqlbiDEkeG7gSBsw6AgU3BPMiDEkcaw4bujYyB0aGVvIGTDtWkgbGnDqm4gdOG7pWMgdHJvbmcga2hpIG5n4bunIG5oxrAgMSB0aW1lIHNlcmllcyBkYXRhIHbhu5tpIHThuqduIHN14bqldCBs4bqleSBt4bqrdSBsw6AgMSBIei4gxJDDonkgbMOgIDEgYsOgaSB0b8OhbiBjaOG7qW5nIG1pbmggbGnDqm4gaOG7hyBuaMOibiBxdeG6oyBy4bqldCBoaeG6v20gZ+G6t3AgdHJvbmcgbmdoacOqbiBj4bupdSB5IGjhu41jLg0KDQpUw6xtIGtp4bq/bSB0csOqbiBt4bqhbmcsIE5oaSBwaMOhdCBoaeG7h24gcmEgYsOgaSB0b8OhbiB0xrDGoW5nIHThu7EgY8OzIHThu5NuIHThuqFpIHRyb25nIG5o4buvbmcgbMSpbmggduG7sWMgaG/DoG4gdG/DoG4ga2jDoWMsIMSRw7MgbMOgIE1hcmtldGluZywgcXXhuqNuZyBjw6FvLCBLaW5oIHThur8gdMOgaSBjaMOtbmguIE3hu5l0IGdp4bqjaSBwaMOhcCBy4bqldCBoaeG7h3UgcXXhuqMgxJHDoyDEkcaw4bujYyBnaeG7m2kgdGhp4buHdSBi4bufaSBCcm9kZXJzZW4gS0ggKG5ow6JuIHZpw6puIGPhu6dhIEdvb2dsZSkgdOG7qyBuxINtIDIwMTUuDQoNCkJyb2RlcnNlbiBLSCwgR2FsbHVzc2VyIEYsIEtvZWhsZXIgSiwgUmVteSBOLCBTY290dCBTTC4gSW5mZXJyaW5nIGNhdXNhbCBpbXBhY3QgdXNpbmcgQmF5ZXNpYW4gc3RydWN0dXJhbCB0aW1lLXNlcmllcyBtb2RlbHMuIEFubmFscyBvZiBBcHBsaWVkIFN0YXRpc3RpY3MsIDIwMTUsIFZvbC4gOSwgTm8uIDEsIDI0Ny0yNzQuIGh0dHA6Ly9yZXNlYXJjaC5nb29nbGUuY29tL3B1YnMvcHViNDE4NTQuaHRtbA0KDQpN4buZdCBjw6FjaCDEkcahbiBnaeG6o24sIHTDoWMgZ2nhuqMgw6FwIGThu6VuZyAxIG3DtCBow6xuaCBj4bqldSB0csO6YyBCYXllcyB0csOqbiBt4buZdCBwaOG6p24gY+G7p2EgVGltZSBzZXJpZXMgdHLGsOG7m2Mga2hpIDEgY2FuIHRoaeG7h3AgYuG6r3QgxJHhuqd1LCBtw7QgaMOsbmggbsOgeSDEkcawYSByYSBk4buxIGLDoW8gduG7gSBkaeG7hW4gdGnhur9uIGPhu6dhIFRpbWUgc2VyaWUgbuG6v3UgbmjGsCBraMO0bmcgY8OzIGNhbiB0aGnhu4dwIG7DoG8geOG6o3kgcmEsIHThu6sgxJHDsyBzbyBzw6FuaCB24bubaSBr4bq/dCBxdeG6oyB0csOqbiB0aOG7sWMgdOG6vyBzYXUgY2FuIHRoaeG7h3AuIEvhur90IHF14bqjIHBow6JuIHTDrWNoIGzDoCAxIHBow6JuIHBo4buRaSBo4bqtdSBuZ2hp4buHbSB24buBIMSR4buZIGzhu5tuIGPhu6dhIGhp4buHdSDhu6luZyBjYW4gdGhp4buHcCBsw6puIGRp4buFbiB0aeG6v24gY+G7p2EgdGltZSBzZXJpZXMgbsOgeS4NCg0KQuG6sW5nIGPDoWNoIG7DoHksIG5nxrDhu51pIHRhIGPDsyB0aOG7gyBraOG6o28gc8OhdCBkb2FuaCBz4buRIGPhu6dhIDEgc+G6o24gcGjhuqltIFNBVSBraGkgYuG6r3QgxJHhuqd1IDEgY2hp4bq/biBk4buLY2ggcXXhuqNuZyBjw6FvIHRyw6puIHRydXnhu4FuIHRow7RuZz8NCg0KIyBCw6BpIHRvw6FuIG3DtCBwaOG7j25nDQoNCk5oaSBxdXnhur90IMSR4buLbmggdGjhu60gbmdoaeG7h20gcGjGsMahbmcgcGjDoXAgbsOgeSB0csOqbiAxIGLDoGkgdG/DoW4gbcO0IHBo4buPbmcgbmjGsCBzYXU6DQoNCkdp4bqjIHPhu60gdGEgbXXhu5FuIGto4bqjbyBzw6F0IGhp4buHdSBxdeG6oyDEkWnhu4F1IHRy4buLIGPhu6dhIDEgbG/huqFpIHRodeG7kWMgxJFp4buBdSB0cuG7iyBoZW4gcGjhur8gcXXhuqNuIHRyw6puIDEgYuG7h25oIG5ow6JuLiBL4bq/dCBxdeG6oyBsw6Aga2jhuqMgbsSDbmcgdGjDtG5nIGtow60gxJHGsOG7o2MgxJFvIGLhurFuZyBsxrB1IGzGsOG7o25nIMSR4buJbmggKFBlYWtmbG93KS4gQ2jhu4kgc+G7kSBQZWFrZmxvdyDEkcaw4bujYyDEkW8gbeG7l2kgbmfDoHkgdsOgbyAxIGdp4budIG5o4bqldCDEkeG7i25oIHRyb25nIDIgdGjDoW5nIHRyxrDhu5tjIHbDoCBzYXUga2hpIMSRaeG7gXUgdHLhu4sgYuG6sW5nIGxv4bqhaSB0aHXhu5FjIG7DoHkuIFRhIGPDsyAxIHRpbWUgc2VyaWVzIDEyMCDEkWnhu4NtLCDEkWnhu4NtIGto4bufaSDEkeG6p3UgxJFp4buBdSB0cuG7iyBu4bqxbSDhu58gbmfDoHkgNjEuDQoNClPhu5EgbGnhu4d1IMSRxrDhu6NjIG3DtCBwaOG7j25nIG5oxrAgc2F1Og0KDQpgYGB7cixtZXNzYWdlID0gRkFMU0Usd2FybmluZz1GQUxTRX0NCmxpYnJhcnkodGlkeXZlcnNlKQ0Kc2V0LnNlZWQoMTIzKQ0KDQpwZWFrZmxvdz1jKHJvdW5kKHJub3JtKDYwLDE2MCwyNSkrYXJpbWEuc2ltKG1vZGVsPWxpc3QoYXI9MC4zKSxuPTYwKSwwKSwNCiAgICAgICAgICAgcm91bmQocm5vcm0oNjAsMzUwLDQwKSthcmltYS5zaW0obW9kZWw9bGlzdChhcj0wLjMpLG49NjApLDApKQ0KDQpkZj1hcy5kYXRhLmZyYW1lKHBlYWtmbG93KSU+JW11dGF0ZSguLFNCUD1yb3VuZChybm9ybSgxMjAsMTAwLDIpK2FyaW1hLnNpbShtb2RlbD1saXN0KGFyPTAuMyksbj0xMjApLDApLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVHJlYXRtZW50PWMocmVwKCJCZWZvcmUiLDU5KSxyZXAoIkFmdGVyIiw2MSkpKQ0KDQpkZiU+JWdncGxvdChhZXMoeD1jKDE6MTIwKSxjb2xvcj1UcmVhdG1lbnQsZ3JvdXA9MSkpKw0KICAgICAgICAgICAgICBnZW9tX3JlY3QoeG1pbj0wLHhtYXg9NjAseW1pbj1taW4oZGYkcGVha2Zsb3dbMTo2MF0pLA0KICAgICAgICAgICAgICAgICAgICAgICAgeW1heD1tYXgoZGYkcGVha2Zsb3dbMTo2MF0pLGZpbGw9InJlZCIsYWxwaGE9MC4wMDUsY29sb3I9TkEpKw0KICAgICAgICAgICAgICBnZW9tX3JlY3QoeG1pbj02MSx4bWF4PTEyMCx5bWluPW1heChkZiRwZWFrZmxvd1sxOjYwXSksDQogICAgICAgICAgICAgICAgICAgICAgICB5bWF4PW1heChkZiRwZWFrZmxvd1s2MToxMjBdKSxmaWxsPSJza3libHVlIixhbHBoYT0wLjAxLGNvbG9yPU5BKSsNCiAgICAgICAgICAgICAgZ2VvbV9wYXRoKGFlcyh5PXBlYWtmbG93KSxsaW5ldHlwZT0xLHNpemU9MSkrDQogICAgICAgICAgICAgIHNjYWxlX3hfY29udGludW91cygiRm9sbG93aW5nIHVwIHRpbWUgKERheSkiKSsNCiAgICAgICAgICAgICAgc2NhbGVfeV9jb250aW51b3VzKCJWYWx1ZSIsYnJlYWtzPWMoMCwxMDAsMTUwLDIwMCwyNTAsMzAwLDM1MCw0MDAsNDUwLDUwMCkpKw0KICAgICAgICAgICAgICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9NjAsc2l6ZT0xLGxpbmV0eXBlPTIpKw0KICAgICAgICAgICAgICB0aGVtZV9idygpK3NjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXM9YygiYmx1ZTMiLCJyZWQzIikpDQoNCmBgYA0KDQpOaGkgdOG6oW8gcmEgdGjDqm0gMSBiaeG6v24gc+G7kSBu4buvYSBsw6AgSHV54bq/dCDDoXAgdMOibSB0aHUsIHbhu5tpIGdp4bqjIMSR4buLbmggYmnhur9uIG7DoHkga2jDtG5nIGPDsyBsacOqbiBxdWFuIGfDrCDEkeG6v24gbG/huqFpIHRodeG7kWMgdHLhu4sgaGVuIG7Ds2kgdHLDqm4uDQoNCiMgUGFja2FnZSBDYXVzYWxJbXBhY3QNCg0KVMOhYyBnaeG6oyBj4bunYSBtZXRob2QgbsOgeSDEkcOjIHThuqFvIHJhIGjhurNuIDEgcGFja2FnZSBy4bqldCDEkcahbiBnaeG6o24sIGThu4Ugc+G7rSBk4bulbmcgbcOgIGLhuqFuIGPDsyB0aOG7gyB04bqjaSB24buBIHThu6sgQ1JBTjogxJHDsyBsw6AgcGFja2FnZSBDYXVzYWxJbXBhY3QNCg0KPGh0dHBzOi8vZ29vZ2xlLmdpdGh1Yi5pby9DYXVzYWxJbXBhY3QvPg0KDQpL4bq/dCBxdeG6oyBjaG8gdGjhuqV5IGxv4bqhaSB0aHXhu5FjIHRy4buLIGhlbiBuw6B5IMSRw6MgbMOgbSBj4bqjaSB0aGnhu4duIGNo4bupYyBuxINuZyB0aMO0bmcga2jDrSBj4bunYSBi4buHbmggbmjDom4gOiBoaeG7h3Ug4bupbmcgdHV54buHdCDEkeG7kWkgdHJ1bmcgYsOsbmggbMOgIDE4NyAoMTgwIOKAkyAxOTQpLCBoaeG7h3Ug4bupbmcgdMawxqFuZyDEkeG7kWk6IHTEg25nIHRydW5nIGLDrG5oICsxMTYlICh04burIDExMSB04bubaSAxMjAlKSwgeMOhYyBzdeG6pXQgaOG6rXUgbmdoaeG7h20gY+G7p2EgcXVhbiBo4buHIG5ow6JuIHF14bqjIGzDoCA5OS45JSAocD0wLjAwMSkNCg0KVMOhYyBnaeG6oyBjw7JuIGN1bmcgY+G6pXAgMSB2xINuIGLhuqNuIGRp4buFbiBk4buLY2ggcuG6pXQgY+G7pSB0aOG7gyBjaG8gY8OhYyBi4bqhbiB04burIGvhur90IHF14bqjLg0KDQpWw6AgMSBow6xuaCB24bq9IHLhuqV0IMSR4bq5cCBjaG8gdGjhuqV5IGRp4buFbiB0aeG6v24gY+G7p2EgYmnhur9uIGvhur90IHF14bqjIHRyb25nIDIga+G7i2NoIGLhuqNuIDogY8OzIHbDoCBraMO0bmcgY8OzIMSRaeG7gXUgdHLhu4suIFRo4bqtdCB0aMO6IHbhu4ssIMSRaeG7gXUgbsOgeSBnaeG7kW5nIG5oxrAgYuG6oW4gY8OzIHRo4buDIGR1IGjDoG5oIHh1ecOqbiB0aOG7nWkgZ2lhbiB2w6AgbmjDrG4gdGjhuqV5IGPhuqMgMiBoaeG7h24gdMaw4bujbmcga2jDoWMgbmhhdSB04bqhaSBjw7luZyAxIHRo4budaSDEkWnhu4NtIOG7nyB0xrDGoW5nIGxhaSBjaG8gY8O5bmcgMSBi4buHbmggbmjDom4gdMO5eSB0aGVvIHF1eeG6v3QgxJHhu4tuaCBj4bunYSBi4bqhbiB0cm9uZyBxdcOhIGto4bupICjEkcOieSBsw6AgxJFp4buBdSBtw6Aga2jDtG5nIGPDsyB0aGnhur90IGvhur8gdGjhu60gbmdoaeG7h20gbMOibSBzw6BuZyBuw6BvIGzDoG0gxJHGsOG7o2MsIMSRxqFuIGdp4bqjbiBsw6AgdGEga2jDtG5nIHRo4buDIHF1YW4gc8OhdCBjw7luZyBsw7pjIGhp4buHdSDhu6luZyBQbGFjZWJvIGzhuqtuIGhp4buHdSDhu6luZyBj4bunYSB0aHXhu5FjIHRyw6puIGPDuW5nIDEgYuG7h25oIG5ow6JuIG3hu5l0IGtoaSBuZ8aw4budaSBuw6B5IMSRw6MgxJHGsOG7o2MgcGjDom4gbmjDs20pLg0KDQpMxrB1IMO9IGzDoCBuZ2/DoGkgaGnhu4d1IHF14bqjIHRydW5nIGLDrG5oLCBjw7JuIGPDsyBj4bqjIGhp4buHdSBxdeG6oyB0w61jaCBsxal5IChkw7luZyB0cm9uZyB0csaw4budbmcgaOG7o3AgYmnhur9uIGvhur90IHF14bqjIGzDoCAxIHThuqduIHN14bqldCBz4buxIGtp4buHbiBoYXkgc+G7kSDEkeG6v20pLg0KDQpgYGB7cixtZXNzYWdlID0gRkFMU0Usd2FybmluZz1GQUxTRX0NCmxpYnJhcnkoQ2F1c2FsSW1wYWN0KQ0KDQptb2RlbDE9Q2F1c2FsSW1wYWN0KGRmWywtM10sIGMoMSw2MCksIGMoNjEsMTIwKSkNCg0Kc3VtbWFyeShtb2RlbDEpDQoNCnBsb3QobW9kZWwxKQ0KICANCnN1bW1hcnkobW9kZWwxLCAicmVwb3J0IikNCg0KYGBgDQoNCk7hur91IHTDsiBtw7IsIGLhuqFuIGPDsyB0aOG7gyB0aOG7rSBk4buxbmcgbcO0IGjDrG5oIGNobyBiaeG6v24gU0JQLCBr4bq/dCBxdeG6oyBz4bq9IMOibSB0w61uaCAhDQoNCiMgS+G6v3QgbHXhuq1uOg0KDQpOZ2/DoGkgY8O0bmcgdmnhu4djIGNodXnDqm4gbcO0biwgY2jDum5nIHRhIG7Dqm4gcXVhbiBzw6F0IHbDoCBo4buNYyBo4buPaSB0aMOqbSB04burIGPDoWMgY2h1ecOqbiBuZ8OgbmgsIGzEqW5oIHbhu7FjIGtow6FjIG5oxrAgS2luaCB04bq/IGjhu41jLCBYw6MgaOG7mWkgaOG7jWMsIFTDom0gbMO9IGjhu41jLiBiw6puIMSRw7MgbmfGsOG7nWkgdGEgY8OzIG5o4buvbmcgcGjGsMahbmcgcGjDoXAgdGjhu5FuZyBrw6ogxJHhu5ljIMSRw6FvIG3DoCBu4bq/dSBtw6xuaCBsxKluaCBo4buZaSDEkcaw4bujYyB0aMOsIGPDsyB0aOG7gyBtYW5nIMOhcCBk4bulbmcgY2hvIGNodXnDqm4gbcO0biBj4bunYSBtw6xuaC4gDQoNCipU4bqhbSBiaeG7h3QgY8OhYyBi4bqhbiB2w6AgaOG6uW4gZ+G6t3AgbOG6oWkuLi4q