Discussions of model comparison in hierarchical models can easily be found, but I have both a theoretical and a couple of practical questions.

My situation is that I have multiple subjects, `s`

, each being measured multiple times in paired measures `x`

and `y`

. I want to regress on the relationship of `x`

and `y`

per subject, but I’m estimating each of these from the multiple measures. That is, model is essentially this:

```
ys ~ N(beta_0 + beta_1*xs, ys_sigma)
y[i] ~ N(ys[s[i]], y_sigma[s[i]])
x[i] ~ N(xs[s[i]], x_sigma[s[i]])
```

I want to compare this model with a version without the `beta_0`

parameter (one that assumes `y=0`

for a subject for whom `x=0`

). Thus, I need to be able to calculate WAIC or LOO at the subject level.

To me, it seems that the log likelihood I want for this model comparison is evaluated for each subject and the likelihood of should include both how likely `y_s`

is given `x_s`

and how likely the `x`

and `y`

are given `y_s`

and `x_s`

.

**First question**: I haven’t seen this anywhere. Is it wrong?

My solution to getting this was to calculate the log likelihood I wanted after the sampling. I will give code snippets below for a simpler version of the problem with only one variable.

**Second question**: This code is really slow. Can it be speeded up?

**Third question**: This code causes `loo`

and `waic`

to return warnings about influential data points. Is there a way to get a better answer?

**Fourth question**: Is it possible to calculate the log likelihood that I want inside the trace sampling so it doesn’t need to be calculated twice?

**Code below**

I made a 1 dimensional version of this problem where multiple subjects each produce multiple `x`

values and my two models either have the grand mean of the subjects given as 0.0 or estimated as a parameter.

Here is random data generation (so that data structure is clear):

```
xs = stats.norm.rvs(loc=0, scale=1, size=num_subjects)
ns = stats.randint.rvs(low=4, high=8, size=num_subjects)
x_sigma = obs_sigma_sigma*stats.truncnorm.rvs(a=0, b=np.infty, size=num_subjects)
data = pd.DataFrame()
for this_s, (this_xs, this_n, this_x_sigma) in enumerate(zip(xs,ns,x_sigma)):
x = stats.norm.rvs(loc=this_xs, scale=this_x_sigma, size=this_n)
data = pd.concat([data, pd.DataFrame({"s": [this_s for xi in x], "x": x})])
```

I fit a hierarchical model which has either a mean at 0.0 or one which is estimated:

```
def simple_model(x,s,use_mu0):
s_idx, s_vals = pd.factorize(s, sort=True)
coords = {"subject": s_vals, "points":np.arange(len(x))}
with pm.Model(coords=coords) as model:
# Data
s_idx = pm.Data("s_idx", s_idx, dims="points", mutable=False)
# Priors
if use_mu0:
mu0 = pm.Normal('mu0', mu=0, sigma=1)
else:
mu0 = pm.Data('mu0', 0.0, mutable=False)
xs_sigma = pm.HalfNormal("xs_sigma", sigma=1.0)
x_sigma = pm.HalfNormal("x_sigma", sigma=0.5, dims="subject")
# Latent variables
xs = pm.Normal('xs', mu=mu0, sigma=xs_sigma, dims="subject")
# Likelihood
obs_x = pm.Normal('obs_x', mu=xs[s_idx], sigma=x_sigma[s_idx], observed=x, dims="points")
return model
```

and sample it:

```
idata = dict()
mdl = dict()
for use_m0 in [False, True]:
mdl[use_m0] = simple_model(data.x, data.s, use_m0)
with mdl[use_m0]:
idata[use_m0] = pm.sample(2000, tune=3000, cores=3, target_accept=0.95, return_inferencedata=True)
```

Until ultimately this is my *very slow* calculation of the log likelihood I think I want:

```
def calc_ll(idata, mu0):
x = idata.observed_data["obs_x"]
s_idx = idata.constant_data["s_idx"]
xs = idata.posterior["xs"]
x_sigma = idata.posterior["x_sigma"]
xs_sigma = idata.posterior["xs_sigma"]
mu0 = mu0.broadcast_like(xs)
xs_sigma = xs_sigma.broadcast_like(xs)
ll_xs = stats.norm(mu0, xs_sigma).logpdf(xs)
ll_x = np.zeros_like(ll_xs)
x_subj = x.groupby(s_idx)
coords = idata.posterior.coords
for c in coords["chain"]:
print(f'Chain: {c.item()}')
for d in coords["draw"]:
for s in range(len(x_subj)):
ll_x[c,d,s] = stats.norm(xs[c, d, s], x_sigma[c, d, s]).logpdf(x_subj[s]).sum()
ll = ll_xs + ll_x
return xr.DataArray(ll, coords=xs.coords, dims=xs.dims)
```

I can use compare to compare this with either WAIC or LOO

```
az.compare(idata, var_name="xs", ic="loo")
```

But I get complaints about quality of the fit of the pareto distribution.

Thanks for your help!

Opher