15 Application: A Quantal Response Equilibrium with discrete types
In this chapter, I would like to show you how to estimate a quantal response equilibrium model with discrete player types. This is in contrast to the Volunteer’s Dilemma chapter, where I implemented the continuous-type models studied in Goeree, Holt, and Smith (2017). The distinction between discrete and continuous types can become a bit blurred when there are a lot of types in a model. For example in auctions, we often assume a continuous distribution of values (i.e. types), but when we implement this in the lab we must have a discrete distribution, albeit with a lot of types. Here I am going to show you a model with just two types, so that you can hopefully see the way we solve for the model, then construct the likelihood.
15.1 Example dataset and models
This dataset is particularly special to me, because it is from my job market paper (Bland 2019)! In it, I was interested in whether choice bracketing was important in games. Choice bracketing is a phenomenon where people who face more than one decision to make partition their decisions into more than one choice set, and then make decisions for each element of these partitions sepatately. To fix ideas, Table 15.1 shows the two games in the baseline treatment that participants played. These games were played simultaneously and against the same opponenet. Payoffs from both games determined the participants’ earnings. However game theory and broad bracketing (the standard, “rational” assumption in economics) assumes that players will understand that the payoffs will be added up, and so will consider these two “games” to actually be one big game, shown in Table 15.2. Note that the actions in this game are the combinations of actions that a player could play in the two games in Table 15.1.
u1<- rbind(
c(10,10),c(100,0)
)
U1<-paste0(u1,", ",u1 |> t()) |> matrix(nrow = 2)
rownames(U1)<-c("A","B")
colnames(U1)<-c("A","B")
u2<- rbind(
c(56,56),c(160,0)
)
U2<-paste0(u2,", ",u2 |> t()) |> matrix(nrow = 2)
rownames(U2)<-c("C","D")
colnames(U2)<-c("C","D")
list(U1,U2) |>
kbl(caption = "The two games in the baseline treatment, narrow presentation") |>
kable_classic(full_width=FALSE) |>
add_header_above(c("Game 1","Game 2"))|
Game 1
|
Game 2
|
||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Ub<-rbind(
c(66,66,66,66),
c(170,10,170,10),
c(156,156,56,56),
c(260,100,160,0)
)
U<-paste0(Ub,", ",t(Ub)) |> matrix(nrow=4)
colnames(U)<-c("AC","AD","BC","BD")
rownames(U)<-c("AC","AD","BC","BD")
U |>
kbl(caption="The $4\\times 4$ game when players add up the payoffs of the two $2\\times 2$ games") |>
kable_classic(full_width=FALSE)| AC | AD | BC | BD | |
|---|---|---|---|---|
| AC | 66, 66 | 66, 170 | 66, 156 | 66, 260 |
| AD | 170, 66 | 10, 10 | 170, 156 | 10, 100 |
| BC | 156, 66 | 156, 170 | 56, 56 | 56, 160 |
| BD | 260, 66 | 100, 10 | 160, 56 | 0, 0 |
My research question was: do people see these as two separate games (i.e. Table 15.1), or do they integrate the payoffs together (i.e. Table 15.2)? I did this with a 3-treatment design that varied payoffs in a way that would not affect players’ actions under assumptions of broad or narrow bracketing. The way I achieved this was that in a second treatment I added 50 points to all the payoffs in Game 1 compared to the baseleine, and in a third treatment I added 50 points to all the payoffs in Game 2. If players bracketed narrowly then there should be no difference in choice frequencies in Game 1 between the baseline and Treatment 3, and in Game 2 between the baseline and Treatment 2. On the other hand, if people broadly bracketed there should be no difference in choice frequencies between Treatments 2 and 3. Unfortunately, I was able to reject both assumptions: the aggregate data didn’t look like narrow bracketing, but it didn’t look like broad bracketing either!
Fortunately, I wasn’t done analyzing me data. One of the structural analyses I did was estimate some Quantal Response Equilibrium models that assumed either broad bracketing, narrow bracketing, or a mixture of both.67 That way, I could ask which model best organized the data, rather than rejecting both assumptions against a generic alternative assumption of “something else”.
All of the quantal response equilibrium models I estimated assume a common choice precision \(\lambda>0\) and common risk-aversion parameter \(r>0\), but differ in their assumptions about how payoffs are calculated. I will start by walking you through the models that assume only broad or narrow bracketing, and then use these to build up to a model where there is a mixture of broad-bracketing and narrow-bracketing participants in the sample.
15.2 A note on replication
What I am about to show you is most certainty not an attempt to replicate the quantal response equilibrium analysis of Bland (2019). Aside from using very different priors for the model’s parameters (which I suspect wouldn’t have changed things too much), I am using a very differ technique to compute quantal response equilibrium. In particular, in my job market paper, I use the empirical payoff technique, which is a computational shortcut that avoids having the solve for the equilibrium at all. However as shown in Bland and Turocy (2025), this technique can produce unreliable results because it can be fitting data to a very different curve than the actual quantal response equilibrium. Instead, I am using the equilibrium correspondence approach, which is what I have been teaching you in the previous quantal response equilibrium chapters. Had I known about the issues of using the empirical payoff approach when I wrote my job market paper, I certainly would not have used it, especially since it is fairly easy to compute the equilibrium in this situation. All of this is to say, don’t expect what follows to be an exact replication of Bland (2019).
15.3 Three models that make different assumptions about bracketing
15.3.1 Broad bracketing only
This is probably the most standard Quantal Response Equilibrium model. Here we are going to derive the system of equations that characterize the equilibrium as a function of \(\lambda\) and \(r\), and then propose a method for solving for the equilibrium mixed strategies iteratively.
First, if \(U^r\) is the payoff matrix of the row player (in Table 15.2) raised element-wise to the power of \(r\), and \(\rho\) is vector of the logged mixed strategy, then we can write the equilibrium conditions in log differences as follows:
\[ \begin{aligned} 0=H_a(\rho;\lambda,r)&\doteq \rho_a-\rho_1-\lambda [U_a^r-U_1^r]\quad \text{for }a\in\{2,3,4\} \end{aligned} \]
and then we need one more constraint to ensure that the probabilities add up to one:
\[ 0=H_1(\rho;\lambda,r)=1-\sum_{a=1}^4\exp(\rho_a) \]
We can then write these in matrix notation as:
\[ H(\rho;\lambda,r)=A\rho-\lambda AU^r\exp(\rho)-B\exp(\rho)+c \]
where:
\[ \begin{aligned} A&=\begin{bmatrix} -1 & 1 & 0 & 0\\ -1 & 0 & 1 & 0\\ -1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 \end{bmatrix},\quad B=\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 1 & 1 & 1 & 1 \end{bmatrix},\quad c=\begin{pmatrix} 0\\ 0 \\ 0 \\ 1 \end{pmatrix} \end{aligned} \]
we can then solve for iterativelythe equilibrium (log) mixed strategy iteratively using Newton’s method:
\[ \rho^{t+1}=\rho^t-J^{-1}(\rho^t;\lambda,r)H(\rho^t;\lambda,r) \]
Where \(J\) is the Jacobian of the system with respect to \(\rho\):
\[ \begin{aligned} J(\rho;\lambda,r)&=A-(\lambda AU^r+B)\mathrm{diag\_matrix}(\exp(\rho)) \end{aligned} \]
For many games this won’t be enough to find an equilibrium, and you should then go and use the full predictor-corrector method described earlier in this book and in Bland and Turocy (2025), but in this instance I found that it worked well.
Once we have the log equilibrium mixed strategy \(\rho\), this is exactly the quantity we need to increment the likelihood: just add \(\rho_a\) to it every time a participant chooses action \(a\).
Here is the Stan program I wrote to estimate this model. One thing to note is that there is a \(\rho\) for every treatment (there were three treatments) becuase the payoffs change between treatments.
functions {
vector QRElp(
matrix U, real lambda,
data matrix A, data matrix B, data vector c,
data real tol, data int iter
) {
// returns the log QRE probabilities
int nactions = dims(U)[1];
// initialization
vector[nactions] rho = rep_vector(log(1.0/nactions),nactions);
for (tt in 1:iter) {
vector[nactions] H = A*rho-(lambda*A*U+B)*exp(rho)+c;
matrix[nactions,nactions] J = A-(lambda*A*U+B)*diag_matrix(exp(rho));
// change in rho
vector[nactions] drho = -J\H;
rho += drho;
//print(max(abs(H)));
//print(rho);
//print("---------------");
if (max(abs(H))<tol) {
break;
}
if (tt==iter && max(abs(H))>tol) {
print("maximum number of iterations reached without reaching tolerance");
print(max(abs(H)));
}
}
return rho;
}
}
data {
// payoff matrices for the 3 treatments
matrix[4,4] Ubroad[3];
int N; // number of observations
int treatment[N];
// coded as 1=AC, 2=AD, 3=BC, 4=BD
int action[N];
vector[2] prior_r;
vector[2] prior_lambda;
real<lower=0> tol;
int<lower=1> iter;
int<lower=0,upper=1> UseData;
int<lower=0,upper=1> contextual;
}
transformed data {
matrix[4,4] A_broad = [
[-1,1,0,0],
[-1,0,1,0],
[-1,0,0,1],
[0,0,0,0]
];
matrix[4,4] B_broad = [
[0,0,0,0],
[0,0,0,0],
[0,0,0,0],
[1,1,1,1]
];
vector[4] c_broad = [0,0,0,1]';
}
parameters {
// risk aversion
real<lower=0> r;
// choice precision
real<lower=0> lambda;
}
transformed parameters {
vector[4] rho[3];
for (tt in 1:3) {
matrix[4,4] U = pow(Ubroad[tt],r);
if (contextual==1) {
// contextual utility normalization
U = U/(U[4,1]-U[4,4]);
}
rho[tt] = QRElp(
U, lambda,
A_broad, B_broad, c_broad,
tol, iter
);
}
}
model {
// likelihood contribution
if (UseData==1) {
for (ii in 1:N) {
target += rho[treatment[ii]][action[ii]];
}
}
// priors
target += lognormal_lpdf(r | prior_r[1],prior_r[2]);
target += lognormal_lpdf(lambda | prior_lambda[1],prior_lambda[2]);
}
generated quantities {
vector[4] sigma[3];
for (tt in 1:3) {
sigma[tt] = exp(rho[tt]);
}
}15.3.2 Narrow bracketing only
Here we need to compute two equilibrium mixed strategies per treatment. One for Game 1 and one for Game 2. The steps for doing this are similar to that for the broad bracketing assumption, except that we have two log mixed strategies to keep track of, which I call \(\rho^1\) and \(\rho^2\). For each game, we can write the system of eqations describing the equilibrium as:
\[ \begin{aligned} H(\rho^g;\lambda,r)&=A\rho^g-\lambda A (U^g)^r\exp(\rho^g)-B\exp(\rho^g)+c,\quad g\in\{1,2\}\\ \text{where: }A&=\begin{bmatrix}-1 & 1\\ 0 & 0\end{bmatrix},\quad B = \begin{bmatrix} 0 & 0\\ 1 & 1 \end{bmatrix},\quad c = \begin{pmatrix}0 \\ 1\end{pmatrix} \end{aligned} \]
So it is basically the same process to solve for the quantal response equilibrium compared to the broad bracketing model.
One additional assumption I had to make here was about dependence between actions in the two games. While broad bracketing implies a probability distribution over all four combinations of actions, narrow bracketing only says something about the marginal strategy profiles. Here I assumed that they were independent, so the log mixed strategy would be (since multiplying probability levels is the same as adding logged probabilities):
\[ \rho = \begin{pmatrix} \rho^1_A+\rho^2_C\\ \rho^1_A+\rho^2_D\\ \rho^1_B+\rho^2_C\\ \rho^1_B+\rho^2_D\\ \end{pmatrix} \] Here is the Stan program I wrote to estimate this model:
functions {
vector QRElp(
matrix U, real lambda,
data matrix A, data matrix B, data vector c,
data real tol, data int iter
) {
// returns the log QRE probabilities
int nactions = dims(U)[1];
// initialization
vector[nactions] rho = rep_vector(log(1.0/nactions),nactions);
for (tt in 1:iter) {
vector[nactions] H = A*rho-(lambda*A*U+B)*exp(rho)+c;
matrix[nactions,nactions] J = A-(lambda*A*U+B)*diag_matrix(exp(rho));
// change in rho
vector[nactions] drho = -J\H;
rho += drho;
//print(max(abs(H)));
//print(rho);
//print("---------------");
if (max(abs(H))<tol) {
break;
}
if (tt==iter && max(abs(H))>tol) {
print("maximum number of iterations reached without reaching tolerance");
print(max(abs(H)));
}
}
return rho;
}
}
data {
// payoff matrices for the 3 treatments
matrix[2,2] Unarrow1[3];
matrix[2,2] Unarrow2[3];
int N; // number of observations
int treatment[N];
// coded as 1=AC, 2=AD, 3=BC, 4=BD
int action[N];
vector[2] prior_r;
vector[2] prior_lambda;
real<lower=0> tol;
int<lower=1> iter;
int<lower=0,upper=1> UseData;
int<lower=0,upper=1> contextual;
}
transformed data {
matrix[2,2] A_narrow = [
[-1,1],
[0,0]
];
matrix[2,2] B_narrow = [
[0,0],
[1,1]
];
vector[2] c_narrow = [0,1]';
}
parameters {
// risk aversion
real<lower=0> r;
// choice precision
real<lower=0> lambda;
}
transformed parameters {
vector[2] rho1[3];
vector[2] rho2[3];
vector[4] rho[3];
for (tt in 1:3) {
matrix[2,2] U1 = pow(Unarrow1[tt],r);
matrix[2,2] U2 = pow(Unarrow2[tt],r);
if (contextual==1) {
// contextual utility normalization
U1 = U1/(U1[2,1]-U1[2,2]);
U2 = U2/(U2[2,1]-U2[2,2]);
}
rho1[tt] = QRElp(
U1, lambda,
A_narrow, B_narrow, c_narrow,
tol, iter
);
rho2[tt] = QRElp(
U2, lambda,
A_narrow, B_narrow, c_narrow,
tol, iter
);
rho[tt] = to_vector(log(exp(rho1[tt])*exp(rho2[tt]')));
}
}
model {
// likelihood contribution
if (UseData==1) {
for (ii in 1:N) {
target += rho[treatment[ii]][action[ii]];
}
}
// priors
target += lognormal_lpdf(r | prior_r[1],prior_r[2]);
target += lognormal_lpdf(lambda | prior_lambda[1],prior_lambda[2]);
}
generated quantities {
vector[2] sigma1[3];
vector[2] sigma2[3];
vector[4] sigma[3];
for (tt in 1:3) {
sigma1[tt] = exp(rho1[tt]);
sigma2[tt] = exp(rho2[tt]);
sigma[tt] = exp(rho[tt]);
}
}15.3.3 A mixture of broad and narrow bracketing
When we have more than one thing that people could be doing, it is possible (and often likely) that all of the things are important. I therefore also estimate a quantal response equilibrium model that assumes fraction \(\phi\in(0,1)\) of participants narrowly bracket, and fraction \(1-\phi\) broadly bracket. That is, we have two types of player in our game. If we weren’t doing quantal response equilibrium, we would be looking to solve a Bayesian Nash equilibrium. For quantal response equilibrium, while conceptually this is not too hard to implement, some care needs to be taken when writing down the equilibrium conditions. For this, we will need a set of equations for the narrow-bracketing type, a set for the broad-bracketing type, and a description of how the types’ mixed strategies relate to the aggregate mixed strategy (when we integrate out the types). By the “aggregate mixed strategy” here I mean the average choice probabilities implied by the mixing probability \(\phi\) and the mixed strategies played by each type. In this kind of game this is the quantity players need to know in order to evaluate their expected utilities (irrespective of whether they are broad or narrow bracketers).
Let \(\rho\) be the log probabilities of the broad and narrow actions:
\[ \rho=\begin{pmatrix} \rho_{AC}^B & \rho_{AD}^B & \rho_{BC}^B & \rho_{BD}^B & \rho_{A}^N & \rho_{B}^N & \rho_{C}^N & \rho_{D}^N \end{pmatrix}^\top \]
We can write the aggregate mixed strategy (in levels, rather than logs) from a broad bracketer’s perspective as:
\[ \begin{aligned} \sigma^{B*}&=\begin{pmatrix} (1-\phi)\exp(\rho_{AC}^B)+\phi\exp(\rho_A^N+\rho_C^N)\\ (1-\phi)\exp(\rho_{AD}^B)+\phi\exp(\rho_A^N+\rho_D^N)\\ (1-\phi)\exp(\rho_{BC}^B)+\phi\exp(\rho_B^N+\rho_C^N)\\ (1-\phi)\exp(\rho_{BD}^B)+\phi\exp(\rho_B^N+\rho_D^N) \end{pmatrix} \end{aligned} \]
In order to evaluate the Jacobian, we need the derivative of this, which is:
\[ \begin{aligned} \frac{\partial \sigma^{B*}}{\partial \rho^\top}&=\begin{bmatrix} (1-\phi)\exp(\rho^B_{AC}) & 0 & 0 & 0 & \phi\exp(\rho^N_A+\rho^N_C) & 0 &\phi\exp(\rho^N_A+\rho^N_C) & 0\\ 0 & (1-\phi)\exp(\rho^B_{AD}) & 0 & 0 & \phi\exp(\rho^N_A+\rho^N_D) & 0 & 0 & \phi\exp(\rho^N_A+\rho^N_D)\\ 0 & 0 & (1-\phi)\exp(\rho^B_{BC}) & 0 & 0 & \phi\exp(\rho^N_B+\rho^N_C)& \phi\exp(\rho^N_B+\rho^N_C)& 0 \\ 0 & 0 & 0 & (1-\phi)\exp(\rho^B_{BD}) & 0 & \phi\exp(\rho^N_B+\rho^N_D)& 0 & \phi\exp(\rho^N_B+\rho^N_D) \end{bmatrix} \end{aligned} \]
Then we can write the aggregate mixed strategy from a narrow bracketer’s perspective as:
\[ \begin{aligned} \sigma^{N*}&=\begin{pmatrix} \sigma^{B*}_{AC}+\sigma^{B*}_{AD}\\ \sigma^{B*}_{BC}+\sigma^{B*}_{BD}\\ \sigma^{B*}_{AC}+\sigma^{B*}_{BC}\\ \sigma^{B*}_{AD}+\sigma^{B*}_{BD} \end{pmatrix} \end{aligned} \]
Then to get the Jacobian we can write:
\[ \begin{aligned} \frac{\partial \sigma^{N*}}{\partial \rho^\top}&=\frac{\partial \sigma^{N*}}{\partial {\sigma^{B*}}^\top}\frac{\partial \sigma^{B*}}{\partial \rho^\top},\quad \text{where }\frac{\partial \sigma^{N*}}{\partial {\sigma^{B*}}^\top} =\begin{bmatrix} 1 & 1 & 0 & 0\\ 0 & 0 & 1& 1\\ 1& 0 & 1 & 0\\ 0 & 1 & 0 & 1 \end{bmatrix} \end{aligned} \]
We can then write the system of equations describing the equilibrium as:
\[ \begin{aligned} H_{1:4}(\rho;r,\phi,\lambda)&=A^B\rho^B-\lambda A^B(U^B)^r\sigma^{B*}-B^B\exp(\rho^B)+c^B\\ H_{5:6}(\rho;r,\phi,\lambda)&=A^N\rho^N_1-\lambda A^N(U_1^N)^r\sigma_1^{N*}-B^N\exp(\rho_1^N)+c^N\\ H_{7:8}(\rho;r,\phi,\lambda)&=A^N\rho^N_2-\lambda A^N(U_2^N)^r\sigma_2^{N*}-B^N\exp(\rho_2^N)+c^N \end{aligned} \]
Where \(A^B\) is the \(A\) matrix defined for the broad-only model, \(A^N\) is the \(A\) matrix defined for the narrow-only model, and so on.
Now that we can set up and solve the system of equations characterizing the quantal response equilibrium mixed strategies, we need to take these mixed strategies and compute the likelihod. Here, becuase participants made more than one decision in the experiment, we could make (at least) two different assumptions about how the likelihood is calculated. These are (1) the dichotomous mixture model assumption, and (2) the toolbox mixture assumption. The dichotomous assumption in this model is that participants are either broad or narrow bracketers, and always make their decisions according to that type. In this case, we compute the log-likelihood contribution for participant \(i\) as follows:
\[ \log p(y_i\mid \rho^B,\rho^N,\phi)=\log\left(\phi\exp\left(y_i^\top\rho^N\right)+(1-\phi)\exp\left(y_i^\top\rho^B\right)\right) \]
where \(y_i\) is a vector of counts of the four combinations of actions.
On the other hand, the toolbox assumption states that participants bracket narrowly for each decision with probability \(\phi\), and so we write the likelihood as follows:
\[ \log p(y_i\mid \rho^B,\rho^N,\phi) =y_i^\top\log\left(\phi\exp(\rho^N)+(1-\phi)\exp(\rho^B)\right) \]
Here is the Stan program I wrote to estimate the dichotomous model:
functions {
vector QRElp(
matrix Ub, matrix U1n, matrix U2n,
real lambda, real mix_n,
data int iter, data real tol
) {
matrix[4,4] Ab = [
[-1,1,0,0],
[-1,0,1,0],
[-1,0,0,1],
[0,0,0,0]
];
matrix[4,4] Bb = [
[0,0,0,0],
[0,0,0,0],
[0,0,0,0],
[1,1,1,1]
];
vector[4] cb = [0,0,0,1]';
matrix[2,2] An = [
[-1,1],
[0,0]
];
matrix[2,2] Bn = [
[0,0],
[1,1]
];
vector[2] cn = [0,1]';
matrix[2,4] Dsigma1n = [
[1,1,0,0],
[0,0,1,1]
];
matrix[2,4] Dsigma2n = [
[1,0,1,0],
[0,1,0,1]
];
// initialization
vector[8] rho;
rho[1:4] = rep_vector(log(0.25),4);
rho[5:8] = rep_vector(log(0.50),4);
// START LOOP HERE-----------------------------------------------------------
for (tt in 1:iter) {
// aggregate mixed strategy, from the broad bracketers' perspective
vector[4] SIGMAb = (1-mix_n)*exp(rho[1:4]);
SIGMAb[1] += mix_n*exp(rho[5]+rho[7]);
SIGMAb[2] += mix_n*exp(rho[5]+rho[8]);
SIGMAb[3] += mix_n*exp(rho[6]+rho[7]);
SIGMAb[4] += mix_n*exp(rho[6]+rho[8]);
// aggregate mixed strategy, from the narrow bracketers' perspective
vector[2] SIGMA1n;
SIGMA1n[1] = SIGMAb[1]+SIGMAb[2];
SIGMA1n[2] = SIGMAb[3]+SIGMAb[3];
vector[2] SIGMA2n;
SIGMA2n[1] = SIGMAb[1]+SIGMAb[3];
SIGMA2n[2] = SIGMAb[2]+SIGMAb[4];
vector[4] rhob = rho[1:4];
vector[2] rho1n = rho[5:6];
vector[2] rho2n = rho[7:8];
vector[8] H;
H[1:4] = Ab*rhob-lambda*Ab*Ub*SIGMAb-Bb*exp(rhob)+cb;
H[5:6] = An*rho1n-lambda*An*U1n*SIGMA1n-Bn*exp(rho1n)+cn;
H[7:8] = An*rho2n-lambda*An*U2n*SIGMA2n-Bn*exp(rho2n)+cn;
matrix[4,8] DSIGMAb = [
[(1-mix_n)*exp(rhob[1]),0, 0, 0, mix_n*exp(rho1n[1]+rho2n[1]),0, mix_n*exp(rho1n[1]+rho2n[1]),0 ],
[0, (1-mix_n)*exp(rho[2]),0, 0, mix_n*exp(rho1n[1]+rho2n[2]),0, 0, mix_n*exp(rho1n[1]+rho2n[2])],
[0, 0, (1-mix_n)*exp(rho[3]),0, 0 , mix_n*exp(rho1n[2]+rho2n[1]),mix_n*exp(rho1n[2]+rho2n[1]),0 ],
[0, 0, 0, (1-mix_n)*exp(rho[4]),0, mix_n*exp(rho1n[2]+rho2n[2]),0, mix_n*exp(rho1n[2]+rho2n[2])]
];
matrix[8,8] J = rep_matrix(0.0,8,8);
J[1:4,1:4] += Ab-Bb*diag_matrix(exp(rho[1:4]));
J[5:6,5:6] += An-Bn*diag_matrix(exp(rho[5:6]));
J[7:8,7:8] += An-Bn*diag_matrix(exp(rho[7:8]));
J[1:4,] += -lambda*Ab*Ub*DSIGMAb;
J[5:6,] += -lambda*An*U1n*Dsigma1n*DSIGMAb;
J[7:8,] += -lambda*An*U2n*Dsigma2n*DSIGMAb;
vector[8] drho = -J\H;
rho += drho;
if (max(abs(H))<tol) {
break;
}
if (tt==iter && max(abs(H))>tol) {
print("Maximum iterations reached without achieveing tolerance");
}
}
return rho;
}
}
data {
// payoff matrices for the 3 treatments
matrix[4,4] Ubroad[3];
matrix[2,2] Unarrow1[3];
matrix[2,2] Unarrow2[3];
int N; // number of observations
int id[N];
int nparticipants;
int treatment[N];
// coded as 1=AC, 2=AD, 3=BC, 4=BD
int action_broad[N];
vector[2] prior_r;
vector[2] prior_lambda;
vector[2] prior_mix;
real<lower=0> tol;
int<lower=1> iter;
int<lower=0,upper=1> UseData;
int<lower=0,upper=1> contextual;
}
transformed data {
vector[4] rho0 = rep_vector(log(1.0/4.0),4);
int action_narrow1[N];
int action_narrow2[N];
for (ii in 1:N) {
action_narrow1[ii] = action_broad[ii]==1 || action_broad[ii]==2 ? 1 : 2;
action_narrow2[ii] = action_broad[ii]==1 || action_broad[ii]==3 ? 1 : 2;
}
}
parameters {
// risk aversion
real<lower=0> r;
// choice precision
real<lower=0> lambda;
real<lower=0,upper=1> mix_narrow;
}
transformed parameters {
vector[8] rho[3];
vector[4] sigma_broad[3];
vector[2] sigma_narrow1[3];
vector[2] sigma_narrow2[3];
for (tt in 1:3) {
matrix[4,4] Ub = pow(Ubroad[tt],r);
matrix[2,2] U1n = pow(Unarrow1[tt],r);
matrix[2,2] U2n = pow(Unarrow1[tt],r);
if (contextual==1) {
// contextual utility normalization
Ub = Ub/(Ub[4,1]-Ub[4,4]);
U1n = U1n/(U1n[2,1]-U1n[2,2]);
U2n = U2n/(U2n[2,1]-U2n[2,2]);
}
rho[tt] = QRElp(
Ub,U1n, U2n,
lambda,mix_narrow,
iter,tol
);
sigma_broad[tt] = exp(rho[tt][1:4]);
sigma_narrow1[tt] = exp(rho[tt][5:6]);
sigma_narrow2[tt] = exp(rho[tt][7:8]);
}
}
model {
// likelihood contribution
if (UseData==1) {
vector[nparticipants] ll_narrow = rep_vector(log(mix_narrow),nparticipants);
vector[nparticipants] ll_broad = rep_vector(log(1-mix_narrow),nparticipants);
for (ii in 1:N) {
ll_broad[id[ii]] += log(sigma_broad[treatment[ii]])[action_broad[ii]];
ll_narrow[id[ii]] += log(sigma_narrow1[treatment[ii]])[action_narrow1[ii]];
ll_narrow[id[ii]] += log(sigma_narrow2[treatment[ii]])[action_narrow2[ii]];
}
for (ii in 1:nparticipants) {
target += log_sum_exp(ll_broad[ii],ll_narrow[ii]);
}
}
// priors
target += lognormal_lpdf(r | prior_r[1],prior_r[2]);
target += lognormal_lpdf(lambda | prior_lambda[1],prior_lambda[2]);
target += beta_lpdf(mix_narrow | prior_mix[1],prior_mix[2]);
}
generated quantities {
vector[4] sigma[3];
for (tt in 1:3) {
sigma[tt] = exp(rho[tt]);
}
}And is the Stan program I wrote to estimate the toolbox model:
functions {
vector QRElp(
matrix Ub, matrix U1n, matrix U2n,
real lambda, real mix_n,
data int iter, data real tol
) {
matrix[4,4] Ab = [
[-1,1,0,0],
[-1,0,1,0],
[-1,0,0,1],
[0,0,0,0]
];
matrix[4,4] Bb = [
[0,0,0,0],
[0,0,0,0],
[0,0,0,0],
[1,1,1,1]
];
vector[4] cb = [0,0,0,1]';
matrix[2,2] An = [
[-1,1],
[0,0]
];
matrix[2,2] Bn = [
[0,0],
[1,1]
];
vector[2] cn = [0,1]';
matrix[2,4] Dsigma1n = [
[1,1,0,0],
[0,0,1,1]
];
matrix[2,4] Dsigma2n = [
[1,0,1,0],
[0,1,0,1]
];
// initialization
vector[8] rho;
rho[1:4] = rep_vector(log(0.25),4);
rho[5:8] = rep_vector(log(0.50),4);
// START LOOP HERE-----------------------------------------------------------
for (tt in 1:iter) {
// aggregate mixed strategy, from the broad bracketers' perspective
vector[4] SIGMAb = (1-mix_n)*exp(rho[1:4]);
SIGMAb[1] += mix_n*exp(rho[5]+rho[7]);
SIGMAb[2] += mix_n*exp(rho[5]+rho[8]);
SIGMAb[3] += mix_n*exp(rho[6]+rho[7]);
SIGMAb[4] += mix_n*exp(rho[6]+rho[8]);
// aggregate mixed strategy, from the narrow bracketers' perspective
vector[2] SIGMA1n;
SIGMA1n[1] = SIGMAb[1]+SIGMAb[2];
SIGMA1n[2] = SIGMAb[3]+SIGMAb[3];
vector[2] SIGMA2n;
SIGMA2n[1] = SIGMAb[1]+SIGMAb[3];
SIGMA2n[2] = SIGMAb[2]+SIGMAb[4];
vector[4] rhob = rho[1:4];
vector[2] rho1n = rho[5:6];
vector[2] rho2n = rho[7:8];
vector[8] H;
H[1:4] = Ab*rhob-lambda*Ab*Ub*SIGMAb-Bb*exp(rhob)+cb;
H[5:6] = An*rho1n-lambda*An*U1n*SIGMA1n-Bn*exp(rho1n)+cn;
H[7:8] = An*rho2n-lambda*An*U2n*SIGMA2n-Bn*exp(rho2n)+cn;
matrix[4,8] DSIGMAb = [
[(1-mix_n)*exp(rhob[1]),0, 0, 0, mix_n*exp(rho1n[1]+rho2n[1]),0, mix_n*exp(rho1n[1]+rho2n[1]),0 ],
[0, (1-mix_n)*exp(rho[2]),0, 0, mix_n*exp(rho1n[1]+rho2n[2]),0, 0, mix_n*exp(rho1n[1]+rho2n[2])],
[0, 0, (1-mix_n)*exp(rho[3]),0, 0 , mix_n*exp(rho1n[2]+rho2n[1]),mix_n*exp(rho1n[2]+rho2n[1]),0 ],
[0, 0, 0, (1-mix_n)*exp(rho[4]),0, mix_n*exp(rho1n[2]+rho2n[2]),0, mix_n*exp(rho1n[2]+rho2n[2])]
];
matrix[8,8] J = rep_matrix(0.0,8,8);
J[1:4,1:4] += Ab-Bb*diag_matrix(exp(rho[1:4]));
J[5:6,5:6] += An-Bn*diag_matrix(exp(rho[5:6]));
J[7:8,7:8] += An-Bn*diag_matrix(exp(rho[7:8]));
J[1:4,] += -lambda*Ab*Ub*DSIGMAb;
J[5:6,] += -lambda*An*U1n*Dsigma1n*DSIGMAb;
J[7:8,] += -lambda*An*U2n*Dsigma2n*DSIGMAb;
vector[8] drho = -J\H;
rho += drho;
if (max(abs(H))<tol) {
break;
}
if (tt==iter && max(abs(H))>tol) {
print("Maximum iterations reached without achieveing tolerance");
}
}
return rho;
}
}
data {
// payoff matrices for the 3 treatments
matrix[4,4] Ubroad[3];
matrix[2,2] Unarrow1[3];
matrix[2,2] Unarrow2[3];
int N; // number of observations
int id[N];
int nparticipants;
int treatment[N];
// coded as 1=AC, 2=AD, 3=BC, 4=BD
int action_broad[N];
vector[2] prior_r;
vector[2] prior_lambda;
vector[2] prior_mix;
real<lower=0> tol;
int<lower=1> iter;
int<lower=0,upper=1> UseData;
int<lower=0,upper=1> contextual;
}
transformed data {
vector[4] rho0 = rep_vector(log(1.0/4.0),4);
int action_narrow1[N];
int action_narrow2[N];
for (ii in 1:N) {
action_narrow1[ii] = action_broad[ii]==1 || action_broad[ii]==2 ? 1 : 2;
action_narrow2[ii] = action_broad[ii]==1 || action_broad[ii]==3 ? 1 : 2;
}
}
parameters {
// risk aversion
real<lower=0> r;
// choice precision
real<lower=0> lambda;
real<lower=0,upper=1> mix_narrow;
}
transformed parameters {
vector[8] rho[3];
vector[4] sigma_broad[3];
vector[2] sigma_narrow1[3];
vector[2] sigma_narrow2[3];
for (tt in 1:3) {
matrix[4,4] Ub = pow(Ubroad[tt],r);
matrix[2,2] U1n = pow(Unarrow1[tt],r);
matrix[2,2] U2n = pow(Unarrow1[tt],r);
if (contextual==1) {
// contextual utility normalization
Ub = Ub/(Ub[4,1]-Ub[4,4]);
U1n = U1n/(U1n[2,1]-U1n[2,2]);
U2n = U2n/(U2n[2,1]-U2n[2,2]);
}
rho[tt] = QRElp(
Ub,U1n, U2n,
lambda,mix_narrow,
iter,tol
);
sigma_broad[tt] = exp(rho[tt][1:4]);
sigma_narrow1[tt] = exp(rho[tt][5:6]);
sigma_narrow2[tt] = exp(rho[tt][7:8]);
}
}
model {
// likelihood contribution
if (UseData==1) {
for (ii in 1:N) {
target += log(
mix_narrow*sigma_narrow1[treatment[ii]][action_narrow1[ii]]*sigma_narrow2[treatment[ii]][action_narrow2[ii]]
+(1.0-mix_narrow)*sigma_broad[treatment[ii]][action_broad[ii]]
);
}
}
// priors
target += lognormal_lpdf(r | prior_r[1],prior_r[2]);
target += lognormal_lpdf(lambda | prior_lambda[1],prior_lambda[2]);
target += beta_lpdf(mix_narrow | prior_mix[1],prior_mix[2]);
}
generated quantities {
vector[4] sigma[3];
for (tt in 1:3) {
sigma[tt] = exp(rho[tt]);
}
}Note that they are exactly the same except for how the likelihood is clacluated.
15.4 Results
I estimate these three models with the following priors:
\[ \begin{aligned} \log r &\sim N(\log(0.5),1)\\ \log\lambda &\sim N(\log(5),1)\\ \phi&\sim \mathrm{Beta}(1,1) \end{aligned} \]
I used data from the final 10 rounds of play in order to hopefully mitigate any learning effects.
Table 15.3 shows the parameter estimates (first three rows) and implied action frequencies (remaining rows) for the four models estimated. The rightmost column of this Table also shows the empirical choice frequencies from the data used to estimate these models. Comparing the “Data” column to the estimated choice frequencies reveals that none the models perform particularly well at organizing these choice frequencies. However model posterior probabilities (assigning equal prior probabilities to all four models) overwhelmingly (one to at least five decimal places) select the dichotomous mixture model. Here we estimate that approximately 46% of participants narrowly bracket the games. If I were forced to select one of the one-type models, the narrow assumption does overwhelmingly better than the broad assumption based on model posterior probability.
fmt<-"%.3f"
ESTIMATES <-rbind(
summary("Code/HowManyGames/Fit_QRE_broad.rds" |> readRDS())$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
model = "Broad"
),
summary("Code/HowManyGames/Fit_QRE_narrow.rds" |> readRDS())$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
model = "Narrow"
),
summary("Code/HowManyGames/Fit_QRE_mix.rds" |> readRDS())$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
model = "Dichotomous"
),
summary("Code/HowManyGames/Fit_QRE_toolbox.rds" |> readRDS())$summary |>
data.frame() |>
rownames_to_column(var = "parameter") |>
mutate(
model = "Toolbox"
)
) |>
mutate(
treatment = parameter |> str_split_i(",",1) |> parse_number(),
action = c("{AC}","{AD}","{BC}","{BD}")[(parameter |> str_split_i(",",2) |> parse_number())]
) |>
filter(parameter!="lp__" & !grepl("rho",parameter) & !grepl("sigma_broad",parameter) & !grepl("sigma_narrow",parameter) & !grepl("sigma1",parameter) & !grepl("sigma2",parameter)) |>
mutate(
parameter = ifelse(is.na(action),parameter,
paste0("$\\sigma^",treatment,"_",action,"$")
)
) |>
mutate(
msd = paste0(sprintf(fmt,mean)," (",sprintf(fmt,sd),")")
)
TAB<-ESTIMATES |>
pivot_wider(
id_cols = parameter,
names_from = model,
values_from = msd
)
D<-"Data/Bland2019HowManyGames.csv" |>
read.csv() |>
arrange(uid,Period) |>
filter(Period>10) |>
mutate(
uid = paste("uid =",uid) |> as.factor() |> as.integer()
) |>
mutate(
BroadActComb = 2*(1-ActionAB)+(1-ActionCD)+1,
treatment = ifelse(
d1==0 & d2==0,1,ifelse(
d1!=0, 2, 3
)
)
) |>
group_by(treatment) |>
mutate(
count = n()
) |>
group_by(treatment,BroadActComb) |>
summarize(
Data = sprintf(fmt,n()/mean(count))
) |>
mutate(
parameter = paste0("$\\sigma^",treatment,"_",c("{AC}","{AD}","{BC}","{BD}")[BroadActComb],"$")
) |>
ungroup() |>
dplyr::select(Data,parameter,-treatment)
TAB<-TAB |>
full_join(D,by="parameter")
TAB[is.na(TAB)]<-""
TAB<-TAB |>
mutate(
order = ifelse(
parameter == "r",1,
ifelse(parameter=="lambda",2,
ifelse(parameter=="mix_narrow",3,4))
)
) |>
arrange(order,parameter) |>
dplyr::select(-order)
TAB |>
kbl(digits=3,caption = "Estimates from the three models. Posterior means with standard deviations in parentheses. The 'Data' column shows the empirical choice frequencies from the data used to estimate the models.") |>
kable_classic(full_width=FALSE) |>
add_header_above(c("","One type"=2,"Two types"=2,""))|
One type
|
Two types
|
||||
|---|---|---|---|---|---|
| parameter | Broad | Narrow | Dichotomous | Toolbox | Data |
| r | 0.474 (1.476) | 0.085 (0.026) | 0.022 (0.010) | 0.021 (0.010) | |
| lambda | 0.149 (0.075) | 2.281 (0.410) | 3.562 (0.330) | 0.935 (0.151) | |
| mix_narrow | 0.457 (0.059) | 0.927 (0.063) | |||
| \(\sigma^1_{AC}\) | 0.250 (0.002) | 0.397 (0.014) | 0.268 (0.003) | 0.258 (0.001) | 0.462 |
| \(\sigma^1_{AD}\) | 0.250 (0.001) | 0.249 (0.007) | 0.267 (0.003) | 0.257 (0.001) | 0.208 |
| \(\sigma^1_{BC}\) | 0.252 (0.002) | 0.217 (0.006) | 0.279 (0.004) | 0.260 (0.001) | 0.142 |
| \(\sigma^1_{BD}\) | 0.248 (0.004) | 0.137 (0.008) | 0.187 (0.004) | 0.225 (0.003) | 0.188 |
| \(\sigma^2_{AC}\) | 0.247 (0.001) | 0.245 (0.006) | 0.232 (0.003) | 0.242 (0.001) | 0.422 |
| \(\sigma^2_{AD}\) | 0.248 (0.001) | 0.402 (0.016) | 0.190 (0.005) | 0.232 (0.003) | 0.320 |
| \(\sigma^2_{BC}\) | 0.252 (0.001) | 0.134 (0.009) | 0.315 (0.007) | 0.268 (0.003) | 0.050 |
| \(\sigma^2_{BD}\) | 0.253 (0.002) | 0.220 (0.004) | 0.263 (0.005) | 0.259 (0.001) | 0.207 |
| \(\sigma^3_{AC}\) | 0.247 (0.001) | 0.309 (0.007) | 0.163 (0.009) | 0.225 (0.005) | 0.406 |
| \(\sigma^3_{AD}\) | 0.248 (0.001) | 0.194 (0.006) | 0.207 (0.003) | 0.236 (0.002) | 0.294 |
| \(\sigma^3_{BC}\) | 0.252 (0.001) | 0.305 (0.006) | 0.286 (0.002) | 0.263 (0.002) | 0.147 |
| \(\sigma^3_{BD}\) | 0.253 (0.002) | 0.192 (0.007) | 0.344 (0.012) | 0.276 (0.005) | 0.153 |
Of some concern here is that the estimated posterior means for \(r\) for all but the broad-only model are not particularly plausible: they are very close to zerom, which implies extreme risk-aversion. One explanation for this could be that risk-aversion isn’t really what is driving behavior in this game, and we might need another model of behavior to better organize the data. Alternatively, it could be that quantal response equilibrium is not a good model for play. In another chapter I will further explore this by taking a learning model to the data (instead of an equilibrium model). This may better respect the heterogeneity in participants’ decision-making processes.
One final thing to note here is that we are evaluating some non-nested models. In particular, broad bracketing cannot generate the same predictions as narrow bracketing and vice versa. On the other hand, both mixture models can generate the same predictions as the broad bracketing model by setting \(\phi=0\), or the narrow bracketing model by setting \(\phi=1\). On a a third hand, if \(\phi\) is strictly in the interior of the unit interval, then the two mixture models will have different implications about the data. The process of evaluating model posterior probabilities, however, is exactly the same whether they are nested or not. This is not the case if we were estimating these models using maximum likelihood, where we might have ended up having to use either AIC or BIC to select a model.
15.5 R code used to estimate these models
library(tidyverse)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(bridgesampling)
D<-"Data/Bland2019HowManyGames.csv" |>
read.csv() |>
arrange(uid,Period) |>
filter(Period>10) |>
mutate(
uid = paste("uid =",uid) |> as.factor() |> as.integer()
) |>
mutate(
BroadActComb = 2*(1-ActionAB)+(1-ActionCD)+1,
treatment = ifelse(
d1==0 & d2==0,1,ifelse(
d1!=0, 2, 3
)
)
)
U0broad<-rbind(
c(66, 66, 66, 66),
c(170,10,170, 10),
c(156,156,56,56),
c(260,100,160,0)
)
U0ABnarrow<-rbind(
c(10,10),
c(100,0)
)
U0CDnarrow<-rbind(
c(56,56),
c(160,0)
)
dStan<-list(
Ubroad = list(
U0broad,
U0broad + 50,
U0broad + 50
),
Unarrow1 = list(
U0ABnarrow,
U0ABnarrow+50,
U0ABnarrow
),
Unarrow2 = list(
U0CDnarrow,
U0CDnarrow,
U0CDnarrow+50
),
N = dim(D)[1],
id = D$uid,
nparticipants = D$uid |> unique() |> length(),
treatment = D$treatment,
action = D$BroadActComb,
action_broad = D$BroadActComb,
prior_r = c(log(0.5),1),
prior_lambda = c(log(5),1),
prior_mix = c(1,1),
tol = 10^-5,
iter = 1000,
UseData = 1,
contextual = 1
)
model_toolbox<-"Code/HowManyGames/QRE_toolbox.stan" |>
stan_model()
Fit_toolbox <- model_toolbox |>
sampling(
#chains=1,iter=100, # debugging settings
data=dStan,seed=42
)
bs_toolbox<-Fit_toolbox |> bridge_sampler()
Fit_toolbox |>
saveRDS("Code/HowManyGames/Fit_QRE_toolbox.rds")
model_mix<-"Code/HowManyGames/QRE_mix.stan" |>
stan_model()
Fit_mix <- model_mix |>
sampling(
#chains=1,iter=100, # debugging settings
data=dStan,seed=42
)
bs_mix<-Fit_mix |> bridge_sampler()
Fit_mix |>
saveRDS("Code/HowManyGames/Fit_QRE_mix.rds")
model_narrow<-"Code/HowManyGames/QRE_narrow.stan" |>
stan_model()
Fit_narrow<-model_narrow |>
sampling(data=dStan,seed=42)
bs_narrow<-Fit_narrow |> bridge_sampler()
Fit_narrow |>
saveRDS("Code/HowManyGames/Fit_QRE_narrow.rds")
model_broad<-"Code/HowManyGames/QRE_broad.stan" |>
stan_model()
Fit_broad <- model_broad |>
sampling(
#iter = 100,chains=1, # debugging settings
data=dStan,seed=42,
control = list(adapt_delta = 0.99)
)
bs_broad<-Fit_broad |> bridge_sampler()
Fit_broad |>
saveRDS("Code/HowManyGames/Fit_QRE_broad.rds")
postprob<-post_prob(bs_mix,bs_toolbox,bs_narrow,bs_broad)
postprob |>
saveRDS("Code/HowManyGames/postprob_QRE.rds")References
I thank an anonymous referee for requesting some QRE analysis in this paper. At the time, QRE was a passing interest of mine, but I hadn’t done any research using it. This was a very helpful nudge!↩︎