# China Party Congress Data # Estimate and Intepret Bayesian Partial Rank Models # # Party Congress 14, Model 1 # # Chris Adolph # Load libraries rm(list = ls()) library(MASS) library(WhatIf) library(tile) library(partialrank) source("setupCF.r") # Controls dorun <- TRUE ignoreranks <- FALSE # Indexes cn <- 14 # Congress number mr <- 1 # Model run number yr <- 1992 # Year iterations <- 10000 # MCMC Iterations burnin <- 5000 # Discarded Iterations thin <- 10 # Skipped Iterations chains <- 3 # MCMC chains ties <- 1 # Eliminate ties up to ci <- 0.95 # CIs for expected value simulations ci2 <- c(0.67,0.95) # CIs for first diff simulations # Boundaries of political bodies politburomin <- 1 politburomax <- 7 centcmtemin <- 8 centcmtemax <- 190 accmin <- 191 accmax <- 320 # Read data data <- read.csv(paste("china",cn,"baselineReplicate.csv",sep=""),header=TRUE) attach(data) # Process covariate data agege63 <- apptage>=63 nonhan <- ethnic>0 female <- 1-gender partyyears <- yr - ptime partyexp <- partyyears/apptage lesshighschool <- (edu==0) highschool <- (edu==1) college <- (edu==2) postgrad <- (edu==3) cntminexp <- e1 cntptyexp <- e2 locgovexp <- e3-e4 cstgovexp <- e4 plaexp <- e5 provleader <- gdp!=0 dengdummy <- dengfi>0 jiangzemindummy <- jiangzeminfi>0 zhaoziyangdummy <- zhaoziyangfi>0 ccp1945 <- (ptime<=1945) ccp1949 <- (ptime<=1949) farmy2 <- rep(0,nrow(data)) for (i in 1:length(farmy2)) if (!is.na(farmy[i])) farmy2[i] <- (farmy[i]==2) escoredata <- cbind(cntminexp,cntptyexp,locgovexp,cstgovexp,plaexp) # Select covariates xstrings <- c("female","nonhan","prince", "apptage","partyexp", "highschool","college","postgrad", "zhaoziyangdummy","dengdummy","jiangzemindummy", "ccp1949", "gdp", "fisrevenue" ) X <- cbind(female,nonhan,prince, apptage,partyexp, highschool,college,postgrad, zhaoziyangdummy,dengdummy,jiangzemindummy, ccp1949, gdp, fisrevenue ) varnames <- c("Female","Ethnic Minority","Prince", "Age","Party Exp", "High School","College","Graduate Degree", "Zhao Ziyang Faction","Deng Faction","Jiang Zemin Fact", "Party by 1949", "GDP Relative Growth", "Revenue Relative Growth" ) # Treat ACC as unknown? if (ignoreranks) { for (i in 1:length(newrank)) { if (!is.na(newrank[i])&&(newrank[i]>(accmin-1))) { newrank[i] <- NA upper[i] <- accmin lower[i] <- accmax } } } # Run latent strength rank model using rank likelihood n <- nrow(X) filename <- paste("china",yr,"m",mr,"run",iterations,"ties",ties,".Rdata",sep="") if (dorun) { result <- partialRankBayes(y=newrank, X=X, lower=lower, upper=upper, ties=ties, meancenter=TRUE, n.chains=chains, n.iter=iterations, n.thin=thin, n.burnin=burnin ) save(result,file=filename) } else { load(filename) } # Save parameter table paramfilename <- paste("param",yr,"m",mr,"run",iterations,"ties",ties,sep="") paramtable(paramfilename,result,varnames) # Setup counterfactuals cfscen <- setupCF(X,xstrings,yr,fd=FALSE) scenarios <- cfscen$scenarios newdata <- cfscen$x xpre <- cfscen$xpre cfscenFD <- setupCF(X,xstrings,yr,fd=TRUE) scenariosFD <- cfscenFD$scenarios newdataFD <- cfscenFD$x xpreFD <- cfscenFD$xpre # Estimate fitted values fittedfile <- paste("pc",cn,"fittedrank_model",mr,"run",iterations,"ties",ties,sep="") fittedrank <- predict(result,ranefs=0) save(fittedrank,file=fittedfile) # Estimate conditional expected values simrankfile <- paste("pc",cn,"simrank_model",mr,"run",iterations,"ties",ties,sep="") simrank3 <- predict(result,newdata,level=ci2) save(simrank3,file=simrankfile) # Estimate conditional expected values with low random effect simranklowerfile <- paste("pc",cn,"simranklower_model",mr,"run",iterations,"ties",ties,sep="") simrank3lower <- predict(result,newdata, level=ci, ranef=mean(result$sims.list$alpha)-sd(as.numeric(result$sims.list$alpha))) save(simrank3lower,file=simranklowerfile) # Estimate conditional expected values with high random effects simrankupperfile <- paste("pc",cn,"simrankupper_model",mr,"run",iterations,"ties",ties,sep="") simrank3upper <- predict(result,newdata, level = ci, ranef=mean(result$sims.list$alpha)+sd(as.numeric(result$sims.list$alpha))) save(simrank3upper,file=simrankupperfile) # Estimate first differences simfdrankfile <- paste("pc",cn,"simfdrank_model",mr,"run",iterations,"ties",ties,sep="") simfdrank3 <- predict(result,newdataFD,xpreFD,level=ci2) save(simfdrank3,file=simfdrankfile) fdtablefile <- paste("efirstdiffranknewpc",cn,"m",mr,"run",iterations,"ties",ties,sep="") firstdiffTable(fdtablefile,simfdrank3,n,scenariosFD) # Print first differences fdres<-cbind(simfdrank3$rank,t(simfdrank3$lower)[,2],t(simfdrank3$upper)[,2]) rownames(fdres) <- scenariosFD print(fdres) print((1-fdres/n)*100-100) fdres2<-cbind(simfdrank3$rank,t(simfdrank3$lower)[,1],t(simfdrank3$upper)[,1]) rownames(fdres2) <- scenariosFD print(fdres2) print((1-fdres2/n)*100-100) # Plot conditional expected values thresholds = c(0,politburomax,centcmtemax,accmax) source("predictRanklik.R") plotranksrl(simrank3, names=scenarios, n=n, interval=c("confidence"), limits=NA, at=NA, thresholds=thresholds, comparerank=cbind(simrank3lower$rank, simrank3upper$rank), outfile=paste("eranksnewpc",cn,"m",mr,"run",iterations,"ties",ties,sep=""), rows=50, width=7, col=c("black","blue","blue"), entryheight=0.05, tiercol=c("gray85","gray70","gray55"), vertlines=(n-simrank3$rank[length(simrank3$rank)])/n*100 ) # Plot fitted values plotranksrl(fittedrank, names=name1, n=n, interval=c("confidence"), limits=NA, at=NA, thresholds=thresholds, comparerank=newrank, comparelower=lower, compareupper=upper, outfile=paste("frankspeoplepc",cn,"m",mr,"run",iterations,"ties",ties,sep=""), maintitle=paste("Fitted Values, ",cn,"th Congress, Model ",mr,sep=""), rows=50, width=7, col="black", lwd=1.5, lwdcompare=0.25, colcompare="red", tiercol=c("gray85","gray70","gray55"), vertlines=(n-simrank3$rank[length(simrank3$rank)])/n*100 ) # Plot random effects by person plotre(result,name1, outfile=paste("ranefxnew_pc",cn,"m",mr,"run",iterations,"ties",ties,sep=""), maintitle=paste(cn,"th Congress, Model ",mr,sep=""), rows=50) ## Estimate all age and party tenure interactions # Loop over party join date and age, construct hyp data agepartyCFs <- agepartyCF(X,xstrings,yr,ptime,apptage) hyppartyexp <- agepartyCFs$hyppartyexp hypptime <- agepartyCFs$hypptime hypage <- agepartyCFs$hypage agepartyCFs <- agepartyCFs$newdata # Check for inclusion in convex hull whatif <- whatif(data=cbind(apptage,ptime), cfact=cbind(hypage,hypptime) ) inhull <- as.logical(whatif$in.hull) hyppartyexp <- hyppartyexp[inhull] hypptime <- hypptime[inhull] hypage <- hypage[inhull] agepartyCFs <- agepartyCFs[inhull,] # Simulate age/party scenarios agepartyfile <- paste("pc",cn,"ageparty_model",mr,"run",iterations,"ties",ties,sep="") simageparty <- predict(result,agepartyCFs,level=ci) save(simageparty,file=agepartyfile) # Make output table ageout <- cbind(simageparty$rank,simageparty$upper,simageparty$lower) pageout <- (n-ageout)/n*100 ageout <- as.data.frame(cbind(hypage,hypptime,hyppartyexp,ageout,pageout)) varsout <- c("age","ptime","partyexp","rank","upper","lower","pctrank","pctupper","pctlower") colnames(ageout) <- varsout write.csv(ageout,file=paste("pc",cn,"agesim_model",mr,"run",iterations,"ties",ties,".csv",sep=""))