###################################### SETUP ########################################################
committee.dat <- read.table("http://jgill.wustl.edu/data/committee.dat",header=TRUE,row.names=1)
attach(committee.dat)
library(MASS)

###################################### MODEL ########################################################
committee.out <- glm.nb(BILLS104 ~ SIZE + SUBS ttee.out <- glm.nb(BILLS104 ~ SIZE + 
	SUBS * (log(STAFF)) + PRESTIGE + BILLS103)
summary.glm(committee.out)
diag.plot(committee.out)
sort(abs( eigen(-glm.vc(committee.out))$values^(-1) ))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.05953  -1.02881  -0.07153   0.53517   1.64052  
 
Coefficients:
                Estimate  Std. Error t value Pr(>|t|)    
(Intercept)     -6.805426  2.546511  -2.672  0.019175 *  
SIZE            -0.028246  0.020926  -1.350  0.200115    
SUBS             1.301595  0.543699   2.394  0.032450 *  
log(STAFF)       3.009711  0.794503   3.788  0.002258 ** 
PRESTIGE        -0.323665  0.441024  -0.734  0.476042    
BILLS103         0.006558  0.001392   4.709  0.000408 ***
SUBS.log(STAFF) -0.323644  0.124888  -2.591  0.022365 *  
---
Signif. codes:  0  `***'  0.001  `**'  0.01  `*'  0.05  `.'  0.1  ` '  1 
 
(Dispersion parameter for Negative Binomial(4.2918) family taken to be 1.494362)
 
    Null deviance: 107.314  on 19  degrees of freedom
Residual deviance:  20.948  on 13  degrees of freedom
AIC: -21130
 
Number of Fisher Scoring iterations: 1
 
committee.se <- sqrt(diag(committee.covmat))
round( cbind(
        committee.out$coefficients, committee.se,
        committee.out$coefficients - committee.se*qt(0.975,committee.out$df.residual),
        committee.out$coefficients + committee.se*qt(0.975,committee.out$df.residual)
),7)

(Intercept)     -6.8054258 2.5465114 -12.3068292 -1.3040224
SIZE            -0.0282463 0.0209262  -0.0734546  0.0169621
SUBS             1.3015945 0.5436985   0.1270053  2.4761838
log(STAFF)       3.0097113 0.7945025   1.2932929  4.7261296
PRESTIGE        -0.3236648 0.4410239  -1.2764391  0.6291095
BILLS103         0.0065575 0.0013924   0.0035494  0.0095656
SUBS.log(STAFF) -0.3236441 0.1248882  -0.5934487 -0.0538394

################################# RESIDUALS GRAPH  ##################################################

z.matrix <- matrix(0,200,200)
for(i in 1:200)  {
        for(j in 1:200)  {
                if(j < 70)    z.matrix[i,j] <- 1
                if(j < 40)    z.matrix[i,j] <- 2
                if(j < 10)    z.matrix[i,j] <- 3
                if(j == 1)    z.matrix[i,j] <- 3.001
                if(j > 130)   z.matrix[i,j] <- 1
                if(j > 160)   z.matrix[i,j] <- 2
                if(j > 190)   z.matrix[i,j] <- 3
                if(j == 200)  z.matrix[i,j] <- 3.001
        }
}

resp <- resid(committee.out,type="response")
pears <- resid(committee.out,type="pearson")
working <- resid(committee.out,type="working")
devs <- resid(committee.out,type="deviance")

r1 <- committee.out$y/max(committee.out$y)
r2 <- committee.out$fitted.values/max(committee.out$fitted.values)
phi.r1 <- cs.string[ceiling(r1*1000)+1,2]
phi.r2 <- cs.string[ceiling(r2*1000)+1,2]
anscombe <- (phi.r1^(2/3) - phi.r2^(2/3))/( (r2^(1/6)) * ((1-r2)^(1/6)) )
anscombe[4] <- devs[4]

par(mfrow=c(4,1),mar=c(0,1,0,1))
plot(pears[order(BILLS104)],type="l",ylim=c(-2,2),xaxt="n")
plot(working[order(BILLS104)],type="l",ylim=c(-2,2),xaxt="n")
plot(devs[order(BILLS104)],type="l",ylim=c(-2,2),xaxt="n")
plot(anscombe[order(BILLS104)],type="l",ylim=c(-2,2),xaxt="n")

postscript("/export/home/jgill/Book.GLM/glm.fig5.ps")
par(mfrow=c(1,1),oma=c(3,3,1,1))
image(seq(0,50,length=200),seq(-2000,2000,length=200),z.matrix,xlim=c(0,50),ylim=c(-2000,2000),
	xaxt="n",yaxt="n",xlab="",ylab="")
points(seq(1,50,length=20),(2000/3)*pears[order(BILLS104)],pch=15)
lines(seq(1,50,length=20),(2000/3)*devs[order(BILLS104)],type="h")
mtext(line=3,side=1,cex=2,"Order of Fitted Outcome Variable")
mtext(line=2,side=2,cex=2,"Residual Effect")
abline(0,0)
dev.off()


########################### COVARIANCE ANALYSIS #####################################################

p1 <- 1:committee.out$rank
committee.covmat <- chol2inv(committee.out$qr$qr[p1, p1, drop = FALSE])*1.494362
committee.se <- sqrt(diag(committee.covmat))[-1]
committee.coef <- committee.out$coef
committee.mat <- cbind(BILLS104, SIZE, SUBS, log(STAFF), PRESTIGE, BILLS103, SUBS*log(STAFF))

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