Logit Choice Set Up and Comparison

[alexkellogg]

Hi all,

I’m having some difficulties with a logit choice model and was hoping to get some help. First, I’m new to stan, so I was hoping to see if I set my model up correctly! Second, I was trying to show that truncating the priors would be bad for the model, but I’m not finding what I expected.

I created synthetic data for this example based on an economic model I’m using. Basically, the data consists of synth_N individuals who each make a choice between a pen and a mug. The valuation of the pen and the mug depends on several factors: 1) the intrinsic value (unknown and to be estimated, but for the synthetic data I fix these in synth_mug_u and synth_pen_u), 2) a parameter \lambda (to be estimated, drawn from a lognormal in the synthetic data), 3) which of the objects they are endowed with (this is observable data), and 4) a treatment condition q (observable). The overall valuations for each individual are computed as pen_utils and mug_utils below. Each individual then has a choice (observed data) which reflects whether the subject traded their endowed good for the alternative.

set.seed(19)
# Synthetic Data To Test Truncated vs Not
synth_N = 400
synth_exchange <- tibble(id = seq_len(synth_N),
                         synth_mug_u = 1,
                         synth_pen_u = rnorm(synth_N, 0.62, 0.1),
                         endowment = sample(c(1,2), synth_N, 0.5),
                         treatment = sample(c(0,1), synth_N, 0.5),
                         lambda = rlnorm(synth_N, 0.17, exp(-0.62)),
                         q =  0.9 * (1 - treatment) + 0.1 * (treatment),
                         pen_utils = (endowment == 2) * synth_pen_u + 
                           (endowment == 1)*
                           (q*synth_pen_u + (1 -q)*synth_mug_u 
                            + q*(1 - q)*(1-lambda)*(synth_mug_u + synth_pen_u)),
                         mug_utils = (endowment == 1) * synth_mug_u + 
                           (endowment == 2)*
                           (q*synth_mug_u + (1 -q)*synth_pen_u 
                            + q*(1 - q)*(1-lambda)*(synth_mug_u + synth_pen_u)),
                         logit_noise = rlogis(synth_N),
                         prefer_mug = mug_utils - pen_utils > logit_noise,
                         exchange = prefer_mug * (endowment ==2) + 
                           (1-prefer_mug)*(endowment == 1)
                         ) %>%
  group_by(endowment, treatment) %>% 
  mutate(group_id = cur_group_id()) %>%
  ungroup() 

synth_dat <- list(
  N = nrow(synth_exchange),
  exchange = synth_exchange$exchange,
  treatment = synth_exchange$treatment,
  q = synth_exchange$q,
  endowment = synth_exchange$endowment,
  id = synth_exchange$id,
  G = nrow(distinct(synth_exchange, group_id)),
  group_id = synth_exchange$group_id
)

With this synthetic data, my goal is to estimate the distribution of \lambda in the population, as well as the relative utility of the two goods (in this example, the distribution of synth_pen_u). I set my model up as a logit choice model, as shown below, where the probability of “exchanging” is a function of the parameters described in the synthetic data. I think that by fixing the mug_utils to 1, both parameters are identified (but could be wrong).

data{
  int N;
  int G;
  int group_id[N];
  int exchange[N];
  int treatment[N];
  int endowment[N];
  real q[N];
}
parameters{
  real<lower=0, upper= 6> lambda;
  real<lower=0, upper = 10> utils_pen;
}

model{
  vector[N] kr_utils_mug;
  vector[N] kr_utils_pen;
  vector[N] p;
  real utils_mug = 1;
  //Priors
  utils_pen ~ normal(utils_mug, 0.25);
  lambda ~ lognormal( 0.168 , exp(-0.62));

  for ( i in 1:N ) {
    kr_utils_mug[i]= utils_mug * (endowment[i] == 1) + 
      (q[i] * utils_mug + (1 - q[i]) * utils_pen + 
      q[i] * (1 - q[i]) * (1 - lambda) * (utils_mug + utils_pen)) * 
      (endowment[i] == 2) ;
     
    kr_utils_pen[i]= utils_pen * (endowment[i] == 2) +
      (q[i] * utils_pen + (1 - q[i])*utils_mug + 
      q[i] * (1 - q[i]) * (1 - lambda) * (utils_mug + utils_pen)) * 
      (endowment[i] == 1);
    // Probability of Exchanging
    p[i] = inv_logit(kr_utils_mug[i] - kr_utils_pen[i]) * (endowment[i] == 2) +
    inv_logit(kr_utils_pen[i] - kr_utils_mug[i]) * (endowment[i] == 1);
  }
  exchange ~ binomial(1, p);
}

generated quantities {
  vector[N] kr_utils_mug;
  vector[N] kr_utils_pen;
  vector[N] p;
  vector[N] log_lik;
  real utils_mug = 1;
  for ( i in 1:N ) {
     kr_utils_mug[i]= utils_mug * (endowment[i] == 1) + 
      (q[i] * utils_mug + (1 - q[i]) * utils_pen + 
      q[i] * (1 - q[i]) * (1 - lambda) * (utils_mug + utils_pen)) * 
      (endowment[i] == 2) ;
      
    kr_utils_pen[i]= utils_pen * (endowment[i] == 2) +
      (q[i] * utils_pen + (1 - q[i])*utils_mug + 
      q[i] * (1 - q[i]) * (1 - lambda) * (utils_mug + utils_pen)) * 
      (endowment[i] == 1);
      
     p[i] = inv_logit(kr_utils_mug[i] - kr_utils_pen[i]) * (endowment[i] == 2) +
     inv_logit(kr_utils_pen[i] - kr_utils_mug[i]) * (endowment[i] == 1);
    
  }
  for(i in 1:N) log_lik[i] = binomial_lpmf(exchange[i] | 1, p[i]);
}

Does this seem like the right way to estimate the \lambda and utils_pen distributions implied by the data? I can also do an aggregated version (groups are defined by endowment and treatment, and I can look at the sum of exchanges vs the total group size), but I don’t see why it would be different.

As for the model comparison – I was hoping to show that truncating \lambda to be larger than 1 (real<lower=1, upper= 6> lambda;) while maintaining the prior would lead to a drop in predictive ability. I do this by estimating the model with and without this change, and then using loo_compare. However, I’m finding that the truncated model may be slightly better (often not more than 1 SE, though). Why would that be, especially given that 1) \lambda < 1 exist in the synthetic data and 2) the non-truncated model is less restrictive, so if having no \lambda <1 did improve performance, the non-truncated model would be able to do that regardless?

Thanks!

Alex

[martinmodrak]

Sorry for not getting to you earlier, your question is relevant and well written.

I am not sure I understand the simulation code completely, but generally it is a good sign if your Stan model matches with a sensible code to simulate data for your problem. In fact, simulating data exactly from a model (starting from priors and moving “up”) is a good way to check that your priors (and potentially your whole model) are sensible - this is known as prior predictive check. Some more advice on how to develop model that is good fit for your data is in the Bayesian workflow preprint.

If I understand things correctly then mug_utils - pen_utils > logit_noise is the same as having rbinom(N_synth, size = 1, prob = inv_logit(mug_utils - pen_utils), (if yes, than this makes sense to me). But my math is far from fresh and confident in the area.

Some things I find weird (without understanding the economics side of things):

• The value 1 - lambda is constrained between negative infinity and 1. Why would the be a sensible multiplier for utility? (once again, maybe this is econ 101 that I am missing).
• I find it strange that the effect of treatment on the utilities (q) is assumed to be known - but maybe that’s available from the experiment design?

A minor notes on coding:

@alexkellogg

inv_logit(kr_utils_mug[i] - kr_utils_pen[i]) * (endowment[i] == 2) + ...

In Stan, those types of expression are likely to be somewhat faster to compute with an explicit if statement - the way you have coded means that Stan will always evaluate the expression and thus have to compute the derivative of inv_logit(kr_utils_mug[i] - kr_utils_pen[i]) even when endowment[i] is 1. (also Stan generally has much lower penalty for loops and branches than R)

binomial(1,p) is equivalent to bernoulli(p) (which could once again save you a few processor cycles).

Best of luck with your model and hope I’ve clarified at least something :-)

[alexkellogg]

Thanks for the feedback!

You’re correct about the binomial specification of choice with the probability of “success” as inv\_logit(mug\_utils - pen\_utils). I’ve since transitioned the model to an aggregated version of the problem, and re-written the if statements explicitly (which is much cleaner, ultimately).

Regarding the constraints on 1-lambda, lambda is taken as the loss aversion parameter which is assumed non-negative in the literature. I did have an upper bound on the parameter block so that 0\leq \lambda \leq 6, since a \lambda>5 is theoretically implausible, but perhaps this is not the right way to model this constraint. That said, the lognormal distribution set as a prior have almost no mass above 5.

The second question regarding the known treatment effect is a good one, and it helped me clarify what I’m hoping to achieve, so thanks! The way I’m thinking about it is that the model for utility is known (I’m assuming people behave according to this KR model of reference-dependence), and the observed choices follow this structural model (with logit noise) on the basis of 2 latent variables: \lambda and the relative utility of pens. Since the treatment is embedded into the model, I know how it interacts with choice behavior. So I can use experimental variation in treatments to help identify the distribution of these latent variables, the \lambda and “consumption utils”.

Thanks again for helping me understand the set up better. I’ve run into a new problem which I think warrants a separate thread as it has to do with a BFMI -Low warning and sample sizes in binomial models. Will edit a link it once I’ve submitted it in case it’s helpful to anyone.

Edit: Thread on recovery in a multi-level set up.