Add missing value support to glm

[fonnesbeck]

When attempting to run a glm, any missing values in the response variable results in an exception, which is appropriate, but the error message is not clear:

TypeError: reduce() of empty sequence with no initial value

For example, in the following model, you get this error when oxygen has one or more missing values.

glm = pm3.glm.glm
​
with pm3.Model() as oxygen_no_hx:

    formula = 'oxygen ~ ' + '+'.join(no_hx_vars)
    glm(formula, rsv_no_hx, family=pm3.glm.families.Binomial())
    oxygen_trace = pm3.sample(2000, pm3.NUTS())

Might be worth checking for that further up the chain and returning a more informative error.

[jsalvatier]

Does the missing values PR #712 automatically add this to GLM or is there more work necessary?

[fonnesbeck]

How would we get a masked array to glm? The data is expected as a DataFrame or 2D ndarray, is it not?

[jsalvatier]

nan values in a DataFrame should do it. Or a 2d masked array I think would also work.

[fonnesbeck]

Really? What about for missing values in the predictor variables? How does it impute without specifying a model?

[jsalvatier]

Good point, it obviously cannot do missing values without specifying a distributions.

[twiecki]

How would we do this manually in pymc3 (missing values in predictors)?

[jsalvatier]

See the Lasso Missing example.

Below mother is a predictor.

p_mother = Beta('p_mother', 1, 1)
mother_imp = Bernoulli('mother_imp', p_mother, observed=masked_values(mother_hs, value=-999))

[twiecki]

Ah, that's neat.

It doesn't sound straight forward to add it though. Internally we create a full (fixed) design-matrix out of the individual predictors which we then dot-multiply with the coefficients. So somehow columns of that design-matrix (those with missing values) would need to be RVs. Perhaps patsy could be convinced to take as input objects instead of numerical values as well. @njsmith, any ideas?

[njsmith]

Sorry, I don't really understand what's going on here well enough to be very useful.

Obviously you can call Patsy to get a dmatrix, and then scan that for nan values and do something if you find any.

Can you give an example of what it would look like to use this imputation functionality as a user, in your dream api?

[twiecki]

patsy.dmatrix('y ~ x1 + x2', data={'x1': [.1, .2, .3], 'x2': theano_vec_of_len_3})

Note that the last object supports indexing, has a dtype etc, but is definitely not a number.

Alternatively we could fill in the nans with the distribution that the user provides. Then patsy would be out of the equation.

[njsmith]

Ah, interesting! Right. Trying to think how this would work... In
principle, if the Theano vector is numeric and you tell patsy not to do
missing value detection, then I guess this could be done straightforwardly
-- use Theano.multiply to assemble interactions, and Theano.column_stack or
whatever the equivalent is to assemble the final design matrix as a Theano
matrix. The tricky bit I think is categorical variables, because in that
case Patsy needs to be able to see the actual values when setting up the
model (well, I guess you could require users to always specify them
explicitly, but...). Is there a way to immediately get a look at the data
when you're initially building the model, like actually compute the Theano
vector value e.g.? And then when running some data through the model to
generate a design matrix, categoricals involve indexing the contrast matrix
with the data matrix, which probably has done equivalent in Theano but I
don't know for sure.
.
The categorical case is really the interesting one, too, because presumably
you'd like to do your imputation in terms of the categorical variables, and
then let Patsy encode the imputed values to design matrix columns.
Postprocessing the design matrix is not so bad for numerics, but for
categoricals it's really ugly. (I guess it's not so pretty for more complex
numericals either -- consider "x1:impute(x2)" or "log(impute(x1))".)
On Jun 6, 2015 10:46 AM, "Thomas Wiecki" notifications@github.com wrote:

patsy.dmatrix('y ~ x1 + x2', data={'x1': [.1, .2, .3], 'x2':
theano_vec_of_len_3})

Note that the last object supports indexing, has a dtype etc, but is
definitely not a number.

Alternatively we could fill in the nans with the distribution that the
user provides. Then patsy would be out of the equation.

—
Reply to this email directly or view it on GitHub
#647 (comment).

[jsalvatier]

@njsmith yup, there's a way to immediately look at the data when you're initially building the model: variable.tag.test_value.

The way I would expect this to work is that patsy builds the design matrix (with the actual data, not theano tensors). Then we figure out which columns have missing values, and we pass those columns to the appropriate distribution (which might be automatically inferred from the kind of data or passed in) which uses the missing values logic and creates the imputed missing values. Then we rebuild the matrix and multiply it with the data. Maybe that has difficulties I don't see.

[njsmith]

@jsalvatier: From the docs it looks like tag.test_value is only available in a special debugging mode?

The more I think about it, the more I think that trying to do imputation on the design matrix directly is a bad idea. I gave examples of categoricals and complicated numerical expressions like x1:impute(x2) or log(impute(x1)) above -- an even nastier one that occurred to me is cases like fitting a polynomial model x1 + I(x1**2), where a single missing value in x1 will produce multiple missing values in the design matrix and -- crucially! -- you will get nonsensical results if you try imputing these values separately; whatever value you impute for the x1 column, you need to impute that same value squared for the I(x1**2) column. Similar issues arise for spline fits, factors that appear in multiple interactions like x1 + x2 + x1:x2, ...

Teaching patsy to pass accept and return theano tensors, and then doing imputation in the data space, is just going to work much much more robustly than trying to do imputation in the design matrix space.

To elaborate on what I'm thinking: Maybe you already know this, but probably it could be emphasized in the docs more, because a lot of people seem to miss it: patsy actually has two different phases in constructing a design matrix, which you can think of as "fit" and "predict". If you pass a formula to dmatrix it does both of these together so you can't see, but they're pretty distinct logically. In the first phase, patsy looks at your data to figure out which factors are numerical and which are categorical, and what the levels are for each categorical variable, and all the other stuff it needs to know in order to figure out what the output columns of your design matrix will be. (Another example: if you write center(x1) then it figures out at this stage what the mean of x1 is, and remembers it.) And at the end of this phase it has locked everything down, and knows exactly what set of operations it will apply to convert any input factors into an output design matrix. This knowledge is represented as a DesignMatrixBuilder object.

Then in the second phase, patsy takes a DesignMatrixBuilder and just applies the operations specified in there to whatever data you pass in, and outputs an actual design matrix.

While the first phase by definition needs to see some representative copy of the actual data, it would be totally reasonable to add a way to convert the second phase into a theano operation, so you end up with a theano tensor that represents the design matrix and which depends on the theano tensors that were passed in as factors.