Multinomial logistic regression on conjoint choice data (original) (raw)
Dear all, I am a new user of PyMC (previously using Stan). I am in the process of converting what I had before with Stan to PyMC. Now I am trying to write the PyMC model to estimate the partworth utilities based on a conjoint choice data.
Specifically, the input data would be a tensor with shape (N, C, K), where N is the total number of tasks; C is the number of alternatives/options within one task; K is the number of attributes of one alternative/option.
The process I am thinking is that:
thetais a matrix such that \theta \in \mathcal{R}^{K}, where K is the number of attribtues in each choice task.Xis a tensor, such that X \in \mathcal{R}^{N, C, K}, where N is the total number of tasks; C is the number of alternatives in each choice task.- For each task x_i\in X where x_i \in \mathcal{R}^{C, K}, we can compute the utility of each alternative as u_i = x_i \theta ; in this case, u_i\in \mathcal{R}^{C}
- Then I can just input u_i into a softmax function p_i = \text{softmax}(u_i) such that \sum p_i = 1
- Finally, these would be used as the probablity in a Multinomial distribution
Based on this, I wrote this following code:
with pm.Model() as multi_logit:
# priors:
theta = pm.Normal(name='theta', mu=0, sigma=1, shape=(1,K))
# use scan to compute xi*θ
U, updates = pytensor.scan(fn=lambda xi: xi.dot(theta),
sequences=X)
# softmax of each row
U_normed = pm.Deterministic('U_normed', pt.special.softmax(U, axis=0))
# distribution of the choice
y_obs = pm.Categorical('y_obs', p=U_normed, observed=y)
where X is the choice data;
array([[[1, 1, 1],
[1, 0, 0]],
[[1, 1, 1],
[0, 0, 0]],
[[1, 1, 0],
[1, 0, 1]],
...,
[[0, 1, 0],
[1, 0, 0]],
[[0, 0, 1],
[0, 1, 1]],
[[0, 0, 0],
[1, 0, 0]]])
and y is the observed choices:
array([0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, ..., 1, 1])
However, I was getting the following errors:
--------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[82], line 12
9 print(U_normed.shape)
11 # sample multinomial
---> 12 y_obs = pm.Categorical('y_obs', p=U_normed, observed=y)
File /opt/homebrew/lib/python3.11/site-packages/pymc/distributions/distribution.py:458, in Discrete.__new__(cls, name, *args, **kwargs)
455 if kwargs.get("transform", None):
456 raise ValueError("Transformations for discrete distributions")
--> 458 return super().__new__(cls, name, *args, **kwargs)
File /opt/homebrew/lib/python3.11/site-packages/pymc/distributions/distribution.py:310, in Distribution.__new__(cls, name, rng, dims, initval, observed, total_size, transform, *args, **kwargs)
307 elif observed is not None:
308 kwargs["shape"] = tuple(observed.shape)
--> 310 rv_out = cls.dist(*args, **kwargs)
312 rv_out = model.register_rv(
313 rv_out,
314 name,
(...)
319 initval=initval,
320 )
322 # add in pretty-printing support
File /opt/homebrew/lib/python3.11/site-packages/pymc/distributions/discrete.py:1117, in Categorical.dist(cls, p, logit_p, **kwargs)
1115 p_ = p_ / at.sum(p_, axis=-1, keepdims=True)
1116 p = at.as_tensor_variable(p_)
-> 1117 return super().dist([p], **kwargs)
File /opt/homebrew/lib/python3.11/site-packages/pymc/distributions/distribution.py:389, in Distribution.dist(cls, dist_params, shape, **kwargs)
387 ndim_supp = cls.rv_op(*dist_params, **kwargs).owner.op.ndim_supp
388 create_size = find_size(shape=shape, size=size, ndim_supp=ndim_supp)
--> 389 rv_out = cls.rv_op(*dist_params, size=create_size, **kwargs)
391 rv_out.logp = _make_nice_attr_error("rv.logp(x)", "pm.logp(rv, x)")
392 rv_out.logcdf = _make_nice_attr_error("rv.logcdf(x)", "pm.logcdf(rv, x)")
File /opt/homebrew/lib/python3.11/site-packages/pytensor/tensor/random/basic.py:1847, in CategoricalRV.__call__(self, p, size, **kwargs)
1828 def __call__(self, p, size=None, **kwargs):
1829 r"""Draw samples from a discrete categorical distribution.
1830
1831 Signature
(...)
1845
1846 """
-> 1847 return super().__call__(p, size=size, **kwargs)
File /opt/homebrew/lib/python3.11/site-packages/pytensor/tensor/random/op.py:289, in RandomVariable.__call__(self, size, name, rng, dtype, *args, **kwargs)
288 def __call__(self, *args, size=None, name=None, rng=None, dtype=None, **kwargs):
--> 289 res = super().__call__(rng, size, dtype, *args, **kwargs)
291 if name is not None:
292 res.name = name
File /opt/homebrew/lib/python3.11/site-packages/pytensor/graph/op.py:295, in Op.__call__(self, *inputs, **kwargs)
253 r"""Construct an `Apply` node using :meth:`Op.make_node` and return its outputs.
254
255 This method is just a wrapper around :meth:`Op.make_node`.
(...)
292
293 """
294 return_list = kwargs.pop("return_list", False)
--> 295 node = self.make_node(*inputs, **kwargs)
297 if config.compute_test_value != "off":
298 compute_test_value(node)
File /opt/homebrew/lib/python3.11/site-packages/pytensor/tensor/random/op.py:334, in RandomVariable.make_node(self, rng, size, dtype, *dist_params)
329 elif not isinstance(rng.type, RandomType):
330 raise TypeError(
331 "The type of rng should be an instance of either RandomGeneratorType or RandomStateType"
332 )
--> 334 shape = self._infer_shape(size, dist_params)
335 _, static_shape = infer_static_shape(shape)
336 dtype = self.dtype or dtype
File /opt/homebrew/lib/python3.11/site-packages/pytensor/tensor/random/op.py:201, in RandomVariable._infer_shape(self, size, dist_params, param_shapes)
199 param_batched_dims = getattr(param, "ndim", 0) - param_ndim_supp
200 if param_batched_dims > size_len:
--> 201 raise ValueError(
202 f"Size length is incompatible with batched dimensions of parameter {i} {param}:\n"
203 f"len(size) = {size_len}, len(batched dims {param}) = {param_batched_dims}. "
204 f"Size length must be 0 or >= {param_batched_dims}"
205 )
207 if self.ndim_supp == 0:
208 return size
ValueError: Size length is incompatible with batched dimensions of parameter 0 U_normed:
len(size) = 1, len(batched dims U_normed) = 2. Size length must be 0 or >= 2
I tested scan with following code, which worked out fine so I was confused. May I ask for help on this? Thank you!
X = pt.tensor3('X')
theta = pt.vector("theta")
results, updates = pytensor.scan(fn=lambda xi: xi.dot(theta),
sequences=X)
compute_utility = pytensor.function(inputs=[X, theta], outputs=results)
# test values
X_1 = np.array([[[1, 0], [0, 1]], [[1, 1], [0, 1]]], dtype=np.int32)
print(X_1)
theta = np.array([-1, 1], dtype=np.int32)
print(theta)
# test outcome
compute_utility(X_1, theta)
[[[1 0]
[0 1]]
[[1 1]
[0 1]]]
[-1 1]
array([[-1., 1.],
[ 0., 1.]])