Hierarchical Switch Point Model

Questions
1

Hello guys,

I’m pretty new to PyMC3 and is planning on using for one of my thesis analysis.
I am trying to find out a hierarchical switch point model for body length measurement of my target animal.
There are a total of 100 individuals that I’ve measured across Age (1-20 years).
The idea behind was that there exist some individual heterogeneity but however governed by the with species distribution.

There is no problem at all with the pooled model, just as the pooled model in the radon example (https://docs.pymc.io/notebooks/multilevel_modeling.html)

There are several problems arising when I built the hierarchical model.

  1. The model couldn’t converge as I can see the traces remain horizontal.

    49
    49695×692 262 KB

  2. Should I be modelling intercept1 and slope1 in one line since they should be “switching” at the same age?
    and if so how do I link the 2 variables together?
    like this?
    [intercept,slope] = pm.math.switch(switchpoint[idx] <= df.Age, [intercept1,slope1][idx], [intercept2, slope2][idx])

  3. I tried using the advi method, but as I read more about it, it seems I shouldn’t have use the advi method since I’m not running any differential calculation?
    However the advi method did converge faster, although I fully understand that it is only an approximation.

with hierarchical_model_advi:
  mean_field = pm.variational.inference.fit(n=10000, method='advi')

Sorry it’s my 1st post and first trial with PyMC3. Much appreciated! Please excuse me if this post should be posted anywhere else.

Cheers,
Derek

### Hierarchical Switch Point
with pm.Model() as hierarchical_model:

    # Model Error
    sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)

    # SwitchPoint Distribution
    mean_sp = pm.TruncatedNormal('mean_sp', 5., sd=5., lower=0.)
    sigma_sp = pm.HalfCauchy('sigma_sp', beta=5)
    
    # Intercept Distribution
    intercept_i1 = pm.Uniform('intercept_i1',lower=10., upper=20.)
    intercept_i2 = pm.Uniform('intercept_i2',lower=10., upper=20.)
    sigma_i = pm.HalfCauchy('sigma_i', 5)
    
    # Slope Distribution
    slope_s1 = pm.Normal('slope_s1', mu=0., sd=1.)
    slope_s2 = pm.Normal('slope_s2', mu=0., sd=1.)
    sigma_s = pm.HalfCauchy('sigma_s', 5)

    # Priors for pre- and post-switch intercepts and slopes
    intercept1 = pm.Normal('intercept1', mu = intercept_i1, sd=sigma_i, shape=len(idx))
    intercept2 = pm.Normal('intercept2', mu = intercept_i2, sd=sigma_i, shape=len(idx))
    slope1 = pm.Normal('slope1', mu=slope_s1, sd=sigma_s,  shape=len(idx))
    slope2 = pm.Normal('slope2', mu=slope_s2, sd=sigma_s,  shape=len(idx))
    switchpoint = pm.TruncatedNormal('switchpoint', mean_sp, sd=sigma_sp, lower=0, upper=10, shape=len(idx))
     
    intercept = pm.math.switch(switchpoint[idx] <= df.Age, intercept1[idx], intercept2[idx])
    slope = pm.math.switch(switchpoint[idx] <= df.Age, slope1[idx], slope2[idx])
    
    # Data likelihood
    likelihood = pm.Normal('likelihood', mu = intercept + slope * df.Age, sd=sigma, observed=df.Body_length)

with hierarchical_model:
step1 = pm.NUTS([intercept, slope])
step2 = pm.NUTS([switchpoint])

hierarchical_trace = pm.sample(5000, step=[step1, step2], progressbar=True, 
                  njobs=4,  tune = 1000,  burn = 1000, thin = 10,
                  nuts_kwargs=dict(target_accept=.95))

pm.summary(hierarchical_trace)

2

I think doing this would fix all of your problems – it would certainly help to better parametrize the model. Rather than modeling two intercepts and two slopes, you have one intercept, one slope, and their respective deltas. Something like this:

base_intercept = pm.Normal('base_int', mu=15, sd=2)
base_slope = pm.Normal('base_slope', mu=0, sd=1)
delta_intercept = pm.Normal('delta_int', mu=0, sd=2)
delta_slope = pm.Normal('delta_slope', mu=0, sd=1)
switchpt = pm.Uniform('switch', 0., 10.)
like = pm.Normal('likelihood', mu = base_intercept + base_slope * df.Age + 
        (df.Age > switchpt) * (delta_intercept + delta_slope * df.Age), 
                              sd=sigma, observed=df.Body_length)
3

TypeError Traceback (most recent call last)
in ()
7 switchpt = pm.Uniform(‘switch’, 0., 10.)
8 like = pm.Normal(‘likelihood’, mu = base_intercept + base_slope * df.Age +
----> 9 (df.Age > switchpt) * (delta_intercept + delta_slope * df.Age),
10 sd=sigma, observed=df.Body_length)
11

/usr/local/lib/python3.6/dist-packages/pandas/core/ops.py in wrapper(self, other, axis)
877
878 with np.errstate(all=‘ignore’):
→ 879 res = na_op(values, other)
880 if is_scalar(res):
881 raise TypeError(‘Could not compare {typ} type with Series’

/usr/local/lib/python3.6/dist-packages/pandas/core/ops.py in na_op(x, y)
816 result = getattr(x, name)(y)
817 if result is NotImplemented:
→ 818 raise TypeError(“invalid type comparison”)
819 except AttributeError:
820 result = op(x, y)

TypeError: invalid type comparison

It didn’t work out at the mu = base_intercept + base_slope * df.Age + (df.Age > switchpt) * (delta_intercept + delta_slope * df.Age); with a invalid type comparison

4

Theano, the computational backend used by pymc3 does not handle pandas series, you have to get the underlying numpy array like this:
df.Age.values.

1 Like