背景
线性回归，主成分回归都做烂了，我之前的案例有很多这些模型，但是一直没写因子分析的回归案例，这个也是传统统计学流行的方法，在金融经济心理学等人文社科用得非常多。这个案例就演示一下python怎么做因子分析。
数据介绍
很早之前的财经新闻的数据了，问卷获取的，SAV格式，要用spss软件打开，但是python肯定也可以读取。
我们直接用python读取SAV格式的数据会，这个数据之前案例是做过线性回归，主成分，随机森林回归。
Python数据分析案例22——财经新闻可信度分析(线性回归，主成分回归，随机森林回归)
为什么本文还是用这个数据？因为这个数据的变量之间有很强的相关性，具有多重共线性，所以很适合用主成分，因子分析等方法。本文就做一点描述性统计分析，因子分析，然后进行因子分析回归。
需要本案例数据和全部代码文件的同学还是可以参考：因子分析
代码实现
首先导入包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
import statsmodels.formula.api as smf
plt.rcParams ['font.sans-serif'] ='SimHei' #显示中文
plt.rcParams ['axes.unicode_minus']=False #显示负号
#sns.set_style("darkgrid",{"font.sans-serif":['KaiTi', 'Arial']})读取数据
# 读取数据清洗后的数据
spss = pd.read_spss('数据2.sav')
#spss选取所需要的变量
data=spss[['微博平台可信','专业性','可信赖性','转发量','微博内容质量','时效性','验证程度','人际信任','投资信息可信度']]
data.head()
做点描述性统计
data.describe() 
取出X和y，y就是投资信息可靠度，也就是数据框的最后一列。
X=data.iloc[:,:-1]
y=data.iloc[:,-1]
可视化
画出变量的核密度图
column = data.columns.tolist() # 列表头
fig = plt.figure(figsize=(8,8), dpi=128) # 指定绘图对象宽度和高度
for i in range(9):
plt.subplot(3,3, i + 1) # 2行3列子图
sns.kdeplot(data=data[column[i]],color='blue',fill= True)
plt.ylabel(column[i], fontsize=16)
plt.tight_layout()
plt.show()
都是类似正态分布，没有极端的异常值。
变量两两之前的散点图
sns.pairplot(data[column],diag_kind='kde')
可以看到变量两两之间都有较强的相关性，基本都是正相关。
画出变量之间的相关性热力图：
#画皮尔逊相关系数热力图
corr = plt.subplots(figsize = (10,10),dpi=128)
corr= sns.heatmap(data[column].corr(),annot=True,square=True)
可以看到越红相关性越高，都很红。
统计学学得还不错的同学可能会说，计算相关性数值没用啊，因为有的相关性系数不一定是显著的，是的，所以我自定义了一个函数，可以计算相关性系数和其显著性。
def calculate_correlation_with_significance(dataframe):
significance_levels = [(0.001, '***'), (0.01, '**'), (0.05, '*'), (1.0, '')]
# 生成一个相关系数和显著性符号的数据框
result_df = pd.DataFrame({
var1: {
var2: ( f"{1.0} " if var1 == var2 else
f" {next(sign for threshold, sign in significance_levels if p < threshold)}{r:.3f}"
)
for var2 in dataframe.columns
for r, p in [stats.pearsonr(dataframe[var1], dataframe[var2])]
}
for var1 in dataframe.columns
})
# 仅保留下三角数据
mask = np.tril(np.ones(result_df.shape), k=0).astype(bool)
result_df = result_df.where(mask, '')
return result_df
merged_df = calculate_correlation_with_significance(data)
merged_df
可以看到，很多X之间都存在很显著很高的相关性，直接线性回归的话可能会受到多重共线性的影响。 而这种X变量之间存在高度相关性的数据适合使用因子分析。
线性回归
我们做因子分析前肯定要做线性回归进行对比
all_columns = "+".join(data.columns[:-1])
print('x是：'+all_columns)
formula = '投资信息可信度~' + all_columns
print('回归方程为：'+formula)
results = smf.ols(formula, data=data).fit()
results.summary()
查看主要的几个回归系数
#：
print(results.summary().tables[1])
可以看到基本都不显著，可能是多重共线性导致的，计算VIF看看
def VIF_calculate(df_all,y_name):
x_cols=df_all.columns.to_list()
x_cols.remove(y_name)
def vif(df_exog,exog_name):
exog_use = list(df_exog.columns)
exog_use.remove(exog_name)
model=smf.ols(f"{exog_name}~{'+'.join(list(exog_use))}",data=df_exog).fit()
return 1./(1.-model.rsquared)
df_vif=pd.DataFrame()
for x in x_cols:
df_vif.loc['VIF',x]=vif(df_all[x_cols],x)
df_vif.loc['tolerance']=1/df_vif.loc['VIF']
df_vif=df_vif.T.sort_values('VIF',ascending=False)
df_vif.loc['mean_vif']=df_vif.mean()
return df_vif
VIF_calculate(data,'投资信息可信度')
VIF虽然没到10，但是每个变量都有5-8多，也很高了，说明模型存在中等程度的多重共线性，这种X变量之间存在高度相关性的数据适合使用因子分析。
因子分析
from sklearn.decomposition import FastICA
from sklearn.decomposition import FactorAnalysis首先我们要选择最优的因子数量
fa = FactorAnalysis(n_components=len(X.columns))
fa.fit(X)
loadings= fa.components_
cumulative_variance_ratio = np.cumsum(np.sum(loadings**2, axis=1) /len(X.columns))
print(cumulative_variance_ratio)
plt.plot(range(1, len(cumulative_variance_ratio) + 1), cumulative_variance_ratio, marker='o')
plt.title('Scree Plot')
plt.xlabel('Factor Number')
plt.ylabel('cumulative_variance_ratio')
plt.show()
就选3个因子吧，3个因子就有基本能解释原数据的所有的方差了。
# 初始化因子分析，设定提取的因子数
n_factors = 3
fa = FactorAnalysis(n_components=n_factors)
factors = fa.fit_transform(X)
# 获取因子载荷矩阵
loadings = fa.components_
print("因子载荷矩阵形状:", loadings.shape)
3就是因子的数量，8就是原变量的变量数量。
# 计算共性方差
communalities = np.sum(loadings.T**2, axis=1)
print("共性方差:", communalities)
共性方差可以通过每个变量的因子载荷的平方和来近似计算。共性方差反映了每个变量的方差中可以被提取的因子解释的部分。
eigenvalues = np.sum(loadings**2, axis=1)
eigenvalues # 因子的特征值，表示每个因子可以解释的方差量
# 计算方差贡献率
variance_ratio = eigenvalues / X.shape[1]
variance_ratio
## 累计计算方差贡献率
cumulative_variance_ratio = np.cumsum(variance_ratio)
cumulative_variance_ratio
因子和原始变量的关系
转化的因子就是原始的变量的线性组合，其系数都在载荷矩阵
# 将载荷矩阵转换为DataFrame以便查看
loadings_df = pd.DataFrame(loadings, index=[f'Factor{i+1}' for i in range(n_factors)],columns=X.columns)
loadings_df
# Visualize FA loadings
fig, ax = plt.subplots(3, 1,figsize=(6,6),dpi=128)
plt.subplots_adjust(hspace=1, wspace=0.5)
for i in range(1, 4):
ax = plt.subplot(3, 1, i)
ax.plot(loadings_df.T['Factor' + str(i)], 'o-')
ax.axhline(0, color='k', linestyle='--', linewidth=1)
ax.set_xticks(range(len(X.columns)))
ax.set_xticklabels(X.columns, rotation=10)
ax.set_title('Factor Loadings for FA' + str(i))
plt.tight_layout()
plt.show()
plt.figure(figsize=(7, 2.4),dpi=128)
# 绘制热力图
sns.heatmap(loadings_df, annot=True, fmt=".3f", cmap="viridis", cbar=True, xticklabels=loadings_df.columns, yticklabels=loadings_df.index)
# 显示图表
plt.title("因子载荷矩阵热力图")
plt.xlabel("变量")
plt.ylabel("因子")
plt.show()
因子含义分析 在因子分析中，每个因子通常代表原始变量的一组线性组合，反映共同的特质或维度。根据你的因子载荷矩阵，我们可以尝试推测和解释每个因子的潜在含义：
Factor 1: 具有较高负载的变量有：微博平台可信、专业性、可信赖性、转发量、微博内容质量、时效性、验证程度、人际信任。 因为大多数负载都为负值且较高，Factor 1可能反映了一个综合性的信任或可靠性维度，表明这些特征共同减少时，可能降低某种信任感。 简单来说就是其他几个变量取值越高，Factor1越低。Factor1这个值越大，说明了不专业不可靠 ，是一个描述不专业程度的因子。
Factor 2: 在转发量、微博内容质量、时效性上有较高的正负载。 Factor 2可能与信息传播与有效性相关，表明其可能与信息的快速传播和内容质量提升有关。
Factor 3: 变量的负载较小，且不明显。 Factor 3的含义可能较弱或特定于某些细微特征，可能代表难以解释的随机效应或特定条件下的特性。
因子回归
我们对转化的因子作为X，进行因子分析回归：
X_factors = fa.fit_transform(X) # 生成因子得分
# 将因子得分转换为DataFrame
X_factors = pd.DataFrame(X_factors, columns=[f'Factor{i+1}' for i in range(n_factors)])
X_factors.shape 
all_columns = "+".join(X_factors.columns)
print('x是：'+all_columns)
formula = '投资信息可信度~' + all_columns
print('回归方程为：'+formula)
results = smf.ols(formula, data=pd.concat([X_factors,y],axis=1)).fit()
results.summary()
print(results.summary().tables[1])
回归分析结果为OLS（最小二乘法）回归，解释了因子对于投资信息可信度的影响：
模型整体性能： R-squared: 0.805 和 Adj. R-squared: 0.777：模型能够解释投资信息可信度80.5%的方差，调整后的R平方为77.7%，表明模型拟合良好。 F-statistic: 28.82 且 Prob (F-statistic): 1.23e-07：F统计量显著，表示整体模型在统计上是显著的。
各因子的回归系数： Intercept（截距）：2.9067，这意味着所有因子为零时，投资信息可信度的基线值大约为2.91。 Factor 1：系数为-0.8143（p值 < 0.001），是高度显著的，表明Factor 1对投资信息可信度有显著的负面影响。这与我们推测的信任或可靠性维度一致，可能因为信任特征的减少而降低信息可信度。 Factor 2：系数为0.3099（p值 = 0.005），表示Factor 2对可信度有显著的正面影响，支持其与信息传播有效性相关的假设。 Factor 3：系数为-0.0185（p值 = 0.884），无显著性，因而对可信度的影响可以忽略。这暗示Factor 3可能不重要或与投资信息可信度无关。
其他统计量： Durbin-Watson: 1.417：接近2，表明残差自相关性较低。 Omnibus 和 Jarque-Bera (JB) Tests：用于检验残差的正态性，结果显示残差分布接近正态。
总结 综合来看，Factor 1显著地降低了投资信息的可信度，可能是因为与信任和专业相关的特征减少。Factor 2则有助于提升信息可信度，与信息传播和质量相关。Factor 3对可信度没有显著影响，可能与随机或其他未观测到的特征有关。整体模型对投资信息可信度的解释能力相当强。
所以当普通的ols因为高度相关性导致的多重共线性而不显著的时候，我们可以使用因子分析，找到有用的因子去进行回归，就可能显著了，而且也具有一定的解释性。
以前的文章可以参考：Python数据分析案例_阡之尘埃的博客
创作不易，看官觉得写得还不错的话点个关注和赞吧，本人会持续更新python数据分析领域的代码文章~(需要定制类似的代码可私信)
