# Sum to zero contraints in ANOVA designs (ABn testing)

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?