I ran an experiment where participants were randomly assigned to a control or treatment condition and in both conditions each participant was presented consecutively with three pairs of products where one product was more hedonic and the other more utilitarian. They had to choose one product in each of the three choice sets (1=hedonic, 0=utilitarian), that is they had to make three choices in total. Every choice set was a different product category.

I combined the three binary answers from every participant to form a "hedonic choice index" and I would like to examine the effect of the condition on this choice index. However, I would also like to check if the effect differs for the product categories.

I was wondering if the right way to do it is to run a factorial ANOVA like this:

df <- df %>% mutate(across(c(1,2,3,4,5), as.factor))

library(rstatix)
anova_test(data = df, choice_index~condition*product_category1+condition*product_category2+condition*product_category3, effect.size = "pes")

Here is my data:

 df <- structure(list(respondent_id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20), 
                      condition = c("control", "treatment", "control", "treatment", "control", "treatment", "control", "treatment", "control", "treatment", "control", "treatment", 
                                    "control", "treatment", "control", "treatment", "control", "treatment", "control", "treatment"), 
                      product_category1 = c(0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0), 
                      product_category2 = c(1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0), 
                      product_category3 = c(0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0),
                      choice_index = c(2, 2, 0, 4, 1, 4, 1, 4, 4, 0, 2, 2, 0, 4, 1, 4, 1, 4, 4, 0)), 
                 class = "data.frame", row.names = c(NA, -20L))