Mostly out of curiosity, I would explore a random network approach.
Suppose a random graph $G$ with a fixed vertex set $v(G) = \{ \text{products} \}$ representing your set of products and a random edge set $E(G) \subseteq v(G) \times v(G)$. Your random edge set can be represented as a random adjacency matrix $A$ where $A_{ij}$ representing a Bernoulli variable for whether the $i$th product was paired with the $j$th product for some unique ID. Since this would be an undirected graph, you would really only need to have distinct parameters for the upper/lower triangle of this adjacency matrix (and exclude the diagonal since products always co-occur with themselves).
Here is a just a toy specification to get you started, but I mostly will defer development advice to Gelman et el. 2020.
$$\alpha_{ij} \sim \text{Exponential}(1)$$
$$\beta_{ij} \sim \text{Exponential}(1)$$
$$p_{ij} \sim \text{Beta}(\alpha_{ij}, \beta_{ij})$$
$$A_{ij} \sim \text{Bernoulli}(p_{ij})$$
You'll want to expand this first approximation to account for other aspects of the problem. For example:
- You may want to use a multilevel/mixed effect/hierarchical modelling approach for repeated customers or similarly accounting for a taxonomy on the types of products.
- There may be further covariance to model between the parameters due to cliques. McElreath 2023 provides some useful guidance on a related problem.
- There could also be some non-stationarity over time that turns this into a time series problem (if you can get dates) because customers sometimes change their habits.
After sampling from the posterior distribution, you can interpret you $p_{ij}$ as a tendency for the $i$th product to co-occur with the $j$th product. For something analogous to correlation, you could compute the posterior distribution of normalized independence gaps using the sampled values of $p_{ij}$ and marginalization (to obtain the product of the marginals $p_i p_j$).
