<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Generalized Linear Models
Step1:  GLM
Step2:  Load the data and add a constant to the exogenous (independent) variables
Step3:  The dependent variable is N by 2 (Success
Step4:  The independent variables include all the other variables described above, as
 well as the interaction terms
Step5:  Fit and summary
Step6:  Quantities of interest
Step7:  First differences
Step8:  The interquartile first difference for the percentage of low income households in a school district is
Step9:  Plots
We extract information that will be used to draw some interesting plots
Step10:  Plot yhat vs y
Step11:  Plot yhat vs. Pearson residuals
Step12:  Histogram of standardized deviance residuals
Step13:  QQ Plot of Deviance Residuals
Step14:  GLM
Step15:  Load the data and add a constant to the exogenous variables
Step16:  Model Fit and summary
Step17:  GLM
Step18:  Fit and summary (artificial data)
<ASSISTANT_TASK:>
Python Code:
%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
<END_TASK>