<jupyter_start><jupyter_text>Reading Data<jupyter_code>data = pd.read_csv('/home/atrides/Desktop/R/statistics_with_Python/07_Regression/Data_Files/Supermodel.dat', sep='\t')
data.head()
fig = plt.figure(figsize=(30,12))


ax1=plt.subplot(131)
coef  =np.polyfit(data['age'], data['salary'],1)
poly1d_fn = np.poly1d(coef) 
_=plt.plot(data['age'],data['salary'],'co' , data['age'], poly1d_fn(data['age']), '--k')
_=ax1.set_title('salary vs age', fontdict={'fontsize': 18, 'fontweight': 'medium'})


ax2=plt.subplot(132)
coef  =np.polyfit(data['years'], data['salary'],1)
poly1d_fn = np.poly1d(coef) 
_=plt.plot(data['years'],data['salary'],'co' , data['years'], poly1d_fn(data['years']), '--k')
_=ax2.set_title('salary vs experience', fontdict={'fontsize': 18, 'fontweight': 'medium'})



ax3 = plt.subplot(133)
coef  =np.polyfit(data['beauty'], data['salary'],1)
poly1d_fn = np.poly1d(coef) 
_=plt.plot(data['beauty'],data['salary'],'co' , data['beauty'], poly1d_fn(data['beauty']), '--k')
_=ax3.set_title('salary vs beauty', fontdict={'fontsize': 18, 'fontweight': 'medium'})<jupyter_output><empty_output><jupyter_text>Fitting Model<jupyter_code>m01 = sm.ols('salary~age', data=data)
res_1 = m01.fit()
res_1.summary()
m02 = sm.ols('salary~age+years', data=data)
res_2 = m02.fit()
res_2.summary()
m03 = sm.ols('salary~age+years+beauty', data=data)
res_3 = m03.fit()
res_3.summary()<jupyter_output><empty_output><jupyter_text>Comparing Models<jupyter_code>from statsmodels.stats.anova import anova_lm

anova_res = anova_lm(res_1, res_2, res_3)
print(anova_res)<jupyter_output>df_resid           ssr  df_diff      ss_diff         F    Pr(>F)
0     229.0  49742.554716      0.0          NaN       NaN       NaN
1     228.0  48555.607734      1.0  1186.946982  5.589655  0.018906
2     227.0  48202.790126      1.0   352.817608  1.661514  0.198711<jupyter_text>It seems that including beauty as a predictor in our model doesn't improves the model much as we can see from aic, bic , r_squared and adjusted r_squared.  Checking For Outlier Influence<jupyter_code>from statsmodels.stats.outliers_influence import OLSInfluence
summary_frame = OLSInfluence(res_2).summary_frame()
summary_frame.head()
summary_frame = summary_frame[['cooks_d','standard_resid', 'student_resid', 'hat_diag' ]]
summary_frame.head()
resid = pd.DataFrame(data['salary'] - res_2.fittedvalues)
resid.columns = ['residual']

dfbeta = pd.DataFrame(pd.DataFrame(OLSInfluence(res_2).dfbeta)[0])
dfbeta.columns = ['dfbeta']

cov_ratio = pd.DataFrame(OLSInfluence(res_2).cov_ratio)
cov_ratio.columns = ['cov_ratio']
df_ = [data, resid, summary_frame, dfbeta]
from functools import reduce
final_summary = reduce(lambda left,right: pd.merge(left,right, left_index=True, right_index=True), df_)
final_summary.head()
large_resid  = final_summary[(final_summary['standard_resid']>=2) | (final_summary['standard_resid']<=-2)]
large_resid = pd.merge(large_resid, cov_ratio, how = 'left', right_index=True, left_index=True)
large_resid
k = 2 #number of predictors
n = len(data)#number of objervations

average_leverage = (k+1)/n
average_leverage
cvr_limit_high = 1+3*average_leverage
cvr_limit_low  = 1-3*average_leverage

(cvr_limit_low, cvr_limit_high)<jupyter_output><empty_output><jupyter_text>doctored data , i.e removing large_resid points from data to see how the output changes<jupyter_code>df_new = data.merge(large_resid, how='left', indicator=True)
df_new = df_new[df_new['_merge'] == 'left_only']
df_new.head()
df_new.reset_index(inplace=True,drop=True)
res_2.summary()
m04 = sm.ols('salary~age+years', data=df_new)
res_4 = m04.fit()
res_4.summary()
fig = plt.figure(figsize=(20,10))


a1=plt.subplot(122)
coef  =np.polyfit(df_new['age'], df_new['salary'],1)
poly1d_fn = np.poly1d(coef) 
_=plt.plot(df_new['age'],df_new['salary'],'co' , df_new['age'], poly1d_fn(df_new['age']), '--k')
_=a1.set_ylim([-10,100])
_=a1.set_title('salary vs age in doctored data', fontdict={'fontsize': 18, 'fontweight': 'medium'})


a2=plt.subplot(121)
coef  =np.polyfit(data['age'], data['salary'],1)
poly1d_fn = np.poly1d(coef) 
_=plt.plot(data['age'],data['salary'],'co' , data['age'], poly1d_fn(data['age']), '--k')
_=a2.set_ylim([-10,100])
_=a2.set_title('salary vs age in original data', fontdict={'fontsize': 18, 'fontweight': 'medium'})<jupyter_output><empty_output><jupyter_text>In this case it seems that the line outliers affect our model to some extent!¶  Testing Assumptions<jupyter_code>from statsmodels.stats.stattools import durbin_watson
durbin_watson(res_2.resid)
from statsmodels.tools.tools import add_constant
from statsmodels.stats.outliers_influence import variance_inflation_factor
df_ = add_constant(data)
df_.drop(['salary'],axis=1, inplace=True)
df_.head()
vif = pd.Series([variance_inflation_factor(df_.values, i) 
               for i in range(1, df_.shape[1])], 
              index=df_.columns[1:])
vif
avg_vif  = np.mean(vif)
avg_vif

fig,ax = plt.subplots(figsize=(6, 4))
ax = plt.hist(final_summary['student_resid'],density=True,bins=30, edgecolor='black', linewidth=1.4)
plt.xlabel('student_resid', fontsize=14)
plt.show()
# performing shapiro-wilk test for checking the normality of errors
st.shapiro(final_summary['standard_resid'])<jupyter_output><empty_output><jupyter_text>error deviates from normality<jupyter_code>prediction = pd.DataFrame(res_2.fittedvalues)
prediction.columns = ['predicted']
prediction.head()
prediction['standarized_prediction']  = (prediction['predicted']-prediction['predicted'].mean())/prediction['predicted'].std()
prediction.head()
_ = sns.scatterplot(x= final_summary['standard_resid'], y = prediction['standarized_prediction'] )
_ = plt.axhline(y=0)
_=pg.qqplot(final_summary['standard_resid'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16422594142259414, 'fraction_numerical': 0.04079497907949791, 'mean_word_length': 6.945983379501385, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '361432', 'license': 'MIT', 'n_tokens_mistral': 2405, 'n_tokens_neox': 2226, 'n_words': 366, 'path': 'Python/statistics_with_Python/07_Regression/Smart_Alex/Task2_Supermodel.ipynb', 'repo_name': 'snehilk1312/AppliedStatistics', 'size': '5736'}