pvalue.glm <- function(modelo=fit.model,fi=NA) {

#
# Descrição e detalhes:
# A saída é idêntica à parte dos coeficientes obtida no summary(). As diferenças são:
# - será plotado os níveis descritivos marginais dos testes de se cada coeficiente é igual a zero (teste de Wald);
# - será informado o parâmetro de forma, que é o inverso do parâmetro de dispersão;
# - se a distribuição for a gama, o parâmetro de forma é obtido por máxima verossimilhança (biblioteca MASS)
#   juntamente com seu erro padrão e é feito um teste Wald de se ele é igual a 1 (distribuição exponencial).
#   Os erros padrões dos outros coeficientes também dependerão da estimativa do parâmetro de forma, por isso
#   não serão iguais aos obtidos  pela saída do summary() do S-Plus,
# - para a distribuição binomial negativa é considerado parâmetro de dispersão igual a 1, mas é incluído na
#   saída o valor de theta e seu erro padrão, mas não é feito nenhum teste.
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Nota: O comando extractAIC da biblioteca MASS retorna valores diferentes entre o S-Plus e o R do AIC e do BIC,
#       mas aparentemente, a diferença entre o AIC e o BIC é sempre a mesma.
#
# Argumentos obrigatórios:
# modelo: deve-se informar o objeto onde está o ajuste de um modelo linear generalizado, caso não seja informado,
# 	  a função procurará o ajuste no objeto fit.model;
#
# Argumentos opcionais:
# fi: permite que se especifica o parâmetro de forma para corrigir os erros padrões e estatísticas. Por ex., a família
#     Aranda-Ordaz é ajustada por quase-verossimilhança apesar da função de variância ser da Binomial. O ajuste de
#     quase-verossimilhança fornece a estimativa do parâmetro de dispersão (=1/forma), mas se quiser fazer a análise
#     para a Binomial, então deve-se tomar o parâmetro de forma igual a 1. O padrão da função é utilizar o parâmetro
#     que fica no ajuste do modelo.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Exemplo:
# pvalue.glm(ajuste)
#

if(class(modelo)[1] != "glm" & class(modelo)[1] != "negbin") {
	stop(paste("\nA classe do objeto deveria ser glm ou negbin !!!\n"))
}

library("MASS")
if(modelo$family[[1]] == "Gamma") {
	shape<-as.numeric(gamma.shape(modelo)$alpha) #Função gamma.shape retorna phi do texto, gamma.shape$alpha=1/gamma.dispersion
	shapest<-gamma.shape(modelo)$SE
} else {
	shape<-as.numeric(1/summary(modelo)$dispersion)
	shapest<-0
}
if(is.na(fi)==F) {
	shape<-fi
	shapest<-0
}

coef<-c(modelo$coef)
sterr<-sqrt(diag(summary(modelo)$cov.unscaled/shape))
zvalue<-coef/sterr
pvalue<-round(2*pnorm(-abs(zvalue)),digits=6)
coef<-c(coef,shape=shape)

if(modelo$family[[1]] == "Gamma") {
	sterr<-c(sterr,shape=shapest)
	zvalue<-c(zvalue,shape=(shape-1)/shapest)
	pvalue<-c(pvalue,round(2*pnorm(-abs(zvalue["shape"])),digits=4))
} else {
	sterr<-c(sterr,shape=0)
	zvalue<-c(zvalue,shape=0)
	pvalue<-c(pvalue,shape=0)
}
if(class(modelo)[1] == "negbin") {
	coef<-c(coef,theta=modelo$theta)
	sterr<-c(sterr,theta=modelo$SE.theta)
	zvalue<-c(zvalue,theta=0)
	pvalue<-c(pvalue,0)
}

n<-length(modelo$y)
saida<-as.data.frame(list("Value"=coef,"Std. Error"=sterr,"z value"=zvalue,"Pr(>|z|)"=pvalue))
print(saida)
cat("\n")
cat("    Null Deviance:",round(summary(modelo)$null.deviance,digits=6),"on",sum(summary(modelo)$df)-1,"degrees of freedom\n")
cat("Residual Deviance:",round(summary(modelo)$deviance,digits=6),"on",summary(modelo)$df[2],"degrees of freedom\n\n")
cat("AIC:",round(extractAIC(modelo,scale=shape,k=2)[2],digits=6),"\n")
cat("BIC:",round(extractAIC(modelo,scale=shape,k=log(n))[2],digits=6),"\n")
}