# 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(stargazer) library("spduration");library("foreign") library("magrittr");library("ggplot2") library("xtable");library("scales") library(survival);library(eha) library(texreg) #source("extractspdur.R") #hacked for TexReg. sort of works ### Table 11.3 ####################################### ## some deets omitted sdat<-read.dta("coalum.dta") sdat[1:10,c(1,2,sample(3:33,5,replace=F),35)] # a random sample of columns, with duration, belgium, and censor cox.fit<-coxph(Surv(time=duration, event=censor12) ~ majority_government + investiture + volatility + polarization + fractionalization + crisis + formation_attempts + opposition_party, data=sdat) texreg((cox.fit)) stargazer(cox.fit,digits=2,star.cutoffs=NA) ###################################################### ### Table 11.4 ####################################### ## need to put in table yourself, here is the data cox.zph(cox.fit) ###################################################### ### Figure 11.2 ###################################### # pdf(file = "fig9.1.pdf",width = 6, height = 6, onefile = TRUE, family = "Times") par(mfrow=c(3,3),pch=19,las=1) plot(cox.zph(cox.fit),pch=19) #dev.off() ###################################################### ### Figure 11.3 ###################################### deets<-coxph.detail(cox.fit) times<-c(0,deets$time) #observed failure times lambda0<- c(0,deets$hazard) #hazard evaluated a sample mean of covariates x0<-c(0)-cox.fit$mean lambda0<-lambda0*exp(t(coef(cox.fit))%*%x0) #baseline hazard x.inv<- c(1,1,median(sdat$volatility), #investiture scenario mean(sdat$polarization), mean(sdat$fractionalization), mean(sdat$crisis), median(sdat$formation_attempts),mean(sdat$opposition_party)) x.ninv<-x.inv x.ninv[2]<-0 #no investiture scenario lambda.inv<-lambda0*exp(t(coef(cox.fit))%*%x.inv) lambda.ninv<-lambda0*exp(t(coef(cox.fit))%*%x.ninv) Sinv<- exp(-cumsum(lambda.inv)) #cumulative survival function Sninv<- exp(-cumsum(lambda.ninv)) #cumulative survival function #pdf("investiture.pdf") plot(times,Sinv, type="l", xlab="months", lwd=2, ylab=expression(hat(S)(t)), bty="n", las=1) lines(times,Sninv, lwd=2, lty=2) text(x=20, y=c(0.2, 0.7), labels=c("investiture", "no investiture")) #dev.off() ###################################################### ### Table 11.5 ####################################### data(bscoup) bscoup$coup <- ifelse(bscoup$coup=="yes", 1, 0) #table(bscoup$coup) bscoup <- add_duration(bscoup, "coup", unitID="countryid", tID="year", freq="year", ongoing = FALSE) weib_model <- spdur( duration ~ milreg + instab + regconf, atrisk ~ couprisk + wealth + milreg + rwar + regconf + samerica + camerica, data=bscoup, silent = TRUE) extract.spdur(weib_model) cph_model<-coxph(Surv(time=t.0, time2=duration, event=coup) ~ milreg + instab + regconf,data=bscoup[rownames(bscoup)%in%rownames(weib_model$mf.dur),]) texreg(cph_model) loglog_model <- spdur( duration ~ milreg + instab + regconf, atrisk ~ couprisk + wealth + milreg + rwar + regconf + samerica + camerica, data=bscoup, distr="loglog", silent = TRUE) extract.spdur(loglog_model) weib_aft <- aftreg(Surv(time=t.0, time2=duration, event=coup) ~ milreg + instab + regconf, dist="weibull", param = "lifeExp", data=bscoup[rownames(bscoup)%in%rownames(weib_model$mf.dur),]) -2*weib_aft$loglik[2]+4*log(4250) #BIC -2*weib_aft$loglik[2]+4*2 #AIC # Note the extraction functions texreg and stargazer don't support these models. # Note Generic functions don't plot the spduration output. # this has to be pieced together, at present. ###################################################### ### Figure 11.5 ##################################### #pdf("hazard-rates.pdf", height=4, width=8) par(mfrow=c(1, 2)) plot(weib_model, type="hazard", main="Weibull",bty="l",las=1) plot(loglog_model, type="hazard", main="Loglog",bty="l",las=1) #dev.off() ##################################################### ### Figure 11.6 ###################################### #pdf("hazard-ex.pdf", height=4, width=8) par(mfrow = c(1, 2)) plot_hazard(loglog_model, main = "A",bty="n",las=1,bty="l") plot_hazard(loglog_model, xvals = c(1, 1, 10, 0.05), zvals = c(1, 7, 8.64, 1, 1, 0.05, 0, 0), main = "B",bty="l",las=1) #dev.off() ###################################################### ### Figure 11.7 ###################################### # Out-of-sample testing rm(list=ls()) data(bscoup) bscoup$coup <- ifelse(bscoup$coup=="yes", 1, 0) coup_train <- bscoup[bscoup$year < 1996, ] coup_train <- add_duration(coup_train, "coup", unitID="countryid", tID="year", freq="year", ongoing = FALSE) coup_test <- add_duration(bscoup, "coup", unitID="countryid", tID="year", freq="year", ongoing = FALSE) coup_test <- coup_test[coup_test$year >= 1996, ] weib_model2 <- spdur( duration ~ milreg + instab + regconf, atrisk ~ couprisk + wealth + milreg + rwar + regconf + samerica + camerica, data = coup_train) loglog_model2 <- spdur( duration ~ milreg + instab + regconf, atrisk ~ couprisk + wealth + milreg + rwar + regconf + samerica + camerica, data = coup_train, distr="loglog") weib2_test_p <- predict(weib_model2, newdata = coup_test) loglog2_test_p <- predict(loglog_model2, newdata = coup_test) # Out-of-sample separation plots # ______________________________ ### Figure 11.7 ###################################### obs_y <- coup_test[complete.cases(coup_test), "coup"] #pdf("graphics/oos-sepplots.pdf", height=4, width=10) par(mfrow=c(2,1),mar=c(2,2,2,2)) separationplot(weib2_test_p[complete.cases(coup_test)], obs_y, newplot = FALSE, lwd1 = 1, lwd2 = 1, heading = "Weibull") separationplot(loglog2_test_p[complete.cases(coup_test)], obs_y, newplot = FALSE, lwd1 = 1, lwd2 = 1, heading = "Loglog") #dev.off() ######################################################