# 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. # # This code is relevant to Chapter 7, but does not produce any tables or figures that # appear therein, because it contains neither figures nor tables. rm(list=ls()); gc() library(foreign); library(pscl) set.seed(123) #### Table 7.1 #### bs<-read.csv("mediation.csv",header=T) bs$pop90<-bs$pop90/100000 # transform to more easily interpret output bs$gdp90<-bs$gdp90/1000 bs<-subset(bs,select=c("medteam","mediat","council","pop90","gdp90","duration")) bs.out<-glm(medteam ~ -1 + gdp90 + council, data=bs,family="poisson") summary(bs.out) #### Table 7.2 #### offset.out <- glm(medteam ~ -1 + gdp90+council + offset(log(duration)), data=bs,family="poisson") summary(offset.out) #### Section 7.2.2 #### bs.out<-glm(medteam ~ -1 + gdp90 + council, data=bs,family="quasipoisson") summary(bs.out) #### Figure 7.1 #### docks <- read.csv("WWFstoppageData102210.csv") #remove NAs for a density plot docks.1 <- docks[which(!is.na(docks$pol.hrs.imp)),] #create summed data docks.summed <- data.frame(matrix(ncol=40,nrow=183)) colnames(docks.summed) <- colnames(docks.1) docks.summed$port <- "summed" docks.summed$yr.qtr <- docks.1$yr.qtr[1:183] docks.summed$year <- docks.1$year[1:183] for(i in 1:nrow(docks.summed)){ docks.summed$pol.hrs.imp[i] <- sum(docks.1$pol.hrs.imp[which(docks.1$yr.qtr == docks.summed$yr.qtr[i] & docks.1$port != "top2")]) } docks.1 <- rbind(docks.1,docks.summed) #plot pdf("fig7.1.pdf", width = 6, height = 6, onefile = TRUE, family = "Times") docks.1$pol.hrs.imp.log <- log(docks.1$pol.hrs.imp + 1) plot(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "summed")],bw=1.377),lwd=2,ylim=c(0,.3), xlim=c(-5,17),axes=F,frame.plot=T,main="") lines(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "syd")],bw=1.377),lty=2) lines(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "mel")],bw=1.377),lty=3) lines(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "kem")],bw=1.377),lty=4,col="blue") lines(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "new")],bw=1.377),col="red",lty=4) lines(density(docks.1$pol.hrs.imp.log[which(docks.1$port == "bri")],bw=1.377),lty=4,col="seagreen") legend('topright',c("Summed","Sydney","Melbourne","Kembla","Newcastle","Brisbane"),col=c("black","black", "black","blue","red","seagreen"),lty=c(1,2,3,4,4,4),lwd=c(2,1,1,1,1,1)) axis(2,at=c(0,.1,.2,.3),labels=c(0.0,0.1,0.2,0.3)) axis(1,at=c(-5,0,5,10,15),labels=c(-5,0,5,10,15)) title("Log quarterly hours lost to political strikes 1951-2001") dev.off() #clean-up rm(docks.1,docks.summed,i) #### Table 7.3 #### #(Still can't recover the estimates in Table 7.3) hurdle.mod <- hurdle(ceiling(pol.hrs.imp) ~ contract.exp.2f + grgdpch + as.factor(accord) + as.factor(labor.govt), dis="negbin",data=docks) zinf.mod <- zeroinfl(ceiling(pol.hrs.imp) ~ contract.exp.2f + grgdpch + as.factor(accord) + as.factor(labor.govt), dis="negbin",data=docks) summary(hurdle.mod) summary(zinf.mod)