统计推断——假设检验——python代码检验两变量相关性，检验线性模型（置换法）

一、检验连变量相关性

简介

在全国家庭增长调查数据集中，新生儿体重和母亲年龄的相关性均为0.07，年龄较大的母亲似乎产下的孩子更重，但是，这种效应是偶然产生的吗？

检验方案

选择Pearson相关性作为检验统计量，此命题，我们使用双侧检验。

原假设：母亲年龄和新生儿体重之间没有相关性。

备择假设：母亲年龄和新生儿体重之间有相关性。

思路：假设母亲年龄和新生儿体重之间没有相关性，那么我们将“母亲年龄”这组数据打散1000次，每次“母亲年龄”的排列都不同，这1000次计算得来的Pearson相关系数应该服从均值为0的某一分布，如果当前情况的测量值偏离均值0太远，完全是小概率事件，则我们完全有理由认为，两变量之间的相关系数不为0，即两变量是有相关性的。

python代码实现

1.导入两变量

live,firsts,others=first.MakeFrames()  #产生数据
live=live.dropna(subset=['agepreg','totalwgt_lb'])  #删除两变量为空值的行
data=live['agepreg'].values,live['totalwgt_lb'].values  #转换为list
xs,ys=data

2、定义函数（协方差，相关系数等）

import collections
import math
#将两个样本中其中一个样本的元素打散。
def RunModel3(xs,ys): 
    xs=np.random.permutation(xs)#与shuffle类似，只是permutation不改变原来的数据结构
    data=xs,ys
    return data

#计算协方差
def Cov(xs,ys):
    xs=np.asarray(xs) #将其他数据类型转换为array数组
    ys=np.asarray(ys)
    mean_x=xs.mean()
    mean_y=ys.mean()
    cov=np.dot(xs-mean_x,ys-mean_y)/len(xs) #计算两数组的內积
    return cov

#计算Pearson相关性
def Corr(xs,ys):
    xs=np.asarray(xs) #将其他数据类型转换为array数组
    ys=np.asarray(ys)
    var_x=xs.var()
    var_y=ys.var()
    corr=Cov(xs,ys)/math.sqrt(var_x*var_y)
    return corr

#计算每次重新划分之后的相关系数的绝对值。
def TestStatistic3(data):
    xs,ys=data
    test_stat=abs(Corr(xs,ys))
    return test_stat

3、计算p值

a=[TestStatistic3(RunModel3(xs,ys)) for _ in range(1000)] #模拟1000次试验
sorted(a,reverse=True)

输出：

b=sum(1 for x in a if x>=Corr(xs,ys)) #计算比当前情况乃至更差情况出现的次数
b/1000

输出：

0.0

实际数据的相关性为0.07。检验计算得到的p值为0，在1000次重复试验中，模拟得到的最大相关性为0.036。因此，虽然观擦到的变量相关性很小，但这种相关性是显著的。

二、检验线性模型

拟合优度

要度量一个线性模型的质量优劣，或拟合优度，最简单的方法之一是计算决定系数，决定系数（）和Pearson（）相关系数之间有一个简单的关系：，我们计算得到。

我们已经知道母亲年龄和新生儿体重的影响很小，几乎没有预测能力。那么，这个明显的关系可能是偶然产生的吗？我们可以使用一下以下方法检验线性拟合的结果。

1、检验方法一

检验决定系数是否为0。

原假设为=0，上面，我们已经对母亲年龄和新生儿体重的相关系数，由于，的单侧检验等效于的双侧检验。我们已经进行了这项检验，得到p<0.001，因此认为母亲年龄和新生儿体重之间的关系是统计显著的。

2、检验方法二

检验斜率是否偶然。

原假设是斜率为0，因此，我们可以使用均值附近的随机变化（残差）建立新生儿体重模型，我们先普及回归方程斜率和截距的公式：

#导入样本
live,firsts,others=first.MakeFrames()
live=live.dropna(subset=['agepreg','totalwgt_lb'])
ages=live['agepreg']
weights=live['totalwgt_lb']


import collections
#根据样本X、Y，计算线性方程的斜率与截距
def LeastSquares(xs,ys):
    mean_x,var_x=xs.mean(),xs.var()  #X的样本的均值、方差
    mean_y=ys.mean()   #Y的均值
    cov=np.dot(xs-mean_x,ys-mean_y)/len(xs) #协方差
    slope=cov/var_x   #斜率
    inter=mean_y-slope*mean_x   #截距
    return inter,slope

#假设斜率为0，则所有的y轴应该都等于y的均值加上偏差。
def RunModel(ys): 
    ybar=ys.mean()   #计算因变量y的均值
    res=ys-ybar  #计算体重和均值的差，作为残差
    weight=ybar+np.random.permutation(res)  #随机残差，与shuffle类似
    return weight

#计算每次重新划分之后的斜率。
def TestStatistic(xs,ys):
    _,slope=LeastSquares(xs,ys)
    return slope

计算1000次对残差随机抽样的斜率结果：

a=[TestStatistic(ages,RunModel(weights)) for _ in range(1000)]
sorted(a,reverse=True)

输出：

其中在斜率为0的原假设下，我们通过抽样随机残差能够获得的斜率最大值为0.0092。

b=sum(1 for x in a if x>=LeastSquares(ages,weights)[1]) #计算比当前情况乃至更差情况出现的次数
b/1000

输出：

0.0

当前数据情况算得的斜率为0.017，原假设条件下，得到当前情况的概率为0，故认为斜率不为0是显著的。

3、检验方法三

绘制原假设（斜率为0）成立的条件下，和抽样分布下，斜率的累积密度分布。

何为抽样分布？

我们将观测到的妊娠数据看成一个整体，从这个整体当中进行有放回的重复抽样。

1000次抽样分布斜率的函数如下：首先需要对原数据进行重复抽样。

#有放回的从df中重新选择行，组成新的df,新df和原df行数相同。
def ResampleRows(df):
    indices=np.random.choice(df.index,len(df),replace=True)
    return df.loc[indices]

#计算斜率和截距
def SampleDistributions(df,iters=1000):   
    t=[]
    for _ in range(iters):
        sample=ResampleRows(df)
        ages=sample['agepreg']
        weights=sample['totalwgt_lb']
        estimates=LeastSquares(ages,weights)
        t.append(estimates)
    inters,slopes=zip(*t)   #t为一个列表，*表示将列表t展开。
    return inters,slopes

#计算1000次抽样的斜率和截距，以列表的形式
_,slopes=SampleDistributions(live)

1000次原解释的斜率分布请参考方法二当中

a=[TestStatistic(ages,RunModel(weights)) for _ in range(1000)]
sorted(a,reverse=True)

计算两种分布下，斜率的累积概率密度

cdf=Cdf(Pmf(Hist(slopes)))   #抽样分布下斜率的累积概率分布
slope,cdf_slope=zip(*cdf.items())   

cdf_h=Cdf(Pmf(Hist(a)))     #原假设条件下的累积概率分布
slope_h,cdf_slope_h=zip(*cdf_h.items())    

画图：

#画图
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
plt.figure(figsize=(6,4.5)) #表示绘制图形的画板尺寸为6*4.5；
plt.plot(slope,cdf_slope,color ='g',label='抽样分布')   
plt.vlines(Pecentile(cdf,0.5),ymin=0,ymax=1) #辅助线，中位值

plt.plot(slope_h,cdf_slope_h,color ='b',label='原假设斜率为0')
plt.vlines(Pecentile(cdf_h,0.5),ymin=0,ymax=1) #辅助线，中位值

plt.legend(loc='upper left')  #添加图例
plt.xlim(-0.03,0.03) 
plt.ylim(0,1) 
plt.xlabel('斜率（磅/年）')
plt.ylabel('CDF')
plt.show()

 输出：

根据以上图片，我们可以用两种方法估计p值：  
1、计算原假设下斜率大于观测斜率的概率。  
2、计算抽样分布中斜率小于0的概率（如果估计斜率为负数，那么应该计算抽样分布中斜率大于0的概率）

由图我们可以看到，抽样分布（绿线）斜率小于0的概率为0，图中不相交。故原假设不成立，斜率不为0是显著的。