This notebook analyzes the results of Bayesian inference for the recommendation lists.
import sys
from pathlib import Path
from textwrap import dedent
import pandas as pd
import numpy as np
import zarr
import seaborn as sns
import matplotlib.pyplot as plt
from plotnine import *
from scipy import stats
from scipy.special import expit, logit, logsumexp
from IPython.display import display, Markdown
import bookgender.datatools as dt
from bookgender.config import rng_seed
from lenskit.util import init_rng
from bookgender.nbutils import *
fig_dir = init_figs('RecModel')
seed = init_rng(rng_seed(), 'RecModelAnalysis')
rng = np.random.default_rng(seed)
seed
datasets = list(dt.datasets.keys())
Load the profile and list data for context:
profiles = pd.read_pickle('data/profile-data.pkl')
rec_lists = pd.read_pickle('data/rec-data.pkl')
rec_lists.head()
rec_lists.groupby('Set')['Total'].count()
Compute the algorithm names:
algo_names = rec_lists.reset_index().groupby('Set')['Algorithm'].apply(lambda x: sorted(x.unique()))
algo_names
And compute rec list length / distinctness stats:
recs = pd.read_parquet('data/study-recs.parquet')
recs.rename(columns={'dataset': 'Set', 'algorithm': 'Algorithm'}, inplace=True)
recs.info()
list_stats = recs.groupby(['Set', 'Algorithm'])['item'].agg(['count', 'nunique'])
list_stats['distfrac'] = list_stats['nunique'] / list_stats['count']
list_stats['distpct'] = list_stats['distfrac'] * 100
list_stats.head()
Let's load the samples!
samples = {}
summary = {}
for ds in datasets:
_zf = zarr.ZipStore(f'data/{ds}/inference/full/samples.zarr', mode='r')
_c = zarr.LRUStoreCache(_zf, 2**30)
samples[ds] = zarr.group(_c)
summary[ds] = pd.read_csv(f'data/{ds}/inference/full-summary.csv', index_col='name')
summary = pd.concat(summary, names=['Set'])
summary.head()
sample_size = len(samples['AZ']['lp__'])
sample_size
Do we have any parameters with troubling $\hat{R}$ values?
summary.sort_values('R_hat', ascending=False).head()
And let's compute LPPD and WAIC to assess model fit:
def ll_stats(ds):
ll_exp = samples[ds]['ll_exp']
ll_var = samples[ds]['ll_var']
lppd = np.sum(ll_exp)
pwaic = np.sum(ll_var)
return pd.Series({'lppd': lppd, 'pWAIC': pwaic, 'WAIC': -2 * (lppd - pwaic)})
pd.Series(datasets).apply(ll_stats).assign(Set=datasets).set_index('Set')
We want some functions for extracting data. Some of our things are per-list samples:
def list_samples(field, mean=False):
def _extract():
for ds, zg in samples.items():
data = zg[field][...].T
data = pd.DataFrame(data, index=rec_lists.loc[ds, :].index)
data.columns.name = 'Sample'
if mean:
data = data.mean(axis=1)
yield ds, data
return pd.concat(dict(_extract()), names=['Set'])
Others are per-algorithm samples:
def algo_samples(field):
def _extract():
for ds, zg in samples.items():
data = zg[field][...].T
names = algo_names[ds]
data = pd.DataFrame(data, index=names)
data.columns.name = 'Sample'
data.index.name = 'Algorithm'
yield ds, data
return pd.concat(dict(_extract()), names=['Set'])
Relable Algorithms:
algo_labels = {
'als': 'ALS',
'bpr': 'BPR',
'item-item': 'II',
'user-user': 'UU'
}
Repeat these helper functions for extracting implicit/explicit results:
def select_implicit(data, reset=True):
if reset:
data = data.reset_index()
implicit = data['Set'].str.endswith('-I')
if 'Algorithm' in data.columns:
implicit |= data['Algorithm'].str.endswith('-imp')
else:
implicit |= data['Set'] == 'AZ'
data = data.loc[implicit].assign(Set=data['Set'].str.replace('-I', ''))
if 'Algorithm' in data.columns:
algos = data['Algorithm'].str.replace('-imp', '').str.replace('wrls', 'als')
algos = algos.astype('category')
algos = algos.cat.rename_categories(algo_labels)
data['Algorithm'] = algos
return data
def select_explicit(data, reset=True):
if reset:
data = data.reset_index()
implicit = data['Set'].str.endswith('-I')
if 'Algorithm' in data.columns:
implicit |= data['Algorithm'].str.endswith('-imp')
data = data[~implicit].assign(Set=data['Set'].str.replace('-E', ''))
if 'Algorithm' in data.columns:
algos = data['Algorithm'].astype('category')
algos = algos.cat.rename_categories(algo_labels)
data['Algorithm'] = algos
return data
Let's plot the distributions of rec list biases. First we need to extract mean biases from the underlying samples, grouped by algorithm family. These are in log-odds space, expit
translates them back:
bias_smooth = algo_samples('biasP').stack()
bias_smooth = expit(bias_smooth)
bias_smooth.head()
Now we need expected new recommendation list proportions. This starts with the bias, plus the recommender's variance; the MCMC sampler outputs this as thetaRP
. We then feed it in to a binomial distribution, with known-item counts sampled from the data, as with user profiles.
bias_pred = algo_samples('thetaRP').stack()
bias_pred.head()
def _samp_obs(s):
known = rec_lists.loc[s.name, 'Known']
ns = rng.choice(known, len(s), replace=True)
ys = rng.binomial(ns, s)
return pd.Series(ys / ns, index=s.index)
bias_pred = bias_pred.groupby('Set').apply(_samp_obs)
bias_pred
def resample(x, n=sample_size):
s = pd.Series(rng.choice(x, n, replace=True))
s.index.name = 'Sample'
return s
Now we need the observed biases. For comparability, these should be the damped biases:
rec_lists['Bias'] = (rec_lists['female'] + 1) / (rec_lists['Known'] + 2)
bias_obs = rec_lists.groupby(['Set', 'Algorithm'])['PropFemale'].apply(resample)
bias_obs.head()
bias_data = pd.concat(dict(
Smoothed=bias_smooth,
Predicted=bias_pred,
Observed=bias_obs
), names=['Mode']).reset_index(name='Value')
bias_data['Mode'] = bias_data['Mode'].astype('category')
bias_data['Mode'].cat.reorder_categories(['Smoothed', 'Predicted', 'Observed'], inplace=True)
bias_data.head()
Let's plot the implicit runs:
grid = sns.FacetGrid(col='Set', row='Algorithm', hue='Mode',
data=select_implicit(bias_data),
sharey=False, aspect=1.2, height=2, margin_titles=True)
grid.map(sns.kdeplot, 'Value', clip=(0,1)).add_legend()
#plt.savefig(fig_dir / 'rec-implicit-dense.pdf')
make_plot(select_implicit(bias_data), aes('Value', color='Mode', linetype='Mode'),
geom_line(stat='density', bw='scott', clip=(0,1)),
scale_color_brewer('qual', 'Dark2'),
facet_grid('Algorithm ~ Set', scales='free_y'),
xlab('Recommender Proportion of Female Authors'), ylab('Density'),
legend_position='top', legend_title=element_blank(),
file='rec-implicit-dense', width=8, height=7)
And the explicit runs:
grid = sns.FacetGrid(col='Set', row='Algorithm', hue='Mode',
data=select_explicit(bias_data),
sharey=False, aspect=1.2, height=2.1)
grid.map(sns.kdeplot, 'Value', clip=(0,1)).add_legend()
# plt.savefig(fig_dir / 'rec-explicit-dense.pdf')
make_plot(select_explicit(bias_data), aes('Value', color='Mode', linetype='Mode'),
geom_line(stat='density', bw='scott', clip=(0,1)),
scale_color_brewer('qual', 'Dark2'),
facet_grid('Algorithm ~ Set', scales='free_y'),
xlab('Recommender Proportion of Female Authors'), ylab('Density'),
legend_position='top', legend_title=element_blank(),
file='rec-explicit-dense', width=8, height=5.5)
params = pd.DataFrame({
'Intercept': algo_samples('recB').mean(axis=1),
'Slope': algo_samples('recS').mean(axis=1),
'Variance': algo_samples('recV').mean(axis=1)
})
params.head()
param_samples = pd.DataFrame({
'Intercept': algo_samples('recB').stack(),
'Slope': algo_samples('recS').stack(),
'Variance': algo_samples('recV').stack()
})
param_samples.head()
def reg_ints(df):
return pd.Series({
'I_mean': df.Intercept.mean(),
'I_lo': df.Intercept.quantile(0.025),
'I_hi': df.Intercept.quantile(0.975),
'S_mean': df.Slope.mean(),
'S_lo': df.Slope.quantile(0.025),
'S_hi': df.Slope.quantile(0.975)
})
What are the implicit parameters?
select_implicit(params).pivot(index='Algorithm', columns='Set', values=['Intercept', 'Slope', 'Variance'])
imp_recB = select_implicit(algo_samples('recB')).set_index(['Set', 'Algorithm'])
imp_recS = select_implicit(algo_samples('recS')).set_index(['Set', 'Algorithm'])
imp_recV = select_implicit(algo_samples('recV')).set_index(['Set', 'Algorithm'])
imp_ds = imp_recB.index.levels[0]
imp_as = imp_recB.index.levels[1]
grid = sns.FacetGrid(col='Algorithm', row='Set', data=select_implicit(param_samples), margin_titles=True, height=2.2)
def _si_kde(*args, **kwargs):
plt.axhline(0, 0, 1, ls='-.', lw=0.5, color='grey')
plt.axvline(1, 0, 1, ls='-.', lw=0.5, color='grey')
sns.kdeplot(*args, levels=5, **kwargs)
grid.map(_si_kde, 'Slope', 'Intercept')
plt.savefig(fig_dir / 'reg-param-implicit.pdf')
imp_cis = select_implicit(param_samples).groupby(['Set', 'Algorithm']).apply(reg_ints)
make_plot(imp_cis.reset_index(),
aes(x='S_mean', xmin='S_lo', xmax='S_hi',
y='I_mean', ymin='I_lo', ymax='I_hi'),
geom_hline(yintercept=0, color='grey'),
geom_vline(xintercept=1, color='grey'),
geom_point(),
geom_errorbar(width=0.02),
geom_errorbarh(height=0.02),
facet_grid('Set ~ Algorithm'))
What are the explicit parameters?
select_explicit(params).pivot(index='Algorithm', columns='Set', values=['Intercept', 'Slope', 'Variance'])
grid = sns.FacetGrid(col='Algorithm', row='Set', data=select_explicit(param_samples), margin_titles=True, height=2.2)
def _si_kde(*args, **kwargs):
plt.axhline(0, 0, 1, ls='-.', lw=0.5, color='grey')
plt.axvline(1, 0, 1, ls='-.', lw=0.5, color='grey')
sns.kdeplot(*args, levels=5, **kwargs)
grid.map(_si_kde, 'Slope', 'Intercept')
plt.savefig(fig_dir / 'reg-param-explicit.pdf')
exp_cis = select_explicit(param_samples).groupby(['Set', 'Algorithm']).apply(reg_ints)
make_plot(exp_cis.reset_index(),
aes(x='S_mean', xmin='S_lo', xmax='S_hi',
y='I_mean', ymin='I_lo', ymax='I_hi'),
geom_hline(yintercept=0, color='grey'),
geom_vline(xintercept=1, color='grey'),
geom_errorbar(width=0.02, color='grey'),
geom_errorbarh(height=0.01, color='grey'),
geom_point(aes(color='Set', shape='Set')),
facet_grid('~ Algorithm'))
Let's try to put these parameter plots into a single display.
param_conf = pd.concat({
'Explicit': exp_cis.reset_index().astype({'Algorithm': 'str'}),
'Implicit': imp_cis.reset_index().astype({'Algorithm': 'str'})
}, names=['Mode'])
make_plot(param_conf.reset_index(),
aes(x='S_mean', xmin='S_lo', xmax='S_hi',
y='I_mean', ymin='I_lo', ymax='I_hi'),
geom_hline(yintercept=0, color='lightsteelblue', linetype='dashed'),
geom_vline(xintercept=1, color='lightsteelblue', linetype='dashed'),
geom_errorbar(width=0.05, color='grey'),
geom_errorbarh(height=0.03, color='grey'),
geom_point(aes(color='Set', shape='Set'), size=2, alpha=0.8),
facet_grid('Mode ~ Algorithm'),
scale_color_brewer('qual', 'Dark2'),
xlab('Slope'), ylab('Intercept'),
file='reg-params-all.pdf', width=8, height=5,
panel_grid=element_blank())
Let's plot those regressions.
We will start with implicit data:
imp_points = pd.merge(
select_implicit(profiles)[['Set', 'user', 'PropFemale']],
select_implicit(rec_lists)[['Set', 'Algorithm', 'user', 'Bias']]
).set_index(['Set', 'Algorithm', 'user'])
imp_points.head()
imp_params = select_implicit(params).set_index(['Set', 'Algorithm'])
xs = np.linspace(0, 1, 201)
imp_curves = np.outer(imp_params['Slope'], logit(xs))
imp_curves = imp_curves + imp_params['Intercept'].values.reshape((len(imp_params), 1))
imp_curves = expit(imp_curves)
imp_curves = pd.DataFrame(imp_curves, index=imp_params.index, columns=xs)
imp_curves.columns.name = 'x'
imp_curves = imp_curves.stack().reset_index(name='y')
imp_curves
make_plot(imp_points.reset_index(),
aes('PropFemale', 'Bias'),
geom_point(alpha=0.2, color='lightskyblue', size=1),
geom_rug(alpha=0.1),
geom_line(aes('x', 'y'), imp_curves, color='crimson'),
facet_grid('Set ~ Algorithm'),
xlab('Profile Proportion of Female Authors'),
ylab('Recomender Proportion of Female Authors'),
file='rec-scatter-imp.png', width=8, height=7, dpi=300)
And the explicit data:
exp_points = pd.merge(
select_explicit(profiles)[['Set', 'user', 'PropFemale']],
select_explicit(rec_lists)[['Set', 'Algorithm', 'user', 'Bias']]
).set_index(['Set', 'Algorithm', 'user'])
exp_points.head()
exp_params = select_explicit(params).set_index(['Set', 'Algorithm'])
xs = np.linspace(0, 1, 201)
exp_curves = np.outer(exp_params['Slope'], logit(xs))
exp_curves = exp_curves + exp_params['Intercept'].values.reshape((len(exp_params), 1))
exp_curves = expit(exp_curves)
exp_curves = pd.DataFrame(exp_curves, index=exp_params.index, columns=xs)
exp_curves.columns.name = 'x'
exp_curves = exp_curves.stack().reset_index(name='y')
exp_curves
make_plot(exp_points.reset_index(),
aes('PropFemale', 'Bias'),
geom_point(alpha=0.2, color='lightskyblue', size=1),
geom_rug(alpha=0.1),
geom_line(aes('x', 'y'), exp_curves, color='crimson'),
facet_grid('Set ~ Algorithm'),
xlab('Profile Proportion of Female Authors'),
ylab('Recomender Proportion of Female Authors'),
file='rec-scatter-exp.png', width=8, height=5.5, dpi=300)