Priors, begginer questions and Titanic Data

Questions v5
1

Hello,

Apologies for the likely common and oft repeated prior quandary. I’ve been reading through the excellent examples, questions here, and articles such as this one by the Stan team: Prior Choice Recommendations · stan-dev/stan Wiki · GitHub

I’m at the point where I just need to play and try things to get towards feeling able to take on real problems. I’m using the titanic data set from Kaggle to begin with: Titanic - Machine Learning from Disaster | Kaggle

The major questions I have are around chosing priors and how best to calculate and display marginal effects. Lets start by trying a few different models. One with super wide priors, one narrower and one attempting to add informed priors.

I tried to structure the data in a way which makes it easier for the model to fit. Ages have been centered on 12 (so that children are negative, and ‘adults’ positive), with a roughly 1.0 standard deviation. Class was changed so 3rd class became -1, 2nd 0 and first 1.

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
import seaborn as sb
import matplotlib.pyplot as plt
import xarray as xr
from sklearn.metrics import classification_report

df=pd.read_csv('titanic_train.csv')
##some extra stuff here to impute ages based off title etc
Fare_Mean=df.Fare.mean()
Fare_Std=df.Fare.std()
Age_Mean=df.AgeImpute_Title.mean()
Age_Std=df.AgeImpute_Title.std()
df['Age_Scaled']=(df.AgeImpute_Title-12)/Age_Std
df['Fare_Scaled']=(df.Fare-Fare_Mean)/Fare_Std
df['Sex_Bin']=np.where(df['Sex']=='male', -0.5, 0.5)
df['Class_Adj']=df.Pclass.replace({1:1, 2:0, 3:-1})

coords={"obs":df.PassengerId}

with pm.Model(coords=coords) as WidePriors:
    Age=pm.Data("Age", df.Age_Scaled,dims='obs')
    Sex=pm.Data("Sex", df.Sex_Bin, dims='obs')
    Class=pm.Data("Class", df.Class_Adj, dims='obs')
    beta0=pm.Normal("Int", mu=0, sigma=100)
    beta1=pm.Normal("Age_Param", mu=0, sigma=100)
    beta2=pm.Normal("Sex_Param", mu=0, sigma=100)
    beta3=pm.Normal("Class_Param", mu=0, sigma=100)
    mu=beta0+beta1*Age+beta2*Sex+beta3*Class
    p=pm.Deterministic("p", pm.invlogit(mu), dims='obs')
    y=pm.Bernoulli("y", p=p, observed=df.Survived, dims='obs')
    idata_Wide=pm.sample()
    pm.sample_posterior_predictive(idata_Wide,extend_inferencedata=True)
    pm.compute_log_likelihood(idata_Wide)

with pm.Model(coords=coords) as VaguePriors:
    Age=pm.Data("Age", df.Age_Scaled,dims='obs')
    Sex=pm.Data("Sex", df.Sex_Bin, dims='obs')
    Class=pm.Data("Class", df.Class_Adj, dims='obs')
    beta0=pm.Normal("Int", mu=0, sigma=1)
    beta1=pm.Normal("Age_Param", mu=0, sigma=1)
    beta2=pm.Normal("Sex_Param", mu=0, sigma=1)
    beta3=pm.Normal("Class_Param", mu=0, sigma=1)
    mu=beta0+beta1*Age+beta2*Sex+beta3*Class
    p=pm.Deterministic("p", pm.invlogit(mu), dims='obs')
    idata_Vague=pm.sample()
    pm.sample_posterior_predictive(idata_Vague,extend_inferencedata=True)
    pm.compute_log_likelihood(idata_Vague)
    
with pm.Model(coords=coords) as VaguePriors:
    Age=pm.Data("Age", df.Age_Scaled,dims='obs')
    Sex=pm.Data("Sex", df.Sex_Bin, dims='obs')
    Class=pm.Data("Class", df.Class_Adj, dims='obs')
    beta0=pm.Normal("Int", mu=0.4, sigma=0.1)
    beta1=pm.Normal("Age_Param", mu=-1.0, sigma=4.0)
    beta2=pm.Normal("Sex_Param", mu=2.45, sigma=0.1)
    beta3=pm.Normal("Class_Param", mu=0.8, sigma=0.2)
    mu=beta0+beta1*Age+beta2*Sex+beta3*Class
    p=pm.Deterministic("p", pm.invlogit(mu), dims='obs')
    y=pm.Bernoulli("y", p=p, observed=df.Survived, dims='obs')
    idata_Informed=pm.sample()
    pm.sample_posterior_predictive(idata_Informed,extend_inferencedata=True)
    pm.compute_log_likelihood(idata_Informed)
    
loo1=az.loo(idata_Wide)
loo2=az.loo(idata_Vague)
loo3=az.loo(idata_Informed)

az.compare({"Wide":loo1, "Vague":loo2, "Informed":loo3})

We can see that the models have the same rough values, but far smaller spread for the informed case, even for the age parameter which we gave vague priors on. Is this because giving informative priors where we can makes it easier for the MCMC algorithm to resolve the more vague ones? If so, does this make it easier to include factors which have overlap with the ones you’ve already included, EG here an interaction term between class and sex, or another term of the fare.

image
image1040×514 14 KB

In order to find prior values I had a look at what is in our data. For the intercept we have around 342 people who survived and 549, so a survival chance of 0.39. So beta0=ln(0.39/(1-0.39))~0.4 seems like a reasonable start point? Similar for age and sex;

df.groupby(['Survived', 'Sex_Bin'])['PassengerId'].count()

Survived  Sex_Bin
0         -0.5       468
           0.5        81
1         -0.5       109
           0.5       233
Name: PassengerId, dtype: int64

Where we can get log odds for the -0.5 group of ~-1.45, and for the 0.5 group 1.0. So beta~2.45 seems reasonable?

This seems sensible for binary or ordinal data with small numbers of levels, but what about for continuous such as Age? From exploratory data analysis we understand that it should most likely be a negative value, but I have no insight beyond that. The chosen prior is pretty excessive, but I’m not sure how I would hone in on a specific value here?

This also feels like it could very much lead to the wrong values. For example if we had a case where a misbalanced group was caused purely by one factor we would be overestimating the intercept if we used the overall group balanced to set our prior.

Finally, is there any functionality or easy way to handle calculating marginal effects. If doing models in R with the brms package theres the handy marginaleffects package which can calculate and do the plotting for me. I can see an easy enough way to do it using set data for EG the Sex variable here, but I cant see a neat way to handle a continuous variable.

2

I think you want to jump up a level of abstraction and think about what the prior and likelihood are doing for you. You usually do not want to use the data to inform the prior other than loosely. You want to use prior knowledge of the problem.

You’re not going to be able to tell anything from just inspecting the posterior intervals. You can make them as small as you want with with tighter and tighter priors. Then you want to do posterior predictive checks to check calibration of within-data fits and cross-validation to test calibration on out-of-data prediction for quantities of interest. The goal is to measure whether the prior combined with the likelihood lead to well calibrated posteriors in the sense of having the correct coverage either at the parameter level or the posterior predictive level for new data.

When the combination of prior and likelihood don’t work well, you have the freedom to change either one. On the likelihood side, maybe you add more covariates, maybe you interact covariates, maybe you transform covariates in other ways (like logit transforming proportions) or swap in robust (Student-t) noise for normal noise in the regression, and so on.

You might find the Bayesian workflow paper (on arXiv) and now the book (Cambridge) useful. It’s the same authors as the ones who wrote the Stan prior choice wiki and they discuss all of this. Unfortunately it’s all using Stan in R, so not ideal. I suspect someone will translate the book to PyMC or NumPyro, so it’ll at least be in Python.