boxcox.ig <- function(modelo=fit.model,var=stop("faltou informar a variavel"),tetas=seq(-2,2,0.1),link="1/mu^2",maxit=20) { # # Descrição e detalhes: # Esta função ajusta o modelo com distribuição normal inversa (ligação log, inversa ou identidade) testando qual o # teta da transformação Box-Cox [(var^teta)/teta, teta!=0 e log(var), teta==0] diminui mais o desvio escalonado. # # A sugestão desse procedimento é feita por McCullagh e Nelder (pág.402, teoria e pág.412, exemplo). # # Por motivos óbvios deve-se analisar se a variável possui valores negativos ou iguais a zero ao se escolhar os tetas. # # O S-Plus só permite ajustar a normal inversa com ligação canônica (1/mu^2). O R permite também com ligações # inverse e log. No entanto, o R tem um bug, pois armazena sempre no objeto ajustado que utilizou a ligação canônica. # Assim, se utilizar uma ligação que não seja a canônica no R, deve-se informar para a função através da opção link. # # Os dados devem estar disponíveis pelo comando attach( ). # # Argumentos obrigatórios: # modelo: deve-se informar o objeto onde está o ajuste do modelo com distribuição normal inversa, sendo que este ajuste # não deve conter a variável que se está testando qual o melhor teta da transformação de Box-Cox; # var: é a variável que será avaliada qual o melhor teta da transformação de Box-Cox. # # Argumentos opcionais: # tetas: permite determinar quais tetas se deseja avaliar o desvio, o padrão é seq(-2,2,0.1); # link: padrão "1/mu^2", que é a ligação canônica (1/mu^2). Se utilizar a ligação log, informe "log" e se utilizar # a ligação inversa, informe "inverse"; # maxit: máximo de iterações permitido para o processo iterativo convergir, o padrão é 20. # # Autor: Frederico Zanqueta Poleto , arquivo disponível em http://www.poleto.com # # Referência: # MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London. # # Exemplo: # boxcox.ig(ajuste,seq(-0.1,0.6,0.01)) # if(class(modelo)[1] != "glm") { stop(paste("\nA classe do objeto deveria ser glm e nao ",class(modelo),"!!!\n")) } if(modelo$family[[1]] != "Inverse Gaussian" & modelo$family[[1]] != "inverse.gaussian") { stop(paste("\nA familia do objeto deveria ser da normal inversa !!!\n")) } if(link!="1/mu^2") { if(is.null(version$language) == T) { stop(paste("\nO S-Plus só aceita a ligação canônica!!!\n")) } else { if(link!="inverse" & link!="log") { stop(paste("\nO R só aceita as ligações canônica (1/mu^2), inversa (1/mu) e log!!!\n")) } } } X <- model.matrix(modelo) y<-modelo$y gl<-nrow(X)-ncol(X)+1 n<-length(tetas) devi<-numeric(n) for(i in 1:n) { if (tetas[i]==0) { h<-log(var) } else { h<-(var^tetas[i]-1)/tetas[i] } X2<-cbind(X,h) if (link == "1/mu^2") { fit <- glm(y ~ X2-1,family=inverse.gaussian,maxit=maxit) } else { if (link == "inverse") { fit <- glm(y ~ X2-1,family=inverse.gaussian(link=inverse),maxit=maxit) } if (link == "log") { fit <- glm(y ~ X2-1,family=inverse.gaussian(link=log),maxit=maxit) } } fi <- gl/sum((resid(modelo,type="response")/(predict(fit,type="response")*sqrt(y)))^2) devi[i]<-deviance(fit)*fi } plot(tetas,devi,xlab="teta", ylab="desvio escalonado",type="l",pch=16,lwd=1) cat("\nDo conjunto\n") print(tetas) cat(", o que minimizou o desvio foi:\n") tetas[order(devi)[1]] }