Python进行常见的气象水文序列统计分析

不可否认的是python的统计功能若于R和SAS，但对于常见的统计分析，python亦可以实现。本文介绍Python中的ggplot2绘图库：plotnine，使用python完成常见的统计描述、分布差异检验、相关分析和回归分析方法。   【学习交流群：853991558】

使用plotnine库自带的mtcars汽车数据集：

选取mtcars的子集df，共32个记录，6个变量

将变量vs、am、gear由数值型连续变量转为字符型分类变量：

 变量分布：

plotnine绘图：

plotnine是python的一个绘图库，模仿了ggplot2的语法和绘图样式，如果熟悉R的ggplot2，那么该库可快速上手。

该库与R的ggplot2使用中主要有以下2点不同：

ggplot()的mapping参数，R中写为aes(x = mpg, y = wt)，plotnine写为aes(x="mpg", y="wt")。

plotnine要求整个绘图语句为一个语句，若中间需要换行，需要使用\连接或者语句首尾加括号。

以下是几个简单的例子：

散点图+回归线：

 

分组：

按照vs分组，分别绘制散点图、回归线及其95%置信区间。

 

分面（xkcd主题）：

按照变量vs和gear进行分面。

 

描述性统计量：

各变量的缺失值。

 连续变量和分类变量的分布信息。

 其他描述性统计量。

 正态分布样本均值的标准差（即标准误）及其95%置信区间。

统计学检验：

正态性检验：

p = 0.09，在0.05的显著性水平下接受原假设，即未发现变量wt不符合正态分布

两独立样本均值t检验：

p = 0.38，未发现方差不齐

 p = 0.001，两总体均值存在显著性差异。

两相关样本均值t检验：

关于独立样本还是相关样本，举个例子：

一项RCT（随机对照实验）想研究药物A和药物B对老年人的降压效果：研究分A组和B组两组，两组各用药14天，分别在基线（入组时）和14天时测血压值。对于A组或B组内部来说，想研究用药前后的血压值是否存在显著性差异，可使用两相关样本t检验（相当于比较前后差值的均值是否显著异于0）；如果想比较A组和B组的降压效果是否存在显著性差异，可先分别求得A组和B组血压前后的变化值，然后对二组的变化值进行独立样本t检验【其他方法如协方差分析，基线血压作为协变量】。

两独立样本分布的秩和检验：

 两相关样本分布的秩和检验：

方差分析：

先拟合回归模型，再对回归模型进行方差分析得到方差分析表
p = 0.0003，说明三组间均值存在显著性差异

 进一步进行多重比较，发现3组和4组、3组和5组间存在显著性差异

卡方检验：

方法使用选择：（n为总样本数，E为每个单元格的理论频数）。

当n >= 40，且E>=5：使用普通卡方检验；当p≈α（α为显著性水平）使用Fisher确切概率法；

当n >= 40，但1<E<5：使用校正卡方检验或者Fisher确切概率法；

当n < 40或者E <= 1：使用Fisher确切概率法。

 由于n=32<40，应直接使用Fisher确切概率法；这里使用校正卡方检验（correction=True）仅为显示其用法。

 Fisher确切概率法：p=0.47，未发现两变量分布存在显著性关联。

相关分析：

pearson相关：

 spearman相关：

 kendall相关：

回归分析：

多重线性回归：

lm.summary()返回大量信息，各变量回归系数为coef。如mpg的回归系数为-0.1191且p<0.05，意味着在保持cyl不变的情况下，mpg每减少1个单位，wt相应减少0.12 (95%CI,-0.18~-0.06)个单位。

 wt的回归拟合值：

含有定性变量：

定性变量写入C()内，reference指定参照水平。

如回归系数0.0738意味着，保持mpg不变的情况下，vs=1组的wt相比于vs=0组的wt的均值高出0.07个单位，但p值不显著（一般情况下先看p值，若p值显著大于设定的显著性水平如p=0.65，则不必再看回归系数）。

加入非线性变换：

对因变量进行变换，此处相当于先对wt取对数获得变量log_wt，再利用log_wt对其它变量回归，mpg的回归系数-0.0495可解释为：mpg每增加1个单位，wt减少 [e^−0.05−1]×100 % ，p值显著。

对自变量进行变换，此处相当于对mpg取平方后获得变量mpg2，再使用mpg2参与回归，回归方程可写为wt^=4.53−0.003×mpg2−0.10×vs。

Logistic回归：

Binomial指定因变量为二项分布，logit指定连接函数为logit: logit(p)=ln(p/(1−p ))。

OR值：

Logistic回归中各变量回归系数经自然指数转化后可解释为OR（Odds Ratio，比值比或优势比），当结局发生率较低时（<10%）可用作RR的估计值反映暴露变量效应大小：如am的回归系数为3.0073，意味着保持mpg不变时，am=1组的结局发生风险是am=0组的e^3.01=20.23倍，p值边缘不显著（p=0.06）。【由于vs的结局发生率40%多，此处使用Logistic回归估计OR值反应暴露变量效应是不合适的，此处可使用Log-binomial模型估计现患比PR】。

泊松回归：

泊松回归的因变量为单位时间或单位空间事件数的发生数，为计数变量（正整数）。

RR值：

泊松回归中各变量回归系数经自然指数转化后可解释为RR（Relative Risk，相对危险度），关联性研究中反映了暴露变量效应大小。如vs的回归系数为-0.0290，意味着保持mpg不变时，vs=1组的结局发生风险是vs= 0组的e^−0.03=0.98倍，但p值不显著（p=0.89）。