rm(list=ls(all=TRUE))
#you will need the following libraries:
library(arm)
library(boot)
library(plotrix) #for calculating standard error of mean



#1) The analysis of the distribution of partisanship rates requires converting data on which individual judges served in a given year to circuit-year information. This is done in the following section

####################################################################
###################################################################
#CIRCUIT JUDGES

#read in dataset with circuit-year as unit of analysis, and indicator for
    #active and senior judges on circuit that is 1 if they are a dem and 0 if gop.

judges <- read.dta("circuit_judges_by_year_party_only.dta", convert.underscore=T)

#for now, keep only 1925-2004, to match up with songer and sunstein data.

judges <- judges[judges$year>=1925 & judges$year<=2004,] #drop 1891 since only 1 judge in many cases
detach()
attach.all(judges, overwrite=T)
circuit.year <- paste(circuit, year)
circuit <- ifelse(circuit == 0, 12, circuit)#recode DC circuit
unique.circuit.year <- unique(circuit.year)
unique.circuit <- unique(circuit)
unique.year <- unique(year)
loop1 <- length(circuit.year)

#29 active judge slots (with many missing in many circ-years)
#24 senior judge slots (with many missing in many circ-years)
    #note that DC circuit is coded as "0"

number.seats <- judgeships.authorized#number of authorized seats by circuit

#transpose data to ease calculations
#active.tranpose <- data.frame(t(judges[1:32]))

#first, get total number of active and senior judges by circuit year
    #j1dem:j29dem appear in columns 4:32 
    #also get number of dems
xx <- 4
yy <- 32
#pull out for easier calculation
active.judges.matrix <- judges[xx:yy]
active.judges.circyear <-  rowSums(ifelse(is.na(active.judges.matrix), 0, 1)) 
active.dems.circyear <-  rowSums(active.judges.matrix, na.rm=T)
active.gop.circyear <- active.judges.circyear - active.dems.circyear
#senior1.dem:senior24.dem appear in columbs 33:56 DOUBLE CHECK THIS !!
rr <- 33
tt <- 56
senior.judges.matrix <- judges[rr:tt]
senior.judges.circyear <-  rowSums(ifelse(is.na(senior.judges.matrix), 0, 1)) 
senior.dems.circyear <-  rowSums(senior.judges.matrix, na.rm=T)
senior.gop.circyear <- senior.judges.circyear - senior.dems.circyear

total.judges.circyear <- active.judges.circyear + senior.judges.circyear
total.dems.circyear <- active.dems.circyear + senior.dems.circyear
total.gop.circyear <- active.gop.circyear + senior.gop.circyear

dem.prop.active.circyear <- active.dems.circyear/active.judges.circyear
dem.prop.senior.circyear <- senior.dems.circyear/senior.judges.circyear
dem.prop.total.circyear <- total.dems.circyear/total.judges.circyear

#build yearly meaures
    #including only active judges
loop <- length(unique.year)
    
number.seats.year <- rep(NA, loop)#number of authorized seats per year
active.judges.year <- rep(NA, loop)#number of active judges per year
active.dems.year <- rep(NA, loop)
active.gop.year <- rep(NA, loop)

    #senior judges
senior.judges.year <- rep(NA, loop)#number of senior judges per year
senior.dems.year <- rep(NA, loop)
senior.gop.year <- rep(NA, loop)
  
    #total combined
total.judges.year <- rep(NA, loop)
total.dems.year <- rep(NA, loop)
total.gop.year <- rep(NA, loop)


for (i in 1:loop){
  keep <- year == unique.year[i]
  number.seats.year[i] <- sum(judgeships.authorized[keep])
  active.judges.year[i] <- sum(active.judges.circyear[keep])
  active.dems.year[i] <-  sum(active.dems.circyear[keep])
  active.gop.year[i] <- sum(active.gop.circyear[keep])
  senior.judges.year[i] <- sum(senior.judges.circyear[keep])
  senior.dems.year[i] <- sum(senior.dems.circyear[keep])
  senior.gop.year[i] <- sum(senior.gop.circyear[keep])
  total.judges.year[i] <- sum(total.judges.circyear[keep])
  total.dems.year[i] <- sum(total.dems.circyear[keep])
  total.gop.year[i] <- sum(total.gop.circyear[keep])
  }

active.dems.per.year <- active.dems.year/(active.dems.year+active.gop.year)
senior.dems.per.year <- senior.dems.year/(senior.dems.year+senior.gop.year)    
total.dems.per.year <- total.dems.year/(total.dems.year+total.gop.year)

year.matrix <- data.frame(year = unique.year,  number.seats.year = number.seats.year ,  active.judges.year = active.judges.year,  
  active.dems.year = active.dems.year,  active.dems.per.year = active.dems.per.year,  
   senior.judges.year = senior.judges.year,  senior.dems.year = senior.dems.year,  senior.dems.per.year = senior.dems.per.year,   
   total.judges.year = total.judges.year,  total.dems.year = total.dems.year,  total.dems.per.year = total.dems.per.year)


#use these commands to export summary info to a dataset.
partisan.composition <- data.frame(year, circuit, number.seats, active.judges.circyear,  
  active.dems.circyear, active.gop.circyear, dem.prop.active.circyear, senior.judges.circyear, dem.prop.senior.circyear, dem.prop.total.circyear, 
  senior.judges.circyear, senior.dems.circyear, senior.gop.circyear, total.judges.circyear, total.dems.circyear,total.gop.circyear )
write.csv(partisan.composition, "partisan.composition.csv")
  

####################################################################
####################################################################
#GET ESTIMATED PROPORTION OF VARIOUS PANEL TYPES
####################################################################
####################################################################

#this analysis uses a condensed version of the Songer database that
    #contains all cases heard by a three-judge panel, including those in
    #which the outcome cannot be be coded as liberal or conservative.
#It differs from the dataset that analyzes votes, which contains only
    #cases where the outcome CAN be coded as liberal or conservative

#the unit in this dataset is the case, with multiple columns for each judge on the panel.
case.data <- read.dta("songer_caselevel_1925_2002_for_mixed_panels.dta", convert.underscore=T)
detach()
attach.all(case.data, overwrite=T)
case <- casenum
year <- year
circuit <- ifelse(circuit == 0, 12, circuit)#recode DC
j1.dem <- ifelse(j1presparty == 100,1,0)#code each judge as 1 if they were appointed by a DEM prz
j2.dem <- ifelse(j2presparty == 100,1,0)
j3.dem <- ifelse(j3presparty == 100,1,0)


keep.data <- data.frame(case, year, circuit, j1.dem, j2.dem, j3.dem)
detach()
attach.all(keep.data, overwrite = T)
keep.data <- keep.data[order(year,circuit, case),] #order by year, circuit, case number
detach()
attach.all(keep.data, overwrite=T)
keep.data$circuit.year <- paste(circuit,year)#create indicator for circuit-year
detach()
attach.all(keep.data, overwrite=T)
unique.circuit.year <- unique(circuit.year)
unique.year <- unique(year)
loop <- length(unique.circuit.year)#shortcut for loops below

#The Songer database is a stratified random sample, with a certain number of cases
#sampled in each-year. To make population-level inferences, it is necessary to weight the data 
#by the number of cases actually hear in a circuit-year
#the procedure for doing this is described in the Songer database codebook
#and is implemented here:

#attach weights to data  create table like this:

#                sample of circuit          universe of circuit      
#circuit     # cases     proportion   # cases   proportion   weight   
#_________________________________________________________________    
#01              15      .1              095     .049       0.49     
#02              15      .1              329     .170       1.70     
#03              15      .1              116     .060       0.60     
#04              15      .1              099     .051       0.51     
#05              15      .1              175     .091       0.91     
#06              15      .1              222     .115       1.15     
#07              15      .1              081     .042       0.42     
#08              15      .1              330     .171       1.71     
#09              15      .1              289     .150       1.50     
#DC              15      .1              196     .101       1.01     
#_______________________________________________________________      
#total          150     1.0              1932    1.0      

#use circuit-year as unit of analysis
songer.weights <- read.dta("songer_weight_data.dta")#
songer.weights$circuit <- ifelse(songer.weights$circuit == 0, 12, songer.weights$circuit)#recode DC circuit
detach()
attach.all(songer.weights)
#note: from 1931 to 1937, dc circuit not included  -- remove from weights data
drop.var <- ifelse(songer.weights$circuit == 12 & songer.weights$year > 1930 & songer.weights$year < 1938, 1,0)
songer.weights <- songer.weights[drop.var == 0,]
detach()
attach.all(songer.weights)
songer.weights <- songer.weights[order(year,circuit),]
detach()
attach.all(songer.weights)
n.cases.circuit.universe <- cases #number of cases in universe, by circuit-year
songer.weights$circuit.year <- paste(circuit,year)#for removing cases below
#since keep.data has fewer circuit-years than songer, use length of keep.data for looping
n.cases.circuit.sample <- rep(NA, loop)#number of cases in a circuit, in a given year
n.cases.year.sample <- rep(NA, loop)#number of cases total, in a given year
sample.proportion <- rep(NA, loop)#propotion of all in-sample cases, for a given year
universe.proportion <- rep(NA, loop)#proportion of all cases in universe, for a given year

#SAMPLE MEASURES:  get number of cases and proportions from actual Songer data, by circuit and year
detach()
attach.all(keep.data, overwrite=T)

#get number of cases, per circuit-year for sample 
for (i in 1:loop){ #loop over circuit year
     n.cases.circuit.sample[i] <- length(case[circuit.year == unique.circuit.year[i]])
            }
            
#first, create vector based on unique-year, then we will attach to circuit-year data.
n.cases.byyear.sample <- rep(NA, length(unique.year))
for(i in 1:length(unique.year)){
    n.cases.byyear.sample[i] <-  sum(length(case[year == unique.year[i]]))
    }
n.cases.byyear.data <- data.frame(n.cases.byyear.sample, unique.year)#for merging to circuit-year level
n.cases.byyear.data$year <- n.cases.byyear.data$unique.year
n.cases.byyear.data <- n.cases.byyear.data[,c(1,3)]
    
#expand to circuit-year    
detach()
attach.all(songer.weights)
songer.weights <- data.frame(cbind(songer.weights, n.cases.circuit.sample))
songer.weights <- merge(songer.weights, n.cases.byyear.data)
#get sample proportion
songer.weights$sample.proportion  <- songer.weights$n.cases.circuit.sample/songer.weights$n.cases.byyear.sample
detach()
attach.all(songer.weights, overwrite=T)

###############################
##UNIVERSE MEASURES

#get total number of cases in each year:
#first, create vector based on unique-year, then we will attach to circuit-year data.
n.cases.byyear.universe <- rep(NA, length(unique.year))
for (i in 1:length(n.cases.byyear.universe)){
    n.cases.byyear.universe[i] <- sum(cases[year==unique.year[i]])
    }
#expand to circuit-year    
universe.year.data <- data.frame(n.cases.byyear.universe, unique.year)#for merging
universe.year.data$year <- universe.year.data$unique.year
universe.year.data <- universe.year.data[,c(1,3)]
songer.weights <- merge(songer.weights, universe.year.data) 
 #get universe proportion
songer.weights$universe.proportion  <- songer.weights$cases/songer.weights$n.cases.byyear.universe
#GET WEIGHTS = DIVIDE UNIVERSE PROPORTION BY SAMPLE PROPORTION
songer.weights$weight <- songer.weights$universe.proportion / songer.weights$sample.proportion
detach()
attach.all(songer.weights, overwrite=T)

#merge weights with case-level data

keep.data <- merge(keep.data, songer.weights)
detach()
attach.all(keep.data, overwrite=T)
keep.data <- keep.data[order(year,circuit, case),]
detach()
attach.all(keep.data, overwrite=T)
songer.data.to.use <- data.frame(year, circuit, circuit.year, j1.dem, j2.dem, j3.dem, weight)

###################################################
#GET AGGREGATE STATISTICS
###################################################
#FIRST, DO SONGER DATA, YEARS 1925-2002
#1) ALL CIRCUITS COMBINED
songer.data.to.use$circuit <- ifelse(circuit == 0, 12,circuit)#recode DC circuit to ease matrix replacement
detach()
attach.all(songer.data.to.use)
unique.year <- unique(year)
detach()
attach.all(songer.data.to.use)

#code party configurations of panel
ddd <- ifelse(j1.dem == 1 & j2.dem == 1 & j3.dem == 1, 1, 0)
rrr   <- ifelse(j1.dem == 0 & j2.dem == 0 & j3.dem == 0, 1, 0)
rrd <-  ifelse(j1.dem + j2.dem + j3.dem == 1, 1,0)
ddr <- ifelse(j1.dem + j2.dem + j3.dem == 2, 1,0)
mixed <- ifelse(ddr == 1 | rrd ==1,1,0)
dem.maj <- ifelse(j1.dem + j2.dem + j3.dem >= 2,1,0)#dem majority panels (ddd|ddr)
gop.maj <- ifelse(j1.dem + j2.dem + j3.dem <= 1,1,0)#gop majority panels (rrr|rrd)
unified <- ifelse(ddd==1 | rrr==1,1,0) #unified panels

#for each year, get estimated number of each type of panels, along with 
    #bootstrapped confidence intervals

#get pooled means, combining all years
mean.rrr.all.songer <- weighted.mean(rrr, weight)
mean.ddd.all.songer <- weighted.mean(ddd, weight)
mean.rrd.all.songer <- weighted.mean(rrd, weight)
mean.ddr.all.songer <- weighted.mean(ddr, weight)
mean.mixed.all.songer <- weighted.mean(mixed, weight)


loop <- length(unique.year)#use shortcut for looping

#set up empty  vectors for yearly calculations
mixed.prop.songer <- rep(NA, loop)
mixed.prop.ci.songer  <- matrix(nrow = loop, ncol =2)
rrr.prop.songer <- rep(NA, loop)
rrr.prop.ci.songer  <- matrix(nrow = loop, ncol =2)
ddd.prop.songer <- rep(NA, loop)
ddd.prop.ci.songer  <- matrix(nrow = loop, ncol =2)
rrd.prop.songer <- rep(NA, loop)
rrd.prop.ci.songer  <- matrix(nrow = loop, ncol =2)
ddr.prop.songer <- rep(NA, loop)
ddr.prop.ci.songer  <- matrix(nrow = loop, ncol =2)
dem.maj.prop.songer <- rep(NA, loop) #number of maj dem panels
dem.maj.prop.ci.songer   <- matrix(nrow = loop, ncol =2)
gop.maj.prop.songer <- rep(NA, loop) #number of maj gop panels
gop.maj.prop.ci.songer   <- matrix(nrow = loop, ncol =2)
unified.prop.songer <- rep(NA, loop)
unified.prop.ci.songer  <- matrix(nrow = loop, ncol =2)

#for bootstrapping
wmean <- function(x,k) { #x is cbind(proportion,weight)
  weighted.mean(x[k,1],x[k,2])
}

for (i in 1:loop){#loop over years
  keep <-  year == unique.year[i]
    mixed.prop.songer[i] <- weighted.mean(mixed[keep], weight[keep])
    mixed.prop.boot <-  boot(cbind(mixed[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    mixed.prop.ci.songer[i,] <- quantile(mixed.prop.boot$t,probs=c(.025,.975))
    rrd.prop.songer[i] <- weighted.mean(rrd[keep], weight[keep])
    rrd.prop.boot <-  boot(cbind(rrd[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    rrd.prop.ci.songer[i,] <- quantile(rrd.prop.boot$t,probs=c(.025,.975))
    ddr.prop.songer[i] <- weighted.mean(ddr[keep], weight[keep])
    ddr.prop.boot <-  boot(cbind(ddr[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    ddr.prop.ci.songer[i,] <- quantile(ddr.prop.boot$t,probs=c(.025,.975))
    rrr.prop.songer[i] <- weighted.mean(rrr[keep], weight[keep])
    rrr.prop.boot <-  boot(cbind(rrr[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    rrr.prop.ci.songer[i,] <- quantile(rrr.prop.boot$t,probs=c(.025,.975))
    ddd.prop.songer[i] <- weighted.mean(ddd[keep], weight[keep])
    ddd.prop.boot <-  boot(cbind(ddd[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    ddd.prop.ci.songer[i,] <- quantile(ddd.prop.boot$t,probs=c(.025,.975))
    dem.maj.prop.songer[i] <- weighted.mean(dem.maj[keep], weight[keep])
    dem.maj.prop.boot <-  boot(cbind(dem.maj[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    dem.maj.prop.ci.songer[i,] <- quantile(dem.maj.prop.boot$t,probs=c(.025,.975))
    gop.maj.prop.songer[i] <- weighted.mean(gop.maj[keep], weight[keep])
    gop.maj.prop.boot <-  boot(cbind(gop.maj[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    gop.maj.prop.ci.songer[i,] <- quantile(gop.maj.prop.boot$t,probs=c(.025,.975))
    unified.prop.songer[i] <- weighted.mean(unified[keep], weight[keep])
    unified.prop.boot <-  boot(cbind(unified[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    unified.prop.ci.songer[i,] <- quantile(unified.prop.boot$t,probs=c(.025,.975))
    }


#2) ESTIMATES FOR EACH CIRCUIT

#create one matrix for each variable, with circuits as columns to store means
#create separate matrices for upper and lower bounds of confidence intervals

mixed.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
mixed.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
mixed.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrr.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrr.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrr.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddd.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddd.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddd.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrd.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrd.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
rrd.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddr.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddr.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
ddr.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
dem.maj.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
dem.maj.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
dem.maj.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
gop.maj.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
gop.maj.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
gop.maj.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
unified.prop.circuit.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
unified.circ.ci.lo.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))
unified.circ.ci.hi.songer <- matrix(NA, nrow = loop, ncol = length(unique(circuit)))

#for confidence, set lower and upper bounds at 0 and 1 using ifelse commands
#also, in some years, cells will be empty because either circuit does not exist (10th and 11th circuit)
    #or no cases entered the data that year (DC).  us ifelse commands to pull out these years as missing

#NOTE: This will take a while to run. To speed up, you can adjust the number of bootstraps:
n.bootstraps <- 1000

for (i in 1:loop){ # loop over year
     keep1 <-  year == unique.year[i]
    for (j in 1:ncol(mixed.prop.circuit.songer)){ #then loop over each circuit within year
      keep2 <- circuit == j
      
      empty <-  is.nan(mean(mixed[keep1 & keep2]))#denote circuit-years with no cases
  
      mixed.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(mixed[keep1 & keep2]))
      if (empty == FALSE)      mixed.prop.boot <-  boot(cbind(mixed[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      mixed.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(mixed.prop.boot$t,probs=c(.025)))
      mixed.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(mixed.prop.boot$t,probs=c(.975)))
 
      ddd.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(ddd[keep1 & keep2]))
      if (empty == FALSE)      ddd.prop.boot <-  boot(cbind(ddd[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      ddd.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(ddd.prop.boot$t,probs=c(.025)))
      ddd.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(ddd.prop.boot$t,probs=c(.975)))       
    
      rrr.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(rrr[keep1 & keep2]))
      if (empty == FALSE)      rrr.prop.boot <-  boot(cbind(rrr[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      rrr.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(rrr.prop.boot$t,probs=c(.025)))
      rrr.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(rrr.prop.boot$t,probs=c(.975)))     
              
      rrd.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(rrd[keep1 & keep2]))
      if (empty == FALSE)      rrd.prop.boot <-  boot(cbind(rrd[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      rrd.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(rrd.prop.boot$t,probs=c(.025)))
      rrd.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(rrd.prop.boot$t,probs=c(.975)))
 
      ddr.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(ddr[keep1 & keep2]))
      if (empty == FALSE)      ddr.prop.boot <-  boot(cbind(ddr[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      ddr.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(ddr.prop.boot$t,probs=c(.025)))
      ddr.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(ddr.prop.boot$t,probs=c(.975)))
 
      dem.maj.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(dem.maj[keep1 & keep2]))
      if (empty == FALSE)      dem.maj.prop.boot <-  boot(cbind(dem.maj[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      dem.maj.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(dem.maj.prop.boot$t,probs=c(.025)))
      dem.maj.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(dem.maj.prop.boot$t,probs=c(.975)))
    
      gop.maj.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(gop.maj[keep1 & keep2]))
      if (empty == FALSE)      gop.maj.prop.boot <-  boot(cbind(gop.maj[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      gop.maj.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(gop.maj.prop.boot$t,probs=c(.025)))
      gop.maj.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(gop.maj.prop.boot$t,probs=c(.975)))
 
      unified.prop.circuit.songer[i,j] <-  ifelse(empty, NaN, mean(unified[keep1 & keep2]))
      if (empty == FALSE)      unified.prop.boot <-  boot(cbind(unified[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      unified.circ.ci.lo.songer[i,j] <-  ifelse(empty, NaN, quantile(unified.prop.boot$t,probs=c(.025)))
      unified.circ.ci.hi.songer[i,j] <-  ifelse(empty, NaN, quantile(unified.prop.boot$t,probs=c(.975)))
       
        }
    }
        
#NOW DO SAME FOR SUNSTEIN DATA
#note that since data is not weighted, can use regular means

#Sunstein case data has > 200 cases for each year from 1995-2004
#grab all years to do validity check and compare to songer for 1995-2002
#use 2003-2004 and add to songer data

sunstein.data <- read.dta("sunstein_case_level.dta", convert.underscore=T)
detach()
attach.all(sunstein.data, overwrite=T)
keep <- year > 1994 & year <2005
year <- year[keep]
circuit <- circuit[keep]
circuit.year <- paste(circuit, year)
j1.dem <- pa1[keep]
j2.dem <- pa2[keep]
j3.dem <- pa3[keep]
weight <- rep(1, length(year[keep]))#sunstein data is not weighted
unique.year <- unique(year)
unique.circuit <- unique(circuit)
sunstein.data.to.use <- data.frame(year, circuit, circuit.year, j1.dem, j2.dem, j3.dem, weight)
detach()
attach.all(sunstein.data.to.use, overwrite=T)
detach()
attach.all(sunstein.data.to.use, overwrite=T)

#code party configurations of panel
ddd <- ifelse(j1.dem == 1 & j2.dem == 1 & j3.dem == 1, 1, 0)
rrr   <- ifelse(j1.dem == 0 & j2.dem == 0 & j3.dem == 0, 1, 0)
rrd <-  ifelse(j1.dem + j2.dem + j3.dem == 1, 1,0)
ddr <- ifelse(j1.dem + j2.dem + j3.dem == 2, 1,0)
mixed <- ifelse(ddr == 1 | rrd ==1,1,0)
dem.maj <- ifelse(j1.dem + j2.dem + j3.dem >= 2,1,0)#dem majority panels (ddd|ddr)
gop.maj <- ifelse(j1.dem + j2.dem + j3.dem <= 1,1,0)#gop majority panels (rrr|rrd)
unified <- ifelse(ddd==1 | rrr==1,1,0) #unified panels

#1) All Circuits Combined

#for each year, get estimated number of each type of panels, along with 
    #standard error of means

loop.sunstein <- length(unique.year)

#set up emptyy  vectors for yearly calculations
mixed.prop.sunstein <- rep(NA, loop.sunstein)
mixed.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)
rrr.prop.sunstein <- rep(NA, loop.sunstein)
rrr.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)
ddd.prop.sunstein <- rep(NA, loop.sunstein)
ddd.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)
rrd.prop.sunstein <- rep(NA, loop.sunstein)
rrd.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)
ddr.prop.sunstein <- rep(NA, loop.sunstein)
ddr.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)
dem.maj.prop.sunstein <- rep(NA, loop.sunstein) #number of maj dem panels
dem.maj.prop.ci.sunstein   <- matrix(nrow = loop.sunstein, ncol =2)
gop.maj.prop.sunstein <- rep(NA, loop.sunstein) #number of maj gop panels
gop.maj.prop.ci.sunstein   <- matrix(nrow = loop.sunstein, ncol =2)
unified.prop.sunstein <- rep(NA, loop.sunstein)#prop of unified panels (DD
unified.prop.ci.sunstein  <- matrix(nrow = loop.sunstein, ncol =2)

for (i in 1:loop.sunstein){#loop.sunstein over years or two-years
  keep <-  year == unique.year[i]
   
    mixed.prop.sunstein[i] <- weighted.mean(mixed[keep], weight[keep])
    mixed.prop.boot <-  boot(cbind(mixed[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    mixed.prop.ci.sunstein[i,] <- quantile(mixed.prop.boot$t,probs=c(.025,.975))
    rrd.prop.sunstein[i] <- weighted.mean(rrd[keep], weight[keep])
    rrd.prop.boot <-  boot(cbind(rrd[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    rrd.prop.ci.sunstein[i,] <- quantile(rrd.prop.boot$t,probs=c(.025,.975))
    ddr.prop.sunstein[i] <- weighted.mean(ddr[keep], weight[keep])
    ddr.prop.boot <-  boot(cbind(ddr[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    ddr.prop.ci.sunstein[i,] <- quantile(ddr.prop.boot$t,probs=c(.025,.975))
    rrr.prop.sunstein[i] <- weighted.mean(rrr[keep], weight[keep])
    rrr.prop.boot <-  boot(cbind(rrr[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    rrr.prop.ci.sunstein[i,] <- quantile(rrr.prop.boot$t,probs=c(.025,.975))
    ddd.prop.sunstein[i] <- weighted.mean(ddd[keep], weight[keep])
    ddd.prop.boot <-  boot(cbind(ddd[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    ddd.prop.ci.sunstein[i,] <- quantile(ddd.prop.boot$t,probs=c(.025,.975))
    dem.maj.prop.sunstein[i] <- weighted.mean(dem.maj[keep], weight[keep])
    dem.maj.prop.boot <-  boot(cbind(dem.maj[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    dem.maj.prop.ci.sunstein[i,] <- quantile(dem.maj.prop.boot$t,probs=c(.025,.975))
    gop.maj.prop.sunstein[i] <- weighted.mean(gop.maj[keep], weight[keep])
    gop.maj.prop.boot <-  boot(cbind(gop.maj[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    gop.maj.prop.ci.sunstein[i,] <- quantile(gop.maj.prop.boot$t,probs=c(.025,.975))
    unified.prop.sunstein[i] <- weighted.mean(unified[keep], weight[keep])
    unified.prop.boot <-  boot(cbind(unified[keep], weight[keep]), wmean, R = 1000, stype = "i", sim = "ordinary")
    unified.prop.ci.sunstein[i,] <- quantile(unified.prop.boot$t,probs=c(.025,.975))
    }


#2)  CIRCUITS
  
    
mixed.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
mixed.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
mixed.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrr.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrr.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrr.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddd.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddd.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddd.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrd.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrd.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
rrd.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddr.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddr.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
ddr.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
dem.maj.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
dem.maj.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
dem.maj.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
gop.maj.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
gop.maj.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
gop.maj.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
unified.prop.circuit.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
unified.circ.ci.lo.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))
unified.circ.ci.hi.sunstein <- matrix(NA, nrow = loop.sunstein, ncol = length(unique(circuit)))

#for confidence, set lower and upper bounds at 0 and 1 using ifelse commands
#also, in some years, cells will be empty because either circuit does not exist (10th and 11th circuit)
    #or no cases entered the data that year (DC).  us ifelse commands to pull out these years as missing
    
for (i in 1:loop.sunstein){ # loop over year
  keep1 <-  year == unique.year[i]
   for (j in 1:ncol(mixed.prop.circuit.sunstein)){ #then loop over each circuit within year
    keep2 <- circuit == j
      
      empty <-  is.nan(mean(mixed[keep1 & keep2]))#denote circuit-years with no cases
  
      mixed.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(mixed[keep1 & keep2]))
      if (empty == FALSE)      mixed.prop.boot <-  boot(cbind(mixed[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      mixed.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(mixed.prop.boot$t,probs=c(.025)))
      mixed.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(mixed.prop.boot$t,probs=c(.975)))
 
      ddd.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(ddd[keep1 & keep2]))
      if (empty == FALSE)      ddd.prop.boot <-  boot(cbind(ddd[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      ddd.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(ddd.prop.boot$t,probs=c(.025)))
      ddd.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(ddd.prop.boot$t,probs=c(.975)))       
    
      rrr.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(rrr[keep1 & keep2]))
      if (empty == FALSE)      rrr.prop.boot <-  boot(cbind(rrr[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      rrr.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(rrr.prop.boot$t,probs=c(.025)))
      rrr.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(rrr.prop.boot$t,probs=c(.975)))     
              
      rrd.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(rrd[keep1 & keep2]))
      if (empty == FALSE)      rrd.prop.boot <-  boot(cbind(rrd[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      rrd.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(rrd.prop.boot$t,probs=c(.025)))
      rrd.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(rrd.prop.boot$t,probs=c(.975)))
 
      ddr.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(ddr[keep1 & keep2]))
      if (empty == FALSE)      ddr.prop.boot <-  boot(cbind(ddr[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      ddr.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(ddr.prop.boot$t,probs=c(.025)))
      ddr.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(ddr.prop.boot$t,probs=c(.975)))
 
      dem.maj.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(dem.maj[keep1 & keep2]))
      if (empty == FALSE)      dem.maj.prop.boot <-  boot(cbind(dem.maj[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      dem.maj.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(dem.maj.prop.boot$t,probs=c(.025)))
      dem.maj.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(dem.maj.prop.boot$t,probs=c(.975)))
    
      gop.maj.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(gop.maj[keep1 & keep2]))
      if (empty == FALSE)      gop.maj.prop.boot <-  boot(cbind(gop.maj[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      gop.maj.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(gop.maj.prop.boot$t,probs=c(.025)))
      gop.maj.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(gop.maj.prop.boot$t,probs=c(.975)))
 
      unified.prop.circuit.sunstein[i,j] <-  ifelse(empty, NaN, mean(unified[keep1 & keep2]))
      if (empty == FALSE)      unified.prop.boot <-  boot(cbind(unified[keep1 & keep2], weight[keep1 & keep2]), wmean, R = n.bootstraps, stype = "i", sim = "ordinary")
      unified.circ.ci.lo.sunstein[i,j] <-  ifelse(empty, NaN, quantile(unified.prop.boot$t,probs=c(.025)))
      unified.circ.ci.hi.sunstein[i,j] <-  ifelse(empty, NaN, quantile(unified.prop.boot$t,probs=c(.975)))
       
        }
    }



#NOTE: THIS IS FOR YEAR-LEVEL ONLY
#FOR AGGREGATE STATISTICS (POOLING ALL CIRCUITS) 
#check to see whether songer and sunstein correlate over years of overlap (1995 to 2002)
    #this is 71:78 in songer vectors, 1:8
keep1 <- c(1:8)
keep2 <- c(71:78)

#take mean absolute difference
options(digits=2)
print( mean (abs(mixed.prop.sunstein[keep1] - mixed.prop.songer[keep2] )))
print( mean (abs(rrd.prop.sunstein[keep1] - rrd.prop.songer[keep2] )))
print( mean (abs(ddr.prop.sunstein[keep1] - ddr.prop.songer[keep2] )))
print( mean (abs(rrr.prop.sunstein[keep1] - rrr.prop.songer[keep2] )))
print( mean (abs(ddd.prop.sunstein[keep1] - ddd.prop.songer[keep2] )))



##########################################################3
#NOW APPEND SUNSTEIN 2003 AND 2004 % TO SONGER, TO FORM SINGLE TIME SERIES FOR ENTIRE JUDICIARY
#that is, use elements 9 and 10 from sunstein
keep <- c(9:10)

#estimates
mixed.prop.all <- c( mixed.prop.songer, mixed.prop.sunstein[keep])
rrr.prop.all <- c( rrr.prop.songer, rrr.prop.sunstein[keep])
ddd.prop.all <- c( ddd.prop.songer, ddd.prop.sunstein[keep])
rrd.prop.all <- c( rrd.prop.songer, rrd.prop.sunstein[keep])
ddr.prop.all <- c( ddr.prop.songer, ddr.prop.sunstein[keep])
dem.maj.prop.all <- c( dem.maj.prop.songer, dem.maj.prop.sunstein[keep])
gop.maj.prop.all <- c( gop.maj.prop.songer, gop.maj.prop.sunstein[keep])
unified.prop.all <- c( unified.prop.songer, unified.prop.sunstein[keep])

#confidence intervals
mixed.prop.ci.all <- rbind( mixed.prop.ci.songer, mixed.prop.ci.sunstein[keep,])
rrr.prop.ci.all <- rbind( rrr.prop.ci.songer, rrr.prop.ci.sunstein[keep,])
ddd.prop.ci.all <- rbind( ddd.prop.ci.songer, ddd.prop.ci.sunstein[keep,])
rrd.prop.ci.all <- rbind( rrd.prop.ci.songer, rrd.prop.ci.sunstein[keep,])
ddr.prop.ci.all <- rbind( ddr.prop.ci.songer, ddr.prop.ci.sunstein[keep,])
dem.maj.prop.ci.all <- rbind( dem.maj.prop.ci.songer, dem.maj.prop.ci.sunstein[keep,])
gop.maj.prop.ci.all <- rbind( gop.maj.prop.ci.songer, gop.maj.prop.ci.sunstein[keep,])
unified.prop.ci.all <- rbind( unified.prop.ci.songer, unified.prop.ci.sunstein[keep,])

############################################
#FOR CIRCUITS, ROW BIND LAST 2 ROWS OF MATRIX TO SONGER DATA FOR 2003 AND 2004

mixed.prop.circuit.all <- rbind(mixed.prop.circuit.songer, mixed.prop.circuit.sunstein[keep,])
mixed.circ.ci.lo.all <- rbind(mixed.circ.ci.lo.songer,mixed.circ.ci.lo.sunstein[keep,])
mixed.circ.ci.hi.all <- rbind(mixed.circ.ci.hi.songer,mixed.circ.ci.hi.sunstein[keep,])

rrr.prop.circuit.all <- rbind(rrr.prop.circuit.songer, rrr.prop.circuit.sunstein[keep,])
rrr.circ.ci.lo.all <- rbind(rrr.circ.ci.lo.songer,rrr.circ.ci.lo.sunstein[keep,])
rrr.circ.ci.hi.all <- rbind(rrr.circ.ci.hi.songer,rrr.circ.ci.hi.sunstein[keep,])

ddd.prop.circuit.all <- rbind(ddd.prop.circuit.songer, ddd.prop.circuit.sunstein[keep,])
ddd.circ.ci.lo.all <- rbind(ddd.circ.ci.lo.songer,ddd.circ.ci.lo.sunstein[keep,])
ddd.circ.ci.hi.all <- rbind(ddd.circ.ci.hi.songer,ddd.circ.ci.hi.sunstein[keep,])

rrd.prop.circuit.all <- rbind(rrd.prop.circuit.songer, rrd.prop.circuit.sunstein[keep,])
rrd.circ.ci.lo.all <- rbind(rrd.circ.ci.lo.songer,rrd.circ.ci.lo.sunstein[keep,])
rrd.circ.ci.hi.all <- rbind(rrd.circ.ci.hi.songer,rrd.circ.ci.hi.sunstein[keep,])

ddr.prop.circuit.all <- rbind(ddr.prop.circuit.songer, ddr.prop.circuit.sunstein[keep,])
ddr.circ.ci.lo.all <- rbind(ddr.circ.ci.lo.songer,ddr.circ.ci.lo.sunstein[keep,])
ddr.circ.ci.hi.all <- rbind(ddr.circ.ci.hi.songer,ddr.circ.ci.hi.sunstein[keep,])

dem.maj.prop.circuit.all <- rbind(dem.maj.prop.circuit.songer, dem.maj.prop.circuit.sunstein[keep,])
dem.maj.circ.ci.lo.all <- rbind(dem.maj.circ.ci.lo.songer,dem.maj.circ.ci.lo.sunstein[keep,])
dem.maj.circ.ci.hi.all <- rbind(dem.maj.circ.ci.hi.songer,dem.maj.circ.ci.hi.sunstein[keep,])

gop.maj.prop.circuit.all <- rbind(gop.maj.prop.circuit.songer, gop.maj.prop.circuit.sunstein[keep,])
gop.maj.circ.ci.lo.all <- rbind(gop.maj.circ.ci.lo.songer,gop.maj.circ.ci.lo.sunstein[keep,])
gop.maj.circ.ci.hi.all <- rbind(gop.maj.circ.ci.hi.songer,gop.maj.circ.ci.hi.sunstein[keep,])

unified.prop.circuit.all <- rbind(unified.prop.circuit.songer, unified.prop.circuit.sunstein[keep,])
unified.circ.ci.lo.all <- rbind(unified.circ.ci.lo.songer,unified.circ.ci.lo.sunstein[keep,])
unified.circ.ci.hi.all <- rbind(unified.circ.ci.hi.songer,unified.circ.ci.hi.sunstein[keep,])



#############
#CREATE GRAPHS
################
#DESCRIPTIVE PLOTS


#############################
#1st plot for paper: 2 stacked time series:
#  a) proportion of dems on judiciary -- prop. of dem majority panels (DDD| DDR)
#  b)  proportion of gop on judiciary -- prop. of gop majority panels (RRR| RRD)   

#Panel Estimates are of length 39 for pooling every 2 years -- extend to align with partisan_distribution_measures
dem.maj.prop.all.graph <- rep(dem.maj.prop.all,each= 2)
duplicate <- seq(1, length(active.dems.per.year-1), by = 2)#designate cells to switch to NA
dem.maj.prop.all.graph[duplicate]<- NA
#set missing values to mean of value above and below it so we can plot lines
for (i in 2:length(dem.maj.prop.all.graph)) {#start with 2 so 1st element remains bland
  if (is.na(dem.maj.prop.all.graph[i])) dem.maj.prop.all.graph[i] <- (dem.maj.prop.all.graph[i-1] +  dem.maj.prop.all.graph[i+1])/2
  }
gop.maj.prop.all.graph <- rep(gop.maj.prop.all,each= 2)
duplicate <- seq(1, length(active.dems.per.year-1), by = 2)#designate cells to switch to NA
gop.maj.prop.all.graph[duplicate]<- NA
#set missing values to mean of value above and below it so we can plot lines
for (i in 2:length(gop.maj.prop.all.graph)) {#start with 2 so 1st element remains bland
  if (is.na(gop.maj.prop.all.graph[i])) gop.maj.prop.all.graph[i] <- (gop.maj.prop.all.graph[i-1] +  gop.maj.prop.all.graph[i+1])/2
  }

#USE THIS FOR 1ST PLOT AND IF USING EACH YEAR
x <- 1925:2004#for years
x.axis.label <- c(1926,seq(1932,2004, by = 8))
x.axis.ticks <- seq(1926, 2004, by = 2)
axis.text.size <- 1.5
text.size <- 1.5


##################
#FIGURE 1
################

pdf("figure1.pdf", height =8, width = 12)   

par(mfrow = c(2,1), mar = c(3,3.5,1,3))
#dems
plot(x,active.dems.per.year, main = "", type = "n", xlab = "", ylab = "",
    xlim = c(min(x), max(x)), ylim = c(0,1), axes = F, xaxs = "i", yaxs = "i")#call empty plot
source("dem_control_shading.R")#add shading
points(x,active.dems.per.year,type ="l")#percent dem among active judges
lines(x,  dem.maj.prop.all, type = "l", lty = 2, col = "blue") #percent of panels with maj dems
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))    
text(1931.3,.65, "Proportion of\nDemocratic\nappointees",cex = text.size)
arrows(1936,.6,1937,.5,length=.1)
text(1944.5,.3, "Proportion of\npanels with a\nDemocratic majority", cex = text.size)
arrows(1938.3,.24,1937.4,.31,length=.1)
box(bty="l")


#gop
plot(x,1-active.dems.per.year, main = "", type = "n", xlab = "", ylab = "",
    xlim = c(min(x), max(x)), ylim = c(0,1), axes = F, xaxs = "i", yaxs = "i")#call empty plot
source("dem_control_shading.R")#add shading
points(x,1-active.dems.per.year,type ="l", col = "dark green")#percent dem among active judges
points(x,  gop.maj.prop.all, type = "l", lty = 2, col = "red") #percent of panels with maj dems
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0)) 
text(1931.4,.4, "Proportion of\nRepublican\nappointees",cex = text.size)
arrows(1932,.6,1933,.7,length=.1)
text(1947.5,.8, "Proportion of\npanels with a\nRepublican majority",cex = text.size)
arrows(1939,.68,1937.5,.68,length=.1)
box(bty="l")

dev.off()
   


##################
#FIGURE 2
################

rrr.color <- "red"
ddd.color <- "blue"
rrd.color <- "purple"
ddr.color <- "darkcyan"
mixed.color <- "brown"

pdf("figure2.pdf", height = 10, width = 8)   

par(mfrow = c(4,1), mar = c(3,3.5,1,3))
plot(x, rrr.prop.all, type = "n", ylim = c(0,1), xlim = c(min(x), max(x)), main ="",
    xlab = "", ylab = "", axes = F, xaxs = "i", yaxs = "i")
source("dem_control_shading.R")#add shading, which is contained in a separate script
#draw plot
points(x, mixed.prop.all, type = "l", lty = 1, col = mixed.color)   # mixed panels
points(x, ddd.prop.all, type = "l", lty = 6, col = ddd.color, lwd = 2)   # ddd panels
points(x, rrr.prop.all, type = "l", lty = 2, col = rrr.color)   # rrd panels
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
box(bty="l")
text(1944, .7, "Mixed", cex = text.size, col =  mixed.color)
text(1936, .35, "RRR", cex = text.size, col =  rrr.color)#labels for each line
text(1947, .34, "DDD ", cex = text.size, col =  ddd.color)#labels for each line
#b)RRD and DDR

plot(x, rrd.prop.all, type = "n", ylim = c(0,1), xlim = c(min(x), max(x)), main ="",
    xlab = "", ylab = "", axes = F, xaxs = "i", yaxs = "i")
source("dem_control_shading.R")#add shading, which is contained in a separate script
#draw plot
points(x, ddr.prop.all, type = "l", lty = 1, col = ddr.color)
points(x, rrd.prop.all, type = "l", lty = 2, col = rrd.color)
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
box(bty="l")
text(1935, .62, "RRD", cex = text.size, col = rrd.color)
text(1954, .56, "DDR", cex = text.size, col = ddr.color)#labels for each line

#c)DDD and DDR

plot(x, ddd.prop.all, type = "n", ylim = c(0,1), xlim = c(min(x), max(x)), main ="",
    xlab = "", ylab = "", axes = F, xaxs = "i", yaxs = "i")
source("dem_control_shading.R")#add shading, which is contained in a separate script
#draw plot
points(x, ddr.prop.all, type = "l", lty = 2, col = ddr.color)
points(x, ddd.prop.all, type = "l", lty = 1, col = ddd.color)
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
box(bty="l")
text(1973, .25, "DDD", cex = text.size, col = ddd.color)
text(1929, .29, "DDR", cex = text.size, col = ddr.color)#labels for each line

#d)RRR and RRD

plot(x, rrd.prop.all, type = "n", ylim = c(0,1), xlim = c(min(x), max(x)), main ="",
    xlab = "", ylab = "", axes = F, xaxs = "i", yaxs = "i")
source("dem_control_shading.R")#add shading, which is contained in a separate script
#draw plot
points(x, rrr.prop.all, type = "l", lty = 1, col = rrr.color)
points(x, rrd.prop.all, type = "l", lty = 2, col = rrd.color)
abline(h=.5, lty =2, col = "light gray")
axis(1, at = x.axis.ticks , labels = NA, mgp = c(2,1.5,0))
axis(1, at = x.axis.label, labels = x.axis.label, cex.axis = axis.text.size)
axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0))
box(bty="l")
text(1935, .62, "RRD", cex = text.size, col = rrd.color)
text(1936, .1, "RRR", cex = text.size, col =  rrr.color)#labels for each line

dev.off()


#Circuit Plots
#for each circuit, plot:
#  a) Proportion of Dems
#  #  c) Propotion of Mixed Panels
#plot proportion of dems and proportion of mixed panels  (with dotted line)
#set up circuit proportions data
detach()
attach.all(partisan.composition, overwrite=T)
unique.year <- unique(year)

#get dem proportion and  probabilities of circuit types by circuit in same format as circuit mixed panels matrix
#use this to do every year

dem.prop.circuit <- matrix(NA, nrow = length(unique.year), ncol = length(unique(circuit)))
for (i in 1:length(unique.year)){ # loop over year
    keep1 <- year == unique.year[i]
        for (j in 1:ncol(dem.prop.circuit)){ #then loop over each circuit within year
    keep2 <- circuit== j
    empty <- (j == 10 & i < 5) | (j==11 & i < 58) | (j == 12 & (i>=7 & i<=13)) #allow empty cells for 10th and 11th circuits before they were created
                                                                                #and years when DC circ had no 3 judge panels
    dem.prop.circuit[i,j] <-  ifelse(empty, NaN,dem.prop.active.circyear[keep1 & keep2])
     }
}

circuit.label <- c("1st", "2nd", "3rd", "4th", "5th","6th","7th","8th","9th","10th","11th", "DC")

#use this for diagonal x-axis labels
axis.text <- 2
axis.diag <- function(axis=1,labels,at,cex, srt=45,tcl=-.7, y = par("usr")[3]-.02) {
 axis(axis,at=at,labels=FALSE,tcl=tcl)
 text(at,y,srt=srt,adj=1,labels=labels,xpd=TRUE, font = 1, cex = axis.text)
}
left.panel.width <- .6 # control width of left-hand panel (see "widths" in layout command)
second.panel.width <- 2.5 #so there's more room to put y-axis labels on 1st col. of plots

pdf("figure3.pdf", height = 12, width = 12)    

#use this for 4x3
layout(rbind(c(1, 2,3), c(4, 5,6), #order so that each circuit first, then y-axis labels, then x-axis lables
    c(7, 8, 9), c(10, 11, 12)),
    widths = c(1,.92,1),#make 2nd column smaller since it has no vertical axes
    heights = c(1,1,1,1.1))
#layout.show(12)
#do 1st 9 circuits 1st (1st 3 rows)
for (i in 1:9){
#set margins based on which column each plot will appear
if (i == 1 | i == 4 | i == 7 | i == 10) par(mar=c(3,3.3,3,2))
if (i == 2 | i == 5 | i ==8 | i == 11)  par(mar=c(3,2,3,2))
if (i == 3 | i == 6 | i == 9 | i == 12) par(mar=c(3,2,3,3.3))

plot(x, mixed.prop.circuit.all[,i], type = "n", xlim = c(1924,2004), ylim = c(0,1), axes = F, xlab = "",
    ylab = "", main = "", xaxs ="i", yaxs ="r")
source("dem_control_shading.R")#add shading
abline(h = .5, col = "light grey", lwd =.7)
points(x, mixed.prop.circuit.all[,i], type = "l", col = "purple")#mixed
points(x, dem.prop.circuit[,i], type = "l", lty =2, col = "blue")#dems control 
#points(x, unified.prop.circuit.all[,i], type = "l", lty =4, col = "dark green")#unified panels

#x-axis
axis(1, at = x.axis.label, label = NA) #add tick marks on x-axis for 1st 3 rows

if (i == 1 | i == 4 | i == 7 | i == 10) axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0)) #y-axis for 1st col
if (i != 1 & i != 4 & i != 7 & i != 10) axis(2, at = c(0,.25,.50,.75,1), 
    label = NA,  mgp = c(2,.5,0), las = 1) #y-axis for 2nd and 3rd col
if (i == 3 | i == 6 | i == 9 | i == 12) axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0)) #y-axis for 1st col #right y-axis for 3rd col
mtext(circuit.label[i], 3, line = .5, cex = 2)
#add labels just to first panel
if (i==1) text(1950,.9, "Proportion\nof Dems", cex = 1.2)
if (i==1) text(1950,.3, "Mixed\nPanels",cex = 1.2)

box()

}
#do last 3 circuits (4th row)
for (i in 10:12){
#set margins based on which column each plot will appear
if (i == 10) par(mar=c(6,3.3,3,2))#adjust for more bottom margin
if (i== 11)  par(mar=c(6,2,3,2))
if (i  == 12) par(mar=c(6,2,3,3.3))
plot(x, mixed.prop.circuit.all[,i], type = "n", xlim = c(1924,2004), ylim = c(0,1), axes = F, xlab = "",
    ylab = "", main = "", xaxs ="i", yaxs ="r")
source("dem_control_shading.R")#add shading
abline(h = .5, col = "light grey", lwd =.7)
points(x, mixed.prop.circuit.all[,i], type = "l", col = "purple")#mixed
points(x, dem.prop.circuit[,i], type = "l", lty =2, col = "blue")#dems control 
#points(x, unified.prop.circuit.all[,i], type = "l", lty =4, col = "dark green")#unified panels

axis(1,at = x.axis.label, label=NA, tcl = -1)#use longer tick marks for bottom row
axis.diag(1, at = x.axis.label, y=-.12, label = x.axis.label, srt=45, cex = 2)
if (i == 1 | i == 4 | i == 7 | i == 10)axis(2, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0)) #y-axis for 1st col
if (i != 1 & i != 4 & i != 7 & i != 10) axis(2, at = c(0,.25,.50,.75,1), 
    label = NA,  mgp = c(2,.5,0), las = 1) #y-axis for 2nd and 3rd col
    mtext(circuit.label[i], 3, line = .5, cex = 2)
if (i == 3 | i == 6 | i == 9 | i == 12) axis(4, at = seq(0,1, by = .25), labels = c(0, ".25", ".5", ".75",1), las =2,  cex.axis = axis.text.size,
    mgp = c(2,.7,0)) #y-axis for 1st col #right y-axis for 3rd col
box()
}

dev.off()


#The rest of the script converts the info presented in the graphs into tables, which appear in Web Appendix B


#set up exports for tables
#1) aggregate data (by year
#years in rows
  #in columns:  % dem appointees, % of panels with dem maj, 
    #RRR, RRD, DDR, DDD, Mixed

rd <- function (x) {#create shortcut for rounding to 2nd decimal place
  round(x,digits = 2)
  }
  
export.matrix <- cbind(
  x, 
  rd(active.dems.per.year),
  paste(rd(dem.maj.prop.ci.all[,1]),rd(dem.maj.prop.all),rd(dem.maj.prop.ci.all[,2]), sep = ", "),
  paste(rd(rrr.prop.ci.all[,1]),rd(rrr.prop.all),rd(rrr.prop.ci.all[,2]), sep = ", "),
  paste(rd(rrd.prop.ci.all[,1]),rd(rrd.prop.all),rd(rrd.prop.ci.all[,2]), sep = ", "),
  paste(rd(ddr.prop.ci.all[,1]),rd(ddr.prop.all),rd(ddr.prop.ci.all[,2]), sep = ", "),
  paste(rd(ddd.prop.ci.all[,1]),rd(ddd.prop.all),rd(ddd.prop.ci.all[,2]), sep = ", "),
  paste(rd(mixed.prop.ci.all[,1]),rd(mixed.prop.all),rd(mixed.prop.ci.all[,2]), sep = ", ")
)

write.csv(export.matrix, "panel_rates_yearly_table.csv")

#############
#CIRCUIT YEAR DATA
#      1ST CIRCUIT     2ND CIRCUIT ......
#YEAR   % Dems  %Mixed % Dems % Mixed 

export.matrix.2 <- cbind(
  x,
  rd(dem.prop.circuit[,1]), 
  paste(rd(mixed.circ.ci.lo.all[,1]), rd(mixed.prop.circuit.all[,1]), rd(mixed.circ.ci.hi.all[,1]), sep =", "),
  rd(dem.prop.circuit[,2]), 
  paste(rd(mixed.circ.ci.lo.all[,2]), rd(mixed.prop.circuit.all[,2]), rd(mixed.circ.ci.hi.all[,2]), sep =", "),
  rd(dem.prop.circuit[,3]), 
  paste(rd(mixed.circ.ci.lo.all[,3]), rd(mixed.prop.circuit.all[,3]), rd(mixed.circ.ci.hi.all[,3]), sep =", "),
  rd(dem.prop.circuit[,4]), 
  paste(rd(mixed.circ.ci.lo.all[,4]), rd(mixed.prop.circuit.all[,4]), rd(mixed.circ.ci.hi.all[,4]), sep =", "),
  rd(dem.prop.circuit[,5]), 
  paste(rd(mixed.circ.ci.lo.all[,5]), rd(mixed.prop.circuit.all[,5]), rd(mixed.circ.ci.hi.all[,5]), sep =", "),
  rd(dem.prop.circuit[,6]), 
  paste(rd(mixed.circ.ci.lo.all[,6]), rd(mixed.prop.circuit.all[,6]), rd(mixed.circ.ci.hi.all[,6]), sep =", "),
  rd(dem.prop.circuit[,7]), 
  paste(rd(mixed.circ.ci.lo.all[,7]), rd(mixed.prop.circuit.all[,7]), rd(mixed.circ.ci.hi.all[,7]), sep =", "),
  rd(dem.prop.circuit[,8]), 
  paste(rd(mixed.circ.ci.lo.all[,8]), rd(mixed.prop.circuit.all[,8]), rd(mixed.circ.ci.hi.all[,8]), sep =", "),
  rd(dem.prop.circuit[,9]), 
  paste(rd(mixed.circ.ci.lo.all[,9]), rd(mixed.prop.circuit.all[,9]), rd(mixed.circ.ci.hi.all[,9]), sep =", "),
  rd(dem.prop.circuit[,10]), 
  paste(rd(mixed.circ.ci.lo.all[,10]), rd(mixed.prop.circuit.all[,10]), rd(mixed.circ.ci.hi.all[,10]), sep =", "),
  rd(dem.prop.circuit[,11]), 
  paste(rd(mixed.circ.ci.lo.all[,11]), rd(mixed.prop.circuit.all[,11]), rd(mixed.circ.ci.hi.all[,11]), sep =", "),
  rd(dem.prop.circuit[,12]), 
  paste(rd(mixed.circ.ci.lo.all[,12]), rd(mixed.prop.circuit.all[,12]), rd(mixed.circ.ci.hi.all[,12]), sep =", ")
  )

write.csv(export.matrix.2, "panel_rate_circuit_table.csv")