#####################################################################################################
# Death Penalty Example, "Generalized Linear Models: A Unified Approach", Jeff Gill.
#####################################################################################################

###################################### SETUP ########################################################
dp.97 <- read.table("http://jgill.wustl.edu/data/dp.data",header=TRUE)

pairs(dp.97)
attach(dp.97)

###################################### MODEL ########################################################
dp.out <- glm(EXECUTIONS ~ INCOME+PERPOVERTY+PERBLACK+log(VC100k96)+SOUTH+PROPDEGREE, family=poisson)
summary(dp.out)

1-pchisq(18.212,10) 	

round(glm.vc(dp.out),4)

glm.summary <- function (in.object, alpha = 0.05) 
{
        lo <- in.object$coefficient - qnorm(1-alpha/2) * sqrt(diag(summary(in.object)$cov.unscaled))
        hi <- in.object$coefficient + qnorm(1-alpha/2) * sqrt(diag(summary(in.object)$cov.unscaled))
        out.mat <- round(cbind(in.object$coefficient, sqrt(diag(glm.vc(in.object))), lo, hi),5)
	dimnames(out.mat)[[2]] <- c("Coefficient","Std. Error",
		paste(1-alpha,"CI Lower"),paste(1-alpha,"CI Upper"))
	out.mat
}

glm.summary(dp.out)

range(abs( eigen(-glm.vc(dp.out))$values^(-1) ))

####################### FIRST DIFFERENCE CALCULATION ################################################
dp.mean.vector.0 <- c(1,apply(dp.97[,2:4],2,mean),mean(log(dp.97[,5])),0,mean(dp.97[,7]))
dp.mean.vector.1 <- c(1,apply(dp.97[,2:4],2,mean),mean(log(dp.97[,5])),1,mean(dp.97[,7]))

E.dp0 <- exp(dp.out$coefficients%*%dp.mean.vector.0) 
E.dp1 <- exp(dp.out$coefficients%*%dp.mean.vector.1) 
E.dp1-E.dp0 

########################### RESIDUALS ANALYSIS ######################################################
resp <- resid(dp.out,type="response")
pears <- resid(dp.out,type="pearson")
working <- resid(dp.out,type="working")
devs <- resid(dp.out,type="deviance")
dp.mat <- as.matrix( cbind(rep(1,nrow(dp.97)), dp.97[,-1]) )
anscombe <- (3/2)*(EXECUTIONS^(2/3)-(exp(dp.mat%*%dp.out$coef))^(2/3))/EXECUTIONS^(1/6)
dp.resid.mat <- cbind(resp,pears,working,devs,anscombe)
dimnames(dp.resid.mat)[[2]] <- c("response","pearson","working","deviance","anscombe")
dimnames(dp.resid.mat)[[1]] <- as.character(dp.97[,1])

##################### LEVELS PLOT BASED ON SOUTHERN #################################################
dp.mat.0 <- cbind(dp.mat[,1:5],rep(0,length=nrow(dp.mat)),dp.mat[,7])
dimnames(dp.mat.0)[[2]] <- names(dp.out$coefficients)
dp.mat.1 <- cbind(dp.mat[,1:5],rep(1,length=nrow(dp.mat)),dp.mat[,7])
dimnames(dp.mat.1)[[2]] <- names(dp.out$coefficients)

par(mfrow=c(3,2),mar=c(6,8,2,2),oma=c(3,1,1,1))
for (i in 2:(ncol(dp.mat.0)-1))  {
	if (i==6) i <- i+1
	ruler <- seq(min(dp.mat.0[,i]),max(dp.mat.0[,i]),length=1000)
	xbeta0 <- exp(dp.out$coefficients[-i]%*%apply(dp.mat.0[,-i],2,mean) + dp.out$coefficients[i]*ruler)
	xbeta1 <- exp(dp.out$coefficients[-i]%*%apply(dp.mat.1[,-i],2,mean) + dp.out$coefficients[i]*ruler)

	plot(ruler,xbeta0,type="l", xlab="",ylab="", ylim=c(min(xbeta0,xbeta1)-2,max(xbeta0,xbeta1)) )
	lines(ruler,xbeta1,type="b")
	mtext(outer=F,side=1,paste("Levels of",dimnames(dp.mat.0)[[2]][i]),cex=1.3,line=5)
	mtext(outer=F,side=2,"Expected Executions",cex=1.3,line=4) 
}
plot(ruler[100:200],rep(ruler[400],101),bty="n",xaxt="n",yaxt="n",xlab="",ylab="",type="l",
	xlim=range(ruler),ylim=range(ruler))
lines(ruler[100:200],rep(ruler[600],101),type="b")
text(ruler[445],ruler[400],"Non-South State",cex=2.8)
text(ruler[390],ruler[600],"South State",cex=2.8)
 
##################### JACKNIFE OUT CASES FOR INDEX PLOT #############################################
coef.vector <- NULL
for (i in 1:length(EXECUTIONS))  {
	dp.temp <- glm(EXECUTIONS[-i] ~ INCOME[-i]+PERPOVERTY[-i]+PERBLACK[-i]+log(VC100k96)[-i]+
		SOUTH[-i]+PROPDEGREE[-i], family=poisson)
	coef.vector <- rbind(coef.vector,dp.temp$coefficients)
}

par(mfrow=c(3,2),mar=c(6,8,2,2),oma=c(3,1,1,1))
for(i in 2:ncol(coef.vector))  {
	plot(coef.vector[,i],type="b",xlab="",ylab="")
	abline(dp.out$coefficients[i],0)
	mtext(outer=F,side=1,"Index Number",cex=1.3,line=4) 
	mtext(outer=F,side=2,dimnames(dp.mat.0)[[2]][i],cex=1.3,line=5)
}

#####################################################################################################