aranda <- function(alpha = 1) {

#
# Descrição e detalhes:
# Função para utilizar no R 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 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 <fred@poleto.com>, 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 {
	care.exp <- function(x, thresh = about36) {
		about36 <-  - log(.Machine$double.eps)
		thresh <- min(thresh, about36)
		x[x > thresh] <- thresh
		x[x < ( - thresh)] <-  - thresh
		exp(x)
	}
	linkfun <- function(mu) { log(((1-mu)^(-alpha)-1)/alpha) }
	linkinv <- function(eta) { 1-(alpha*care.exp(eta)+1)^(-1/alpha) }
	mu.eta <- function(eta) { care.exp(eta)*(alpha*care.exp(eta)+1)^(-1/alpha-1) }
        valideta <- function(eta) TRUE
	variance <- function(mu) mu * (1 - mu)
	validmu <- function(mu) all(mu > 0) && all(mu < 1)
	dev.resids <- function(y, mu, wt) {
			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)
			wt*devi
		}
	aic <- function(y, n, mu, wt, dev) {
		m <- if (any(n > 1)) 
			n
		else wt
		-2 * sum(ifelse(m > 0, (wt/m), 0) * dbinom(round(m * y), round(m), mu, log = TRUE))
	}
	initialize = expression( {
		if (NCOL(y) == 1) {
			if (is.factor(y)) y <- y != levels(y)[1]
			n <- rep(1, nobs)
			if (any(y < 0 | y > 1)) stop("y values must be 0 <= y <= 1")
			mustart <- (weights * y + 0.5)/(weights + 1)
			m <- weights * y
			if (any(abs(m - round(m)) > 0.001)) warning("non-integer #successes in a binomial glm!")
		} else if (NCOL(y) == 2) {
			if (any(abs(y - round(y)) > 0.001)) warning("non-integer counts in a binomial glm!")
			n <- y[, 1] + y[, 2]
			y <- ifelse(n == 0, 0, y[, 1]/n)
			weights <- weights * n
			mustart <- (n * y + 0.5)/(n + 1)
		} else stop(paste("For the binomial family, y must be", 
			"a vector of 0 and 1's or a 2 column", "matrix where col 1 is no. successes", 
			"and col 2 is no. failures"))
	} )
	structure(list(family = "Aranda", link = alpha, linkfun = linkfun, linkinv = linkinv,
		variance = variance, dev.resids = dev.resids, aic = aic, mu.eta = mu.eta,
		initialize = initialize, validmu = validmu, valideta = valideta), class = "family")
}

}