escore<-function(modelo=fit.model,x0=NULL) { # # Descrição e detalhes: # Esta função calcula o valor da estatística Escore para um teste com hipóteses lineares do tipo H0: B=0. # A expressão utilizada pode ser encontrada em Paula (2003, págs.31-32). # # Os dados devem estar disponíveis pelo comando attach( ). # # Argumentos obrigatórios: # modelo: ajuste do modelo sob H0; # x0: matriz das variáveis correspondentes aos parâmetros utilizados na hipótese nula, obrigatoriamente deve ter o mesmo # número de linhas da matriz modelo do ajuste. # # A saída será o valor da estatística Escore (qui-quadrado) 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: # escore(modelo,cbind(x3,x4)) # 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(x0)) { stop(paste("\nPor favor, informe a matriz x0 (se tiver duvidas veja explicacoes dentro da funcao) !!!\n")) } if(is.matrix(x0)==F) { x0<-as.matrix(x0) } n<-nrow(x0) k<-ncol(x0) X1<-x0 X2<-model.matrix(modelo) if(n!=nrow(X2)) { stop(paste("\nO numero de linhas da matriz x0 deve ser igual ao numero de linhas da matriz modelo do ajuste !!!\n")) } if(class(modelo)[1]=="lm") { W<-diag(n) fi<-1/(summary(modelo)$sigma^2) } else { W<-diag(modelo$weights) 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) } } rp<-resid(modelo,type="pearson")*sqrt(fi) C<-solve( t(X2)%*%W%*%X2 ) %*% t(X2)%*%W%*%X1 R<-X1-X2%*%C x2<-solve(t(R)%*%W%*%R) x2<-t(rp)%*%sqrt(W)%*%X1%*%x2%*%t(X1)%*%sqrt(W)%*%rp x2<-as.numeric(x2) pvalue<-1-pchisq(x2,k) list(x2=x2,df=k,pvalue=pvalue) }