Enkel lineærregresjon med leieindeks-data (ISTX1003)¶
Denne notatboken ble brukt på Zoom-forelsning nr 1, 21.10.2020, se også forelesningsnotater fra timen på https://www.math.ntnu.no/emner/IST100x/ISTx1003/Zoom1InClass20201021.pdf
Spørsmål: hva er sammenhengen mellom leiepris for leiligheter i München og areal for leiligheten?
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import statsmodels.formula.api as smf
import statsmodels.api as sm
from statsmodels.distributions.empirical_distribution import ECDF
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)# fordi en av plottefunksjonene må oppdateres av utvikler og kommer per i dag med lang FutureWarning som vi ikke vil se på
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "last_expr"
Steg 1: Bli kjent med data¶
Hva ser vi etter når vi ser på oppsummeringsmål og plott av data?
In [2]:
filurl="http://www.uni-goettingen.de/de/document/download/64c29c1b1fccb142cfa8f29a942a9e05.raw/rent99.raw"
df=pd.read_csv(filurl," ") # bruker mellomrom mellom kolonner
print(df.shape)
df.sort_values(by=['area'],inplace=True) # Sorterer datasettet mhp area (for plotting)
print(df.head()) # Merk at den originale indekseringen av rader bevares
print(df[['rent','area']].describe())
print(df[['rent','area']].corr())
I denne presentasjonene skal vi bare se på rent som respons og area som forklaringsvariabel.
Legg merke til oppsummerende mål: minste leilighet er 20 m$^2$ og største er 160 m$^2$. Leien går fra 40.5 til 1843 Euro pr måned.
Vi plotter disse hver for seg og mot hverandre.
In [3]:
plt.subplot(1, 2, 1)
sns.distplot(df['rent'],
kde = False, hist_kws=dict(edgecolor="k", linewidth=2))
plt.subplot(1, 2, 2)
sns.boxplot(y='rent',data=df)
plt.tight_layout()
plt.show
Out[3]:
In [4]:
plt.subplot(1, 2, 1)
sns.distplot(df['area'],
kde = False, hist_kws=dict(edgecolor="k", linewidth=2))
plt.subplot(1, 2, 2)
sns.boxplot(y='area',data=df)
plt.tight_layout()
plt.show
Out[4]:
In [5]:
sns.relplot(x='area', y='rent',data = df)
plt.show()
In [6]:
sns.pairplot(df, vars = ['rent','area'],
diag_kind = 'kde',
plot_kws = {'alpha': 0.9, 's': 20})
#plot_kws=dict(alpha=0.4)
plt.show()
Steg 2: Spesifiser en matematisk modell¶
Hvordan bruker vi en såkalt modellformel for å spesifisere regresjonsmodellen?
In [7]:
formel='rent~area'
Steg 3: Initialiser og tilpass modellen¶
In [8]:
modell = smf.ols(formel,data=df)
resultat = modell.fit()
Steg 4: Presenter resultater fra den tilpassede modellen¶
- Hva er ligningen for den estimerte regresjonsmodellen? Hva betyr de ulike elementene?
- Hvordan kan vi tolke den estimerte verdien til skjæringspunktet $\hat{\beta}_0$, og til stigningstallet $\hat{\beta}_1$?
- Hva er den beste måten å kommunisere usikkerhet i stigningstallet (og skjæringspunktet)?
- Hva betyr det at en forklaringsvariabel er signifikant?
In [9]:
print(resultat.summary())
In [12]:
predres = resultat.get_prediction().summary_frame(alpha = 0.05) # Estimert forventningsverdi, prediksjon og intervaller
df_pred = pd.concat([df[["area","rent"]], predres], axis=1)
sns.relplot(x="area", y="rent",data = df_pred, label="Observasjoner")
sns.lineplot(x='area', y='mean', data = df_pred, label = "Tilpasset regresjonslinje")
sns.lineplot(x='area', y='mean_ci_lower', data = df_pred, color = "lightblue",label = "95% Konfidensintervall \n for regresjonslinjen")
sns.lineplot(x='area', y='mean_ci_upper', data = df_pred, color = "lightblue")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.ylabel("rent")
plt.show()
Steg 5: Evaluer om modellen passer til dataene¶
- Hva er modellantagelsene våre?
- Hvordan sjekker vi at modellantagelsene er oppfylt med plott?
- Hvordan evaluerer vi godheten til modellen?
In [13]:
df_resid = pd.DataFrame({"area": df["area"], "Residual": resultat.resid} )
sns.relplot(x='area', y="Residual",data = df_resid,marker="x")
plt.plot([2.5, 8], [0, 0], color = "grey",linewidth=2, linestyle = "--")
plt.show()
sns.scatterplot(resultat.fittedvalues, resultat.resid)
plt.ylabel("Residual")
plt.xlabel("Predikert verdi")
plt.show()
# Lage kvantil-kvantil-plott for residualene
sm.qqplot(resultat.resid,line='45',fit=True)
plt.ylabel("Kvantiler i residualene")
plt.xlabel("Kvantiler i normalfordelingen")
plt.show()
sns.distplot(resultat.resid,
kde = False, hist_kws=dict(edgecolor="k", linewidth=2))
plt.show()
# versjon 11 av seaborn har også ecdfplot
df_resid['Observasjonsnummer'] = df_resid.index
sns.relplot(x='Observasjonsnummer', y="Residual",
data = df_resid,marker="x", kind="line")
Out[13]:
In [ ]:
## Andel forklart variasjon $R^2$
In [22]:
print(resultat.rsquared_adj)
print(resultat.rsquared)
SSE=3082*resultat.mse_resid
SST=3082*np.var(df['rent'])
print(SSE)
print(SST)
1-SSE/SST
Out[22]:
Steg 6: Prediksjoner¶
In [21]:
sns.relplot(x='area', y='rent',data = df_pred, label="Observasjoner")
sns.lineplot(x='area', y='mean', data = df_pred, label = "Tilpasset regresjonslinje")
sns.lineplot(x='area', y='mean_ci_lower', data = df_pred, color = "lightblue",label = "95% Konfidensintervall \n for regresjonslinjen")
sns.lineplot(x='area', y='mean_ci_upper', data = df_pred, color = "lightblue")
sns.lineplot(x='area', y='obs_ci_lower', data = df_pred, color = "black",label = "95% Prediksjonsintervall \n nye leiligheter")
sns.lineplot(x='area', y='obs_ci_upper', data = df_pred, color = "black")
plt.ylabel("rent")
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()
In [ ]:
# Kan vi gjøre det bedre?
In [ ]:
df=pd.concat([df_in['Result'],
df_in['DBF.1']-df_in['DBF.2'],
df_in['ACE.1']-df_in['ACE.2'],
df_in['UFE.1']-df_in['UFE.2']],axis=1)
df.columns=['y','dobbeldiff','acediff','upressetdiff']
df.dropna(inplace=True)
print(df.describe())