Confidence interval for descriptors

As @bienvenu suggested:

But for the report we should also provide "confidence intervals" for the descriptors of the extended models, or things like this. I think the main question a biologist can ask themselves are:

• what classes are sensitive to heterogeneities -> [See #119 ]
• how reliable are the descriptors

Here is a draft proposal for a function that do just that by performing n random expansions:

model = mpm.examples.orcinus_orca
classes = [0,1,3]
descriptors= ['R0', 'mean_age_population', 'lifetable_entropy']
intraclass_heterogeneity_ci(model, descriptors, classes, number)

Text report:

Model: orcinus_orca
Expanding Class 0, 50 times
	 R0                   2.013 ; 90% CI: (2.013 ; 2.014)
	 mean_age_population  26.22 ; 90% CI: (26.22 ; 26.22)
	 lifetable_entropy    0.8277 ; 90% CI: (0.8277 ; 0.8277)
Expanding Class 1, 50 times
	 R0                   2.013 ; 90% CI: (2.004 ; 2.022)
	 mean_age_population  26.22 ; 90% CI: (26.12 ; 26.46)
	 lifetable_entropy    0.8277 ; 90% CI: (0.8267 ; 0.8295)
Expanding Class 3, 50 times
	 R0                   2.013 ; 90% CI: (2.013 ; 2.013)
	 mean_age_population  26.22 ; 90% CI: (26.06 ; 26.47)
	 lifetable_entropy    0.8277 ; 90% CI: (0.8096 ; 0.8627)

Graph Report:

Here is the code:

Click to expand the code

model = mpm.examples.orcinus_orca
classes = [0,1,3]
descriptors=['R0', 'mean_age_population', 'lifetable_entropy']
number = 50
if descriptors is None:
    descriptors = ["R0",
                    "T_R0",
                           "T_a",
                       "mu1",
                       "life_expectancy",
                       "mean_age_population",
                       "damping_ratio",
                       "second_order_period",
                       "lifetable_entropy"]
if classes == 'each':
    classes = list(range(model.dim)) 
elif classes is None:
    classes = ['all'] 

report = {"expanded_models":{}, "data":{}}

for c in classes:
    if c == 'all':
        z = np.zeros(m.dim, dtype=int) + 2
    else:
        z = np.zeros(m.dim, dtype=int) + 1
        z[c] += 1
    if c not in report['expanded_models']:
        report['expanded_models'][c] = []
        report['data'][c] = {}

    for r in range(number):
        try:
            report['expanded_models'][c].append(mpm.expand.random_expansion(model, z))
        except mpm.expand.ExpansionNotFound as ex:
            print(z, ex)
        else:
            for f in descriptors:
                if f not in report['data'][c]:
                    report['data'][c][f] = []
                report['data'][c][f].append(mpm.expand._evaluate(report['expanded_models'][c][-1], f))

# Text Report 

print(f"Model: {model.metadata['ModelName']}")
for i,c in enumerate(classes):
    print(f"Expanding Class {c}, {number} times")
    for j,f in enumerate(descriptors):
        qlow = np.quantile(report['data'][c][f], q)
        qhigh = np.quantile(report['data'][c][f], 1-q)
        val = mpm.expand._evaluate(model, f)
        print(f'\t {f:20} {val:4.4} ; {1-2*q:2.0%} CI: ({qlow:4.4} ; {qhigh:4.4})')

# Graph Report 

M = len(classes)
N = len(descriptors)
fig, ax = plt.subplots(N,M, figsize=(M*5,N*3))
q = 0.05

def hist(ax, c, f, xlabel=True, ylabel=True):
        b = ax.hist(report['data'][c][f], alpha=0.5)
        qlow = np.quantile(report['data'][c][f], q)
        qhigh = np.quantile(report['data'][c][f], 1-q)
        val = mpm.expand._evaluate(model, f)
        if ylabel:
            ax.set(ylabel=f'Frequency of\n {f} values')
        ax.vlines(val, 0, np.max(b[0])+1, color='k', )
        ax.hlines((np.max(b[0])+1)/2, qlow, qhigh, color='k', ls='--')
        ax.text(0.95, 0.95,f'Model: {val:4.4}\n {1-2*q:2.0%} CI: ({qlow:4.4} ; {qhigh:4.4})',
                 horizontalalignment='right',
         verticalalignment='top',
         transform = ax.transAxes)

if N==1: 
    for i,c in enumerate(classes):
        try:
            hist(ax[i], c, f)
        except Exception as ex: 
            print(ex)
        
if M==1: 
    for j,f in enumerate(descriptors):
        try:
            hist(ax[j], c, f)
        except Exception as ex: 
            print(ex)
else:
    for i,c in enumerate(classes):
        for j,f in enumerate(descriptors):
            try:
                hist(ax[j,i], c, f, ylabel=(i==0))
                if j == 0:
                    ax[j,i].set(title=f"Expanding Class {c}, {number} times")
            except Exception as ex: 
                print(ex)