Using this tutorial as a guide: Multi-output Gaussian Processes: Coregionalization models using Hamadard product — PyMC example gallery - I have implemented this method to analyse common effects reported across longitudinal studies. Following this and plotting as per the original authors suggestion I get separate estimates of the trajectory of the param for each included study - Im wondering; how can I combine these outputs to estimate the common effect between studies?
Code is as follows:
def get_icm(input_dim, kernel, W=None, kappa=None, B=None, active_dims=None):
“”"
This function generates an ICM kernel from an input kernel and a Coregion kernel.
“”"
coreg = pm.gp.cov.Coregion(input_dim=input_dim, W=W, kappa=kappa, B=B, active_dims=active_dims)
icm_cov = kernel * coreg # Use Hadamard Product for separate inputs
return icm_cov
with pm.Model() as model_contin:
# Priors
ell = pm.Gamma(“ell”, alpha=2, beta=0.5)
eta = pm.Gamma(“eta”, alpha=3, beta=1)
kernel = eta**2 * pm.gp.cov.ExpQuad(input_dim=2, ls=ell, active_dims=[0])
sigma = pm.HalfNormal(“sigma”, sigma=3)
# Get the ICM kernel
W = pm.Normal("W", mu=0, sigma=3, shape=(n_outputs, 2), initval=np.random.randn(n_outputs, 2))
kappa = pm.Gamma("kappa", alpha=1.5, beta=1, shape=n_outputs)
B = pm.Deterministic("B", pt.dot(W, W.T) + pt.diag(kappa))
cov_icm = get_icm(input_dim=2, kernel=kernel, B=B, active_dims=[1])
# Define a Multi-output GP
mogp = pm.gp.Marginal(cov_func=cov_icm)
y_ = mogp.marginal_likelihood("f", X, Y, sigma=sigma)
gp.conditional was then used to estimate new x values however this returns each constitutent study individually.
Apologies if the ability to combine these preferably with shrinking is obvious, however its not obvious to me at the time of writing.