aranda <- function(alpha = 1) { # # Descrição e detalhes: # Função para utilizar no S-Plus a ligação de Aranda-Ordaz no modelo binomial da função glm. # # 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. # # Se necessitar extrair o valor do alpha do ajuste utilize as.numeric(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 opcionais: # alpha: deve-se informar alpha ligação, caso não seja especificado a função irá assumir alpha=1. # # 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: # glm(y~x,family=aranda(0.5)) # if(alpha == 0) { binomial(link = cloglog) } else { aranda.lnk <- list( names = alpha, link = substitute( function(mu, al = .Alpha) { log(((1-mu)^(-al)-1)/al) }, list(.Alpha = alpha) ), inverse = substitute( function(eta, al = .Alpha) { 1-(al*care.exp(eta)+1)^(-1/al) }, list(.Alpha = alpha) ), deriv = substitute( function(mu, al = .Alpha) { (al*(1-mu)^(-al-1))/((1-mu)^(-al)-1) }, list(.Alpha = alpha) ), initialize = expression( { if(is.matrix(y)) { if(dim(y)[2]>2) { stop("Only binomial response matrices (2 columns)") } n<-drop(y%*%c(1,1)) y<-y[,1] } else { if(is.category(y)) { y<-y!=levels(y)[1] } else { y<-as.vector(y) } n<-rep(1,length(y)) } w<-w*n n[n==0]<-1 y<-y/n mu<-y+(0.5-y)/n } ) ) aranda.var <- list(names = "mu*(1-mu)", variance = substitute( function(mu) { mu*(1-mu) } ), deviance = substitute( function(mu, y, w, residuals = F) { devy<-y nz<-y!=0 devy[nz]<-y[nz]*log(y[nz]) nz<-(1-y)!=0 devy[nz]<-devy[nz]+(1-y[nz])*log(1-y[nz]) devmu<-y*log(mu)+(1-y)*log(1-mu) if(any(small <- mu*(1-mu) < .Machine$double.eps)) { warning("fitted values close to 0 or 1") smu<-mu[small] sy<-y[small] smu<-ifelse(smu < .Machine$double.eps, .Machine$double.eps,smu) onemsmu<-ifelse((1-smu) < .Machine$double.eps, .Machine$double.eps,1-smu) devmu[small]<-sy*log(smu)+(1-sy)*log(onemsmu) } devi<-2*(devy-devmu) if(residuals) { sign(y-mu)*sqrt(abs(devi)*w) } else { sum(w*devi) } } ) ) make.family("Aranda",link=aranda.lnk,variance=aranda.var) } }