R Notebook

《白话统计学》读书笔记

第三章：关于统计资料的思考

计数资料和分类资料的区别

计数资料：其实就是对于某个事件进行计数的资料，例如对于咳嗽次数的计数。而分类资料则是一个事件的两种类别，比如男女。所以两者是不同的数据类型。
统计分析上来看，

由于计数资料，一般来说服从Poisson分布，所以在回归分析的时候使用Poisson回归或者负二项回归。两个回归之间的区别在于，Poisson一般用于个体之间相互独立的情形，而负二项则用于个体之间不独立的情形，比如说咳嗽是相互传染的，那么分析的时候需要用到负二项。另外如果条件方差超过条件均值时也是推荐使用负二项的。我们使用robust包中的Breslow癫痫数据进行常见的Poisson回归分析演示。我们就遭受轻微或严重间歇性癫痫的病人的年龄和癫痫发病数收集了数据，包含病人被随机分配到药物组或者安慰剂组前八周和随机分配后八周两种情况。响应变量为sumY(随机化后八周内癫痫发病数)，预测变量为治疗条件(Trt)、年龄(Age)和前八周内的基础癫痫发病数(Base)。

library(robust)
library(tidyverse)
data("breslow.dat")
smallData <- breslow.dat[,c(6:8,10)]
summary(smallData)

      Base             Age               Trt          sumY       
 Min.   :  6.00   Min.   :18.00   placebo  :28   Min.   :  0.00  
 1st Qu.: 12.00   1st Qu.:23.00   progabide:31   1st Qu.: 11.50  
 Median : 22.00   Median :28.00                  Median : 16.00  
 Mean   : 31.22   Mean   :28.34                  Mean   : 33.05  
 3rd Qu.: 41.00   3rd Qu.:32.00                  3rd Qu.: 36.00  
 Max.   :151.00   Max.   :42.00                  Max.   :302.00

####查看响应变量的基本分布
hist(smallData$sumY, breaks = 20)

####Poisson回归分析
model <- glm(sumY ~ Base + Age + Trt, data = smallData, family = poisson())
summary(model)


Call:
glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = smallData)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-6.0569  -2.0433  -0.9397   0.7929  11.0061  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***
Base          0.0226517  0.0005093  44.476  < 2e-16 ***
Age           0.0227401  0.0040240   5.651 1.59e-08 ***
Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2122.73  on 58  degrees of freedom
Residual deviance:  559.44  on 55  degrees of freedom
AIC: 850.71

Number of Fisher Scoring iterations: 5

####检验过度离势
qcc::qcc.overdispersion.test(smallData$sumY, type = "poisson")

                   
Overdispersion test Obs.Var/Theor.Var Statistic p-value
       poisson data          62.87013  3646.468       0

过度离势检验<0.05的时候。表示存在过度离势。因此需要把family参数转换为:quasipoisson

###负二项回归
library(MASS)
modle <- glm.nb(sumY ~ Base + Age + Trt, data = smallData)
summary(modle)


Call:
glm.nb(formula = sumY ~ Base + Age + Trt, data = smallData, init.theta = 3.361505307, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.3755  -0.7517  -0.1916   0.2817   2.6103  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.931832   0.408372   4.731 2.24e-06 ***
Base          0.027736   0.002817   9.846  < 2e-16 ***
Age           0.016815   0.012535   1.341     0.18    
Trtprogabide -0.193518   0.154271  -1.254     0.21    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(3.3615) family taken to be 1)

    Null deviance: 181.856  on 58  degrees of freedom
Residual deviance:  63.667  on 55  degrees of freedom
AIC: 476.46

Number of Fisher Scoring iterations: 1

              Theta:  3.362 
          Std. Err.:  0.710 

 2 x log-likelihood:  -466.457

而分类资料一般选用Logistic回归。对于分类资料而言，如果是二分类资料使用基本的Logistic回归即可。如果是多分类的则需要看分组之间是否有序，适当的选择有序/无序Logistic回归即可。

###简单的二分类Logistic回归
data(anesthetic, package = "DAAG")
head(anesthetic)

fit_logit <- glm(nomove ~ conc, family = binomial(), data = anesthetic)
summary(fit_logit)


Call:
glm(formula = nomove ~ conc, family = binomial(), data = anesthetic)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.76666  -0.74407   0.03413   0.68666   2.06900  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   -6.469      2.418  -2.675  0.00748 **
conc           5.567      2.044   2.724  0.00645 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 41.455  on 29  degrees of freedom
Residual deviance: 27.754  on 28  degrees of freedom
AIC: 31.754

Number of Fisher Scoring iterations: 5

多分类的话，如果分类是无序的则使用nnet::multinom即可。

###无序多分类
require(foreign)
require(nnet)
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
ml$prog2 <- relevel(ml$prog, ref = "academic")
head(ml)

###无序logit
test <- multinom(prog2 ~ ses + write, data = ml)

Error in multinom(prog2 ~ ses + write, data = ml) : 参数没有用(data = ml)

如果分类是有序的，则可是使用MASS::polr

require(foreign)
require(MASS)
require(Hmisc)
###读取数据
dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
####有序logit回归
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)
####定义结果返回函数
order_logires <- function(model, n){
    ctable <- coef(summary(model))
    p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
    ci <- confint(model)
    if(n == 1){
        res <- c(exp(c(OR = coef(model), ci)), p = p[1])
    }else{
        res <- cbind(exp(cbind(OR = coef(model), ci)), p = p[1:n])
    }
    return(res)
}
###返回OR, 95%CI, P
order_logires(m, 3)

Waiting for profiling to be done...

              OR     2.5 %   97.5 %            p
pared  2.8510579 1.6958376 4.817114 8.087072e-05
public 0.9429088 0.5208954 1.680579 8.435464e-01
gpa    1.8513972 1.1136247 3.098490 1.811594e-02

计数资料可否采用连续性变量来分析

由于计数资料也是有很多分类的，那么在看其回归的时候，能否使用线性回归。这个一般来说是不能的。线性回归需要满足的条件之一是正态。计数资料因为一般没有负数所以很有可能是偏态的。所以是不能是有的。
不过如果计数资料取值离0很远。同时呈正态，如果研究者不介意预测结果出现小数点或者负数则可以使用。

分类资料分析统计方法的选择

对于分类资料，如果要观察其等级关系的结果的时候，最好使用秩和检验(Wilcoxon)。如果只是观察各组之间有没有差异则最好使用卡方检验

寻找cutoff的多种方法

根据专业知识来决定

基于专业知识划分是最好的cutoff的方式。模型计算出来的东西。必定不具备一定的说服力。

利用广义客家模型结合专业划分

当我们在选择cut-off的时候，没有办法完全利用专业知识来划分的时候。但是如果有一定的参考，再结合专业知识的话就可以进行划分的时候，那么就可以用到广义可加模型。相较于一般的线性模型，广义可加模型绘制出来的曲线不一定是线性的。这个模型主要可以用来探索自变量和因变量的关系。 - R语言中可是使用mgcv::gam来进行模型预测

library(mgcv)
dat <- read.csv("dat1.csv")
head(dat)

fit <- gam(hyper ~ s(age), data = dat, family = binomial())
plot(fit)
abline(h = 0)

由例子可见。可以吧年龄三分类。

利用ROC曲线

如果有一个明确的二分类结局的话，可以使用ROC曲线来寻找最佳的cutoff。 - R里面可是使用pROC包进行结果的统计

利用最大选择秩统计量

如果当结局变量函数时间变量，或者说是分类变量的时候.可是使用最大选择秩统计量来进行最佳cutoff的选择。其主要思想也是把所有可能的分组都计算一遍。然后寻找最佳的结果。这种方式最多用于在生存分析的时候寻找最佳的cutoff - R中可以通过maxstat来实现。

library(maxstat)
set.seed(29)
x <- sort(runif(20))
y <- c(rnorm(10), rnorm(10, 2))
mydata <- data.frame(cbind(x,y))
head(mydata)

mod <- maxstat.test(y ~ x, data=mydata, smethod="Wilcoxon", pmethod="HL",
                    minprop=0.25, maxprop=0.75, alpha=0.05)
plot(mod)

使用分类树进行划分

使用分类树，类似于决策曲线一样的来说明如果进行分组。

library(rpart)
fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis)
plot(fit)
text(fit, use.n = TRUE)

聚类分析

之前的方法都有一个条件即必须有一个明确的确定的结局。这样根据结局对自变量进行划分，通常将这种情况称为有监督的。但是如果我们没有结局变量的时候，如果要讲样本进行划分的话，这样就需要用到无监督的聚类分析。

常用的方式发层次法。K-means法等等。

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMg44CK55m96K+d57uf6K6h5a2m44CL6K+75Lmm56yU6K6wCgojIyMg56ys5LiJ56ug77ya5YWz5LqO57uf6K6h6LWE5paZ55qE5oCd6ICDCgojIyMjIOiuoeaVsOi1hOaWmeWSjOWIhuexu+i1hOaWmeeahOWMuuWIqwoKLSAqKuiuoeaVsOi1hOaWmSoq77ya5YW25a6e5bCx5piv5a+55LqO5p+Q5Liq5LqL5Lu26L+b6KGM6K6h5pWw55qE6LWE5paZ77yM5L6L5aaC5a+55LqO5ZKz5Ze95qyh5pWw55qE6K6h5pWw44CC6ICMKirliIbnsbvotYTmlpkqKuWImeaYr+S4gOS4quS6i+S7tueahOS4pOenjeexu+WIq++8jOavlOWmgueUt+Wls+OAguaJgOS7peS4pOiAheaYr+S4jeWQjOeahOaVsOaNruexu+Wei+OAggoKLSDnu5/orqHliIbmnpDkuIrmnaXnnIvvvIwKMS4g55Sx5LqOKirorqHmlbDotYTmlpkqKu+8jOS4gOiIrOadpeivtOacjeS7jioqUG9pc3NvbuWIhuW4gyoq77yM5omA5Lul5Zyo5Zue5b2S5YiG5p6Q55qE5pe25YCZ5L2/55SoKipQb2lzc29u5Zue5b2S5oiW6ICF6LSf5LqM6aG55Zue5b2SKirjgILkuKTkuKrlm57lvZLkuYvpl7TnmoTljLrliKvlnKjkuo7vvIwqKlBvaXNzb24qKuS4gOiIrOeUqOS6juS4quS9k+S5i+mXtOebuOS6kueLrOeri+eahOaDheW9ou+8jOiAjCoq6LSf5LqM6aG5KirliJnnlKjkuo7kuKrkvZPkuYvpl7TkuI3ni6znq4vnmoTmg4XlvaLvvIzmr5TlpoLor7TlkrPll73mmK/nm7jkupLkvKDmn5PnmoTvvIzpgqPkuYjliIbmnpDnmoTml7blgJnpnIDopoHnlKjliLAqKui0n+S6jOmhuSoq44CC5Y+m5aSW5aaC5p6cKirmnaHku7bmlrnlt67otoXov4fmnaHku7blnYflgLzml7YqKuS5n+aYr+aOqOiNkOS9v+eUqCoq6LSf5LqM6aG5KirnmoTjgIIKICAgIOaIkeS7rOS9v+eUqGByb2J1c3Rg5YyF5Lit55qEYEJyZXNsb3dg55mr55er5pWw5o2u6L+b6KGM5bi46KeB55qEUG9pc3NvbuWbnuW9kuWIhuaekOa8lOekuuOAguaIkeS7rOWwsemBreWPl+i9u+W+ruaIluS4pemHjemXtOath+aAp+eZq+eXq+eahOeXheS6uueahOW5tOm+hOWSjOeZq+eXq+WPkeeXheaVsOaUtumbhuS6huaVsOaNru+8jOWMheWQq+eXheS6uuiiq+maj+acuuWIhumFjeWIsOiNr+eJqee7hOaIluiAheWuieaFsOWJgue7hOWJjeWFq+WRqOWSjOmaj+acuuWIhumFjeWQjuWFq+WRqOS4pOenjeaDheWGteOAgioq5ZON5bqU5Y+Y6YePKirkuLpzdW1ZKOmaj+acuuWMluWQjuWFq+WRqOWGheeZq+eXq+WPkeeXheaVsCnvvIwqKumihOa1i+WPmOmHjyoq5Li65rK755aX5p2h5Lu2KFRydCnjgIHlubTpvoQoQWdlKeWSjOWJjeWFq+WRqOWGheeahOWfuuehgOeZq+eXq+WPkeeXheaVsChCYXNlKeOAggpgYGB7cn0KbGlicmFyeShyb2J1c3QpCmxpYnJhcnkodGlkeXZlcnNlKQpkYXRhKCJicmVzbG93LmRhdCIpCnNtYWxsRGF0YSA8LSBicmVzbG93LmRhdFssYyg2OjgsMTApXQpzdW1tYXJ5KHNtYWxsRGF0YSkKIyMjI+afpeeci+WTjeW6lOWPmOmHj+eahOWfuuacrOWIhuW4gwpoaXN0KHNtYWxsRGF0YSRzdW1ZLCBicmVha3MgPSAyMCkKIyMjI1BvaXNzb27lm57lvZLliIbmnpAKbW9kZWwgPC0gZ2xtKHN1bVkgfiBCYXNlICsgQWdlICsgVHJ0LCBkYXRhID0gc21hbGxEYXRhLCBmYW1pbHkgPSBwb2lzc29uKCkpCnN1bW1hcnkobW9kZWwpCiMjIyPmo4Dpqozov4fluqbnprvlir8KcWNjOjpxY2Mub3ZlcmRpc3BlcnNpb24udGVzdChzbWFsbERhdGEkc3VtWSwgdHlwZSA9ICJwb2lzc29uIikKYGBgCj4g6L+H5bqm56a75Yq/5qOA6aqMPDAuMDXnmoTml7blgJnjgILooajnpLrlrZjlnKjov4fluqbnprvlir/jgILlm6DmraTpnIDopoHmiopgZmFtaWx5YOWPguaVsOi9rOaNouS4ujpgcXVhc2lwb2lzc29uYCAKCmBgYHtyfQojIyPotJ/kuozpobnlm57lvZIKbGlicmFyeShNQVNTKQptb2RsZSA8LSBnbG0ubmIoc3VtWSB+IEJhc2UgKyBBZ2UgKyBUcnQsIGRhdGEgPSBzbWFsbERhdGEpCnN1bW1hcnkobW9kbGUpCmBgYAoKMi4g6ICMKirliIbnsbvotYTmlpkqKuS4gOiIrOmAieeUqCoqTG9naXN0aWPlm57lvZIqKuOAguWvueS6juWIhuexu+i1hOaWmeiAjOiogO+8jOWmguaenOaYr+S6jOWIhuexu+i1hOaWmeS9v+eUqOWfuuacrOeahCoqTG9naXN0aWPlm57lvZIqKuWNs+WPr+OAguWmguaenOaYr+WkmuWIhuexu+eahOWImemcgOimgeeci+WIhue7hOS5i+mXtOaYr+WQpuacieW6j++8jOmAguW9k+eahOmAieaLqSoq5pyJ5bqPL+aXoOW6j0xvZ2lzdGlj5Zue5b2SKirljbPlj6/jgIIKCmBgYHtyfQojIyPnroDljZXnmoTkuozliIbnsbtMb2dpc3RpY+WbnuW9kgpkYXRhKGFuZXN0aGV0aWMsIHBhY2thZ2UgPSAiREFBRyIpCmhlYWQoYW5lc3RoZXRpYykKZml0X2xvZ2l0IDwtIGdsbShub21vdmUgfiBjb25jLCBmYW1pbHkgPSBiaW5vbWlhbCgpLCBkYXRhID0gYW5lc3RoZXRpYykKc3VtbWFyeShmaXRfbG9naXQpCmBgYAoKPiDlpJrliIbnsbvnmoTor53vvIzlpoLmnpzliIbnsbvmmK8qKuaXoOW6jyoq55qE5YiZ5L2/55SoYG5uZXQ6Om11bHRpbm9tYOWNs+WPr+OAggoKYGBge3J9CiMjI+aXoOW6j+WkmuWIhuexuwpyZXF1aXJlKGZvcmVpZ24pCnJlcXVpcmUobm5ldCkKbWwgPC0gcmVhZC5kdGEoImh0dHBzOi8vc3RhdHMuaWRyZS51Y2xhLmVkdS9zdGF0L2RhdGEvaHNiZGVtby5kdGEiKQptbCRwcm9nMiA8LSByZWxldmVsKG1sJHByb2csIHJlZiA9ICJhY2FkZW1pYyIpCmhlYWQobWwpCiMjI+aXoOW6j2xvZ2l0CnRlc3QgPC0gbXVsdGlub20ocHJvZzIgfiBzZXMgKyB3cml0ZSwgZGF0YSA9IG1sKQpzdW1tYXJ5KHRlc3QpCiMjI+iuoeeul1DlgLwKeiA8LSBzdW1tYXJ5KHRlc3QpJGNvZWZmaWNpZW50cy9zdW1tYXJ5KHRlc3QpJHN0YW5kYXJkLmVycm9ycwpwIDwtICgxIC0gcG5vcm0oYWJzKHopLCAwLCAxKSkgKiAyCnAKYGBgCgo+IOWmguaenOWIhuexu+aYryoq5pyJ5bqPKirnmoTvvIzliJnlj6/mmK/kvb/nlKhgTUFTUzo6cG9scmAKCmBgYHtyfQpyZXF1aXJlKGZvcmVpZ24pCnJlcXVpcmUoTUFTUykKcmVxdWlyZShIbWlzYykKIyMj6K+75Y+W5pWw5o2uCmRhdCA8LSByZWFkLmR0YSgiaHR0cHM6Ly9zdGF0cy5pZHJlLnVjbGEuZWR1L3N0YXQvZGF0YS9vbG9naXQuZHRhIikKIyMjI+acieW6j2xvZ2l05Zue5b2SCm0gPC0gcG9scihhcHBseSB+IHBhcmVkICsgcHVibGljICsgZ3BhLCBkYXRhID0gZGF0LCBIZXNzPVRSVUUpCiMjIyPlrprkuYnnu5Pmnpzov5Tlm57lh73mlbAKb3JkZXJfbG9naXJlcyA8LSBmdW5jdGlvbihtb2RlbCwgbil7CiAgICBjdGFibGUgPC0gY29lZihzdW1tYXJ5KG1vZGVsKSkKICAgIHAgPC0gcG5vcm0oYWJzKGN0YWJsZVssICJ0IHZhbHVlIl0pLCBsb3dlci50YWlsID0gRkFMU0UpICogMgogICAgY2kgPC0gY29uZmludChtb2RlbCkKICAgIGlmKG4gPT0gMSl7CiAgICAgICAgcmVzIDwtIGMoZXhwKGMoT1IgPSBjb2VmKG1vZGVsKSwgY2kpKSwgcCA9IHBbMV0pCiAgICB9ZWxzZXsKICAgICAgICByZXMgPC0gY2JpbmQoZXhwKGNiaW5kKE9SID0gY29lZihtb2RlbCksIGNpKSksIHAgPSBwWzE6bl0pCiAgICB9CiAgICByZXR1cm4ocmVzKQp9CiMjI+i/lOWbnk9SLCA5NSVDSSwgUApvcmRlcl9sb2dpcmVzKG0sIDMpCmBgYAoKIyMjIyDorqHmlbDotYTmlpnlj6/lkKbph4fnlKjov57nu63mgKflj5jph4/mnaXliIbmnpAKCi0g55Sx5LqO6K6h5pWw6LWE5paZ5Lmf5piv5pyJ5b6I5aSa5YiG57G755qE77yM6YKj5LmI5Zyo55yL5YW25Zue5b2S55qE5pe25YCZ77yM6IO95ZCm5L2/55So57q/5oCn5Zue5b2S44CC6L+Z5Liq5LiA6Iis5p2l6K+05piv5LiN6IO955qE44CC57q/5oCn5Zue5b2S6ZyA6KaB5ruh6Laz55qE5p2h5Lu25LmL5LiA5pivKirmraPmgIEqKuOAguiuoeaVsOi1hOaWmeWboOS4uuS4gOiIrOayoeaciei0n+aVsOaJgOS7peW+iOacieWPr+iDveaYryoq5YGP5oCBKirnmoTjgILmiYDku6XmmK/kuI3og73mmK/mnInnmoTjgIIKLSDkuI3ov4flpoLmnpzorqHmlbDotYTmlpnlj5blgLznprsw5b6I6L+c44CC5ZCM5pe25ZGIKirmraPmgIEqKu+8jOWmguaenOeglOeptuiAheS4jeS7i+aEj+mihOa1i+e7k+aenOWHuueOsCoq5bCP5pWw54K5KirmiJbogIUqKui0n+aVsCoq5YiZ5Y+v5Lul5L2/55So44CCCgojIyMjIOWIhuexu+i1hOaWmeWIhuaekOe7n+iuoeaWueazleeahOmAieaLqQoK5a+55LqO5YiG57G76LWE5paZ77yM5aaC5p6c6KaB6KeC5a+f5YW2562J57qn5YWz57O755qE57uT5p6c55qE5pe25YCZ77yM5pyA5aW95L2/55SoKirnp6nlkozmo4DpqowoV2lsY294b24pKirjgILlpoLmnpzlj6rmmK/op4Llr5/lkITnu4TkuYvpl7TmnInmsqHmnInlt67lvILliJnmnIDlpb3kvb/nlKgqKuWNoeaWueajgOmqjCoqCgojIyMjIOWvu+aJvmN1dG9mZueahOWkmuenjeaWueazlQoKIyMjIyMg5qC55o2u5LiT5Lia55+l6K+G5p2l5Yaz5a6aCgrln7rkuo7kuJPkuJrnn6Xor4bliJLliIbmmK/mnIDlpb3nmoRjdXRvZmbnmoTmlrnlvI/jgILmqKHlnovorqHnrpflh7rmnaXnmoTkuJzopb/jgILlv4XlrprkuI3lhbflpIfkuIDlrprnmoTor7TmnI3lipvjgIIKCiMjIyMjIOWIqeeUqOW5v+S5ieWuouWutuaooeWei+e7k+WQiOS4k+S4muWIkuWIhgoK5b2T5oiR5Lus5Zyo6YCJ5oupY3V0LW9mZueahOaXtuWAme+8jOayoeacieWKnuazleWujOWFqOWIqeeUqOS4k+S4muefpeivhuadpeWIkuWIhueahOaXtuWAmeOAguS9huaYr+WmguaenOacieS4gOWumueahOWPguiAg++8jOWGjee7k+WQiOS4k+S4muefpeivhueahOivneWwseWPr+S7pei/m+ihjOWIkuWIhueahOaXtuWAme+8jOmCo+S5iOWwseWPr+S7peeUqOWIsCoq5bm/5LmJ5Y+v5Yqg5qih5Z6LKirjgILnm7jovoPkuo7kuIDoiKznmoQqKue6v+aAp+aooeWeiyoq77yMKirlub/kuYnlj6/liqDmqKHlnosqKue7mOWItuWHuuadpeeahOabsue6v+S4jeS4gOWumuaYr+e6v+aAp+eahOOAgui/meS4quaooeWei+S4u+imgeWPr+S7peeUqOadpeaOoue0ouiHquWPmOmHj+WSjOWboOWPmOmHj+eahOWFs+ezu+OAggotIFLor63oqIDkuK3lj6/mmK/kvb/nlKhgbWdjdjo6Z2FtYOadpei/m+ihjOaooeWei+mihOa1iwoKYGBge3J9CmxpYnJhcnkobWdjdikKZGF0IDwtIHJlYWQuY3N2KCJkYXQxLmNzdiIpCmhlYWQoZGF0KQpmaXQgPC0gZ2FtKGh5cGVyIH4gcyhhZ2UpLCBkYXRhID0gZGF0LCBmYW1pbHkgPSBiaW5vbWlhbCgpKQpwbG90KGZpdCkKYWJsaW5lKGggPSAwKQpgYGAKCj4g55Sx5L6L5a2Q5Y+v6KeB44CC5Y+v5Lul5ZCn5bm06b6E5LiJ5YiG57G744CCCgojIyMjIyDliKnnlKhST0Pmm7Lnur8KCuWmguaenOacieS4gOS4quaYjuehrueahOS6jOWIhuexu+e7k+WxgOeahOivne+8jOWPr+S7peS9v+eUqFJPQ+absue6v+adpeWvu+aJvuacgOS9s+eahGN1dG9mZuOAggotIFLph4zpnaLlj6/mmK/kvb/nlKhgcFJPQ2DljIXov5vooYznu5PmnpznmoTnu5/orqEKCiMjIyMjIOWIqeeUqOacgOWkp+mAieaLqeenqee7n+iuoemHjwoK5aaC5p6c5b2T57uT5bGA5Y+Y6YeP5Ye95pWw5pe26Ze05Y+Y6YeP77yM5oiW6ICF6K+05piv5YiG57G75Y+Y6YeP55qE5pe25YCZLuWPr+aYr+S9v+eUqCoq5pyA5aSn6YCJ5oup56ep57uf6K6h6YePKirmnaXov5vooYzmnIDkvbNjdXRvZmbnmoTpgInmi6njgILlhbbkuLvopoHmgJ3mg7PkuZ/mmK/miormiYDmnInlj6/og73nmoTliIbnu4Tpg73orqHnrpfkuIDpgY3jgILnhLblkI7lr7vmib7mnIDkvbPnmoTnu5PmnpzjgILov5nnp43mlrnlvI/mnIDlpJrnlKjkuo7lnKgqKueUn+WtmOWIhuaekCoq55qE5pe25YCZ5a+75om+5pyA5L2z55qEY3V0b2ZmCi0gUuS4reWPr+S7pemAmui/h2BtYXhzdGF0YOadpeWunueOsOOAggoKYGBge3J9CmxpYnJhcnkobWF4c3RhdCkKc2V0LnNlZWQoMjkpCnggPC0gc29ydChydW5pZigyMCkpCnkgPC0gYyhybm9ybSgxMCksIHJub3JtKDEwLCAyKSkKbXlkYXRhIDwtIGRhdGEuZnJhbWUoY2JpbmQoeCx5KSkKaGVhZChteWRhdGEpCm1vZCA8LSBtYXhzdGF0LnRlc3QoeSB+IHgsIGRhdGE9bXlkYXRhLCBzbWV0aG9kPSJXaWxjb3hvbiIsIHBtZXRob2Q9IkhMIiwKICAgICAgICAgICAgICAgICAgICBtaW5wcm9wPTAuMjUsIG1heHByb3A9MC43NSwgYWxwaGE9MC4wNSkKcGxvdChtb2QpCmBgYAoKIyMjIyMg5L2/55So5YiG57G75qCR6L+b6KGM5YiS5YiGCgrkvb/nlKjliIbnsbvmoJHvvIznsbvkvLzkuo7lhrPnrZbmm7Lnur/kuIDmoLfnmoTmnaXor7TmmI7lpoLmnpzov5vooYzliIbnu4TjgIIKCmBgYHtyfQpsaWJyYXJ5KHJwYXJ0KQpmaXQgPC0gcnBhcnQoS3lwaG9zaXMgfiBBZ2UgKyBOdW1iZXIgKyBTdGFydCwgZGF0YSA9IGt5cGhvc2lzKQpwbG90KGZpdCkKdGV4dChmaXQsIHVzZS5uID0gVFJVRSkKYGBgCgojIyMjIyDogZrnsbvliIbmnpAKCuS5i+WJjeeahOaWueazlemDveacieS4gOS4quadoeS7tuWNs+W/hemhu+acieS4gOS4quaYjuehrueahOehruWumueahOe7k+WxgOOAgui/meagt+agueaNrue7k+WxgOWvueiHquWPmOmHj+i/m+ihjOWIkuWIhu+8jOmAmuW4uOWwhui/meenjeaDheWGteensOS4uioq5pyJ55uR552jKirnmoTjgILkvYbmmK/lpoLmnpzmiJHku6zmsqHmnInnu5PlsYDlj5jph4/nmoTml7blgJnvvIzlpoLmnpzopoHorrLmoLfmnKzov5vooYzliJLliIbnmoTor53vvIzov5nmoLflsLHpnIDopoHnlKjliLAqKuaXoOebkeedoyoq55qE6IGa57G75YiG5p6Q44CCCgotIOW4uOeUqOeahOaWueW8j+WPkeWxguasoeazleOAgkstbWVhbnPms5XnrYnnrYnjgII=