gaussian process with multiple outputs

[Jon Zelner]

Hi Folks - 

I have two-part question about using the multivariate gaussian process functions in stan. I hunted around for some examples and couldn't find any, so I thought I'd put this out here. The first is a general modeling question about when it is appropriate to use this approach and the second is a simple implementation question using multi_gp in stan (pg 387 of the stan 2.6 guide).

Modeling question: I have a spatial dataset of age-specific disease prevalence data. The data are each associated with a health care center in a large city. I want to model spatial variation in the age-prevalence relationship, using a GP based on the distances between health centers for spatial smoothing. In addition, there are spatially varying treatment effects that I want to include in the model as well. So, my model has 3 parameters: intercept, slope for age and the treatment effect. I expect all of these to vary similarly in space, so that neighborhoods with a similar intercept should definitely have a similar age effect and likely similar treatment effects. My question is, is this a situation in which it is ok/advisable to use a single kernel matrix to model all three of these coefficients, or at least the two that I'm sure should vary together?  Also, if I'm including the noise term on the diagonal of the kernel matrix, is it OK to not include an unstructured random effect term for each of these?

My implementation question related to how to specify the kernel matrix and vector of inverse scale parameters to sample all of these coefficients at once. My plan is to use the typical approach for a squared-exponential covariance function, allowing for noise on the diagonal. The inverse scales will be specified by a vector of length three, with the first element (representing the intercept) constrained to 1, so that the fitted noise term represents the unstructured variation in the area-level intercepts, and the other terms in the vector of inverse scales control the noise and covariances appropriately for the age & treatment coefficients. Is this correct?

Thanks!

-Jon

[Andrew Gelman]

Hi, I’m ccing Seth F. as this seems up his alley!

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

There are a few ways to model correlations in spatial dependence.

Firstly, if you assume that the effects vary spatially with similar

smoothing but not necessarily similar trends (for example, the

effects have the same length scale but they might not all be

positive in the same regions) then you can use a covariance 

function with common hyperparameters.  If that’s too restrictive

then you can make the hyperparameters hierarchical across

the sites.

On the other hand if the effects are expected to have similar

smoothing and similar trends then you can consider a

hierarchical GP, for example

y ~ GP(0, X)

y_age ~ GP(y, X)

y_prevalence ~ GP(y, X)

You can even go double hierarchical and have a hierarchical

GP with hierarchical hyperparameters.  Again, it really depends

on what prior relationships you want to encode in the model.

I’m not sure diagonal noise will completely replace random

effects at each site — especially if there may be nonspatially-

corrected effects.  I’d start with just the GP and then add

random effects to see if you get better PPCs.

[Herra Huu]

One possible solution is to use multi-output covariance functions (review article of the subject: http://arxiv.org/pdf/1106.6251). The basic idea is quite simple: say we have outputs X and Y, then we would jointly model them with a single GP, ie.  

(X,Y) ~ GP(0, K), where covariance function

K = [b11 * K_xx, b12 * K_xy,
       b21 * K_yx, b22 * K_yy]

and B = [b11, b12; b21, b22] is symmetric positive definite matrix capturing correlations between X and Y. For example, if B is the identity matrix, then they are independent etc. Unfortunately, the computational burden is quite large, this doesn't really scale up nicely. How many health centers there is in the dataset?

[Seth Flaxman]

Seconding everything Michael and Herra said. Since it can be tricky to
get even simple versions of what you describe working, I'd start as
simple as possible--no spatial model--and build from there. (Maybe
you've already done this?) Certainly email the list with code if you
get stuck.

I have to say that the existence of multi_gp is news to me; are there
examples somewhere of it being used? (It's weird to me that the mean
is fixed at 0...)

The model you're describing is interesting--do I understand correctly
that you have a linear model prevalence = intercept + b1 * age + b2 *
treatment and you want intercept, b1, and b2 to be spatially varying?
I was going to suggest above that you stick with a Gaussian
observation model if possible as that seems to help a lot with speed,
but now I'm not actual clear on what your observation model is, since
the GP is a prior on the coefficients, right? So are you using a
Binomial distribution for the prevalence?

An alternative you might consider if there are non-linearities in the
relationship between age and prevalence is that a GP would also be
appropriate there as well, and you could even do everything in one
model:

f(age, location, treatment) ~ GP(0, K_location * (1 + K_age + K_treatment))

Or more interpretable:

f(age, location, treatment) ~ GP(b0 + b1 * age + b2 * treatment,
K_location * (1 + K_age + K_treatment))

In this model f(.) is the log-prevalence. As noted above a Gaussian
observation model helps a lot with speed, so if you can accept the
approximation, I've had success using a "pseudo-likelihood" approach
with case counts k_i out of a population n_i with \hat p_i =
(k_i+1)/(n_i+2):

k_i/n_i ~ GP(invlogit(b0 + b1 * age_i), K_age + diag(\hat p_i (1-\hat
p_i) / n_i))

Seth

On Tue, May 12, 2015 at 9:48 AM, Herra Huu <heikki.uu...@gmail.com> wrote:
> One possible solution is to use multi-output covariance functions (review
> article of the subject: http://arxiv.org/pdf/1106.6251). The basic idea is

> quite simple: say we have inputs X and Y, then we would jointly model them

> with a single GP, ie.
>
> (X,Y) ~ GP(0, K), where covariance function
>

> K = [b11 * K_x, b12 * K_y,
> b21 * K_x, b22 * K_y]
>
> and B = [b11, b12; b21, b22] is positive definite matrix capturing

> correlations between X and Y. For example, if B is the identity matrix, then
> they are independent etc. Unfortunately, the computational burden is quite
> large, this doesn't really scale up nicely. How many health centers there is
> in the dataset?
>

[Jon Zelner]

Hi Everyone - 

Thanks for all of the helpful replies! Just to clarify, it's a binomial observation model. For now, I'm sticking with a single kernel matrix and focusing on just a spatially varying intercept and treatment effect and leaving the age effect fixed across units. I did try a version without the spatial smoothing and it was very noisy and PPCs were way off; these are pretty small units, some of which have only a few observations. They're all completely contiguous and form subsets of neighborhoods, so it seems like keeping the spatial info in there is worthwhile, even though it is indeed something of a pain.

The 'treatment' in this case is a household exposure of varying intensity, so the model also includes covariates for household-level attributes; which is a total bummer b/c otherwise the gaussian approximation seems perfect for this purpose!

Just a quick question for future reference on combining GPs with possibly differing dimensionalities: If K_age is a 30 x 30 kernel matrix with its own hyperparameters and K_space is a 100 x 100 matrix with its own hyperparameters, is the covariance between age i and j at locations k and l just the sum of K_age[i,j] and K_space[k,l]? Seems like it, just wanted to be sure.

Thanks!

-Jon

[Seth Flaxman]

> so the model also includes covariates for household-level attributes

But can't you put these into the mean function and still use the
Gaussian approximation? e.g.

f(age,location) ~ GP(invlogit(b0 + b1 * age + b2 * location + b3 *
household covariate 1 + b4 * household covariate 2), K)

> is the covariance between age i and j at locations k and l just the sum of K_age[i,j] and K_space[k,l]?

You could try an additive model, which is what you're describing, but
you could also try a multiplicative model instead. The additive model
is equivalent to two separate GPs:

f ~ GP(0, K1 + K2) is equivalent to f ~ GP(0, K1) + GP(0,K2)

But the multiplicative model might be better. Both should eventually
allow for some speedups--I haven't done much with speeding up additive
GP models but maybe someone else has code to share?--though it may
take some changes to core Stan; I'm working on code for the
multiplicative case (using Kronecker products) at the moment, will
hopefully share with the list soon.

Quick caution--if you have a multiplicative model:

k1(x,x') = v1 * exp( - (|x-x'|/l1)^2)
k2(y,y') = v2 * exp( - (|y-y'|/l2)^2)

Then you'll want to get rid of v2 for identifiability:

k1(x,x')*k2(y,y') = v1 * exp( - (|x-x'|/l1)^2) * exp( - (|y-y'|/l2)^2)

Seth

[Jon Zelner]

But can't you put these into the mean function and still use the
Gaussian approximation? e.g.

f(age,location) ~ GP(invlogit(b0 + b1 * age + b2 * location + b3 *
household covariate 1 + b4 * household covariate 2), K)

My concern about this is that the gaussian approximation would break down b/c some of the cell sizes will be very small, possibly with only one or two individuals in a given age x location x covariate combo. But I'm not sure if I'm misinterpreting the usage here b/c I've only worked with the zero-mean GPs in Stan. 

[Michael Betancourt]

I’d actually go the other way and think about using the GP to

model the slopes and intercepts as in a GLM.  So, for example,

you have binomial data which you could write in Stan as

y ~ binomial( inv_logit( f(covariates) ) );

y ~ binomial_logit( f(covariates) );

Usually we take f to be linear,

f = X * beta + alpha

with the beta, alpha hierarchical over groups.  But you could

also treat f as GP or just think about modeling beta and alpha

as GPs,

f = X * beta(spatial location) + alpha(spatial location);

alpha ~ GP(0, Sigma(hyperparams));

beta ~ GP(0, Sigma(hyperparams);

// add hyperprior here

And then if you had multiple effects you could do something like

f = X * (beta1(spatial location) + beta2(spatial location)) 

     + alpha1(spatial location) + alpha2(spatial location);

alpha1 ~ GP(y_alpha, Sigma(hyperparameters);

alpha2 ~ GP(y_alpha, Sigma(hyperparameters);

y_alpha ~ GP(0, Sigma);

// hyper priors

beta1 ~ GP(y_beta, Sigma(hyperparameters);

beta2 ~ GP(y_beta, Sigma(hyperparameters);

y_beta ~ GP(0, Sigma);

// hyper priors

This assumes that the mean effects and the group variations

have the same length scales, but you could easily generalize

this with different hyperparameters or hierarchical hyperparameters,

etc...