3
$\begingroup$

I've been trying to gain a better understanding of factor potentials in PYMC. In reading this article by Cam Davidson-Pilon on Yhat, I got confused about how observed nodes are understood by PYMC.

Are they fixed in PYMC when observed=True? Or are they random like a stochastic node?

The factor potential is specified as:

@mc.potential
def censorfactor(obs=obs): 
    if np.any((obs + birth < 10)[lifetime_.mask] ): 
        return -100000
    else:
        return 0

The variable birth and lifetime_ are fixed numpy arrays. Hence, it seems, the only changeable variable is obs. But I thought that was fixed in PYMC because it is defined by:

obs = mc.Weibull( 'obs', alpha, beta, value = lifetime_, observed = True )

If obs is fixed then the potential function will return a log probability of -100000 every time.

I suspect I'm missing a fundamental point of MCMC here. Appreciate if someone can help me understand.

UPDATE

I did some additional research and it appears that observed nodes can indeed be (partly) stochastic. If you pass a masked array as observed data to a stochastic node, each of the masked values will be individually represented as a stochastic node.

Chris Fonnesbeck's tutorial from this year's Scipy conference covers this topic at the end of this IPython notebook. Here he describes how the posterior predictive distribution is used to estimate each of the missing data points:

$$p(\tilde{y}|y) = \int p(\tilde{y}|\theta) f(\theta|y) d\theta$$

The ability to impute missing data conditional on the data and the model parameters (simultaneously) a powerful technique. Plus we get traces, distributions, HDI etc for each missing value.

$\endgroup$
3
$\begingroup$

Quite right. Looking back, there are a few mistakes in that code. My recommendation: don't look at it. I'd like to go back and update it. Luckily, the code is in gists, so I can edit it directly from there.

Sorry about the confusion! To answer the question: observed nodes are fixed!

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I did some additional research. It appears that observed nodes can be stochastic when a masked array is passed as the observed value. It treats the masked observations as unknown stochastic nodes to be estimated as part of the model. $\endgroup$ – Mark Regan Sep 7 '14 at 18:59
  • $\begingroup$ Keen to hear your thoughts @cam-davidson-pilon: I'm trying to understand how the bayesian approach of imputing right censored observations compares to Aalen's additive or Cox Regression model (your lifelines library is awesome btw!). It seems the bayesian approach would naturally incorporate covariates of the data into the estimation of masked/missing observations. But the bayesian approach has no notion of hazard rate etc. What is your perspective when comparing lifelines (Aalens/Cox/KM) to the bayesian approach discussed in this article (ignoring mistakes etc). $\endgroup$ – Mark Regan Sep 7 '14 at 19:32
  • $\begingroup$ Ah, that's right! That's little-known PyMC thing. Thanks for bringing that back to my attention. Personally, I've moved away from Bayesian survival analysis for three reasons: i) computational difficulties - this post goes into them, and it can get worse. ii) I rarely need a point-estimate distribution when I perform survival analysis - I'm mostly interested in kaplan-meier curves. iii) I like the non-parametric nature of kaplan-meier estimates, and the semi-parametric form of Aalens. Bayesian SA requires a few too many strong assumptions. $\endgroup$ – Cam.Davidson.Pilon Sep 8 '14 at 0:22
  • $\begingroup$ Very interesting. Thanks for the response, BMH book and lifelines library. Really great learning resources! $\endgroup$ – Mark Regan Sep 8 '14 at 15:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.