https://stats.stackexchange.com/q/592434/69508

592434_69508.py

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import shapiro, probplot
import statsmodels.api as sm

data = sm.datasets.scotland.load()
data.exog = sm.add_constant(data.exog)

links = [sm.families.Gaussian(), sm.families.Gamma()]
link_names = [link.__str__().split(' ')[0].split('.')[-1] for link in links]
models = [sm.GLM(data.endog, data.exog, family=link) for link in links]
results = [model.fit() for model in models]
aics = [result.aic for result in results]
resid_dfs = [result.df_resid for result in results]

d = {'link':link_names, 'df':resid_dfs, 'AIC':aics}
pooled_results = pd.DataFrame(d)
pooled_results = pooled_results.set_index('link')

print(pooled_results)

print(shapiro(data.endog))

fig = plt.figure(figsize=(8,6), dpi=72, facecolor="white")
axes = plt.subplot(111)
heights, positions, patches = axes.hist(data.endog, facecolor='lightgray', edgecolor='k')

axes.spines['right'].set_color('none')
axes.spines['top'].set_color('none')
axes.xaxis.set_ticks_position('bottom')

axes.spines['bottom'].set_position(('axes', -0.05))
axes.yaxis.set_ticks_position('left')
axes.spines['left'].set_position(('axes', -0.05))

axes.set_xlim([np.floor(positions.min()), np.ceil(positions.max())])
axes.set_ylim([0,70])
axes.xaxis.grid(False)
axes.yaxis.grid(False)
axes.set_ylabel('Frequency')
axes.set_xlabel('Y')
fig.tight_layout()
plt.show()

(osm, osr), (slope, intercept, r) = probplot(data.endog)
plt.scatter(osm, osr, facecolors='none', edgecolors='k')
plt.plot(osm, osm *  slope + intercept, color='k')
plt.title('Normal Q-Q Plot')
plt.xlabel('Theoretical Quantiles')
plt.ylabel('Sample Quantiles')
plt.show()