gee.bino <- function(modelo=fit.model,trace=T,maxit=20,epsilon=0.001) {

#
# Descrição e detalhes:
# Esta programa ajusta o modelo logístico de quase-verossimilhança para dados correlacionados supondo correlacao rr
# intra-classe (ou intra-grupo). A teoria por trás desse ajuste pode ser vista em Paula (2003, págs.158-162) e em
# McCullagh e Nelder (1989, págs.124-128).
#
# Para obter um ajuste com mais opções utilize algum pacote que ajusta GEE.
# 
# Atenção: 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 (sob independência).
#
# Argumentos opcionais:
# trace: como padrão a função irá mostrar as estimativas de cada iteração, se trace=F ela não mostrará;
# maxit: máximo de iterações permitido para o processo iterativo convergir, o padrão é 20;
# epsilon: é a diferença absoluta avaliada entre os coeficientes de duas iterações consecutivas. Se a diferença for menor
#	   ou igual a epsilon para todos os coeficientes, então considera-se que o processo convergiu, padrão=0.001.
#
# A função retorna o ajuste do GEE logístico (fit) e a estimativa da correlação intra-classe (rho).
# Veja também a função corr.bino que testa se rho=0.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# MCCULLAGH, P. e NELDER, J. A. (1989). Generalized Linear Models. 2ª ed. Chapman and Hall, London.
# 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:
# gee.bino(ajuste)
#

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") {
	stop(paste("\nA familia do objeto deveria ser Binomial e nao ",modelo$family[[1]],"!!!\n"))
}
if(modelo$family["link"] != "Logit: log(mu/(1 - mu))" & modelo$family["link"] != "logit") {
	stop(paste("\nA ligacao deve ser a logito !!!\n"))
}

X <- model.matrix(modelo)
if(is.null(modelo$prior.weights)) {
	ntot<-rep(1,n)
} else {
	ntot<-modelo$prior.weights
}
yt <- round(modelo$y*ntot,0)
fr <- predict(modelo,type="response")
rho <- corr.bino(modelo)$rho
i <- 1
converged <- F
while(i <= maxit & converged == F) {
	if (is.null(version$language) == T) {
		#No S-Plus, a opção start é para entrar com o preditor linear
		pm<-predict(modelo)
	} else {
		#No R, a opção start é para entrar com os coeficientes
		pm<-coef(modelo)
	}
	antigo<-coef(modelo)
	modelo <- glm(cbind(yt,ntot-yt) ~ X - 1, family=binomial,start=pm,maxit=1,weights=1/(1+(ntot-1)*rho))
	fr <- fitted(modelo)
	rho <- corr.bino(modelo)$rho
	if (trace == T) {
		cat(paste("\nIter",i,"\n"))
		print(modelo$coef)
	}
	if (any(abs(modelo$coef-antigo)>epsilon)) {
		i <- i+1
	} else {
		converged <- T
	}
}
cat("\n")
if (trace == T) {
	cat(paste("\nConverged",converged,"\n"))
}
list(fit=modelo,rho=rho,converged=converged,iter=i)
}