%matplotlib inline

import numpy as np
import statsmodels.api as sm
from scipy import stats
from matplotlib import pyplot as plt

plt.rc("figure", figsize=(16,8))
plt.rc("font", size=14)

"""
Explanation: Generalized Linear Models
End of explanation
"""


print(sm.datasets.star98.NOTE)

"""
Explanation: GLM: Binomial response data
Load Star98 data
In this example, we use the Star98 dataset which was taken with permission
 from Jeff Gill (2000) Generalized linear models: A unified approach. Codebook
 information can be obtained by typing:
End of explanation
"""


data = sm.datasets.star98.load(as_pandas=False)
data.exog = sm.add_constant(data.exog, prepend=False)

"""
Explanation: Load the data and add a constant to the exogenous (independent) variables:
End of explanation
"""


print(data.endog[:5,:])

"""
Explanation: The dependent variable is N by 2 (Success: NABOVE, Failure: NBELOW):
End of explanation
"""


print(data.exog[:2,:])

"""
Explanation: The independent variables include all the other variables described above, as
 well as the interaction terms:
End of explanation
"""


glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())
res = glm_binom.fit()
print(res.summary())

"""
Explanation: Fit and summary
End of explanation
"""


print('Total number of trials:',  data.endog[0].sum())
print('Parameters: ', res.params)
print('T-values: ', res.tvalues)

"""
Explanation: Quantities of interest
End of explanation
"""


means = data.exog.mean(axis=0)
means25 = means.copy()
means25[0] = stats.scoreatpercentile(data.exog[:,0], 25)
means75 = means.copy()
means75[0] = lowinc_75per = stats.scoreatpercentile(data.exog[:,0], 75)
resp_25 = res.predict(means25)
resp_75 = res.predict(means75)
diff = resp_75 - resp_25

"""
Explanation: First differences: We hold all explanatory variables constant at their means and manipulate the percentage of low income households to assess its impact on the response variables:
End of explanation
"""


print("%2.4f%%" % (diff*100))

"""
Explanation: The interquartile first difference for the percentage of low income households in a school district is:
End of explanation
"""


nobs = res.nobs
y = data.endog[:,0]/data.endog.sum(1)
yhat = res.mu

"""
Explanation: Plots
We extract information that will be used to draw some interesting plots:
End of explanation
"""


from statsmodels.graphics.api import abline_plot

fig, ax = plt.subplots()
ax.scatter(yhat, y)
line_fit = sm.OLS(y, sm.add_constant(yhat, prepend=True)).fit()
abline_plot(model_results=line_fit, ax=ax)


ax.set_title('Model Fit Plot')
ax.set_ylabel('Observed values')
ax.set_xlabel('Fitted values');

"""
Explanation: Plot yhat vs y:
End of explanation
"""


fig, ax = plt.subplots()

ax.scatter(yhat, res.resid_pearson)
ax.hlines(0, 0, 1)
ax.set_xlim(0, 1)
ax.set_title('Residual Dependence Plot')
ax.set_ylabel('Pearson Residuals')
ax.set_xlabel('Fitted values')

"""
Explanation: Plot yhat vs. Pearson residuals:
End of explanation
"""


from scipy import stats

fig, ax = plt.subplots()

resid = res.resid_deviance.copy()
resid_std = stats.zscore(resid)
ax.hist(resid_std, bins=25)
ax.set_title('Histogram of standardized deviance residuals');

"""
Explanation: Histogram of standardized deviance residuals:
End of explanation
"""


from statsmodels import graphics
graphics.gofplots.qqplot(resid, line='r')

"""
Explanation: QQ Plot of Deviance Residuals:
End of explanation
"""


print(sm.datasets.scotland.DESCRLONG)

"""
Explanation: GLM: Gamma for proportional count response
Load Scottish Parliament Voting data
In the example above, we printed the NOTE attribute to learn about the
 Star98 dataset. statsmodels datasets ships with other useful information. For
 example:
End of explanation
"""


data2 = sm.datasets.scotland.load()
data2.exog = sm.add_constant(data2.exog, prepend=False)
print(data2.exog[:5,:])
print(data2.endog[:5])

"""
Explanation: Load the data and add a constant to the exogenous variables:
End of explanation
"""


glm_gamma = sm.GLM(data2.endog, data2.exog, family=sm.families.Gamma(sm.families.links.log()))
glm_results = glm_gamma.fit()
print(glm_results.summary())

"""
Explanation: Model Fit and summary
End of explanation
"""


nobs2 = 100
x = np.arange(nobs2)
np.random.seed(54321)
X = np.column_stack((x,x**2))
X = sm.add_constant(X, prepend=False)
lny = np.exp(-(.03*x + .0001*x**2 - 1.0)) + .001 * np.random.rand(nobs2)

"""
Explanation: GLM: Gaussian distribution with a noncanonical link
Artificial data
End of explanation
"""


gauss_log = sm.GLM(lny, X, family=sm.families.Gaussian(sm.families.links.log()))
gauss_log_results = gauss_log.fit()
print(gauss_log_results.summary())

"""
Explanation: Fit and summary (artificial data)
End of explanation
"""


