GitModel 假设检验下|Task02多元数值向量的假设检验

多元数值向量的假设检验

#必要的包
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from IPython.display import display
1
2
3
4
5
6

3.1 均值向量的检验（多元数值数据中的组间均值向量假设检验）

如何通过样本数据对其所在多元总体的均值向量的性质进行检验

均值向量检验的内容架构与均值检验架构是一样的，主要内容:

1. 比较一组数据的总体均值向量与一个固定向量是否相等的检验
2. 比较两组数据的总体均值向量之间是否相等的检验
3. 比较两组以上的多组数据的总体均值向量之间是否相等的检验

在多元领域中，Hotelling T2检验，该检验在变量个数等于1时退化成t检验

3.1.1单组样本均值向量假定的检验

单组样本均值向量检验与单组样本的均值检验类似
假定的对象从一个数值变成一个向量（可以使用单样本hotelling T2检验）
$\mathbf{\mu }$是加粗的，代表向量

设$x_1,...,x_n$为p元正态分布的一个样本，对于单样本hotelling T2检验，做如下两个假设：
$$H_0:\mathbf{\mu }={\mathbf{\mu }}_0\leftrightarrow H_1:\mathbf{\mu }\ne {\mathbf{\mu }}_0$$

检验统计量为:
$$Teststatistics=\frac{n-p}{p\left(n-1\right)}{T}^2,T^2=n{\left(\bar{\mathbf{x}}-{\mathbf{\mu }}_0\right)}^{\prime }{\mathbf{S}}^{-1}\left(\bar{\mathbf{x}}-{\mathbf{\mu }}_0\right)$$
$\mathbf{S}$样本的无偏协方差矩阵，检验统计量服从F分布
$$Teststatistics\sim F(p,n-p)$$

p值计算公式为
$pvalue=P(F>Teststatistics)$

尽管备择假设的符号是≠，但这里p值并不是双边p值，多元检验中p值形式的判断并不能沿用一元均值检验的判断方法。

---

Exp5:我们想研究某地区农村2岁男婴的发育状况是否达到国家参考标准，
于是抽取了6名男婴并测量了他们的身高$x_1$、胸围$x_2$、上半臀围$x_3$。
已知这三个指标的国家参考标准均值为${\mathbf{\mu }}_0={\left(90,58,16\right)}^{\prime }$
6名男婴的数据如下所示，问能否认为该地区农村男婴的发育状况达到国家标准？

由于scipy包中并没有Hotelling T2的api，因此我们需要使用上述的检验统计量及其分布自己编写假设检验。

代码 hotelling T2检验

#引入数据
data=pd.read_excel('./data/ex5.xlsx')
#构造一个函数，帮助我们进行单样本Hotelling T2检验

def multi_checkmean(data:pd.DataFrame,mean:np;ndarray,confidence=0.0.5):
	'''
	data:待比较的矩阵，格式最好为pd.DataFrame()，以列为变量，以行为样本
	mean:假定的均值向量
	
	#计算检验统计量
	S=np.cov(data.T)
	S_inv=np.linalg.inv(S)
	x_bar=np.mean(data).values.T
	#样本的均值向量
	
	n=len(data)
	p=len(mean)
	T2=len(data)*np.dot((x_bar-mean).T,np.dot(S_inv,x_bar-mean))
	Test_statistics=T2*(len(data)-len(mean))/(len(mean)*(len(data)-1))
	
	#计算p值
	pvalue=f.sf(Test_statistics,len(mean),len(data)-len(mean))
	
	#输出样本的均值向量
	print('样本均值向量为:{}'.format(x_bar))
	print('目标均值向量为:{}'.format(mean))
	
	#比较p值和显著性水平

	if pvalue<confidence:
		print('在显著性水平{0:}下，样本均值向量不等于假定的均值向量。（p={1:.4f}）'.format(confidence,pvalue))
	else:
		print('在显著性水平(0:)下，样本均值向量等于假定的均值向量.(p={1:.4f})'.format(confidence.pvalue))
	return pvalue	

#use
mean=np.array([90,58,16]).T#均值向量
multi_checkmean(data, mean)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

3.1.2两组样本之间的均值向量相等性检验

与单变量的双样本均值检验一样，双样本均值向量的检验同样分为组别间独立与组别成对两种情况，使用的方法是一样的

组别间独立

ExP6

我们想知道3周岁的男女婴儿的身体指标是否存在显著差异，于是各抽取了男女各15名婴儿，并测量了他们的身高$X_1$与体重$X_2$，问：在显著性水平0.05下能否认为两者存在显著差异？

由于两个样本性别不同，因此我们可以认为两个样本相互独立。

设两个样本${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_{n1}$和${\mathbf{y}}_1,\cdots ,{\mathbf{y}}_{n2}$分别来自两个$p$元正态分布，我们做如下两个假设：

$$H_0:{\mathbf{\mu }}_1={\mathbf{\mu }}_2\leftrightarrow H_1:{\mathbf{\mu }}_1\ne {\mathbf{\mu }}_2$$

检验统计量为：
$$\begin{gathered} Teststatistics=\frac{n_1+n_2-p-1}{p\left(n_1+n_2-2\right)}{T}^2, \\ {T}^2=\frac{n_1{n}_2}{n_1+n_2}(\overline{\mathbf{x}}-\overline{\mathbf{y}}{)}^{{}^{\prime }}{\mathbf{S}}_p^{-1}(\overline{\mathbf{x}}-\overline{\mathbf{y}}), \\ {\mathbf{S}}_p=\frac{\left(n_1-1\right){\mathbf{S}}_1+\left(n_2-1\right){\mathbf{S}}_2}{n_1+n_2-2} \end{gathered}$$

检验统计量服从F分布：
$$Teststatistics\sim F\left(p,n_1+n_2-p-1\right)$$
p值形式为：
$$pvalue=P\left(F>Teststatistics\right)$$

#读取数据
data=pd.read_excel('./data/ex5.xlsx')
#将两个独立样本分开成两组
group_boy=data[data['gender']=='b'].drop('gender',axis=1)
group_girl=data[data['gender']=='g'].drop('gender', axis=1)

def multi_unpaired_data(group1:pd.DataFrame,group2:pd.DataFrame,confidence=0.05):
#计算检验统计量
	n1=len(group1)
	n2=len(group2)
	p=np.shape(group1)[1]#变量维度
	mean1=np.mean(group1).values.T
	mean2=np.mean(group2).values.T
	S1=np.cov(group1.T)
	S2=np.cov(group2.T)
	Sp=((n1-1)*S1+(n2=1)*S2)/(n1*n2-2)
	T2=n1*n2*(mean1-mean2).T@np.linalg.inv(Sp)@(mean1-mean2)
	#矩阵相乘a@b
	Test_statistic(n1+n2-p-1)*T2/(p*(n1+n2-2))
	
	#计算p值
	from scipy.stats import f
	pvalue=f.sf(Test_statistics,p,n1+n2-p-1)

	#比较p值和显著性水平	
	if pvalue<confidence:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(confidence,pvalue))
    else:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(confidence,pvalue))
    return pvalue
#use

multi_unpaired_data(group_boy,group_girl)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

成对组别

多元均值向量检验中，成对检验的原理与一元均值检验的成对检验原理是相似的，本质上做两个均值向量之差与零向量的之间的单样本均值向量检验。

设$\left({\mathbf{x}}_i,{\mathbf{y}}_i\right),i=1,\cdots ,n$为服从$p$元正态分布的成对数据，令他们相减得：
$${\mathbf{d}}_i={\mathbf{x}}_i-{\mathbf{y}}_i$$
检验统计量为：
$$Teststatistics=\frac{n-p}{p\left(n-1\right)}{T}^2,T^2=n{\bar{\mathbf{d}}}^{\prime }{\mathbf{S}}_{\mathbf{d}}^{-1}\bar{\mathbf{d}}$$
检验统计量服从分布：
$$Teststatistics\sim F\left(p,n-p\right)$$
p值计算公式为
$$pvalue=P\left(F>Teststatistics\right)$$

def multi_paried_data(group1:pd.DataFrame,group2:pd.DataFrame,confidence=0.05):
#计算检验统计量
n=len(group1)
p=np.shape(group1)[1]
mean1=np.mean(group1).values.T
mean2=np.mean(group2).values.T
d=group1-group2
d_mean=mean1-mean2
Sd=np.cov(d.T)
Sd_inv=np.linalg.inv(Sd)#求逆
T2=n*d_mean.T@Sd_inv@d_mean
Test_statistics=(n-p)*T2/(p*(n-1))

#计算p值
from scipy.stats import f
pvalue=f.sf(Test_statistics,p.n-p)

#比较p值与显著性水平
if pvalue<confidence:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(confidence,pvalue))
else:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(confidence,pvalue))
return pvalue
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

3.1.3 多元方差分析-多组样本间的均值向量相等性比较

多元方差分析与一元方差分析研究的内容是相似的，只是前者研究的是多个总体均值向量的相等性，而后者是均值的相等性。多元方差分析的检验统计量被称为威尔克斯(Wilks)$\Lambda$统计量，该统计量服从的分布为$\Lambda$分布。这个分布的分位点值可通过查F分布表得,转换复杂

statsmodels.multivariate.manova模块有多元方差分析(MANOVA)的api

Example.8

某医生测量了16个正常人的早晨中3个小时各小时的低频心电频值谱（LF）与高频心电频值谱（HF），数据在下面给出，问：这16个人的两个指标表现在3次测量中有无显著的差别？

我们简单对这个例子进行分析：由于问题考察的变量有两个：LF与HF，因此这是一个多元问题；又由于样本有三个：早晨的3次测量，因此这是一个多组均值比较问题。结合以上分析，该问题是一个（单因素）多元方差分析问题。

#读取数据
data=pd.read_excel('./data/ex8.xlsx')
group_1=data[data['Time']==1].drop('Time',axis=1).reset_index(drop=True)
group_2=data[data['Time']==2].drop('Time',axis=1).reset_index(drop=True)
group_3=data[data['Time']==3].drop('Time',axis=1).reset_index(drop=True)
from statsmodels.multivariable.manova import MANOVA
#多元方差分析(MANOVA)的api
model=MANOVA.from_formula('LF + HF ~ Time',data=data).mv_test()
#formula参数中填写计算公式 ~前自变量，，~因素变量
print(models.results['Times']['stat'])
1
2
3
4
5
6
7
8
9
10

3.2 从多元到一元——寻找显著的变量

• 多元均值向量的检验：不同分析对象的某个宏观性质（如城市企业的发展状况）是否存在显著区别，如果我们拒绝了原假设，便可以说明它们确实存在显著区别。
• 一元的均值向量检验：是什么变量导致了这种显著区别

一元数值与多元数值向量的常用假设检验在进行实际分析中的功能是互补的。

Example.9

为了研究3种销售方式对商品销售额的影响，我们选择了4种商品按照这3种销售方式进行销售。这四种商品的销售额分别为x1,x2,x3,x4。数据如下所示，请用多元/一元假设检验方法对该问题进行研究分析。

该例中有三个样本，因此我们会使用方差分析进行研究；销售额方面，本例用了四个商品的销售额去衡量，也就是说变量有4个，因此在分析的第一步我们可以使用多元方差分析，分析三种销售方式在“宏观”上是否会给商品销售额带来影响。

1. 多元方差分析

data=pd.read_excel('./data/ex9.xlsx')
# 先做多元方差分析
model=MANOVA.from_formula('x1 + x2 + x3 + x4 ~ sale', data=data).mv_test()
print(model.results['sale']['stat'])
1
2
3
4

p值为0.002807，显然三种销售方式的均值向量不全部相等。分析到这一步肯定是不够的，我们可以进一步思考：这三种销售方式的显著差异究竟是由哪些商品引起的呢？这个问题可以转化成——这三种销售方式会对哪种商品的销售产生显著的差异，为了研究这个问题，我们可以分别对变量x1,x2,x3,x4做一元方差分析。
2. 单变量方差检验

# 分别检验

## 先分组
group_1=data[data['sale']==1].drop('sale',axis=1).reset_index(drop=True)
group_2=data[data['sale']==2].drop('sale',axis=1).reset_index(drop=True)
group_3=data[data['sale']==3].drop('sale',axis=1).reset_index(drop=True)

## 分别做单变量方差检验
print(stats.f_oneway(group_1.x1.values,group_2.x1.values,group_3.x1.values))
print(stats.f_oneway(group_1.x2.values,group_2.x2.values,group_3.x2.values))
print(stats.f_oneway(group_1.x3.values,group_2.x3.values,group_3.x3.values))
print(stats.f_oneway(group_1.x4.values,group_2.x4.values,group_3.x4.values))

1
2
3
4
5
6
7
8
9
10
11
12
13

3. 剔除变量，再次进行多元方差分析
   x2x3的p值很大，因此我们可以认为这三种销售方式对于商品2/3而言不会产生显著的影响；x1、x4的p值小于0.05，且x4的p值尤其小，我们猜测，可能是商品4销售额均值的差异引起了均值向量的差异。对于该问题，我们可以剔除x4，对x1x3再做一次多元方差分析。

# 去掉变量x4，再做一次多元方差分析
model=MANOVA.from_formula('x1 + x2 + x3 ~ sale', data=data).mv_test()
print(model.results['sale']['stat'])
1
2
3

去除了x4后，多元方差分析检验不显著，因此说明对于商品1/2/3而言，三种销售方式的总体均值向量之间没有显著差异；此外，尽管对x1单独进行方差分析是显著的，可是x1x3一起做多元方差分析时却不显著，这说明商品1对三种销售方式的差异无明显影响

作业

为研究东、中、西部各省市规模以上的企业发展状况，我们收集了各城市企业的主要经济指标，包括：总资产贡献率、资产负债率、流动资产周转次数、工业成本费用利润率、产品销售率。我们用变量“类别”定义了各类城市，其中1为东部城市；2为中部城市；3为西部城市。数据文件为homework2.xlsx。假设显著性水平为𝛼，问：

对三个类别的城市进行均值向量间的两两比较，查看结果

对三个类别的城市同时进行均值向量间的比较，查看结果

承接问题2，你认为哪些变量导致了三个类别城市均值向量的差异？说出你的理由。

import pandas as pd
import numpy as np
from scipy.stats import f
import matplotlib.pyplot as plt
from IPython.display import display
data=pd.read_excel('./data/homework2.xlsx')
#数据分类
group_1=data[data['类别']==1].drop(['类别', '地区'], axis=1)#东部地区
group_2=data[data['类别']==2].drop(['类别', '地区'], axis=1)#中部地区
group_3=data[data['类别']==3].drop(['类别', '地区'], axis=1)#西部地区


############
# 对三个类别的城市进行均值向量间的两两比较，查看结果
def multi_unparied_data(group1:pd.DataFrame,group2:pd.DataFrame,confidence=0.05):
    # 计算检验统计量
    n1=len(group1)
    n2=len(group2)
    p=np.shape(group1)[1] # 变量维度
    mean1=np.mean(group1).values.T
    mean2=np.mean(group2).values.T
    S1=np.cov(group1.T)
    S2=np.cov(group2.T)
    Sp=((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)
    T2=n1*n2*(mean1-mean2).T@np.linalg.inv(Sp)@(mean1-mean2)/(n1+n2)
    Test_statistics=(n1+n2-p-1)*T2/(p*(n1+n2-2))

    # 计算p值
    pvalue=f.sf(Test_statistics,p,n1+n2-p-1)

    # 比较p值与显著性水平
    if pvalue<confidence:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(confidence,pvalue))
    else:
        print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(confidence,pvalue))
    return pvalue

#三个类别属于不同地区的城市，采用组间均值相等性检验
print('东部地区和中部地区的比较',multi_unparied_data(group_1,group_2))
print('西部地区和中部地区的比较',multi_unparied_data(group_2,group_3))
print('西部地区和东部地区的比较',multi_unparied_data(group_1,group_3))
#可以看出，东部和中部的均值向量相等，其他组别之间的都不相等

#对三个类别的城市同时进行均值向量间的比较，查看结果
#多个类别进行均值向量比较应该采用多元向量分析
from statsmodels.multivariate.manova import MANOVA
model=MANOVA.from_formula('总资产贡献率 + 资产负债率 + 流动资产周转次数 + 工业成本费用利润率 + 产品销售率 ~ 类别', data=data).mv_test()
print(model.results['类别']['stat'])

##承接问题2，你认为哪些变量导致了三个类别城市均值向量的差异？说出你的理由。
## 进行一元方差分析
#提取出所有的列表头
col=group_1.columns.tolist()
def onevar_analysis(g1:pd.DataFrame,g2:pd.DataFrame,g3:pd.DataFrame,col,confidence=0.05):
    for i in col:
        print('一元方差分析-列：{}'.format(i))
        ans=stats.f_oneway(g1[i].values,g2[i].values,g3[i].values)
        print(ans)
        pvalue=ans.pvalue
        if pvalue<confidence:
            print('在显著性水平{0:}下，两组样本所在总体的均值向量不相等。(p={1:.4f})'.format(confidence,pvalue))
        else :      
            print('在显著性水平{0:}下，两组样本所在总体的均值向量相等。(p={1:.4f})'.format(confidence,pvalue))
 
onevar_analysis(group_1,group_2,group_3,col)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

#移除最小的一列：流动资产周转次数
# col.remove('流动资产周转次数')
model=MANOVA.from_formula('总资产贡献率 + 资产负债率 + 工业成本费用利润率 + 产品销售率 ~ 类别', data=data).mv_test()
print(model.results['类别']['stat'])
#p=0.002995
#剔除第二小的产品销售率
model=MANOVA.from_formula('总资产贡献率 + 资产负债率 + 工业成本费用利润率  ~ 类别', data=data).mv_test()
print(model.results['类别']['stat'])
#p=0.012114
#剔除第三小的产品销售率
model=MANOVA.from_formula('资产负债率 + 工业成本费用利润率  ~ 类别', data=data).mv_test()
print(model.results['类别']['stat'])
#p=0.01364
1
2
3
4
5
6
7
8
9
10
11
12
13

可以认为资产负债率 + 工业成本费用利润率是在单变量方差分析中是显著的