aranda.bino <- function(modelo=fit.model,alphas=seq(0,1.5,0.01),maxit=20) { # # Descrição e detalhes: # Este programa ajusta o modelo com distribuição Binomial e ligação Aranda-Ordaz, estimando o alpha para o menor # desvio obtido. # # Essa ligação é eta=log{[(1-mu)^(-alpha)-1]/alpha}. # Casos particulares: se alpha=1, a ligação é a logito e se alpha=0, a ligação é a complementar log-log. # # Para se ter uma idéia das curvas estimadas da ligação de Aranda-Ordaz com alphas 0 (cloglog), 0.5, 1 (logito), 2 e 4 # e também compará-las com a ligação probito, use os comandos abaixo: #eta<-seq(-3,3,0.01) #w1<-pnorm(eta) #probito #w2<-1-(4*exp(eta)+1)^(-1/4) #aranda(4) #w3<-1-(2*exp(eta)+1)^(-1/2) #aranda(2) #w4<-exp(eta)/(1+exp(eta)) #logito=aranda(1) #w5<-1-(0.5*exp(eta)+1)^(-1/0.5) #aranda(0.5) #w6<-1-exp(-exp(eta)) #cloglog=aranda(0) #plot(cbind(eta,eta,eta,eta,eta,eta),cbind(w1,w2,w3,w4,w5,w6),type="n",xlab="eta",ylab="mu") #lines(eta,w1,col=1,lty=1,lwd=3) #lines(eta,w2,col=2,lty=1,lwd=3) #lines(eta,w3,col=3,lty=1,lwd=3) #lines(eta,w4,col=4,lty=1,lwd=3) #lines(eta,w5,col=5,lty=1,lwd=3) #lines(eta,w6,col=6,lty=1,lwd=3) #abline(0.5,0,lty=2) #legend(x=-3,y=1,legend=c("probito","aranda(4)","aranda(2)","logito","aranda(0.5)","cloglog"),col=c(1,2,3,4,5,6),lty=c(1,1,1,1,1,1),lwd=c(3,3,3,3,3,3)) # # Se necessitar extrair o valor do alpha do ajuste utilize ajuste$family[[2]] # # Nota: o R e o S-Plus ajustam essa ligação como um modelo de quase-verossimilhança, isto é, eles estimam o parâmetro # de dispersão. Assim, se quiser trabalhar com a distribuição Binomial, é necessário considerar o dispersão=1. # Para isso, utilize summary(ajuste,dispersion=1) ou a função pvalue.glm(ajuste,fi=1) # # 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 Binomial (qualquer ligação), caso # não seja informado, a função procurará o ajuste no objeto fit.model; # # Argumentos opcionais: # alphas: permite determinar quais alphas se deseja avaliar o desvio, o padrão é seq(0,1.5,0.01); # 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: # 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: # aranda.bino(ajuste,seq(-0.2,1.2,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]] != "Binomial" & modelo$family[[1]] != "binomial" & modelo$family[[1]] != "Aranda") { stop(paste("\nA familia do objeto deveria ser Binomial ou Aranda-Ordaz e nao ",modelo$family[[1]],"!!!\n")) } X <- model.matrix(modelo) n <- nrow(X) if(is.null(modelo$prior.weights)) { ntot<-rep(1,n) } else { ntot<-modelo$prior.weights } yt <- round(modelo$y*ntot,0) na<-length(alphas) devi<-numeric(na) for(i in 1:na) { modelo<-glm(cbind(yt,ntot-yt) ~ X - 1,family=aranda(alphas[i]),maxit=maxit) devi[i]<-deviance(modelo) } plot(alphas,devi,xlab="alpha", ylab="desvio",type="l",pch=16,lwd=1) cat("\nDo conjunto\n") print(alphas) cat(", o que minimizou o desvio foi:\n") alphas[order(devi)[1]] }