Sum of two custom likelihoods

Questions
1

I am trying to define custom likelihoods with density dist and add them up but it complains that the RVs cannot be added using a + operation (I also tried to use theano.tensor.add and it did not work). I would appreciate it if someone could help me understand how to compute this with pymc3.

def create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=[], status_col='status'):
    with pm.Model() as model:
        BoundedNormal = pm.Bound(pm.Normal, lower=-30, upper=0)
        BoundedNormal2 = pm.Bound(pm.Normal, lower=0, upper=0.5)
        BoundedNormal3 = pm.Bound(pm.Normal, lower=-2, upper=2)
        BoundedHalfNormal = pm.Bound(pm.HalfNormal, lower=0, upper=2)
        a0 = BoundedNormal('a0', mu=0, sd=30)
        a_sigma = BoundedHalfNormal('a_sigma', sd=10)
        a_Sex = BoundedNormal3('a_Sex',
                               mu=0,
                               sd=a_sigma
                              )
        b0 = BoundedNormal2('b0', sd=10)
        b_sigma = BoundedHalfNormal('b_sigma', sd=10)
        b_Sex = BoundedNormal3('b_Sex',
                               mu=0,
                               sd=b_sigma
                              )
        mu_a = a0 + a_Sex*events['Sex']
        mu_b = b0 + b_Sex*events['Sex']
        alpha_com, beta_com = [], []
        for curr_com in comorbidities:
            alpha_com.append(BoundedNormal3('a_{}'.format(curr_com), mu=0, sd=a_sigma))
            beta_com.append(BoundedNormal3('b_{}'.format(curr_com), mu=0, sd=b_sigma))
        num_events = trajectories_contributions.shape[2]
        print(num_events)
        for i in range(num_events):
            logalpha = mu_a + trajectories_contributions[:, 1, i] * alpha_com[0]
            beta = mu_b + trajectories_contributions[:, 1, i] * beta_com[0]
            logalpha += beta*trajectories_age[:, i]
            if i == 0 and num_events != 1:
                S = pm.DensityDist('S_{}', survival,
                                   observed={'event': np.zeros((events.shape[0],)),
                                             'logalpha': logalpha,
                                             'beta': beta,
                                             'time': trajectories_time[:, i]
                                             }
                                   )
            elif i == 0 and i == (num_events - 1):
                S = pm.DensityDist('S_{}'.format(i), 
                                   survival,
                                   observed={'event': events[status_col],
                                             'logalpha': logalpha,
                                             'beta': beta,
                                             'time': trajectories_time[:, i]
                                             }
                                    )
            elif i == (num_events - 1):
                print(1)
                t = trajectories_time[:, i] - trajectories_time[:, i-1]
                S1 = pm.DensityDist('S_{}'.format(i), 
                                    survival,
                                    observed={'event': events[status_col].values,
                                              'logalpha': logalpha,
                                              'beta': beta,
                                              'time': t
                                             }
                                    )
            else:
                t = trajectories_time[:, i] - trajectories_time[:, i-1]
                S1 = pm.DensityDist('S_{}'.format(i), 
                                                  survival,
                                                  observed={'event': events[status_col].values,
                                                                     'logalpha': logalpha,
                                                                      'beta': beta,
                                                                       'time': t
                                                                    }
                                                )

        S_tot = S + S1
                
    return model

def survival(event, logalpha, beta, time):
    log_h = logalpha + beta*time
    H = (np.exp(logalpha)/beta) * (np.exp(beta*time) - 1)
    return ((event * log_h) - H).sum()


model = create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=['B'], status_col='status')

Leads to the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-49-e05f6466f5f3> in <module>
----> 1 model = create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities=['B'], status_col='status')

<ipython-input-48-2db965ebc090> in create_model_comorbidity(events, trajectories_time, trajectories_contributions, trajectories_age, comorbidities, status_col)
     71                                                   )
     72                          )
---> 73         S_tot = S + S1
     74 
     75     return model

TypeError: unsupported operand type(s) for +: 'MultiObservedRV' and 'MultiObservedRV'

Thanks in advance,
Anurag

2

What’s strange to me is that S1 may not even be defined (if num_events == 1) then it looks like only S is defined?

In general, if you have conditionally independent observations like this, where

P[D|\theta] = P[D_1|\theta]P[D_2|\theta]P[D_3|\theta]\dots P[D_k|\theta]

you can use

  for dat_idx in data:
    pm.DensityDist('D%d' % dat_idx, my_lik_funcs[dat_idx], observed=my_args[dat_idx])

and the likelihoods are added up by default.

So: assuming there’s no problem with S1 not being defined in your script; just drop “S_tot = S + S1” and it will work just fine.

3

Thanks @charti. That was good to confirm/know.

Regarding your other question: Yes, this is survival with time dependent hazard functions (step function changes). So, they survive for some time under one hazard function and then they have a change of hazard function and may or may not have an event (i.e., they survived until final time T without undergoing the event or they have an event at t < T). In addition, there is stochasticity associated with the exact time each individual undergoes a change in the hazard function. I probablt could have probably used pm.switch as well but this was another way to code it up and I went with this.