# China Party Congress Data # Estimate and Intepret Bayesian Partial Rank Models # # Party Congress 12, 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 <- 12 # Congress number mr <- 1 # Model run number yr <- 1982 # 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 <- 9 centcmtemin <- 10 centcmtemax <- 272 accmin <- 273 accmax <- 410 # 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 maodummy <- maofi>0 dengdummy <- dengfi>0 huyaobangdummy <- huyaobangfi>0 zhaoziyangdummy <- zhaoziyangfi>0 chenyundummy <- chenyunfi>0 yejianyingdummy <- yejianyingf>0 ccp1935 <- ptime<=1935 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", "maodummy","dengdummy","huyaobangdummy","chenyundummy","yejianyingdummy","farmy2", "ccp1935","ccp1945", "lmarch", "gdp", "fisrevenue" ) X <- cbind(female,nonhan,prince, apptage,partyexp, highschool,college,postgrad, maodummy,dengdummy,huyaobangdummy,chenyundummy,yejianyingdummy,farmy2, ccp1935,ccp1945, lmarch, gdp, fisrevenue ) varnames <- c("Female","Ethnic Minority","Prince", "Age","Party Exp", "High School","College","Graduate Degree", "Mao Faction","Deng Faction","Hu Yaobang Fact","Chen Yun Fact","Ye Jianying Fact","2nd Field Army", "Party by 1935","Party by 1945", "Long Marcher", "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=""))