wald<-function(modelo=fit.model,C=NULL,b0=NULL) {

#
# Descrição e detalhes:
# Esta função calcula o valor da estatística Wald para um teste com hipóteses lineares do tipo H0: CB=b0, em que C é a
# matriz que determina o teste e b0 um vetor opcional de valores (padrão é um vetor de 0).
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Argumentos obrigatórios:
# modelo: ajuste do modelo sob H1 (irrestrito);
# C: a matrix deve ser de dimensão k por p, considerando p o número de parâmetros do modelo e k o número de igualdades
#    a serem testadas;
#
# Argumentos opcionais:
# b0: opcionalmente pode-se informar um vetor em que se quer testar se CB é igual. Caso nada seja informado, a função
#     automaticamente estara testando CB=0, ou seja, será utilizado b0 como um vetor de k zeros.
#
# A saída será o valor da estatística de Wald (qui-quadrado) juntamente com os graus de liberdade e o nível descritivo.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referência:
# PAULA, G. A. (2003). Modelos de Regressão com apoio computacional. IME-USP, São Paulo. [Não publicado,
#    disponível em http://www.ime.usp.br/~giapaula/Book.pdf]
#
# Exemplos:
# wald(modelo,matrix(c(0,1,0,-1, 0,0,1,-1),nrow=2,ncol=4,byrow=T))
# wald(modelo,cbind(matrix(0,nrow=4,ncol=8),diag(4)))
#

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

if(is.null(C)) {
	stop(paste("\nPor favor, informe a matriz C (se tiver duvidas veja explicacoes dentro da funcao) !!!\n"))
}

if(is.matrix(C)==F) {
	C<-t(as.matrix(C))
}

k<-nrow(C)
p<-ncol(C)

if(p!=length(modelo$coef)) {
	stop(paste("\nO numero de colunas da matriz C deve ser igual ao numero de parametros obrigatoriamente !!!\n"))
}

if(is.null(b0)) {
	b0<-numeric(k)
} else {
	if(length(b0)!=k) {
		stop(paste("\nTamanho de b0 deve ser o mesmo do numero de linhas de C !!!\n"))
	}
}

if(class(modelo)[1]=="lm") {
	fi<-1/(summary(modelo)$sigma^2)
} else {
	if(modelo$family[[1]] == "Gamma") {
		library("MASS")
		fi<-as.numeric(gamma.shape(modelo)$alpha) #Função gamma.shape retorna phi do texto, gamma.shape$alpha=1/gamma.dispersion
	} else {
		fi<-as.numeric(1/summary(modelo)$dispersion)
	}
}

bb<-C%*%modelo$coef-b0
varbb<-C%*%(summary(modelo)$cov.unscaled/fi)%*%t(C)

x2<-as.numeric(t(bb)%*%solve(varbb)%*%bb)

pvalue<-1-pchisq(x2,k)

list(x2=x2,df=k,pvalue=pvalue)
}