OOS predictions of covariate with "out-of-range" missing values

[PaoloAndrich]

Hi all,

apologies for the rather obscure title for this topic. Part of the problem resides in the fact that I do not know exactly how to frame the issue that I am trying to tackle.
The (partially simplified) setup is overall pretty simple:

• I have a data set of measurements as a function of time (yhe value of the dependent variable should decrease with time).
• The time information is missing for a large proportion of the measurements. More importantly, I know that time is available only for measurements collected at early times. The measurements with missing covariate could instead refer to times that are longer than those available in the dataset.
• The goal would be to learn the relationship between the dependent variable and the covariate but also to impute the missing time values (I can reach the first goal without using the missing values although of course that comes with a loss of information).

The first thing I tried (I was a bit sceptical that it would work) was to develop a model that included the missing values in the typical PyMC3 way (using a normal random variable for the time with mean and standard deviation defined by hyperpriors and the masked time array as the observed data). Here’s the code for the model (where train is the dataset for training):

days = train["age_in_days"].values
days_missing = np.ma.masked_array(days, np.isnan(days))
road_cat_idx = train["road_cat"].cat.codes.values 
n_cat = len(np.unique(road_cat_idx))
deg_obs = train["qual_idx"].values
days_max = days_missing.max()
days_std = (days_missing/days_max) 
cat_names = train["road_cat"].cat.categories.values

coords = {
    "categories": cat_names,
    "obs_id": days_std,
}

with pm.Model(coords=coords) as betabinmod:
    
    #Hyperpriors
    a_mu = pm.Normal("a_mu", mu=-3, sd=0.1)
    a_sd = pm.HalfNormal("a_sd", sd=0.3)
    b_sd = pm.HalfNormal("b_sd", sd=1)
    phi_mu = pm.Exponential("phi_mu", lam=0.1)
    phi_sd = pm.HalfNormal("phi_sd", sd=0.1)
    age_mu = pm.Normal("age_mu", mu=1.7, sd=1.7)
    age_sd = pm.HalfNormal("age_sd", sd=1.7)

    #Priors
    a = pm.Normal("a", mu=0, sd=1, dims="categories")
    b = pm.HalfNormal("b", sd=2, dims="categories")            
    phi = pm.HalfNormal("phi", sd=0.1, dims="categories")
    age = pm.Normal("age", mu=age_mu, sd=age_sd,  observed=days_std, dims="obs_id")

    theta = pm.Deterministic("theta", var=(phi_mu + phi[road_cat_idx]*phi_sd) + 2, dims="obs_id")
    pbar = pm.Deterministic("pbar",
                                            pm.math.invlogit((a_mu + a[road_cat_idx]*a_sd) 
                                            + (b[road_cat_idx]*b_sd)*age),
                                            dims="obs_id"
                                            )
    
    # Likelihood
    road_qual = pm.BetaBinomial("road_qual",
                                                    n=200,
                                                    alpha=pbar * theta,
                                                    beta=(1 - pbar)*theta,
                                                    observed=deg_obs,
                                                    dims="obs_id"
                                                    )

I should mention that I can’t seem to be able to sample the posterior if I used a shared theano variable to replace the covariate. If I do that, when I start sampling it sets the initial value for age to NaN and PyMC3 does not return the imputation message warning. Any thoughts about why this is the case?
Regardless, with this approach the model seems to be unable to understand that the lower values of the dependent variable associated to missing time values should refer to larger times and instead assumes that they all fall within the same range of the observed values. Essentially the model squeezes all the data in the range of time values that I observe and thinks the dependent variable drops faster.

I guess this is not surprising given that I have no information for larger times but I was wondering if there could be any way to do this (stringent prior specifications?).

A different (and I guess less proper) approach I wanted to try was instead to split my two goals and to first use only the data for which the time information is available to sample the posterior. This works well and returns the expected behaviour of the dependent variable. Secondly, I would like to use a test dataset with missing time values to do inference and obtain samples for the time variable. I’m struggling to conceptualize if this is something the model should even be able to do. I implemented this approach following hints from previous discussions (shared variables sizes seem to not get updated correctly):

OOS predictions with missing input values

I’m having trouble generating out-of-sample predictions for a simple linear model with missing covariate data. I understand how to fit a model with missing covariate data using a numpy masked array: import numpy as np import pymc3 as pm # generate data with missing values X_train = np.random.randn(4,2) y_train = np.random.randn(4) mask = [[0, 0], [1, 0], [0, 1], [0, 0],] X_train = np.ma.masked_array(X_train, mask=mask) # fit model with pm.Model() as model: alpha = pm.Normal('alpha'…

Posterior Predictive Sampling in PyMC3 by Luciano Paz

Talk Abstract PyMC3 is great for inferring parameter values in a model given some observations, but sometimes we also want to generate random samples from the model as predictions given what we already inferred from the observed data. This kind of sampling is called posterior predictive sampling, and it can be very hard. The typical problems that show up are related to shape mismatches in hierarchical models, latent categorical values that aren’t correctly re-sampled or changing the shape of th…

and in particular I ended up constructing a model factory that I used to separately run the “training” step and the “imputation” step (please disregard the different dimension names).

def betabinom_model_factory(coords, days_std, days_max, categories, obs):
    with pm.Model(coords=coords) as model:
        days_shared = pm.Data("days_shared", days_std, dims="days_array")
        cat_shared = pm.Data("cat_shared", categories, dims="cat_array")
        obs_qual = pm.Data("obs_qual", obs, dims="obs_array")
        
        #Hyperpriors
        a_mu = pm.Normal("a_mu", mu=-3, sd=0.1)
        a_sd = pm.HalfNormal("a_sd", sd=0.3)
        b_sd = pm.HalfNormal("b_sd", sd=1)
        phi_mu = pm.Exponential("phi_mu", lam=0.1)
        phi_sd = pm.HalfNormal("phi_sd", sd=0.1)
        age_mu = pm.Normal("age_mu", mu=1.7, sd=1.7)
        age_sd = pm.HalfNormal("age_sd", sd=1.7)
        
        #Priors
        a = pm.Normal("a", mu=0, sd=1, dims="categories")
        b = pm.HalfNormal("b", sd=2, dims="categories")            
        phi = pm.HalfNormal("phi", sd=0.1, dims="categories")
        age = pm.Normal("age", mu=age_mu, sd=age_sd, observed=days_shared, dims="days_array")
        
        theta = pm.Deterministic("theta", var=(phi_mu + phi[cat_shared]*phi_sd) + 2, dims="days_array")
        pbar = pm.Deterministic("pbar",
                                                pm.math.invlogit((a_mu + a[cat_shared]*a_sd) 
                                                + (b[cat_shared]*b_sd)*age),
                                               dims="days_array"
                                               )
        
        #Likelihood
        road_qual = pm.BetaBinomial("road_qual",
                                                        n=200,
                                                       alpha=pbar * theta,
                                                       beta=(1 - pbar)*theta,
                                                       observed=obs_qual,
                                                       dims="obs_array"
                                                       )
    return model

This “works” in the sense that the sampling runs smoothly but it returns essentially the same value for age for all the samples (but different for each category). This is the case also when I use an approach more similar to that suggested by @junpenglao. Also, with this approach PyMC3 does not create an age_missing variable given that there are no missing values in the sampling of the posterior distribution.

Is what I am trying to do feasible?
I am still learning about PyMC3 capabilities and about Bayesian statistical approaches in general so apologies for the lack of clarity and for the likely blunders!

Thanks for your help!

[PaoloAndrich]

Hi all,

apologies for asking again for help. I was wondering if anyone might have any input particularly for what concerns the approach that follows the steps:

1. Train the model on a complete dataset
2. Generate predictions from the posterior distribution for entries with missing values

If anyone has thoughts on how best to approach this problem it would be super helpful.

Thanks,
Paolo