《应用商务统计分析》第三章 协方差分析

本章旨在通过教学评估数据分析的案例，介绍协方差分析统计回归模型的原理和应用。

目录

一、目的

二、数据来源和相关说明

三、数据清洗

四、描述性分析

五、模型分析

六、模型选择与预测

---

一、目的

1. 找出影响最终教学评估成绩的因素； 
2. 提出一个合理的绩效考核标准。

二、数据来源和相关说明

1.数据来源：北京大学光华管理学院的教学评估记录

import os
import pandas as pd
filePath = r'E:\CH3'
fileName = r'teaching.csv'
df_raw =  pd.read_csv(open(os.path.join(filePath, fileName)))

2.数据信息：时间范围：2002~2004年，数据量：340

print(df_raw.info())
print(df_raw.head())

3.变量信息

(1) 自变量/解释性变量：

• 教员职称(title)：助理教授、副教授、正教授教员；
• 性别(gender)：女、男；
• 学生类别(student)：本科生、MBA、研究生；
• 年份(year)：2002、2003、2004；
• 学期(semester)：春季、秋季；
• 学生人数(size)：3~136；（连续型）

# 离散变量
str_cols = ['title', 'gender', 'student', 'year', 'semester']
for col in str_cols:
    print(df_raw[col].value_counts().sort_index())
    print('-'*10)
# 连续变量
num_cols = ['size']
print(df_raw[num_cols].describe().T)

(2) 因变量：

• 课程得分(score)：2.56~5；（连续型）

y_col = 'score'
print(df_raw[[y_col]].describe().T)

三、数据清洗

原始数据已清洗好。

#原始数据已经清洗好，直接拷贝
df_clean = df_raw.copy() 

四、描述性分析

目的在于对数据进行初阶认识，并形成初步观点，但后续需对它们予以严格的分析和检验。

1.连续型变量

（1）自变量中只含有一个连续变量——班级人数(size)，首先通过散点图探索其与因变量间关系；

df_clean.plot(x='size', y='score', kind='scatter')

（2）观察散点图，发现数据噪音很大，数据规律不明显，因此对其进行了初步离散化降噪处理（20步长的等距分箱），生成了新变量——班级规模(group)；

df_clean['group'] = np.ceil(df_clean['size']/20) #等距分箱

（3）再次观察箱线图，以及各水平的样本量，对上述分箱划分点进行合并优化。

df_clean.boxplot(by='group', column='score', figsize=(8, 5))#箱线图

print(df_clean['group'].value_counts().sort_index()) #各水平数据量

#1和2之间差异明显，之后无太大差异，6和7由于数据小导致异常，因此最终以size=20为阈值进行2分箱
df_clean['group'] = df_clean['group'].map(lambda x: 1 if x<2 else 0)

2.离散型变量

通过箱线图探索以下离散型变量与因变量间的关系：

（1）教员职称(title): 随着教员职称提高（助教，副教授，正教授），平均成绩（中位数）依次增高。

（2）性别(gender): 性别间差异不大。

（3）学生类别(student): 研究生给出的平均成绩要高于MBA和本科生的，但MBA和本科生间的差异不大。

（4）年份(year): 成绩逐年增加。

（5）学期(semester): 春季学期的成绩要高于秋季。

（6）班级规模(group): 小班的成绩高于大班的（班级人数>20）。

for i in df_clean.columns:
    if i not in ['size', 'score']:
        df_clean.boxplot(by=i, column='score')

五、模型分析

df_model = df_clean.copy() #建模数据隔离

1.单因素可加模型

简单出发，先考虑班级规模(group)、班级人数(size)，假定模型为：

其中a0、a1为截距项，b0为斜率项，e为误差，a1为班级规模效应(表现在截距上)。

#建立模型
model = smf.ols('score ~ C(group) + size + 1', data=df_model).fit()
print(model.summary())
print('Residual standard error:  %.4f' % model.mse_resid**0.5)

#可视化预测结果
df_model['score_model'] = model.predict()
ax = df_model.plot(x='size', y='score', kind='scatter', legend='score')
df_model.plot(x='size', y='score_model', kind='scatter', legend='score_model', ax=ax, color='r')

2.单因素交互作用模型

考虑班级规模(group)和班级人数(size)的交互作用，假定模型为：

其中a0、a1为截距项，b0、b1为斜率项，e为误差，a1、b1为班级规模效应(表现在截距和斜率上)。

#建立模型
model = smf.ols('score ~ C(group) + size + C(group)*size + 1', data=df_model).fit()
print(model.summary())
print('Residual standard error:  %.4f' % model.mse_resid**0.5)

#可视化预测结果
df_model['score_model'] = model.predict()
ax = df_model.plot(x='size', y='score', kind='scatter', legend='score')
df_model.plot(x='size', y='score_model', kind='scatter', legend='score_model', ax=ax, color='r')

3.多因素协方差分析

（1）全模型

考虑所有自变量，以及班级规模(group)和班级人数(size)的交互作用，假定模型为：

教学评估成绩 = 教员职称 + 教员性别 + 学生类别 + 年份 + 学期 + 班级规模*学生人数 + e

#建立模型
model = smf.ols('score ~ C(title) + C(gender) + C(student) + C(year) + C(semester) + C(group)*size + 1', data=df_model).fit()
#方差分析
print(sm.stats.anova_lm(model, typ=3)) 
print('AIC: %.1f, BIC: %.1f' % (model.aic, model.bic))

（2）剔除不显著变量

• 手动剔除在10%水平下不显著的特征：教员性别(gender)、学期(semester)；
• 虽然size主效应估计量小，且不显著，但其与group交互显著，所以保留；

#手动剔除在10%水平下不显著的特征：gender、semester
model = smf.ols('score ~ C(title) +  C(student) + C(year) + C(group)*size + 1', data=df_model).fit()
#方差分析
print(sm.stats.anova_lm(model, typ=3)) 

print(model.summary())
print('Residual standard error:  %.4f' % model.mse_resid**0.5)

（3）模型诊断：确保模型结果可靠性

#获取残差
df_model['score_model'] = model.predict(df_model)
df_model['resid'] = model.resid
#残差：方差齐性
df_model.plot(x='score_model', y='resid', kind='scatter', grid=True)
#残差：正态性
sm.qqplot(df_model['resid'], line='s')
plt.show()

六、模型选择与预测

1.模型选择

使用向前逐步回归，基于AIC/BIC进行模型的自变量集选择。（指标越小越好）

#使用向前逐步回归（aic/bic为指标）进行自动选择
def step(baseFormula, data, scoreTyp='aic'):
    f = baseFormula.replace(' ', '')
    response = f.split('~')[0]
    remaining = [i for i in f.split('~')[1].split('+') if i != '1']
    
    selected = []
    current_best_score, best_score = np.PINF, np.PINF # 初始值

    k = 1
    while remaining and current_best_score == best_score:
        print('\n[INFO] No. %d'%k)
        k += 1
        scores_with_candidates = []
        for candidate in remaining:
            formula = "{} ~ {} + 1".format(response, ' + '.join(selected + [candidate]))
            score = eval('smf.ols(formula, data).fit().' + scoreTyp)
            scores_with_candidates.append((score, candidate))
            
        scores_with_candidates.sort()
        current_best_score, current_best_candidate = scores_with_candidates[0]
        print('Add %s, %.1f' % (current_best_candidate, current_best_score))
        
        if current_best_score <= best_score:
            remaining.remove(current_best_candidate)
            selected.append(current_best_candidate)
            best_score = current_best_score
        else:
            break    
      
    formula = "{} ~ {} + 1".format(response,' + '.join(selected))
    model = smf.ols(formula, data).fit()
    print('\n[INFO] FINAL')  
    print('%s, %.1f' % (formula, eval('model.' + scoreTyp)))
    
    return model
#基础模型
baseFormula = 'score ~ C(title) + C(gender) + C(student) + C(year) + C(semester) + C(group)*size + 1'
#基于AIC筛选的逐步回归
model = step(baseFormula, df_model, scoreTyp='aic')
print(sm.stats.anova_lm(model, typ=3)) #方差分析
print('AIC: %.1f, BIC: %.1f' % (model.aic, model.bic))

2.预测

df_model['score_model'] = model.predict(df_model)
df_model['score'].hist()
df_model['score_model'].hist()