#source("E:\\windows\\Unicamp\\Disciplinas\\1_semestre_2014\\ME 714\\Programas\\diag_pois.txt")

#---------------------------------------------------------------#
# Para rodar este programa  deixe no objeto   fit.model a saída 
# do  ajuste  da  regressão com erros Poisson.  Deixe  os dados 
# disponíveis  através do comando attach(...). Depois use o 
# comando source(...) no S-Plus ou R para  executar o programa. 
# A sequência de comandos é a seguinte:
#
#        > fit.model <- ajuste
#        > attach(dados)
#        > source("diag_pois")
#
# A saída terá quatro gráficos: de pontos de alacanca, de pontos
# influentes  e  dois de resíduos. Para identificar os pontos
# que  mais  se destacam usar o comando identify(...). Se por
# exemplo se destacam três pontos no plot(fitted(fit.model),h,...), 
# após esse comando coloque
#     
#        > identify(fitted(fit.model),h,n=3)
#
# O mesmo pode ser feito nos demais gráficos. Nos gráficos de 
# resíduos foram colocados os limites ylim=c(a-1,b+1), onde a
# é o menor valor e b o maior valor para o resíduo.Para voltar
# a  ter apenas um gráfico por tela faça o seguinte:
#
#        > par(mfrow=c(1,1))
# 
#--------------------------------------------------------------#
#
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
w <- fit.model$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
ts <- resid(fit.model,type="pearson")/sqrt(1-h)
td <- resid(fit.model,type="deviance")/sqrt(1-h)
di <- (h/(1-h))*(ts^2)
par(mfrow=c(2,2))
a <- min(td)
b <- max(td)
plot(td,xlab="Indice", ylab="Residuo Componente do Desvio",
ylim=c(a-1,b+1), pch=16)
abline(2,0,lty=2)
abline(-2,0,lty=2)
#identify(td, n=1)
#title(sub="(c)")
#
fited = fitted(fit.model)
plot(fited ,td,xlab="valor ajustado (média)", ylab="Residuo Componente do Desvio",ylim=c(a-1,b+1), pch=16)
abline(2,0,lty=2)
abline(-2,0,lty=2)
eta <- predict(fit.model)
z <- eta + resid(fit.model, type="pearson")/sqrt(w)
plot(predict(fit.model),z,xlab="Preditor Linear", 
ylab="Variavel z", pch=16)
lines(smooth.spline(predict(fit.model), z, df=2))
#title(sub="(d)")
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
w <- fit.model$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- resid(fit.model,type="deviance")/sqrt((1-h))
e <- matrix(0,n,100)
#
for(i in 1:100){
nresp <- rpois(n, fitted(fit.model))
fit <- glm(nresp ~ X, family=poisson)
w <- fit$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
#
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:n){
  eo <- sort(e[i,])
e1[i] <- (eo[2]+eo[3])/2
e2[i] <- (eo[97]+eo[98])/2}
#
med <- apply(e,1,mean)
faixa <- range(td,e1,e2)
#par(pty="s")
qqnorm(td,xlab="Percentil da N(0,1)",
ylab="Componente do Desvio", ylim=faixa, pch=16)
par(new=T)
#
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2)
#---------------------------------------------------------------#                      

#---------------------------------------------------------------#