# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python (bayes)
#     language: python
#     name: bayes
# ---

# # Linear Regression (single variable) demo
#
# The dataset used is the Leinhardt Infant Mortality dataset from the R CAR package. Here we get the same data from the [Rdatasets repository](https://vincentarelbundock.github.io/Rdatasets/).
#
# The [Documentation for Leinhardt dataset](https://vincentarelbundock.github.io/Rdatasets/doc/carData/Leinhardt.html) can also be found in this repository.
#
# ## Problem Definition
#
# We want to predict the infant mortality response variable by using one or more explanatory variables available in the Leinhardt dataset.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
# %matplotlib inline

import warnings
warnings.filterwarnings("ignore")

# ## Load and explore data

try:
    leinhardt_df = pd.read_csv("Leinhardt.csv")
except:
    leinhardt_df = pd.read_csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Leinhardt.csv")
    leinhardt_df.to_csv("Leinhardt.csv")

leinhardt_df.head()

# ## Investigate correlation between pairs of variables

leinhardt_df["region"].unique()

leinhardt_df["oil"].unique()

# +
leinhardt_df.loc[leinhardt_df["region"] == "Asia", "region_i"] = 0
leinhardt_df.loc[leinhardt_df["region"] == "Europe", "region_i"] = 1
leinhardt_df.loc[leinhardt_df["region"] == "Americas" "region_i"] = 2
leinhardt_df.loc[leinhardt_df["region"] == "Africa", "region_i"] = 3

leinhardt_df.loc[leinhardt_df["oil"] == 'no', 'oil_i'] = 0
leinhardt_df.loc[leinhardt_df["oil"] == 'yes', 'oil_i'] = 1

_ = pd.plotting.scatter_matrix(leinhardt_df, diagonal="hist", figsize=(10,10))
# -

# ## Investigate correlation between Infant Mortality and Income
#
# From the scatter plots, there is negative correlation, but correlation doesn't appear to be linear. Thus Linear Regression on these variables is not appropriate.
#
# However, both distributions are right skewed, which is a strong hint that they might become linearly related on a log scale. On adding these log scaled columns, we see a strong linear relationship.

plt.scatter(leinhardt_df["income"], leinhardt_df["infant"])
plt.xlabel("income")
plt.ylabel("infant")
_ = plt.show()

plt.hist(leinhardt_df["infant"])
plt.xlabel("infant")
plt.ylabel("counts")
_ = plt.show()

plt.hist(leinhardt_df["income"])
plt.xlabel("income")
plt.ylabel("counts")
_ = plt.show()

leinhardt_df["log_income"] = np.log(leinhardt_df["income"])
leinhardt_df["log_infant"] = np.log(leinhardt_df["infant"])

plt.scatter(leinhardt_df["log_income"], leinhardt_df["log_infant"])
plt.xlabel("log_income")
plt.ylabel("log_infant")
_ = plt.show()

# ## Linear Regression (single variable) on log scale
#
# Fitting a OLS model against the log scaled data results in NaN results for the coefficients. This is because of missing values in the input or caused during the log transformation. 
#
# The R `lm` call is more forgiving and will automatically deal with this and report that there are 4 missing values and that the results might be incorrect.

leinhardt_df_na = leinhardt_df[(leinhardt_df["log_infant"].isna()) | 
                               (leinhardt_df["log_income"].isna())]
leinhardt_df_na.head()

leinhardt_df.dropna(subset=["log_income", "log_infant"], inplace=True)

# ## Linear Regression (single variable) on log scale after NaN removed
#
# Following information can be seen from the result of the linear regression:
#
# * intercept is 7.14 with standard error 0.317
# * coefficient of income is -0.5 with standard error 0.05
# * both values are large compared to their standard error (deviation) of the posterior, so they appear to be very statistically significant (P>|t|).
# * R-squared tells us how much of the variance is explained by the linear model.

# +
x = leinhardt_df["log_income"].values
X = np.vstack((x, np.ones(len(x)))).T

y = leinhardt_df["log_infant"].values

results = sm.OLS(y, X).fit()
results.summary()
# -



StarcoderdataJupyterScriptsDedupFiltered