# This code is to accompany Maximum Likelihood Methods Strategies for Social Science, # Michael D. Ward and John S. Ahlquist, 2018, Cambridge University Press, # ISBN 9781107185821. The code is presented as is, but is designed to replicate all of # the material in the book, chapter by chapter. Any comments, suggestions, or problems can # be reported to the authors: michael.don.ward@gmail.com and John Ahlquist jahlquist@ucsd.edu. # These programs were designed to run on a variety of platforms, but currently uses # platform x86_64-apple-darwin15.6.0 # arch x86_64 # os darwin15.6.0 # system x86_64, darwin15.6.0 # version.string R version 3.4.4 (2018-03-15) # nickname Someone to Lean On # We used Rstudio, Version 1.1.447 # This code requires installation of the loaded libraries. rm(list=ls()); gc() library(foreign); library(MASS); library(rms); library(stargazer);library(RColorBrewer) set.seed(123) ### Figure 8.1 ####################################### x1<-rlogis(10000,0,1) x2<-rlogis(10000,2,1) plot(density(x1,adjust=2), xlim=c(-6,6), axes=F, main="", xlab = "Y*", ylab="", frame.plot=F) lines(density(x2, adjust=2), lty=2) segments(x0=c(-2,-.9, .7,2.4), y0=0, y1=dlogis(c(-2,-.9, .7,2.4)), col=gray(.5)) xlabs<-c(expression(tau[SD-D]), expression(tau[D-DK]), expression(tau[DK-A]), expression(tau[A-SA])) axis(1, at=c(-2,-.9, .7,2.4), col=grey(.05), labels=xlabs) curvelabs<-c(expression(bold(x)[i]^T~~beta), expression(bold(x)[j]^T~~beta)) text(c(0,2), c(.22,.22), label=curvelabs) text(3.5,.01,"strongly \n agree", cex=.8) text(1.3,.01,"agree", cex=.8) ###################################################### # Moving on anes<-read.csv("anes_pilot_2016.csv", header=TRUE) anes$qtime<-apply(anes[,239:374], 1,sum)/60 #minutes to take the survey anes$age<-2016-anes$birthyr use<-c("caseid", "weight", "age", "qtime", "faminc", "pid7", "educ", "race", "givefut", "follow", "incgap20", "getahead", "disc_m", "disc_b","lcbo", "lazyb") anes.sub<-anes[,colnames(anes)%in%use] anes.sub$faminc[anes.sub$faminc>16]<-NA anes.sub$lcbo[anes.sub$lcbo>7]<-NA anes.sub$disc_b[anes.sub$disc_b>7]<-NA anes.sub$lazyb[anes.sub$lazyb>7]<-NA anes.sub$pid7[anes.sub$pid7>7]<-NA anes.sub$incgap20[anes.sub$incgap20>5]<-NA anes.sub$disc_m[anes.sub$disc_m>7]<-NA anes.sub$white<-anes.sub$race==1 anes.sub$getahead[anes.sub$getahead>7]<-NA names(anes.sub)<-c("caseid", "weight", "follow.politics", "political.donation", "lazy.black", "Obama.LR", "black.disc", "men.disc", "ineq.up", "get.ahead", "race", "education", "income", "pid","survey.time", "age", "white" ) anes.sub<-na.omit(anes.sub) anes.sub$blacks.discrimination<- -1*anes.sub$black.disc + 6 #bigger numbers => more descrimination anes.sub$blacks.lazy<- -1*anes.sub$lazy.black+6 #bigger numbers => more lazy ol.r<-polr.out <- polr(as.ordered(blacks.discrimination)~ blacks.lazy, data=anes.sub, method="logistic", Hess=T) #homework! polr.out <- polr(as.ordered(Obama.LR)~ follow.politics + pid + age + education +income + white, data=anes.sub, method="logistic", Hess=T) summary(polr.out) beta <- coef(polr.out) tau <- polr.out$zeta attach(anes.sub) #create predicted probabilities X <- cbind(median(follow.politics), min(pid):max(pid), median(age), median(education), median(income), TRUE) p1 <- plogis(tau[1] - X %*% beta) p2 <- plogis(tau[2] - X %*% beta) - plogis(tau[1] - X %*% beta) p3 <- plogis(tau[3] - X %*% beta) - plogis(tau[2] - X %*% beta) p4 <- plogis(tau[4] - X %*% beta) - plogis(tau[3] - X %*% beta) p5 <- plogis(tau[5] - X %*% beta) - plogis(tau[4] - X %*% beta) p6 <- plogis(tau[6] - X %*% beta) - plogis(tau[5] - X %*% beta) p7 <- 1.0 - plogis(tau[6] - X %*% beta) ### Figure 8.2 ####################################### gsc<-rev(brewer.pal(7,"Greys")) #pdf(file = "ANESord.pdf") plot(0:6, p1, type="b", col=gsc[1], ylim=c(0,1),las=1,ylab=expression(hat(pi)), xlab="",lwd=2,bty="n",xaxt="n") lines(0:6, p2, col=gsc[2],lwd=2, type="b") lines(0:6, p3, col=gsc[3],lwd=2, type="b") lines(0:6, p4, col=gsc[4],lwd=2, type="b") lines(0:6, p5, col=gsc[5],lwd=2, type="b") lines(0:6, p6, col=gsc[6],lwd=2, type="b") lines(0:6, p7, col=gsc[7],lwd=2, type="b") legend(0, 1, legend=c("Pr(extremely liberal)", "Pr(liberal)", "Pr(slightly liberal)", "Pr(moderate)", "Pr(slightly conservative)", "Pr(conservative)", "Pr(extremely conservative)"), col=gsc, lty=1,lwd=2) Axis(x = 0:6, at = 0:6, side=1, labels = c("Strong Dem","WD","ID","Ind","IR","WR","Strong Rep")) #dev.off() ###################################################### ### Figure 8.3 ####################################### X.wd <- cbind(median(follow.politics), 1, #PID=1 weak democrat median(age), median(education), median(income), TRUE) X.ind <- cbind(median(follow.politics), 3, #PID=3 independent median(age), median(education), median(income), TRUE) draws<-mvrnorm(1000, c(coef(polr.out),polr.out$zeta), solve(polr.out$Hessian)) #1000 draws from posterior; note inclusion of cutpoints B<-draws[,1:length(coef(polr.out))] Taus<-draws[,(length(coef(polr.out))+1):ncol(draws)] pi.lib.wd<- plogis(Taus[,2] - B%*%t(X.wd)) - plogis(Taus[,1] - B%*%t(X.wd)) pi.lib.ind <- plogis(Taus[,2] - B%*%t(X.ind)) - plogis(Taus[,1] - B%*%t(X.ind)) pi.mod.wd<- plogis(Taus[,4] - B%*%t(X.wd)) - plogis(Taus[,3] - B%*%t(X.wd)) pi.mod.ind <- plogis(Taus[,4] - B%*%t(X.ind)) - plogis(Taus[,3] - B%*%t(X.ind)) fd.lib<- pi.lib.ind - pi.lib.wd fd.mod<- pi.mod.ind - pi.mod.wd #pdf(file = "ANESinterp.pdf") plot(density(fd.mod, adjust=1.5), xlim=c(-0.2,0.2),ylim=c(0,50), xlab="Change in predicted probability", bty="n", col=1, yaxt="n", lwd=2, main="", ylab="") lines(density(fd.lib, adjust=1.5), col=grey(0.5), lwd=2, lty=2) text(x=0.06, y=40, labels="Pr(’liberal’ | PID=’independent’) - \n Pr(’liberal’ | PID=’weak dem’)",cex=.8) text(x=-.08, y=33, labels="Pr(’moderate’ | PID=’independent’) - \n Pr(’moderate’ | PID=’weak dem’)",cex=.8) #dev.off() ###################################################### ### Table 8.3 ######################################## krain<-read.dta("isq05.dta",convert.dates = TRUE, convert.factors = TRUE, missing.type = FALSE, convert.underscore=TRUE, warn.missing.labels=TRUE) cbind(colnames(krain[c(5,10,4,18,19,17,7,11,8,12)]),c("Cold War Dummy variable (0=coldwar, 1=post-cold war", "Krain's ethnic fractionalization score","Duration of geno/politicide (1=1st year of g/p, 2=2nd yr of g/p,)", "Presence of contiguous intervener in previous year, dummy variable", "Neutral or anti-perpetrator intervention dummy variable, lagged", "Previous year's magnitude of severity of geno/politicide (0-5)","Magnitude or severity of geno/politicide (0-5)", "Polity IV regime type variable (democ - autocratic, -10 to 10)", "state failure dummy variable (1=state failure, 0=not", "Degree of economic marginalization within the global economy")) ###################################################### ### Table 8.4 ######################################## ols.fit <- glm(magnitud ~ intrvlag+icntglag+maglag+genyr+stfl+regtype+ethkrain+marg+coldwar,data=krain) polr.fit<-polr(as.ordered(magnitud) ~ intrvlag+icntglag+maglag+genyr+stfl+regtype+ethkrain+marg+coldwar, data=krain, method="logistic", Hess=T) stargazer(ols.fit, polr.fit, covariate.labels=c("Intervention", "Contiguity", "Genocide severity", "Genocide duration", "State failures", "Regime type", "Ethnic fractionalization", "Economic marginalization", "Cold War"), star.cutoffs=NA, digits = 2, no.space=T, ord.intercepts=T) # now let's make factors of the factor variables krain$magnitud<-as.ordered(krain$magnitud) krain$icntglag<-as.factor(krain$icntglag) krain$intrvlag<-as.factor(krain$intrvlag) krain$maglag<-as.factor(krain$maglag) krain$stfl<-as.factor(krain$stfl) ###################################################### ### Table 5.5 ######################################## z.out<-polr(magnitud ~ intrvlag+icntglag+maglag+genyr+stfl+regtype+ethkrain+marg+coldwar, data=krain, method="logistic", Hess=T) summary(z.out) ###################################################### ### Table 8.5 ######################################## cuts <- seq(0,5,.5) cuts.table <- data.frame(matrix(ncol=3,nrow=10)) for(i in 1:nrow(cuts.table)){ cuts.table[i,1] <- paste("Yi",cuts[i],sep=",") cuts.table[i,2] <- paste(cuts[1:i],collapse=",") cuts.table[i,3] <- paste(cuts[(i+1):11],collapse=",") } colnames(cuts.table) <- c("Y","Y = 1","Y = 0") cuts.table ###################################################### ### Figure 8.4 ####################################### krain$magnitud<-as.numeric(krain$magnitud) #makes looping easier below #pdf("fig84.pdf") par(mfrow=c(4,3),las=1) plot.xmean.ordinaly(magnitud ~ intrvlag + icntglag + maglag + genyr + stfl + regtype + ethkrain + marg + coldwar, data=krain,pch=19) #dev.off() ###################################################### ### Figure 8.5 ####################################### krain$stfl <- as.numeric(krain$stfl) krain$maglag <- as.numeric(krain$maglag) krain$icntglag <- as.numeric(krain$icntglag) krain$intrvlag <- as.numeric(krain$intrvlag) Y <- krain$magnitud-1 ncut <- length(unique(Y)) - 1 k <- 9 # total no. of coefficients less intercept Coef <- matrix(NA, ncol=k, nrow=ncut, dimnames=list(paste(">=", levels(as.factor(krain$magnitud))[-1],sep=""), NULL)) SEs<-matrix(NA, ncol=k, nrow=ncut, dimnames=list(paste(">=", levels(as.factor(krain$magnitud))[-1],sep=""), NULL)) for(m in 1:ncut){ f <- lrm(Y >= m ~ intrvlag + icntglag + maglag + genyr + stfl + regtype + ethkrain + marg + coldwar, data=krain) Coef[m,] <- coef(f)[-1] SEs[m,]<-sqrt(diag(f$var)[-1]) } colnames(Coef) <- names(coef(f))[-1] round(Coef, 3) #pdf("fig8.5.pdf") par(mfrow=c(4,3),las=1,pch=19) for(i in 1:k){ plot(1:ncut, Coef[,i], type="b", lty=2, xlab="magnitud", ylab=colnames(Coef)[i], ylim=c( (min(Coef[,i])-2*mean(SEs[,i])) , (max(Coef[,i])+2*mean(SEs[,i])))) abline(h=0, lty=1, lwd=.75, col=grey(.5)) } #dev.off() ###################################################### ### Figure 8.6 ####################################### dta<-read.csv("politymortality.csv") using<-subset(dta, YEAR==2000) names(using) attach(using) #pdf("infantmort.pdf") par(bty="l",las=1) plot(as.factor(XCONST),log(MORTALIT), xlab="POLITY IV executive constraints", ylab="log child mortality", bty="n") #dev.off() ###################################################### ### Table 8.6 ######################################## fit.1<-glm(log(MORTALIT) ~ XCONST) fit.2<-glm(log(MORTALIT) ~ as.factor(XCONST)) fit.3<-glm(log(MORTALIT)~ C(as.factor(XCONST),contr.sdif)) stargazer(fit.1, fit.2, fit.3, star.cutoffs=NA, digits = 2, no.space=T) ######################################################