<jupyter_start><jupyter_text>Multi-level Models in Keras PlaygroundLinear Mixed effects models, also known as hiearchical linear models, also known as multi-level models, are powerful linear ensemble modeling tools that can do both regression and classification tasks for many structured data sets. This notebook describes what a multi-level model is, and how to implement a model using the neural network library keras. Model outputs are compared to multi-level models available in the statsmodels package.A comparison among:  [StatsModels](https://github.com/statsmodels/statsmodels)  [Keras](https://github.com/fchollet/keras) with Tensorflow backendFor brevity, this tutorial will ignore cross-validation and hold out data as tools for model assessment. A very brief introduction to multi-level modelsMulti-level models account for different levels within a data set. Levels are groupings of data that apply across several observations. For example, a classic data set (simulated below), is the math achievement versus SES status for students who attend catholic schools versus public schools. The first level is the student's SES, whereas the second level is their attended school. Multi-level models can account for fixed effects (i.e., the variance does not change within groups), and random effects (i.e., the variance is distributed across groups). Multi-level models are linear models. In the case of the catholic school data set the equation that predicts student math achievement is:This equation is wrong lol$$ math\_achievement = \alpha_{01} + \beta_{01} * SES + \beta_{02} * catholic\_school $$Takes the general form of:$$ Y_{ij} = \beta_{0j} + \beta_{1j}X+{ij} + r_{ij} $$$$ \beta_{0j} = \gamma_{00} + \gamma_{01}W_j + u_{0j} $$$$ \beta_{1j} = \gamma_{10} + \gamma_{11}W_j + u_{1j} $$And the more specific form:$$ Y_{math achievement, i} = \beta{0,school} + \beta_{1,school}X_{i, SES} + \beta_{2,school}X_{i,school} + r_{i,j} $$$$ \beta_{0, school} = \gamma_{00} + \gamma_{01}W_{school} + u_{0, school} $$$$ \beta_{1, school} = \gamma_{00} + \gamma_{11}W_{school} + u_{0, school} $$ Where| variable                               | description                                                                                      ||----------------------------------------|--------------------------------------------------------------------------------------------------|| $i=1,2,3...$                           | the student indicator, i.e., the student ID.                                                     || $j=catholic,public$                    | the school group indicator                                                                       || $\beta_{0, school}$,$\beta_{1,school}$ | level-1 coefficients. in this case the SES status and the categorical school belonging variable. || $\gamma_{00}...\gamma{11}$             | level-2 coefficients. Also known as fixed effects.                                               || $X_{ij}$                               | Level-1 variable. SES, etc.                                                                      || $W_{ij}$                               | Level-2 predictor. school belonging, etc.                                                        || $r_{ij}$                               | Level-1 random effect.                                                                           || $u_{0,j},u_{1,j}...$                   | Level-2 random effects.                                                                          |These equations in summation for different intercepts and coefficients for models depending on whether a student attended a catholic school or public school.<jupyter_code>import numpy as np
import statsmodels.formula.api as smf
from patsy import dmatrices
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.metrics import r2_score

%load_ext watermark<jupyter_output><empty_output><jupyter_text>`watermark` prints the versions of libraries, pythons, and computer hardware in case this matters for your use of this notebook. It notes that `keras` is not installed. This is because `keras` comes with a `tensorflow` installation and thus does not need to be installed again.<jupyter_code>%watermark -v -m -p numpy,pandas,statsmodels,tensorflow,keras,matplotlib<jupyter_output>CPython 3.6.9
IPython 7.8.0

numpy 1.16.4
pandas 0.25.1
statsmodels 0.10.1
tensorflow 1.12.0
keras not installed
matplotlib 3.1.1

compiler   : GCC 7.3.0
system     : Linux
release    : 4.15.0-36-generic
machine    : x86_64
processor  : x86_64
CPU cores  : 20
interpreter: 64bit<jupyter_text>Prepare data and quick visualizationFor the first comparison, we use the california housing data set and create a categorical variable called `HouseAgeGroup`. This is due to the data analysis showing there is probably several groups of house ages based on periods of housing expansion. `sklearn` has many built in data sets that include data descriptions. The california housing data set has no strong correlations between data points and makes for a good first attempt at modeling data.<jupyter_code>from sklearn.datasets import fetch_california_housing as get_data
data = get_data()
df = pd.DataFrame(data['data'], columns=data['feature_names'])
# popGroup is high=1; low=0
df['popGroup'] = df.Population.apply(lambda x: 1 if x >= 1000 else 0).astype('category')
def house_age_group(age):
    if age >= 52:
        return 3
    elif age >= 30 and age < 52:
        return 2
    elif age < 30:
        return 1
df['HouseAgeGroup'] = df.HouseAge.apply(lambda x: house_age_group(x))
df.head()
print(data['DESCR'])
fig, ax = plt.subplots(1, 5, figsize=(15, 5), sharey=True)

for n, c in enumerate(df.columns[1:6]):
    df.plot.scatter(ax=ax[n], x=c, y='MedInc', alpha=0.05)
df.HouseAge.hist(bins=np.arange(1, 53, 1))
fig, ax = plt.subplots(figsize=(9, 7))
z = df.corr('spearman')
cbar = ax.pcolormesh(z, cmap='seismic', vmin=-0.7, vmax=0.7)
fig.colorbar(cbar, label='Spearman Correlation')
ax.set_xticks(np.arange(0, 8, 1)+0.5)
ax.set_yticks(np.arange(0, 8, 1)+0.5)
ax.set_xticklabels(z.columns)
ax.set_yticklabels(z.index)
for n,mc in enumerate(z.values):
    for i,m in enumerate(mc):
        ax.text(n+0.3, i+0.35, str(round(m,2)), color='black', fontsize=12)<jupyter_output><empty_output><jupyter_text>Using StatsModels to perform a linear mixed model of reaction time`statsmodels` can use `R`-like formulas to define fixed effects equations. However it uses the groups argument instead of the `|` within the equation to declare random effects. It is common practice to have a variable that has random effects also have fixed effects. This is because [random effects without fixed effects implies that the variable has no average effect.](https://stats.stackexchange.com/questions/173159/can-a-variable-be-both-random-and-fixed-effect-at-the-same-time-in-a-mixed-effec)<jupyter_code># https://www.statsmodels.org/stable/mixed_linear.html
formula = "MedInc ~ AveRooms + AveBedrms + AveRooms*AveBedrms + C(HouseAgeGroup)"
md  = smf.mixedlm(formula, df, groups=df['HouseAgeGroup'])
mdf = md.fit()
print(mdf.summary())

fe_params = pd.DataFrame(mdf.fe_params,columns=['LMM'])
random_effects = pd.DataFrame(mdf.random_effects)
random_effects = random_effects.transpose()
random_effects = random_effects.rename(index=str, columns={'groups': 'LMM'})
random_effects<jupyter_output>Mixed Linear Model Regression Results
======================================================================
Model:                 MixedLM     Dependent Variable:     MedInc     
No. Observations:      20640       Method:                 REML       
No. Groups:            3           Scale:                  1.7434     
Min. group size:       1273        Likelihood:             -35043.1383
Max. group size:       10669       Converged:              Yes        
Mean group size:       6880.0                                         
----------------------------------------------------------------------
                      Coef.    Std.Err.      z     P>|z| [0.025 0.975]
----------------------------------------------------------------------
Intercept              3.729        1.321    2.823 0.005  1.140  6.318
C(HouseAgeGroup)[T.2] -0.057        1.867   -0.031 0.976 -3.717  3.603
C(HouseAgeGroup)[T.3]  0.201        1.868    0.108 0.914 -3.460  3.862
AveRooms               [...]<jupyter_text>In this case it seems like these groups are not so important (high or low population). However the other features in the model seem to be able to explain at least half the variance in the median incomes of home owners in southern california.<jupyter_code>ypred = mdf.predict(df)

fig, ax = plt.subplots()
ax.scatter(df['MedInc'], ypred, alpha=0.05)
ax.set_ylim(0, 10)
ax.set_ylabel('Predicted', fontsize=15)
ax.set_xlabel('Actual', fontsize=15)
ax.plot([0, 10], [0, 10], color='red')
print('R2 score:', r2_score(df['MedInc'], ypred))<jupyter_output>R2 score: 0.517097818228871<jupyter_text>creating a design matrix from a statsmodels formula`statsmodels` can accept `pandas` dataframes directly as input with the defined groups. `keras` can not. Thus we need to create a [design matrix](https://en.wikipedia.org/wiki/Design_matrix) directly for training the `keras` model.<jupyter_code>Y, X = dmatrices(formula, data=df, return_type='matrix')
Terms = X.design_info.column_names
_, Z = dmatrices('MedInc ~ -1 + C(HouseAgeGroup)', data=df, return_type='matrix')
X = np.asarray(X) # fixed effect
Z = np.asarray(Z) # mixed effect
Y = np.asarray(Y).flatten()
nfixed = np.shape(X)
nrandm = np.shape(Z)<jupyter_output><empty_output><jupyter_text>Using Keras`keras` is a library that allows the construction of neural networks. Neural networks at the most basic level are linear combinations of variables, that is linear models. They can include a lot more sophistication but at their core, they are no different than any other model that is based on linear combinations of variables. Thus, `keras` provides a modular and verbose way to construct multi-level models.<jupyter_code>import tensorflow.keras as keras
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Add, Dense
from tensorflow.keras import backend as K
from tensorflow.keras.callbacks import TensorBoard
K.clear_session()
nb_epoch = 500

fixedpred = np.argmax(X,axis=1)
randmpred = np.argmax(Z,axis=1)

Xinput      = Input(batch_shape=(None, nfixed[1]-1), name='level_1_variables')
fixed_keras = Dense(1, input_dim=nfixed[1]-1, name = 'fixedEffect')(Xinput)

Zinput     = Input(batch_shape=(None, nrandm[1]), name='level_2_variables')
randm_keras = Dense(1, input_dim=nrandm[1], use_bias=None, name = 'randomEffect')(Zinput)

merged = keras.layers.add([fixed_keras, randm_keras])

model = Model([Xinput, Zinput], merged)

model.compile(loss='mean_squared_error', optimizer='adam')

# train the model
model.fit([X[:,1:], Z], Y.flatten(), 
          epochs=nb_epoch, 
          batch_size=100, 
          verbose=0,
          shuffle=True,
          )

Ypredict = model.predict([X[:,1:], Z])

betakeras = np.hstack((model.get_weights()[1], model.get_weights()[0].flatten()))
bkeras    = model.get_weights()[2].flatten()
from tensorflow.keras.utils import plot_model

pm = plot_model(model,
    to_file='model.png',
    show_shapes=True,
    show_layer_names=True,
    rankdir='TB')
from IPython.display import display, Image
display(Image(filename='model.png'))
fe_params['Keras'] = pd.Series(betakeras, index=fe_params.index)
random_effects['Keras'] = pd.Series(bkeras, index=random_effects.index)
fe_params
fig, ax = plt.subplots(figsize=(5, 10))

yticks = np.arange(fe_params.shape[0])
ax.plot(fe_params, yticks)
ax.set_yticks(yticks)
ax.set_yticklabels(labels=fe_params.index, rotation=0)

ax.legend(['LMM', 'Keras'], fontsize=15)
ax.set_xlabel('Coefficient value', fontsize=15)
random_effects
fig, ax = plt.subplots(figsize=(10, 5))

random_effects.reset_index().plot(ax=ax)
# ax.set_xticks(np.arange(0, 20, 1))
ax.set_title('random effects', fontsize=15)
fig, ax = plt.subplots(figsize=(9, 7))
ax.scatter(ypred, Ypredict, alpha=0.5, label='model comparison')
ax.plot([-100, 100], [-100, 100], label='perfect match', color='red')
ax.set_ylabel('Keras', fontsize=15)
ax.set_xlabel('statsmodels', fontsize=15)
ax.set_ylim(-20, 50)
ax.set_xlim(-25, 80)
ax.legend(fontsize=15, title='Median Income')
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ypred - Ypredict.flatten(), marker='o', linewidth=0)
ax.set_ylabel('statsmodels(y) - keras(y)', fontsize=15)
print('R2 score of model comparison:', r2_score(ypred, Ypredict))<jupyter_output>R2 score of model comparison: 0.9694711175635172<jupyter_text>Catholic School SimulationThe catholic school data set is a classic hierarchical data set used in education research to justify multi-level models (Bryk and Raudenbush, 2002). The data set is typically a 2D comparison of math achievement (typically described as a math test score) versus socio-economic status (SES) of students who attend catholic or public schools. The catholic school students perform better than the public school students thus justifying the need for a linear model to have two intercepts and slopes depending on the group they belong to.In this case, I have simulated the catholic school data set.<jupyter_code>num_samples = 1000

# The desired mean values of the sample.
mu = np.array([5.0, 10.0])

# The desired covariance matrix.
r = np.array([
        [  3.40, 5.75],
        [ 5.75,  5.50]
    ])

# Generate the random samples.
y = np.random.multivariate_normal(mu, r, size=num_samples)

catholic_data = pd.DataFrame({'SES':y[:,0], 'math_score':y[:,1]})
catholic_data['catholic_student'] = [1 if n>0.5 else 0 for n in np.random.random(num_samples)]
catholic_data['math_score'] = catholic_data.apply(lambda x: x['math_score']*3 if x['catholic_student']==1 else x['math_score'], axis=1)
catholic_data['math_score'] = catholic_data['math_score']/catholic_data['math_score'].max()
catholic_data['SES'] = catholic_data['SES'].apply(lambda x: (x - catholic_data['SES'].mean())/catholic_data['SES'].std())
catholic_data['colors'] = catholic_data['catholic_student'].apply(lambda x: 'green' if x==1 else 'purple')
catholic_data.describe()
fig, ax = plt.subplots(figsize=(9, 7))
catholic_data.plot.scatter(x='SES', y='math_score', color=catholic_data['colors'], alpha=0.5, ax=ax, s=55)
ax.set_ylabel('Math Achievement', fontsize=15)
ax.set_xlabel('Socio-Economic Status', fontsize=15)<jupyter_output><empty_output><jupyter_text>statsmodels catholic data<jupyter_code># https://www.statsmodels.org/stable/mixed_linear.html

# random effects should be fixed effects unless you want to imply the average effect of the random effect is 0
# https://stats.stackexchange.com/questions/173159/can-a-variable-be-both-random-and-fixed-effect-at-the-same-time-in-a-mixed-effec

formula = 'math_score ~ SES + SES * C(catholic_student)'
md  = smf.mixedlm(formula, catholic_data, groups=catholic_data['catholic_student'])
mdf = md.fit()
print(mdf.summary())

fe_params = pd.DataFrame(mdf.fe_params,columns=['LMM'])
random_effects = pd.DataFrame(mdf.random_effects)
random_effects = random_effects.transpose()
random_effects = random_effects.rename(index=str, columns={'groups': 'LMM'})
random_effects

mdf.random_effects
ypred = mdf.predict(catholic_data)

fig, ax = plt.subplots()
ax.scatter(catholic_data['math_score'], ypred, alpha=0.5)
ax.set_ylabel('Predicted', fontsize=15)
ax.set_xlabel('Actual', fontsize=15)
ax.plot([0, 1], [0, 1], color='red')<jupyter_output><empty_output><jupyter_text>keras catholic data<jupyter_code>Y, X = dmatrices(formula, data=catholic_data, return_type='matrix')
Terms = X.design_info.column_names
_, Z = dmatrices('math_score ~ -1 + C(catholic_student)', data=catholic_data, return_type='matrix')
X = np.asarray(X) # fixed effect
Z = np.asarray(Z) # mixed effect
Y = np.asarray(Y).flatten()
nfixed = np.shape(X)
nrandm = np.shape(Z)
K.clear_session()
nb_epoch = 500

fixedpred = np.argmax(X,axis=1)
randmpred = np.argmax(Z,axis=1)

Xinput = Input(batch_shape=(None, nfixed[1]-1), name='individualEffects')
fixed_keras = Dense(1, input_dim=nfixed[1]-1, name='fixedEffect')(Xinput)

Zinput = Input(batch_shape=(None, nrandm[1]), name='schoolEffects')
randm_keras = Dense(1, input_dim=nrandm[1], use_bias=None, name='randomEffect')(Zinput)

merged = keras.layers.add([fixed_keras, randm_keras])

model = Model([Xinput, Zinput], merged)

model.compile(loss='mean_squared_error', optimizer='adam')

# train the model
model.fit([X[:,1:], Z], Y.flatten(), 
          epochs=nb_epoch, 
          batch_size=100, 
          verbose=0,
          shuffle=True,
          )

Ypredict = model.predict([X[:,1:], Z])

betakeras = np.hstack((model.get_weights()[1], model.get_weights()[0].flatten()))
bkeras = model.get_weights()[2].flatten()
from tensorflow.keras.utils import plot_model

pm = plot_model(model,
    to_file='model.png',
    show_shapes=True,
    show_layer_names=True,
    rankdir='TB')
from IPython.display import display, Image
display(Image(filename='model.png'))
fe_params['Keras'] = pd.Series(betakeras, index=fe_params.index)
random_effects['Keras'] = pd.Series(bkeras, index=random_effects.index)
fe_params
fig, ax = plt.subplots(figsize=(5, 10))

yticks = np.arange(fe_params.shape[0])
ax.plot(fe_params, yticks)
ax.set_yticks(yticks)
ax.set_yticklabels(labels=fe_params.index, rotation=0)

ax.legend(['LMM', 'Keras'], fontsize=15)
ax.set_xlabel('Coefficient value', fontsize=15)
random_effects<jupyter_output><empty_output><jupyter_text>It is important to note from before that, although the model fitting (i.e., regression coefficients) are not the same across different approach, the model prediction is highly similar (at least it pass the eyeball test).   statsmodels compared to keras predicted outputAs we can see below, the predicted math achievement value for each student is identical between the `statsmodels` multi-level model and the `keras` model.<jupyter_code>fig, (ax, ax2) = plt.subplots(1, 2, figsize=(13, 6))

ax.scatter(ypred, Ypredict, alpha=0.5, label='model comparison')
ax.plot([0, 1], [0, 1], color='red', linewidth=2)
ax.set_ylabel('Keras', fontsize=15)
ax.set_xlabel('statsmodels', fontsize=15)
ax.legend(fontsize=15, title='catholic school\nmath achievement')
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)

ax2.scatter(Y, ypred, label='statsmodels', alpha=0.5, s=100)
ax2.scatter(Y, Ypredict, label='keras', marker='x', color='black', alpha=0.5)
ax2.plot([0, 1], [0, 1], color='red', linewidth=2, label='perfect prediction')

ax2.legend(fontsize=15, title='catholic school\nmath achievement')
ax2.set_ylabel('Predicted', fontsize=15)
ax2.set_xlabel('Actual', fontsize=15) 

fig.tight_layout()
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ypred - Ypredict.flatten(), marker='o', linewidth=0)
ax.set_ylabel('statsmodels(y) - keras(y)', fontsize=15)
print('R2 score of model comparison:', r2_score(ypred, Ypredict))<jupyter_output>R2 score of model comparison: 0.9999643491763599
{'fraction_non_alphanumeric': 0.14219580198853174, 'fraction_numerical': 0.03051186280183071, 'mean_word_length': 3.8138769308685743, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 45, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '16145', 'license': 'MIT', 'n_tokens_mistral': 6496, 'n_tokens_neox': 6060, 'n_words': 1657, 'path': 'ipython-notebooks/Making_multilevel_models_in_keras.ipynb', 'repo_name': 'mnky9800n/python-plotting', 'size': '19009'}