I am trying to extend a model from the paper Default Bayes Factors for ANOVA designs from a single experiment, to multiple experiments. The motivation is to understand if there is something to be learnt across experiments.
The model is defined as follows:
\beta \sim \mathrm{Logistic}\left(0, 1\right) \\
\alpha^{*} \sim MVNormal(0, \Sigma) \\
\hat{p} = \beta + XQ_k \alpha^*\\
\log{\left( \frac{p}{1-p} \right)} = \hat{p} \\
y \sim Binomial\left(n, p\right)
The matrix Q_k is constructed as follows: we want the differences from the overall mean (in log odds scale) to each of our experiment arms to be constrained to sum to 0; this is to make up for the over-parameterized model we introduce by including an overall mean and a distance to each of the experiment arms. We can do this by projecting our design matrix into a space with one fewer dimension using a contrast matrix; that is an orthonormal set of contrasts that identify k -1 difference parameters in our model. After fitting our model under this parameterization, we can then project the results back to our original parameter space as the estimates and all tuples of posterior samples will sum to 0. We’ll assign \alpha to represent the k -length vector of fixed effects from the experimental arms from an overall mean \beta and define \alpha^* = Q_k^T \alpha , where \alpha is the k -1 vector representation of the effects. We can construct $Q_k$as follows:
$$\Omega_k = I_k - \frac{1}{k}J_k $$
where I_k is a k -sized identity matrix and J_k is a k-sized square matrix of 1s. We can then apply eigendecomposition to \Omega_k to obtain a k(k -1) matrix Q_k containing the k-1 eigenvectors corresponding to the non-zero eigenvalues, such that:
$$ \Omega_k = Q_k I_k Q_k^T $$
This has been coded up as
INPUT = {
"variation_id": ["A", "B", "C"],
"control_id": ["A", "A", "A"],
"conversions": [1000, 1100, 1200],
"participants": [10000, 10000, 10000],
}
data = pd.DataFrame(INPUT)
X = np.eye(number_of_variations)
COORDS = {"variations": np.arange(number_of_variations), "hidden": np.arange(number_of_variations-1)}
Sigma_k = np.eye(number_of_variations) - np.ones((number_of_variations, number_of_variations)) / number_of_variations
V, Q = np.linalg.eig(Sigma_k)
def extract_required_columns(Q, k):
"""
Extracts the first and the last k-2 columns from the eigenvector matrix Q.
Parameters:
- Q: numpy.ndarray, the eigenvector matrix from eigendecomposition.
- k: int, the original number of variations or the size of Q.
Returns:
- Q_k: numpy.ndarray, modified Q with the required columns.
"""
if k < 3:
return Q
columns_to_select = [0] + list(range(2, k))
Q_k = Q[:, columns_to_select]
return Q_k
Q_k = Q_k = extract_required_columns(Q, number_of_variations)
with pm.Model(coords=COORDS) as contrast_model:
beta = pm.Logistic("beta", 0, 1)
alpha_star = pm.MvNormal("alpha_star", mu=np.zeros(number_of_variations-1), cov=np.eye(number_of_variations-1), dims=("hidden"))
alpha = pm.Deterministic("alpha", pm.math.dot(Q_k, alpha_star.T), dims=("variations"))
p_hat = pm.Deterministic("p_hat", beta + pm.math.dot(X, alpha), dims=("variations"))
p = pm.Deterministic("p", pm.math.invlogit(p_hat), dims=("variations"))
y = pm.Binomial("y", n=data['participants'], p=p, observed=data['conversions'], dims=("variations"))
trace = pm.sample(1000, tune=n_tune, chains=4, return_inferencedata=True, random_seed=seed)
Now I want to expand this model to have an experiment dimension in the coords, but I can’t get the model to work. This is what I’ve got so far:
X = np.eye(number_of_variations)
COORDS = {"variations": np.arange(number_of_variations), "expt": np.arange(number_of_experiments), "hidden": np.arange(number_of_variations-1)}
Sigma_k = np.eye(number_of_variations) - np.ones((number_of_variations, number_of_variations)) / number_of_variations
V, Q = np.linalg.eig(Sigma_k)
def extract_required_columns(Q, k):
"""
Extracts the first and the last k-2 columns from the eigenvector matrix Q.
Parameters:
- Q: numpy.ndarray, the eigenvector matrix from eigendecomposition.
- k: int, the original number of variations or the size of Q.
Returns:
- Q_k: numpy.ndarray, modified Q with the required columns.
"""
if k < 3:
return Q
columns_to_select = [0] + list(range(2, k))
Q_k = Q[:, columns_to_select]
return Q_k
Q_k = Q_k = extract_required_columns(Q, number_of_variations)
with pm.Model(coords=COORDS) as contrast_model:
beta = pm.Logistic("beta", 0, 1, dims=("expt"))
alpha_star = pm.MvNormal("alpha_star", mu=np.zeros(number_of_variations-1), cov=np.eye(number_of_variations-1), dims=("expt", "hidden"))
alpha = pm.Deterministic("alpha", pm.math.dot(Q_k, alpha_star.T), dims=("expt", "variations"))
p_hat = pm.Deterministic("p_hat", beta + pm.math.dot(X, alpha), dims=("expt", "variations"))
p = pm.Deterministic("p", pm.math.invlogit(p_hat), dims=("expt", "variations"))
y = pm.Binomial("y", n=trials.T, p=p, observed=successes.T, dims=("expt", "variations"))
contrast_trace_multi_exp = pm.sample(1000, tune=n_tune, chains=4, return_inferencedata=True, random_seed=seed)
This model samples but gives the following error (which is indicative of something fundamental I think I’ve misunderstood):
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/var/folders/1y/pyzwj_b548377277blj2d6xm0000gp/T/ipykernel_59305/972904906.py in <module>
32 y = pm.Binomial("y", n=[data['participants']], p=p, observed=[data['conversions']], dims=("expt", "variations"))
33
---> 34 contrast_trace_multi_exp = pm.sample(1000, tune=n_tune, chains=4, return_inferencedata=True, random_seed=seed)
~/opt/anaconda3/lib/python3.9/site-packages/pymc/sampling/mcmc.py in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
787 # Packaging, validating and returning the result was extracted
788 # into a function to make it easier to test and refactor.
--> 789 return _sample_return(
790 run=run,
791 traces=traces,
~/opt/anaconda3/lib/python3.9/site-packages/pymc/sampling/mcmc.py in _sample_return(run, traces, tune, t_sampling, discard_tuned_samples, compute_convergence_checks, return_inferencedata, keep_warning_stat, idata_kwargs, model)
855 ikwargs: Dict[str, Any] = dict(model=model, save_warmup=not discard_tuned_samples)
856 ikwargs.update(idata_kwargs)
--> 857 idata = pm.to_inference_data(mtrace, **ikwargs)
858
859 if compute_convergence_checks:
~/opt/anaconda3/lib/python3.9/site-packages/pymc/backends/arviz.py in to_inference_data(trace, prior, posterior_predictive, log_likelihood, coords, dims, sample_dims, model, save_warmup, include_transformed)
510 return trace
511
--> 512 return InferenceDataConverter(
513 trace=trace,
514 prior=prior,
~/opt/anaconda3/lib/python3.9/site-packages/pymc/backends/arviz.py in to_inference_data(self)
429 """
430 id_dict = {
--> 431 "posterior": self.posterior_to_xarray(),
432 "sample_stats": self.sample_stats_to_xarray(),
433 "posterior_predictive": self.posterior_predictive_to_xarray(),
~/opt/anaconda3/lib/python3.9/site-packages/arviz/data/base.py in wrapped(cls)
63 if all((getattr(cls, prop_i) is None for prop_i in prop)):
64 return None
---> 65 return func(cls)
66
67 return wrapped
~/opt/anaconda3/lib/python3.9/site-packages/pymc/backends/arviz.py in posterior_to_xarray(self)
279 )
280 return (
--> 281 dict_to_dataset(
282 data,
283 library=pymc,
~/opt/anaconda3/lib/python3.9/site-packages/arviz/data/base.py in dict_to_dataset(data, attrs, library, coords, dims, default_dims, index_origin, skip_event_dims)
304 dims = {}
305
--> 306 data_vars = {
307 key: numpy_to_data_array(
308 values,
~/opt/anaconda3/lib/python3.9/site-packages/arviz/data/base.py in <dictcomp>(.0)
305
306 data_vars = {
--> 307 key: numpy_to_data_array(
308 values,
309 var_name=key,
~/opt/anaconda3/lib/python3.9/site-packages/arviz/data/base.py in numpy_to_data_array(ary, var_name, coords, dims, default_dims, index_origin, skip_event_dims)
253 # filter coords based on the dims
254 coords = {key: xr.IndexVariable((key,), data=np.asarray(coords[key])) for key in dims}
--> 255 return xr.DataArray(ary, coords=coords, dims=dims)
256
257
~/opt/anaconda3/lib/python3.9/site-packages/xarray/core/dataarray.py in __init__(self, data, coords, dims, name, attrs, indexes, fastpath)
446 data = _check_data_shape(data, coords, dims)
447 data = as_compatible_data(data)
--> 448 coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
449 variable = Variable(dims, data, attrs, fastpath=True)
450
~/opt/anaconda3/lib/python3.9/site-packages/xarray/core/dataarray.py in _infer_coords_and_dims(shape, coords, dims)
193 new_coords[dim] = var.to_index_variable()
194
--> 195 _check_coords_dims(shape, new_coords, dims)
196
197 return new_coords, dims
~/opt/anaconda3/lib/python3.9/site-packages/xarray/core/dataarray.py in _check_coords_dims(shape, coords, dims)
131 for d, s in v.sizes.items():
132 if s != sizes[d]:
--> 133 raise ValueError(
134 f"conflicting sizes for dimension {d!r}: "
135 f"length {sizes[d]} on the data but length {s} on "
ValueError: conflicting sizes for dimension 'expt': length 3 on the data but length 1 on coordinate 'expt'
Has anyone encountered something like this before?