回归分析与模型诊断——作业

作业1：含二次项/对数项模型的讨论

我们想要探究婴儿出生的体重与何种因素相关，数据集为bwght2.dta，本次习题所使用的变量解释如下：
因变量：
· bwght：婴儿出生体重

自变量：
· npvis：母亲产前检查次数
· mage：母亲年龄

使用python进行实操并回答以下问题

（1）：使用OLS估计方程
$$\log (bwght)={\beta }_0+{\beta }_{1}npvis+{\beta }_{2}npvi{s}^2+u$$
输出报告表，并回答：自变量npvis的二次项是否显著？自变量npvis是否对因变量有显著影响？

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
from IPython.display import display
import statsmodels.api as sm

# 加载数据
bwght2=pd.read_stata('./bwght2.dta')
bwght2.head()

bwght2_lm=sm.formula.ols('np.log(bwght)~npvis + I(npvis**2)',data=bwght2).fit()
print(bwght2_lm.summary())

自变量npvis的二次项显著，自变量npvis对因变量有显著影响。

（2）：基于（1）的方程，我们认为最大化log(bwght)的产前检查次数npvis约为24，其理论依据是什么？
对（1）的npvis进行求导，并令导数等于0，可得
$${\beta }_1+2{\beta }_{2}npvis=0$$
代入可得
$$npvis=23.625$$
所以最大化log(bwght)的产前检查次数npvis约为24。

（3）：按照这个模型的结果，在24次产前检查后婴儿出生体重会下降，这是为什么？你认为这有实际意义吗？这蕴含了一个含二次项变量模型的常见陷阱，请仔细思考！

这种结果是因为遗漏变量导致的。注意到R方只有0.021，模型的拟合效果很差，应该是有很多相关变量没有被纳入。可能存在某种变量，同时影响产检次数和婴儿出生体重。

（4）：在模型中加入母亲年龄变量及其二次形式。回答：保持npvis不变，母亲在什么生育年龄时，孩子出生体重最大？大于这个年龄时，孩子出生体重下降，这是否具有实际意义呢？请结合问题（3）思考这一问题。

bwght2_lm2=sm.formula.ols('np.log(bwght)~npvis + I(npvis**2) + mage + I(mage**2)',data=bwght2).fit()
print(bwght2_lm2.summary())

保持npvis不变，母亲在生育年龄47.46（0.0009/(9.481e-06*2)）时，孩子出生体重最大。这并没有实际意义。因为这个模型并不正确，特别是mage的一次项、二次项的系数都不显著，并且模型的R方也很小。说明有很多的遗漏变量，影响到所考虑的变量。

（5）：（4）中的模型能否解释log(gwght)大部分变异？
不能，因为R方非常小，说明了这一点。

作业2：异方差模型的讨论

紧接着example13（主管人数与工人人数间关系的问题），我们对异方差问题进行更进一步的探讨。

使用python进行实操并回答以下问题

（1）：在假设方差形式为$\mathrm{Var}(u\mid x)={\sigma }^2{x}^2$并进行wls估计后，比较wls估计与ols估计的残差图，回答：异方差消除了吗？

# 载入数据集
data=pd.read_table('./P176.txt')
data.head()

# 查看ols估计的残差与X的散点图
fig=plt.figure(figsize=(13,6))
ax1=fig.add_subplot(1,2,1)
data_lm=sm.formula.ols('Y~X',data=data).fit()
plt.scatter(data.X,data_lm.resid,axes=ax1)
ax1.set_xlabel('X')
ax1.set_ylabel('resid_ols')
ax1.set_title('resid_ols | X')

# 查看wls估计的残差与X的散点图
ax2=fig.add_subplot(1,2,2)
data_lm_wls=sm.formula.wls('Y~X',weights=1/data.X**2,data=data).fit()
plt.scatter(data.X,data_lm_wls.resid,axes=ax2)
ax2.set_xlabel('X')
ax2.set_ylabel('resid_wls')
ax2.set_title('resid_wls | X')

使用wls估计，异方差并没有消除。

（2）：使用FGLS估计对该模型进行重新估计，观察残差图并回答：异方差消除了吗？

data['resid']=data_lm.resid

# 进行辅助回归
data_lm_log=sm.formula.ols('np.log(resid**2)~X',data=data).fit()
h_hat=np.exp(data_lm_log.fittedvalues)

# 进行wls检验
data_lm_fgls=sm.formula.wls('Y~X',weights=1/h_hat,data=data).fit()

# 查看ols估计的残差与X的散点图
fig=plt.figure(figsize=(13,6))
ax1=fig.add_subplot(1,2,1)
data_lm=sm.formula.ols('Y~X',data=data).fit()
plt.scatter(data.X,data_lm.resid,axes=ax1)
ax1.set_xlabel('X')
ax1.set_ylabel('resid_ols')
ax1.set_title('resid_ols | X')

# 查看fgls估计的残差与X的散点图
ax2=fig.add_subplot(1,2,2)
plt.scatter(data.X,data_lm_fgls.resid,axes=ax2)
ax2.set_xlabel('X')
ax2.set_ylabel('resid_fgls')
ax2.set_title('resid_fgls | X')

使用FGLS估计，异方差并没有消除。

（3）：画出log(Y)与X的散点图，观察方差的状况，说说你的发现；根据散点图的情况，请大胆假设一个你认为正确的模型。

# 直接看log(Y)与X的散点图
fig=plt.figure()
plt.scatter(data.X,np.log(data.Y))
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_title('log(Y) | X')

发现log(Y)在X上异方差情况不明显。
假设模型是
$$\log (Y)={\beta }_0+{\beta }_{1}X+{\beta }_2{X}^2+u$$

data_lm_logy=sm.formula.ols('np.log(Y)~X+I(X**2)',data=data).fit()
print(data_lm_logy.summary())

各变量系数均显著，R方也很大。

（4）：考虑新模型
$$\log (Y)={\beta }_0+{\beta }_{1}X+{\beta }_2{X}^2+u$$
使用ols估计该模型，并画出残差散点图，说说你的发现

估计结果如上。

# 查看ols估计的残差与X的散点图
fig=plt.figure(figsize=(13,6))
ax1=fig.add_subplot(1,2,1)
data_lm=sm.formula.ols('Y~X',data=data).fit()
plt.scatter(data.X,data_lm.resid,axes=ax1)
ax1.set_xlabel('X')
ax1.set_ylabel('resid_ols')
ax1.set_title('resid_ols | X')

# 查看ols估计log(Y)的残差与X的散点图
ax2=fig.add_subplot(1,2,2)
data_lm_logy=sm.formula.ols('np.log(Y)~X+I(X**2)',data=data).fit()
plt.scatter(data.X,data_lm_logy.resid,axes=ax2)
ax2.set_xlabel('X')
ax2.set_ylabel('resid_log(Y)_ols')
ax2.set_title('resid_log(Y) | X')

发现对因变量取log再回归后，异方差的现象消失了。

（5）：综合以上四个问题，谈谈你对纠正模型异方差的见解。

异方差模型有可能是因为模型设计不合理造成的。在这种情况下，即使采用wls或者FGLS估计，也不一定能够消除异方差现象。当出现残差的方差随着自变量变大而变大的现象时，可以考虑对因变量取对数后再回归，可以很好地缓解异方差程度。

感谢Datawhale
感谢Gitmodel
参考文献
【教程地址】https://github.com/Git-Model/Modeling-Universe/tree/main/Data%20Analysis%20and%20Statistical%20Modeling
【备用gitee地址】https://gitee.com/mr-yinbob/Modeling-Universe/tree/main/Data%20Analysis%20and%20Statistical%20Modeling