# How to model a sequence of binary choices: Generalized Estimating Equations (gee)?

For the sake of example, let’s say this is a "costumer research study" with a 3 by 3 factorial design, in which I am studying whether customers make purchases or not depending on various factors.

I have three independent variables: store decor (A), customer service (B), and price (C).

Each variable has three levels:

• A store decor: high level decor, medium level decor, no decor
• B customer service: friendly interaction, minimal interaction, no interaction
• C price: high price, medium price, low price

The binary outcome (x) is either making a purchase (1) or not (0).

This gives me a total of 27 conditions.

I have run the study on five participants using a virtual simulation in which they "visit" 27 different stores, one after the other - one store for each combination of the three independent variables.

These are the results for the five participants (I am using r):

A <- c(0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2)
B <- c(0,0,0,1,1,1,2,2,2,0,0,0,1,1,1,2,2,2,0,0,0,1,1,1,2,2,2)
C <- c(0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2,0,1,2)
1 <- c(1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0)
2 <- c(0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0)
3 <- c(0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,0,0)
4 <- c(0,0,0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0)
5 <- c(0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0)
data <- data.frame(A, B, C, x)


My question is, what is the best way to analyze this data?

I am interested in the effect of A, B and C on whether the customer makes a purchase or not.

I have found the following example online about the use of Generalized Estimating Equations (gee). The example uses the following dataset:

Initially I thought this dataset looked similar to mine if I considered "stores" in place of "visits".

But then I realized that this example would be good to study, say, the influence of the customer's gender on the purchase. My variables, instead, are about the store, so each store is characterized by a different combination of the IVs.

Should I "flip" the data then? Could I apply the same methodology?

It seems to me that you have a multi-level study with 2 random grouping factors: participant and store (presuming you are not only interested in the 5 participants and the 27 stores included in your study, but rather view the 5 participants as representative of a larger set of participants and the 27 stores as representative of a larger set of stores).

The 2 random grouping factors are fully crossed, since each participant "visits" each store.

For each participant-by-store combination, you obtain a single value of your binary outcome variable, which will tell you whether the participant will make the purchase (1) or will not make the purchase (0).

One possible way to analyze the data from your participants would be via a mixed effects binary logistic regression model which will include (at a minimum) a random participant effect and a random store effect. To that end, you would first have to store your data in long format, along these lines:

Participant  Store  Outcome  StoreDecor  CustomerService  Price
1         1       0         High         Friendly     Low
.         .       1         Medium       Minimal      High
.         .       etc.
.         .
1        27
.         .
.         .
.         .
5         1
.         .
.         .
.         .
5        27


To fit the a mixed effects binary logistic regression model to the data specified in this format (Data), you could use the glmer() function from the lme4 package:

library(lme4)

Data$$Participant <- factor(Data$$Participant)

Data$$Store <- factor(Data$$Store)

model <- glmer(Outcome ~ StoreDecor*CustomerService*Price +  # fixed effects
(1|Participant) + (1|Store), # random effects
family = binomial("logit"),
data = Data)


This model allows for a 3-way interaction between your predictor variables and uses dummy variable coding for capturing the fixed effects of these variables. However, note that you might want to use deviation coding in your setting, as explained here: http://www.mypolyuweb.hk/~sjpolit/coding_schemes.html and here: enter link description here

The above model will produce subject-specific fixed effects for your predictor variables (e.g., effects which refer to a 'typical' participant visiting a 'typical' store).

If you are interested in population-level fixed effects for your predictor variables, you could fit your model using the brm() function in the brms package. However, that will put you in a Bayesian setting.

• thanks a lot for your comprehensive explanation and for taking the time to write down the R code. I have run a pilot study in which I asked participants to go through the 27 scenarios 3 times (to gather some data quickly) but now I realize that may not be the best way to do it since now I have a repeated measure situation. What do you think? Also what's the best way to do a power analysis for such a model?
– Emy
Jul 15, 2019 at 23:50