rv.ig <- function(fit0,fit1){

#
# Descrição e detalhes:
# Esta função calcula o valor da estatística da razão de verossimilhanças para testar dois modelos encaixados com
# distribuição normal inversa.
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Argumentos obrigatórios:
# fit0: ajuste do modelo sob H0;
# fit1: ajuste do modelo sob H1.
#
# A saída será o valor da estatística da razão de verossimilhanças 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]
#
# Exemplo:
# rv.ig(ajuste.sobH0,ajuste.sobH1)
#

if( (class(fit0)[1] != "glm") | (class(fit1)[1] != "glm") ) {
	stop(paste("\nAs classes dos dois objetos deveriam ser glm !!!\n"))
}
if(fit0$family[[1]] != "Inverse Gaussian" & fit0$family[[1]] != "inverse.gaussian") {
	stop(paste("\nA familia do objeto deveria ser da normal inversa !!!\n"))
}
if(fit1$family[[1]] != "Inverse Gaussian" & fit1$family[[1]] != "inverse.gaussian") {
	stop(paste("\nA familia do objeto deveria ser da normal inversa !!!\n"))
}

y<-fit0$y

if ( sum(y != fit1$y) != 0 ) {
	stop(paste("\nOs dois modelos nao possuem a mesma variavel resposta!!!\n"))
}

m0 <- fitted(fit0)
m1 <- fitted(fit1)

f0 <- 1/summary(fit0)$dispersion
f1 <- 1/summary(fit1)$dispersion

n <- nrow(model.matrix(fit0))

rv <- 2*sum( f1*((1/m1)-(y/(2*(m1^2))))-f0*((1/m0)-(y/(2*(m0^2)))) + ((log(2*pi*(y^3)/f0)+(f0/y))/2)-((log(2*pi*(y^3)/f1)+(f1/y))/2) )

df <- summary(fit0)$df[2]-summary(fit1)$df[2]

list(rv=rv,df=df,pvalue=1-pchisq(rv,df))
}