本作业内容来自Datawhale和GitModel，学习链接：https://github.com/Git-Model/Modeling-Universe/tree/main/Data-Story
感谢课程开发者的付出与贡献！
1 问题描述
为研究东、中、西部各省市规模以上的企业发展状况，我们收集了各城市企业的主要经济指标，包括：总资产贡献率、资产负债率、流动资产周转次数、工业成本费用利润率、产品销售率。我们用变量“类别”定义了各类城市，其中1为东部城市；2为中部城市；3为西部城市。数据文件为homework2.xlsx。假设显著性水平为𝛼，问：
- 对三个类别的城市进行均值向量间的两两比较，查看结果
- 对三个类别的城市同时进行均值向量间的比较，查看结果
- 承接问题2，你认为哪些变量导致了三个类别城市均值向量的差异？说出你的理由
2 问题解决
# 加载必要的包
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from IPython.display import display
# 加载数据
data = pd.read_excel('./data/homework2.xlsx')
data
group1 = data[data['类别']==1].drop(['类别','地区'],axis=1).reset_index(drop=True)
group2 = data[data['类别']==2].drop(['类别','地区'],axis=1).reset_index(drop=True)
group3 = data[data['类别']==3].drop(['类别','地区'],axis=1).reset_index(drop=True)
group1
# Q1：对三个类别的城市进行均值向量间的两两比较，查看结果
# 分析思路：多元数值数据 > 双样本均值向量的检验 > 组别间独立（分为东、中、西，从地理位置上看独立）
def multi_unparied_data(group1:pd.DataFrame,group2:pd.DataFrame,confidence=0.01):
# 计算检验统计量
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值
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 None
print("东部地区与中部地区的两组样本之间的均值向量相等性检验结果为：")
multi_unparied_data(group1,group2)
print("--------------------")
print("中部地区与西部地区的两组样本之间的均值向量相等性检验结果为：")
multi_unparied_data(group2,group3)
print("--------------------")
print("东部地区与西部地区的两组样本之间的均值向量相等性检验结果为：")
multi_unparied_data(group1,group3)
东部地区与中部地区的两组样本之间的均值向量相等性检验结果为：
在显著性水平0.01下，两组样本所在总体的均值向量相等。(p=0.2793) 企业发展相似。
--------------------
中部地区与西部地区的两组样本之间的均值向量相等性检验结果为：
在显著性水平0.01下，两组样本所在总体的均值向量相等。(p=0.0470) 企业发展相似。
--------------------
东部地区与西部地区的两组样本之间的均值向量相等性检验结果为：
在显著性水平0.01下，两组样本所在总体的均值向量不相等。(p=0.0097) 企业发展不相似。
# Q2：对三个类别的城市同时进行均值向量间的比较，查看结果
# 分析思路：多元数值数据 > 多元方差分析
from statsmodels.multivariate.manova import MANOVA
model=MANOVA.from_formula('总资产贡献率 + 资产负债率 + 流动资产周转次数 + 工业成本费用利润率 + 产品销售率 ~ 类别', data=data).mv_test()
# 在''中填入公式，其中~左侧填入自变量名称，~右侧填入因素名称
print(model.results['类别']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.535029 5 25.0 4.345284 0.005542
Pillai's trace 0.464971 5.0 25.0 4.345284 0.005542
Hotelling-Lawley trace 0.869057 5 25.0 4.345284 0.005542
Roy's greatest root 0.869057 5 25 4.345284 0.005542
因为p<0.01，所以接受原假设，不同类别的城市企业发展状况相似
# Q3：承接问题2，你认为哪些变量导致了三个类别城市均值向量的差异？说出你的理由
# 这部分作业主要参照参考链接中的代码来完成，自己之前写的内容有很多纰漏
data = pd.read_excel('./data/homework2.xlsx')
data.drop(columns='地区',inplace=True)
CityDict = {'T': '总资产贡献率', # Total Asset Contribution Rate
'A': '资产负债率', # Assets And Liabilities
'C': '流动资产周转次数', # Current Asset Turnover
'I': '工业成本费用利润率', # Industrial CostExpense Profit Margin
'P': '产品销售率' # Product Sales Rate
}
columns = ['category', 'T', 'A', 'C', 'I', 'P']
data.columns = columns
group1 = data[data['category']==1].drop('category', axis=1).reset_index(drop=True)
group2 = data[data['category']==2].drop('category', axis=1).reset_index(drop=True)
group3 = data[data['category']==3].drop('category', axis=1).reset_index(drop=True)
# 注意要先做方差齐性检验，再根据置信水平，选择相对应的一元方差分析（单因素方差分析）方法
for key in ['T', 'A', 'C', 'I', 'P']:
_, levene = stats.levene(group1[key].values, group2[key].values, group3[key].values)
if levene < 0.01:
print('-------------------------------')
print('警告: 方差齐性检验的p值小于0.01: p={}，方差分析结果在小样本下可能不准确'.format(round(levene, 4)))
print('三个类别城市企业的经济指标 ', CityDict[key], ' 样本不服从正态性样本')
# 单因素方差分析-假设样本不服从正态分布
print(stats.mstats.kruskalwallis(group1[key].values, group2[key].values, group3[key].values))
print('-------------------------------')
else:
print('三个类别城市企业的经济指标 ', CityDict[key], ' 样本服从正态性样本')
# 单因素方差分析-假设样本服从正态分布
print(stats.f_oneway(group1[key].values, group2[key].values, group3[key].values))
三个类别城市企业的经济指标 总资产贡献率 样本服从正态性样本
F_onewayResult(statistic=3.5776324463321805, pvalue=0.04133778417613981)
三个类别城市企业的经济指标 资产负债率 样本服从正态性样本
F_onewayResult(statistic=0.7703999845168195, pvalue=0.47239023585585066)
三个类别城市企业的经济指标 流动资产周转次数 样本服从正态性样本
F_onewayResult(statistic=6.484356111904944, pvalue=0.0048515449323179)
-------------------------------
警告: 方差齐性检验的p值小于0.01: p=0.0011，方差分析结果在小样本下可能不准确
三个类别城市企业的经济指标 工业成本费用利润率 样本不服从正态性样本
KruskalResult(statistic=1.8717467008797541, pvalue=0.39224315017025624)
-------------------------------
三个类别城市企业的经济指标 产品销售率 样本服从正态性样本
F_onewayResult(statistic=4.847034218142545, pvalue=0.015573643638731122)
可以看出资产负债率和工业成本费用利润率的p值较高，可以认为这两种指标不会对企业发展状况产生显著影响；
总资产贡献率和资产和产品销售率的p值略大于0.01；
流动资产周转次数的p值远小于0.01，因此设想流转资产周转次数的因素影响较大，剔除该因素，进行多元方差分析
# 去掉流转资产周转次数
new_model = MANOVA.from_formula('T + A + I + P ~ category', data=data).mv_test()
print(new_model.results['category']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.551704 4 26.0 5.281678 0.002995
Pillai's trace 0.448296 4.0 26.0 5.281678 0.002995
Hotelling-Lawley trace 0.812566 4 26.0 5.281678 0.002995
Roy's greatest root 0.812566 4 26 5.281678 0.002995
p值小于0.01，无法证明流动资产周转次数的差异引起均值向量的最大差异，去掉多个变量进行分析
# 去掉流转资产周转次数和总资产贡献率
new_model = MANOVA.from_formula('A + I + P ~ category', data=data).mv_test()
print(new_model.results['category']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.590459 3 27.0 6.242381 0.002327
Pillai's trace 0.409541 3.0 27.0 6.242381 0.002327
Hotelling-Lawley trace 0.693598 3 27.0 6.242381 0.002327
Roy's greatest root 0.693598 3 27 6.242381 0.002327
# 去掉流转资产周转次数和产品销售率
new_model = MANOVA.from_formula('T + A + I ~ category', data=data).mv_test()
print(new_model.results['category']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.671812 3 27.0 4.396615 0.012114
Pillai's trace 0.328188 3.0 27.0 4.396615 0.012114
Hotelling-Lawley trace 0.488513 3 27.0 4.396615 0.012114
Roy's greatest root 0.488513 3 27 4.396615 0.012114
# 去掉总资产贡献率和产品销售率
new_model = MANOVA.from_formula('A + C + I ~ category', data=data).mv_test()
print(new_model.results['category']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.715144 3 27.0 3.584874 0.026592
Pillai's trace 0.284856 3.0 27.0 3.584874 0.026592
Hotelling-Lawley trace 0.398319 3 27.0 3.584874 0.026592
Roy's greatest root 0.398319 3 27 3.584874 0.026592
# 去掉总资产贡献率、流转资产周转次数和产品销售率
new_model = MANOVA.from_formula('A + I ~ category', data=data).mv_test()
print(new_model.results['category']['stat'])
Value Num DF Den DF F Value Pr > F
Wilks' lambda 0.735822 2 28.0 5.026337 0.01364
Pillai's trace 0.264178 2.0 28.0 5.026337 0.01364
Hotelling-Lawley trace 0.359024 2 28.0 5.026337 0.01364
Roy's greatest root 0.359024 2 28 5.026337 0.01364
第三问按照小林博主的代码把结果做出来了，但是不知道要怎么分析结果……
参考链接：小林博主