Model comparison

Comparing closed capture-recapture models in PyMC

Author

Affiliations

Philip T. Patton

Marine Mammal Research Program

Hawaiʻi Institute of Marine Biology

Published

November 25, 2024

In this notebook, I demonstrate an approach to model selection in PyMC. To do so, follow the lead of King and Brooks (2008), although not nearly as elegantly. They demonstrate an approach to model selection for a typical suite of closed capture-recapture models. These include the effects of behavior \(b\), time \(t,\) and individual heterogeneity \(h\) on capture probabilities \(p\). The eight models considered here are combinations of the three: \(M_{0},\) \(M_{t},\) \(M_{b},\) \(M_{tb},\) \(M_{h},\) \(M_{th},\) \(M_{bh}\). The full model, \(M_{tbh}\), is
\[ \begin{equation} \text{logit} \; p_{it} = \mu + \alpha_t + \beta x_{it} + \gamma_i, \end{equation} \] where \(\mu\) is the average catchability, \(\alpha_t\) is the effect of each occasion on catchability, \(\beta\) is the behavioral effect, \(x_{it}\) indicates whether the individual has been previously caught, and \(\gamma_i\) is the individual random effect such that \(\gamma_i \sim \text{Normal}(0,\sigma)\). Formulating the model this way makes the other models nested subsets of the full model.

Like King and Brooks (2008), I use the the Moray Firth bottlenose dolphin data as a motivating example. Wilson, Hammond, and Thompson (1999) detected \(n=56\) dolphins over the course of \(T=17\) boat surveys between May and September 1992. They generated the capture-recapture histories by way of photo-identification, which is near and dear to my heart (and my dissertation).

# libraries 
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt 
import seaborn as sns
from pymc.distributions.dist_math import binomln, logpow

plt.rcParams['figure.dpi'] = 600
plt.style.use('fivethirtyeight')
plt.rcParams['axes.facecolor'] = 'white'
plt.rcParams['figure.facecolor'] = 'white'
pal = sns.color_palette("Set2")
sns.set_palette(pal)

# hyperparameters 
SEED = 808
RNG = np.random.default_rng(SEED)

def augment_history(history):
    '''Augment a capture history with all-zero histories.'''
    
    animals_captured, T = history.shape

    # create M - n all zero histories
    zero_history_count = M - animals_captured
    zero_history = np.zeros((zero_history_count, T))

    # tack those on to the capture history
    augmented = np.row_stack((history, zero_history))

    return augmented 

def get_behavior_covariate(history):
    
    # note the occasion when each individual was first seen
    first_seen = (history != 0).argmax(axis=1)
    
    # create the covariate for the behavior effect
    behavior_covariate = np.zeros_like(history)
    for i, f in enumerate(first_seen):
        behavior_covariate[i, (f + 1):] = 1

    return behavior_covariate

def get_occasion_covariate(history):

    _, T = history.shape
    l = []
    for t in range(T):
        oc = np.zeros_like(history)
        oc[:, t] = 1
        l.append(oc)

    return np.stack(l, axis=2)

def sim_N(idata):
    
    psi_samps = az.extract(idata).psi.values
    p_samps = az.extract(idata).p.values
    not_p = (1 - p_samps)
    
    if p_samps.ndim == 1:
        p_included = psi_samps * (not_p) ** T 
        number_undetected = RNG.binomial(M - n, p_included)

    elif p_samps.ndim == 3:
        p_included = psi_samps * not_p.prod(axis=1)
        number_undetected = RNG.binomial(1, p_included).sum(axis=0)

    N = n + number_undetected
    return N

# convert the dolphin capture history from '1001001' to array
dolphin = np.loadtxt('firth.txt', dtype=str)
dolphin = np.array([list(map(int, d)) for d in dolphin])

# augment the capture history with all zero histories
n, T = dolphin.shape
M = 500
dolphin_augmented = augment_history(dolphin)

# covariates for t and b
occasion_covariate = get_occasion_covariate(dolphin_augmented)
behavior_covariate = get_behavior_covariate(dolphin_augmented)

The discovery curve, the number of unique dolphins encountered as a function of the total number of dolphins encountered, may be flattening. This suggests that, at this point in the study, Wilson, Hammond, and Thompson (1999) may have encountered many of the unique individuals in the population.

# how many dolphins have been seen?
total_seen = dolphin.sum(axis=0).cumsum()

# how many new dolphins have been seen?
first_seen = (dolphin != 0).argmax(axis=1)
newbies = [sum(first_seen == t) for t in range(T)]
total_newbies = np.cumsum(newbies)

fig, ax = plt.subplots(figsize=(5, 3.5))
ax.plot(total_seen, total_newbies)
ax.fill_between(total_seen, total_newbies, alpha=0.2)
ax.set_title('Discovery curve')
ax.set_xlabel('Total dolphins')
ax.set_ylabel('Unique dolphins')
plt.show()

Figure 1: Discovery curve for the Moray Firth bottlenose dolphin surveys (Wilson, Hammond, and Thompson 1999).

Training each model

This notebook looks messier than the others, in that I train several models with little commentary along the way. In practice, it would probably be better to wrap these up into a function or a class. To complete the model, I used the following priors, \[ \begin{align} \psi &\sim \text{Uniform}(0, 1)\\ \mu &\sim \text{Logistic}(0, 1) \\ \alpha_t &\sim \text{Normal}(0, \sigma_{\alpha}) \\ \beta &\sim \text{Normal}(0, \sigma_{\beta}) \\ \gamma_i &\sim \text{Normal}(0, \sigma_{\gamma}) \\ \sigma_{\alpha} &\sim \text{InverseGamma}(4, 3) \\ \sigma_{\beta} &\sim \text{InverseGamma}(4, 3) \\ \sigma_{\gamma} &\sim \text{InverseGamma}(4, 3), \end{align} \] which were also used by King and Brooks (2008). Although note that I used an informative \(\text{Beta}(1, 5)\) prior for \(\psi\) in the full model (see below). I use the same logp seen in the occupancy and closed capture-recapture notebooks, which accounts for row-level zero-inflation. Unlike other notebooks, I did not look at the summaries or the trace plots unless the sampler indicated that it had issues during training.

Throughout the notebook, I use the nutpie sampler within PyMC. Nutpie is a NUTS sampler written in Rust, and is often faster than PyMC. Also, I have tweaked the sampling keyword arguments for each model, since they are a little finicky.

def logp(value, n, p, psi):
    
    binom = binomln(n, value) + logpow(p, value) + logpow(1 - p, n - value)
    bin_sum = pm.math.sum(binom, axis=1)
    bin_exp = pm.math.exp(bin_sum)

    res = pm.math.switch(
        value.sum(axis=1) > 0,
        bin_exp * psi,
        bin_exp * psi + (1 - psi)
    )
    
    return pm.math.log(res)

with pm.Model() as m0:

    # Priors
    # inclusion
    psi = pm.Uniform('psi', 0, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # Linear model
    mu_matrix = (np.ones((T, M)) * mu).T
    p = pm.Deterministic('p', pm.math.invlogit(mu_matrix))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(m0)

Figure 2: Visual representation of model \(M_{0}\).

with m0:
    m0_idata = pm.sample(nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for now

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	1.11	3
2000	1.13	3
2000	1.14	3
2000	1.13	7

with pm.Model() as mt:

    # Priors
    # inclusion
    psi = pm.Uniform('psi', 0, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # time effect
    sigma_alpha = pm.InverseGamma('sigma_alpha', 4, 3)
    alpha = pm.Normal('alpha', 0, tau=1/sigma_alpha, shape=T)

    # Linear model
    nu = mu + pm.math.dot(occasion_covariate, alpha)
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mt)

Figure 3: Visual representation of model \(M_t\).

with mt:
    mt_idata = pm.sample(nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 30 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.63	7
2000	0.63	7
2000	0.65	7
2000	0.61	7

with pm.Model() as mb:

    # Priors
    # inclusion
    psi = pm.Uniform('psi', 0, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)
    
    # behavior effect
    sigma_beta = pm.InverseGamma('sigma_beta', 4, 3)
    beta = pm.Normal('beta', 0, tau=1/sigma_beta)

    # Linear model
    nu = mu + behavior_covariate * beta 
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mb)

Figure 4: Visual representation of model \(M_b\).

with mb:
    mb_idata = pm.sample(nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for now

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.80	3
2000	0.84	3
2000	0.79	7
2000	0.83	3

with pm.Model() as mtb:

    # Priors
    # inclusion
    psi = pm.Uniform('psi', 0, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # time effect
    sigma_alpha = pm.InverseGamma('sigma_alpha', 4, 3)
    alpha = pm.Normal('alpha', 0, tau=1/sigma_alpha, shape=T)

    # behavior effect
    sigma_beta = pm.InverseGamma('sigma_beta', 4, 3)
    beta = pm.Normal('beta', 0, tau=1/sigma_beta)

    # Linear model
    nu = mu + pm.math.dot(occasion_covariate, alpha) + behavior_covariate * beta
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mtb)

Figure 5: Visual representation of model \(M_{tb}\).

with mtb:
    mtb_idata = pm.sample(nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 37 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
2000	0.56	7
2000	0.55	7
2000	0.57	7
2000	0.54	7

with pm.Model() as mh:

    # Priors
    # inclusion
    psi = pm.Uniform('psi', 0, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # individual effect
    sigma_gamma = pm.InverseGamma('sigma_gamma', 4, 3)
    gamma = pm.Normal('gamma', 0, tau=1/sigma_gamma, shape=M)

    # Linear model
    individual_effect = (np.ones((T, M)) * gamma).T
    nu = mu + individual_effect
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mh)

Figure 6: Visual representation of model \(M_h\).

with mh:
    mh_idata = pm.sample(3000, target_accept=0.99, nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 34 seconds

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
4000	0.08	63
4000	0.07	31
4000	0.09	31
4000	0.08	63

az.summary(mh_idata, var_names=['psi', 'mu', 'sigma_gamma'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
psi	0.185	0.039	0.122	0.257	0.002	0.001	589.0	912.0	1.01
mu	-2.788	0.329	-3.432	-2.241	0.018	0.013	376.0	572.0	1.01
sigma_gamma	0.775	0.351	0.291	1.453	0.024	0.017	213.0	419.0	1.03

az.plot_trace(mh_idata, figsize=(8, 6), var_names=['psi', 'mu', 'sigma_gamma']);

with pm.Model() as mth:

    # Priors
    # inclusion
    psi = pm.Beta('psi', 1, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # time effect
    sigma_alpha = pm.InverseGamma('sigma_alpha', 4, 3)
    alpha = pm.Normal('alpha', 0, tau=1/sigma_alpha, shape=T)

    # individual effect
    sigma_gamma = pm.InverseGamma('sigma_gamma', 4, 3)
    gamma = pm.Normal('gamma', 0, tau=1/sigma_gamma, shape=M)

    # Linear model
    individual_effect = (np.ones((T, M)) * gamma).T
    nu = mu + pm.math.dot(occasion_covariate, alpha) + individual_effect
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mth)

Figure 7: Visual representation of model \(M_{th}\).

with mth:
    mth_idata = pm.sample(draws=3000, target_accept=0.95, nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 3 minutes

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
4000	0.16	31
4000	0.17	31
4000	0.17	31
4000	0.14	15

az.summary(mth_idata, var_names=['psi', 'mu', 'sigma_alpha', 'sigma_gamma', 'alpha'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
psi	0.182	0.041	0.120	0.251	0.002	0.001	543.0	1020.0	1.01
mu	-3.058	0.437	-3.873	-2.308	0.021	0.015	450.0	951.0	1.00
sigma_alpha	0.923	0.356	0.381	1.577	0.004	0.003	7989.0	7805.0	1.00
sigma_gamma	0.899	0.425	0.293	1.644	0.028	0.020	212.0	512.0	1.01
alpha[0]	-1.095	0.601	-2.230	0.033	0.006	0.005	10515.0	8318.0	1.00
alpha[1]	0.570	0.403	-0.178	1.336	0.005	0.004	5931.0	7290.0	1.00
alpha[2]	-0.802	0.562	-1.829	0.296	0.006	0.005	10454.0	7723.0	1.00
alpha[3]	0.570	0.409	-0.199	1.342	0.005	0.004	5787.0	7480.0	1.00
alpha[4]	0.448	0.411	-0.325	1.207	0.006	0.004	5561.0	7430.0	1.00
alpha[5]	0.782	0.403	0.026	1.536	0.006	0.004	5355.0	6868.0	1.00
alpha[6]	0.177	0.441	-0.664	0.980	0.005	0.004	6445.0	8379.0	1.00
alpha[7]	0.024	0.456	-0.845	0.870	0.006	0.004	5890.0	6913.0	1.00
alpha[8]	-0.804	0.557	-1.832	0.235	0.006	0.005	8067.0	7704.0	1.00
alpha[9]	0.025	0.457	-0.841	0.860	0.005	0.004	7612.0	7826.0	1.00
alpha[10]	1.153	0.377	0.429	1.842	0.006	0.004	4699.0	6330.0	1.00
alpha[11]	-0.332	0.494	-1.230	0.608	0.006	0.004	7986.0	7561.0	1.00
alpha[12]	-1.109	0.609	-2.270	0.000	0.007	0.006	9351.0	6756.0	1.00
alpha[13]	-0.147	0.470	-1.047	0.713	0.006	0.004	7172.0	7546.0	1.00
alpha[14]	-1.106	0.608	-2.253	-0.003	0.006	0.005	9200.0	8052.0	1.00
alpha[15]	1.608	0.360	0.960	2.326	0.005	0.004	4561.0	7334.0	1.00
alpha[16]	-0.817	0.561	-1.884	0.225	0.006	0.005	9633.0	6540.0	1.00

az.plot_trace(mth_idata, figsize=(8, 10),
              var_names=['psi', 'mu', 'sigma_alpha', 'sigma_gamma', 'alpha']);

Figure 8: Trace plots for model \(M_{th}\).

with pm.Model() as mbh:

    # Priors
    # inclusion
    psi = pm.Beta('psi', 1, 1)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # behavior effect
    sigma_beta = pm.InverseGamma('sigma_beta', 4, 3)
    beta = pm.Normal('beta', 0, tau=1/sigma_beta)
    
    # individual effect
    sigma_gamma = pm.InverseGamma('sigma_gamma', 4, 3)
    gamma = pm.Normal('gamma', 0, tau=1/sigma_gamma, shape=M)

    # Linear model
    individual_effect = (np.ones((T, M)) * gamma).T
    nu = mu + behavior_covariate * beta + individual_effect
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mbh)

Figure 9: Visual representation of model \(M_{bh}\).

with mbh:
    mbh_idata = pm.sample(draws=3000, target_accept=0.95, nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for 2 minutes

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
4000	0.15	31
4000	0.15	31
4000	0.15	31
4000	0.15	31

az.summary(mbh_idata, var_names=['psi', 'mu', 'beta', 'sigma_gamma'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
psi	0.591	0.188	0.293	0.958	0.009	0.006	454.0	1087.0	1.01
mu	-3.628	0.560	-4.679	-2.624	0.031	0.022	333.0	766.0	1.01
beta	-1.560	0.252	-2.015	-1.061	0.005	0.003	2754.0	4288.0	1.00
sigma_gamma	2.300	0.740	1.052	3.723	0.043	0.030	291.0	612.0	1.01

az.plot_trace(mbh_idata, figsize=(8, 10),
              var_names=['psi', 'mu', 'beta', 'sigma_beta', 'sigma_gamma']);

with pm.Model() as mtbh:

    # Priors
    # inclusion
    psi = pm.Beta('psi', 1, 5)  

    # mean catchability 
    mu = pm.Logistic('mu', 0, 1)

    # time effect
    sigma_alpha = pm.InverseGamma('sigma_alpha', 4, 3)
    alpha = pm.Normal('alpha', 0, tau=1/sigma_alpha, shape=T)

    # behavior effect
    sigma_beta = pm.InverseGamma('sigma_beta', 4, 3)
    beta = pm.Normal('beta', 0, tau=1/sigma_beta)

    # individual effect
    sigma_gamma = pm.InverseGamma('sigma_gamma', 4, 3)
    gamma = pm.Normal('gamma', 0, tau=1/sigma_gamma, shape=M)

    # Linear model
    individual_effect = (np.ones((T, M)) * gamma).T
    nu = mu + pm.math.dot(occasion_covariate, alpha) + behavior_covariate * beta + individual_effect
    p = pm.Deterministic('p', pm.math.invlogit(nu))

    # Likelihood 
    pm.CustomDist(
        'y',
        1,
        p,
        psi,
        logp=logp,
        observed=dolphin_augmented
    )
    
pm.model_to_graphviz(mtbh)

Figure 10: Visual representation of model \(M_{tbh}\).

with mtbh:
    mtbh_idata = pm.sample(draws=2000, nuts_sampler='nutpie')

Sampler Progress

Total Chains: 4

Active Chains: 0

Finished Chains: 4

Sampling for a minute

Estimated Time to Completion: now

Draws	Step Size	Gradients/Draw
3000	0.30	15
3000	0.29	15
3000	0.29	15
3000	0.29	15

az.summary(mtbh_idata, 
           var_names=['psi', 'mu', 'alpha', 'beta', 'sigma_alpha', 'sigma_beta', 'sigma_gamma'])

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
psi	0.712	0.100	0.528	0.897	0.002	0.001	2553.0	4215.0	1.00
mu	-3.060	0.479	-3.931	-2.133	0.019	0.013	645.0	1313.0	1.00
alpha[0]	-3.133	0.729	-4.568	-1.855	0.015	0.010	2591.0	3674.0	1.00
alpha[1]	-0.675	0.525	-1.640	0.313	0.012	0.008	1992.0	3155.0	1.00
alpha[2]	-1.701	0.658	-2.897	-0.440	0.012	0.008	3312.0	4018.0	1.00
alpha[3]	-0.055	0.491	-0.991	0.865	0.010	0.007	2333.0	4005.0	1.00
alpha[4]	0.211	0.483	-0.710	1.095	0.009	0.006	2818.0	4354.0	1.00
alpha[5]	0.834	0.459	-0.045	1.663	0.009	0.006	2529.0	3670.0	1.00
alpha[6]	0.376	0.505	-0.551	1.353	0.009	0.006	3085.0	3796.0	1.00
alpha[7]	0.372	0.530	-0.659	1.338	0.009	0.006	3420.0	4718.0	1.00
alpha[8]	-0.563	0.664	-1.796	0.674	0.009	0.007	5469.0	4714.0	1.00
alpha[9]	0.467	0.520	-0.552	1.402	0.009	0.006	3226.0	4124.0	1.00
alpha[10]	1.657	0.435	0.887	2.508	0.010	0.007	2000.0	3439.0	1.00
alpha[11]	0.137	0.580	-0.996	1.179	0.009	0.007	4033.0	4856.0	1.00
alpha[12]	-0.892	0.758	-2.323	0.551	0.010	0.009	6313.0	4271.0	1.00
alpha[13]	0.384	0.539	-0.639	1.384	0.009	0.006	3577.0	4661.0	1.00
alpha[14]	-0.874	0.760	-2.373	0.476	0.010	0.009	6301.0	4217.0	1.00
alpha[15]	2.237	0.419	1.483	3.042	0.010	0.007	1762.0	3599.0	1.00
alpha[16]	-0.264	0.675	-1.542	0.994	0.009	0.008	5795.0	5445.0	1.00
beta	-3.293	0.393	-4.026	-2.550	0.009	0.007	1739.0	3061.0	1.00
sigma_alpha	1.567	0.601	0.668	2.669	0.009	0.007	4221.0	4797.0	1.00
sigma_beta	2.418	1.549	0.654	5.047	0.021	0.017	6836.0	5326.0	1.00
sigma_gamma	2.707	0.679	1.540	3.949	0.041	0.029	270.0	674.0	1.02

az.plot_trace(mtbh_idata, figsize=(8,14),
           var_names=['psi', 'mu', 'alpha', 'beta', 'sigma_alpha', 'sigma_beta', 'sigma_gamma']);

Figure 11: Trace plots for the model with \(M_{tbh}\).

The trace plots and summary statistics show convergence issues for many of the individual heterogeneity models. The variance parameter, \(\sigma_{\gamma},\) seems to sample poorly. Further, models with both behavioral and individual effects lead to extremely large estimates of \(\psi\). This appears to happen regardless of the size of the data augmentation \(M.\)

Note that I upped the target_accept value for some models. This slows the sampler, but lowers the risk of divergence.

Model comparison

Next, I select a model for inference using an approximation of leave-one-out (loo) cross-validation (Vehtari, Gelman, and Gabry 2017). This approximation can be calculated using PyMC. To do so, I calculate the log-likelihood for each model, which is added to the InferenceData object. This makes it possible to compare the models using loo and az.compare.

with m0:
    pm.compute_log_likelihood(m0_idata)

with mt:
    pm.compute_log_likelihood(mt_idata)

with mb:
    pm.compute_log_likelihood(mb_idata)

with mtb:
    pm.compute_log_likelihood(mtb_idata)

with mh:
    pm.compute_log_likelihood(mh_idata)

with mth:
    pm.compute_log_likelihood(mth_idata)

with mbh:
    pm.compute_log_likelihood(mbh_idata)

with mtbh:
    pm.compute_log_likelihood(mtbh_idata)

model_dict = {"m0": m0_idata, "mt": mt_idata, "mb": mb_idata, 
              "mtb": mtb_idata, "mh": mh_idata, "mth": mth_idata, 
              "mbh": mbh_idata, "mtbh": mtbh_idata}

comparison = az.compare(model_dict)

/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

The comparison tools notes issues with several of the models, suggesting a lack of robustness. Inspection of the comparison table shows that the struggling models all include the individual effect \(h.\) A more thorough analysis would consider reparameterizing the model, e.g., through the non-centered parameterization. In lieu of that, I simply discard the models that fail this test and re-do the comparison with the passing models.

comparison.round(2)

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
mtbh	0	-431.31	87.25	0.00	0.98	57.35	0.00	True	log
mth	1	-489.44	39.53	58.13	0.00	58.98	10.98	True	log
mtb	2	-493.57	22.01	62.27	0.00	62.09	10.93	True	log
mt	3	-493.85	14.33	62.54	0.00	59.55	12.28	False	log
mbh	4	-494.30	71.50	63.00	0.00	60.74	11.12	True	log
mh	5	-518.38	25.46	87.07	0.02	61.46	15.14	True	log
mb	6	-521.40	5.22	90.09	0.00	63.10	15.03	False	log
m0	7	-522.51	2.45	91.20	0.00	62.09	15.92	False	log

good_dict = {"m0": m0_idata, "mt": mt_idata, "mb": mb_idata, "mtb": mtb_idata}
good_comparison = az.compare(good_dict)
good_comparison.round(2)

/Users/philtpatton/miniforge3/envs/pymc/lib/python3.12/site-packages/arviz/stats/stats.py:795: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
mtb	0	-493.57	22.01	0.00	0.59	62.09	0.00	True	log
mt	1	-493.85	14.33	0.28	0.41	59.55	6.67	False	log
mb	2	-521.40	5.22	27.83	0.00	63.10	8.54	False	log
m0	3	-522.51	2.45	28.93	0.00	62.09	10.26	False	log

az.plot_compare(good_comparison, figsize=(5, 4));

Figure 12: Differences in the ELPD criteria, calculated using loo, for each model (Vehtari, Gelman, and Gabry 2017).

The comparison shows that all of the model weight belongs to two models: \(M_t\) and \(M_{tb}.\)

Model averaged predictions

Finally, we can use the model weights to simulate a weighted posterior of \(N.\) To do so, I take a weighted sample of each of the posteriors of \(N,\) with the weight dictated by the comparison tool.

posteriors = [sim_N(good_dict[model]) for model in good_dict]
weights = [good_comparison.loc[model].weight for model in good_dict]
sample_count = len(posteriors[0])

l = []
for w, p in zip(weights, posteriors):
    weighted_sample = RNG.choice(p, size=int(w * sample_count))
    l.append(weighted_sample)

weighted_posterior = np.concatenate(l)

fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(7, 6), sharex=True, sharey=True, tight_layout=True)

pal = sns.color_palette("Set2")

# labs = [k for k in good_dict.keys()]
labs = [r'$M_{0}$', r'$M_{t}$', r'$M_{b}$', r'$M_{tb}$']
for i, p in enumerate(posteriors):
    ax0.hist(p, color=pal[i], edgecolor='white', bins=60, alpha=0.6, label=labs[i])

ax0.set_title(r'Posteriors of $N$')
# ax1.set_title(r'Weighted posterior')

ax0.set_xlim((53, 150))
ax0.legend()

ax0.set_ylabel('Number of samples')
ax1.set_ylabel('Number of samples')

ax1.hist(weighted_posterior, edgecolor='white', bins=60, alpha=0.9, color=pal[6], label='Weighted')
ax1.legend()

plt.show()

Figure 13: Posteriors of \(N\) from the four models under consideration (top panel), with the model averaged posterior (bottom panel).

We can also look at the posterior densities of \(p\) from Model \(M_t,\) the second most weighted model.

p_samps = az.extract(mt_idata).p.mean(axis=0)

fig, ax = plt.subplots(figsize=(6, 4))

a = 0.4
# ax[0].set_title("Poisson")
pal = sns.color_palette('viridis', T)
for t in range(T):
    label_idx = t % 2
    if label_idx == 0:
        az.plot_dist(p_samps[t], ax=ax, color=pal[t], label=f'$t_{{{t}}}$',
                     plot_kwargs={'linewidth':3, 'alpha': a})
    else:
        az.plot_dist(p_samps[t], ax=ax, color=pal[t],
                     plot_kwargs={'linewidth':3, 'alpha': a})

ax.set_title(r'Posterior densities of $p$ from $M_t$')
ax.set_xlabel(r'$p$')

plt.show()

Figure 14: Posteriors of \(p\) from model \(M_t\)

This notebook demonstrates a simple way to compare models using leave one out cross-validation (loo) and a classic example from capture-recapture. This is just one way, however, to perform model comparison using PyMC. Perhaps a more effective solution for this problem would be placing a shrinkage prior on the \(\sigma\) parameters.

%load_ext watermark

%watermark -n -u -v -iv -w

Last updated: Mon Nov 25 2024

Python implementation: CPython
Python version       : 3.12.7
IPython version      : 8.29.0

numpy     : 1.26.4
pandas    : 2.2.3
seaborn   : 0.13.2
pymc      : 5.18.2
matplotlib: 3.9.2
arviz     : 0.20.0

Watermark: 2.5.0

References

King, Ruth, and SP2526632 Brooks. 2008. “On the Bayesian Estimation of a Closed Population Size in the Presence of Heterogeneity and Model Uncertainty.” Biometrics 64 (3): 816–24.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27: 1413–32.

Wilson, Ben, Philip S Hammond, and Paul M Thompson. 1999. “Estimating Size and Assessing Trends in a Coastal Bottlenose Dolphin Population.” Ecological Applications 9 (1): 288–300.