Supreme Court justices move left over time (more than they move right)

Mark Fisher, Atlanta Fed
and
Nicholas Georgakopoulos, Indiana University

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Genesis of the project

Nicholas asked a question at the Mathematica Stack Exchange Q&A community in March 2022 about “Markov Chain Monte Carlo Estimation”

I answered, saying (1) he was confused and (2) it was not possible to straighten him out (given the constraints of the forum)

He thanked me and provided his email just in case I changed my mind

I used his email to look him up

Based on what I learned I decided he probably wasn't a hopeless case, so I contacted him

And here we are!

Intro/Abstract

Martin and Quinn (2002) apply Item Response Theory (IRT) to Supreme Court votes in order to infer the “location” of justices on a simple left/right spectrum.

We modify their model, adopting a non-conjugate prior for the case parameters and allowing justices' locations to vary with each change in the court's composition (instead of with each term as MQ do).

We examine the joint distribution of justices' locations in their first and last compositions, showing that the tendency to move left over the course of their tenures is greater than the tendency to move right.

Although this finding is not new, we provide interesting illustrations via conditional distributions.

The model is simple

The only inputs are the justices' votes: (Subscript[v, jk])=1 indicates that justice j voted with the majority on case k and (Subscript[v, jk])=0 indicates a vote with the minority.

Let (Subscript[v, jk])∼Bernoulli((Subscript[π, jkt])), where (Subscript[π, jkt])=Φ((((Subscript[θ, jt])-(Subscript[γ, k]))/(Subscript[δ, k]))), Φ(·) is the CDF of the standard normal distribution, (Subscript[θ, jt]) is the location of justice j at "time" t, and ((Subscript[γ, k]),(Subscript[δ, k])) are the case parameters.

The posterior distribution is computed using a Gibbs sampler. We use a Metropolis step to draw the case parameters conditional on the locations. We discretize the location space and adopt a hidden Markov model wherein the locations become the hidden states. We use the forward-backward algorithm to draw the locations conditional on the case parameters.

The data

The data are contained in

"SCDB_2021_01_justiceCentered_Citation.csv"

which can be found at

"supremecourtdatabase.org"

The database covers the 75 terms from 1946 through 2020 inclusive, during which there are 35 compositions (the composition changes when a justice retires and/or a new justice joins the court).

There are 39 justices in the dataset. We allow each justice’s location to change with each composition in which the justice is a member. There are 307 justice-locations.

The votes (i.e., the data)

Given justice j and case k, the data (votes) are coded as follows:

(Subscript[v, jk])=({{{, ({{1, vote with majority}, {0, vote with minority}})}})

We drop cases with either unanimous decisions (the votes of which contain no information about the relative left/right locations) or tied decisions (for which there is no majority). This leaves

40,780 votes on 4,653 cases (about 8.8 votes/case)

The model: Item Response Theory (IRT)

(Subscript[v, jk])∼Bernoulli((Subscript[π, jkt]))

where

(Subscript[π, jkt])=Φ((((Subscript[θ, jt])-(Subscript[γ, k]))/(Subscript[δ, k])))

Φ(·) is the CDF of the standard normal distribution

(Subscript[θ, jt]) is the location of justice j at "time" t (i.e., composition t)

((Subscript[γ, k]),(Subscript[δ, k])) are the case parameters (the cut point and (the inverse of) the discrimination parameter)

Notes: 
(1) We do not use information as to whether a case is liberal or conservative.  Therefore there is an arbitrary degree of freedom as to what “left” and “right” mean. This can be removed by requiring (for example) Justice Thomas to be on the “right” or (equivalently) Justice Sotomayor to be on the “left”.   

(2) There is no “absolute zero.” There are only relative locations for the justices.

To complete the setup we need priors for the 307 justice locations and the 9,306 case parameters.

Bayes Rule

all locations: θ

all case parameters: η

all “parameters” (i.e., unknowns): ϕ=(θ,η)

data: y

joint distribution: p(y,ϕ)=p(y|ϕ)p(ϕ)=p(ϕ|y)p(y)

Bayes rule: p(ϕ|y)=((p(y|ϕ)p(ϕ))/(p(y)))

“Likelihood” p(y|ϕ), “prior” p(ϕ), and “posterior” p(ϕ|y)

Once we observe the data they are fixed and p(y) is a number

So Bayes rule is often expressed as a proportionality: p(ϕ|y)∝p(y|ϕ)p(ϕ)

Gibbs sampler

Joint posterior distribution: p(ϕ|y)=p(θ,η|y)

The complete set of full conditional distributions,  p(θ|y,η) and p(η|y,θ), completely characterizes the joint distribution

Gibbs sampler: Markov chain constructed by alternately sampling from each of the conditional distributions (MCMC)

The case parameters are conditionally independent: p(η|y,θ)=(∏_k)p((Subscript[η, k])|y,θ) where (Subscript[η, k]) is the case parameters for case k

The locations are conditionally independent: p(θ|y,η)=(∏_j)p((Subscript[θ, j])|y,η) where (Subscript[θ, j]) is the locations for justice j

Programs dedicated to Bayesian inference

bugs (Bayesian Updating via Gibbs Sampling)

jags (Just Another Gibbs Sampler)

stan (Stanislaw Ulam, invented Monte Carlo approach along with von Neumann and Metropolis)

Moore’s law

In the 20 years since Martin and Quinn (2002) was published computing power (per dollar in real terms) has increased by (10^3) or (10^4)

We can do far more in a given amount of time than they could

They made choices about the prior for the unknowns that delivered expressions that could be evaluated quickly

We don’t have to, and we don’t

Case parameters

Ours: (Subscript[π, jkt])=Φ(((Subscript[θ, jt])-(Subscript[γ, k]))/(Subscript[δ, k]))

Versus theirs: (Subscript[π, jkt])=Φ((Subscript[α, k])+(Subscript[β, k])(Subscript[θ, jt]))

For each case, we have (at most) 9 votes

So the prior for the case parameters plays an important role

MQ’s prior for β (the discrimination parameter) is normal and centered on zero

But if β=0 then then the case provides no information about locations

Our prior for δ=1/β is normal and centered on zero

(we actually do something a bit more involved)

As a consequence: We don’t get to use the expressions that can be evaluated quickly

Plot[{PDF[NormalDistribution[],x],PDF[TransformedDistribution[(1/\[FormalX]),\[FormalX]\[Distributed]NormalDistribution[0,1]],x]}//Evaluate,{x,-5,5},GridLinesStyle->Black,GridLines->{{0},{0}},PlotLegends->{"x","1/x"},PlotStyle->Thick,ImageSize->Large,BaseStyle->20]

Case parameters

Martin and Quinn (MQ) adopt a parameterization and prior for the case parameters that produces conditionally conjugate posterior distributions that allows for draws from standard distributions with their sampler. This was important 20 years ago when computation was much slower. Now, however, it is feasible to explore alternative parameterizations and priors.

Comparison of parameterizations:
	α + β θ=(θ-γ)/δ

MQ put independent standard normal priors on α and β. Their parameterization along with this prior allows them to make standard linear regression updates. Note, γ=-α/β, and so according to their prior γ has a Cauchy distribution.

We model γ and δ in a highly dependent way and involves hyper-parameters. The the posterior distribution we obtain has much in common with our prior.

ListPlot[Transpose[{alphaMean,betaMean}],ScalingFunctions->{"SignedLog","SignedLog"},GridLines->{{0},{0}},GridLinesStyle->Black,PlotRange->All,PlotStyle->Directive[Opacity[.1],PointSize[Medium]],FrameLabel->{"α","β"},ImageSize->Large,BaseStyle->20,PlotLabel->Column[{"Joint postererior distribution for ((Subscript[α, k]),(Subscript[β, k]))","implied by our prior on ((Subscript[γ, k]),(Subscript[δ, k]))"},Center]]

Locations: Setup

We discretize the location space into 701 possible states:

$thetaRange=Range[-4,3,.01];

The prior for a justice’s location/state in their first composition is a  discretized normal distribution:

$bigG0= 
PDF[
NormalDistribution[0,.2],
$thetaRange
]//Chop//Normalize[#,Total]&

The transition matrix (701×701) is computed from the vector of states and a standard deviation (that is itself unknown and must be inferred):

Normalize[Exp[(((($thetaRange-#))^2)/(-2(σ^2)))],Total]&/@$thetaRange

Each justice has a likelihood matrix (Subscript[L, j]) with a row for each of the justice’s compositions and 701 columns. The likelihood values in the rows depend on the cases in each composition which vary in number across compositions.

(The distribution that determines the relation between the "emissions" and "hidden states" is not time-homogeneous and not equivalent to the third argument in HiddenMarkovProcess.)

Forward-backward algorithm

We use the forward-backward algorithm to sample from the posterior distribution for the locations (conditional on the case parameters).

The forward-backward algorithm can be used to randomly sample from a HiddenMarkovProcess object.

Note the second With is there to work around a bug in the compiled version of FoldList that has been present since Version 3.0 released in 1996.

forwardBackward=Compile[
{{L,_Real,2},{Q,_Real,2},{G0,_Real,1},{thetaRange,_Real,1}},
With[{
Grev=FoldList[#/Total[#]&[#1.Q*#2]&,G0,L]//Rest//Reverse,
range=Range[Length[thetaRange]],
Qtr=Transpose[Q]
},
With[{start=RandomChoice[First[Grev]->range]},
thetaRange[[
FoldList[
RandomChoice[Qtr[[#1]]*#2->range]&,
start,
Rest[Grev]
]//Reverse
]]
]
]
]

Locations: Results

Does it appear that justices tend to move left? (It looks that way to me, but how can we characterize this idea more formally?)

ListLinePlot[courtsThetaMean,PlotLegends->Placed[LineLegend[Values[justiceNames]],Below],PlotRange->{-3.3,3.3},GridLines->{{},{0}},GridLinesStyle->Black,ImageSize->Large,BaseStyle->20,FrameLabel->{"composition","location"},PlotStyle->Directive[Thickness[.008],Opacity[.8]]]

ListLinePlot[meanLocationsByTerm,PlotLegends->Placed[LineLegend[Values[justiceNames]],Below],PlotRange->{-3.3,3.3},GridLines->{{},{0}},GridLinesStyle->Black,ImageSize->Large,BaseStyle->20,FrameLabel->{"term","location"},PlotStyle->Directive[Thickness[.008],Opacity[.8]]]

A sense of the uncertainty: sampling 50 paths at random:

Show[ListLinePlot[byJusticeToByTerm[#],PlotStyle->Opacity[.1]]&/@RandomChoice[thetaDraws,50],ImageSize->Large,BaseStyle->20,FrameLabel->{"term","location"}]

First and last compositions

Only include the 35 justices present in two or more compositions.

We examine the joint distribution of the location of justices in their first and last compositions, which we refer to as (first, last) location pairs.

Consider the following selection process: Choose a justice at random (using equal weights) and then choose a (first, last) pair at random from the posterior distribution of location paths for that justice. This process is equivalent to combining all justices' (first, last) pairs and then drawing a pair at random.

We examine the joint distribution in two ways, both of which rely on DataDistribution as a container. Recall the locations are discretized for convenience.

We can use EmpiricalDistribution to construct a discrete joint distribution of (first, last) pairs, which we can then use to compute various conditional probabilities.

Alternatively, we can use SmoothKernelDistribution to compute a continuous joint distribution, which provides us with a better basis for visualizations.

Discrete distribution

(In the table “you” means “a justice”.)

Dataset[assoc,ItemStyle->Large,HeaderStyle->Large]

Continuous distribution

Plot3D[PDF[skd,{first,last}],{first,-4,4},{last,-4,4},PlotRange->All,PlotPoints->50,PlotLabel->"Joint Distribution",AxesLabel->{"first court","last court"},SphericalRegion->True,ImagePadding->30,ImageSize->Large,BaseStyle->20]

cp0=ContourPlot[PDF[skd,{first,last}],{first,-4,4},{last,-4,4},PlotRange->{0,All},PlotPoints->50,PlotLabel->"Joint Distribution",FrameLabel->{"first court","last court"},Contours->12,ContourShading->None,GridLines->{{0},{0}},GridLinesStyle->Black,Prolog->{Gray,Dashed,Line[{{-5,-5},{5,5}}]},ImageSize->Large,BaseStyle->20]

Show[cp0,
ListPlot[DeleteCases[AssociationThread[Values[justiceNames]->(#[[{1,-1}]]&/@thetaMean)],{x_,x_}]],PlotRange->{{-2,2.5},{-3,2}},PlotLabel->"Joint Distribution with justices' means shown"]

Through[{Mean,StandardDeviation}[#]]&/@margDists//PaddedForm[#,{3,2}]&

Plot[PDF[#,x]&/@margDists//Evaluate,{x,-4,4},Filling->0,GridLines->{{0},{0}},GridLinesStyle->Black,PlotLegends->Placed[{"first court","last court"},{.85,.7}],FrameLabel->{"location","probability density"},PlotLabel->"Marginal Distributions",ImageSize->Large,BaseStyle->20]

Framed@Labeled[Manipulate[
Plot[{condPDF1[-first][last],condPDF1[first][last]},{last,-3.5,3.5},Filling->0,PlotRange->{All,2.1},PlotLabel->"Distribution of the last court given the first court",FrameLabel->{"last court","probability density"},
GridLines->{{{-first,Directive[Thickness[Large],ColorData[97][1]]},{first,Directive[Thickness[Large],ColorData[97][2]]},{0,Dashed}},{}},GridLinesStyle->Black,Epilog->{Text[Style[PaddedForm[N[first],{3,2}],"Text"],{first,1.7},Background->White],Text[Style[PaddedForm[N[-first],{3,2}],"Text"],{-first,1.7},Background->White]},ImageSize->Large,BaseStyle->20],
{{first,.6},0,2.5},SaveDefinitions->True,LabelStyle->Large],Style[Column[{"Two symmetric conditioning values shown","to facilitate comparison of the corresponding distributions"},Center],"Text",20],Top]

Framed@Labeled[Manipulate[
Plot[{condPDF2[-last][first],condPDF2[last][first]},{first,-3.5,3.5},Filling->0,PlotRange->{All,2.1},PlotLabel->"Distribution of the first court given the last court",FrameLabel->{"first court","probability density"},
GridLines->{{{-last,Directive[Thickness[Large],ColorData[97][1]]},{last,Directive[Thickness[Large],ColorData[97][2]]},{0,Dashed}},{}},GridLinesStyle->Black,Epilog->{Text[Style[PaddedForm[N[last],{3,2}],"Text"],{last,1.7},Background->White],Text[Style[PaddedForm[N[-last],{3,2}],"Text"],{-last,1.7},Background->White]},ImageSize->Large,BaseStyle->20],
{{last,2},0,3},SaveDefinitions->True,LabelStyle->Large],Style[Column[{"Two symmetric conditioning values shown","to facilitate comparison of the corresponding distributions"},Center],"Text",20],Top]

The end

Necessary? (initialization) stuff

Note that tab is about 77MB. It is computed from the locations.

tab=;
Clear[first,last];
ed=EmpiricalDistribution[tab];
probFun[args_]:=PaddedForm[Probability[args,{first,last}\[Distributed]ed],{3,2}]

skd=SmoothKernelDistribution[tab,"Oversmooth",InterpolationPoints->200];
margDists=Table[MarginalDistribution[skd,j],{j,2}];condPDF1[first_][last_]=((PDF[skd,{first,last}])/(PDF[margDists[[1]],first]));condPDF2[last_][first_]=((PDF[skd,{first,last}])/(PDF[margDists[[2]],last]));

argList=;
tagLines=;
assoc=AssociationThread[Column/@Transpose[{tagLines,("\t"<>ToString[#])&/@argList}]->probFun/@argList];