rv.gama <- 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 gama. # # 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 , 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.gama(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]] != "Gamma") | (fit1$family[[1]] != "Gamma") ) { stop(paste("\nAs familias dos objetos deveriam ser Gamma !!!\n")) } if(fit0$family[[2]] != fit1$family[[2]]) { stop(paste("\nAs funcoes de ligacao dos objetos deveriam ser as mesmas !!!\n")) } y<-fit0$y if ( sum(y != fit1$y) != 0 ) { stop(paste("\nOs dois modelos nao possuem a mesma variavel resposta!!!\n")) } t0 <- sum(log(y/fitted(fit0)) - y/fitted(fit0)) t1 <- sum(log(y/fitted(fit1)) - y/fitted(fit1)) library("MASS") f0 <- 1/gamma.dispersion(fit0) #Função gamma.shape retorna phi do texto, gamma.shape$alpha=1/gamma.dispersion f1 <- 1/gamma.dispersion(fit1) d0 <- f0*log(f0) - log(gamma(f0)) d1 <- f1*log(f1) - log(gamma(f1)) n <- nrow(model.matrix(fit0)) rv <- 2*(f1*t1 - f0*t0 + n*(d1 - d0)) df <- summary(fit0)$df[2]-summary(fit1)$df[2] list(rv=rv,df=df,pvalue=1-pchisq(rv,df)) }