ligsl.vars <- function(modelo=fit.model,lambda=seq(-2,3,0.1),sep=F,maxit=20) {

#
# Descrição e detalhes:
# A função irá obter o gráfico com desvios escalonados obtidos de ajustes de quase-verossimilhança para as funções
# de variância constante, mu, mu^2 e mu^3, e uma ligação potência/Box-Cox (mu^lambda).
# É uma forma de checagem informal contra má especificação das funções de ligação e variância.
# McCullagh e Nelder (1989, pág.400), fazem o mesmo plotando o desvio contra a função de variância na forma de potência.
# Como estamos variando mais a função de ligação do que a de variância, o foco principal é da função de ligação.
# Usaremos a família quasi e ligação power.
#
# Notas: Como as vezes os valores escolhidos para o lambda da função de ligação quando combinados com certas funções de
#	 variância podem travar o S-Plus, após cada lambda será impressa uma mensagem. Deve-se tomar muito cuidado, pois
#	 as vezes mesmo quando os ajustes aparentemente estão convergindo, podem haver instabilidades em que os
#	 coeficientes do modelo são enormes e o desvio também, nesse caso o parâmetro de forma tenta compensar o
#	 desvio alto e fica muito pequeno, o que faz com que o desvio escalonado seja próximo de zero.
#	 O R apresentou mais problemas do que o S-Plus, mas em ambos esses problemas podem ocorrer.
#
# Os dados devem estar disponíveis pelo comando attach( ).
#
# Argumentos obrigatórios:
# modelo: deve-se informar o objeto onde está o ajuste de um modelo, caso não seja informado, a função procurará
# 	  o ajuste no objeto fit.model;
# 
# Argumentos opcionais:
# lambda: permite que se especifique para quais lambdas deve-se ajustar a ligação. O padrão é seq(-2,3,0.1);
# sep: se for F (False) irá desenhar apenas um gráfico com uma linha para cada função de variância e se for T (True)
#      irá desenhar um gráfico para cada função de variância. O padrão é F;
# maxit: essa opção é utilizada nos ajustes e indica o máximo de iterações permitidas nos ajustes. O padrão é maxit=20.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, 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.
#
# Exemplos:
# ligsl.vars(ajuste)
# ligsl.vars(ajuste,lambda=seq(1.5,2.5,0.01),sep=T,maxit=50)
#

if(class(modelo)[1]=="lm" | class(modelo)[1]=="glm") {

} else {
	stop(paste("\nA classe do objeto deveria ser lm ou glm !!!\n"))
}

X <- model.matrix(modelo)
if(class(modelo)[1]=="glm") {
	Y <- modelo$y
} else {
	Y <- predict(modelo)+resid(modelo)
}

devi<-matrix(0,length(lambda),2)
for(i in 1:length(lambda)) {
	x1<-as.data.frame(summary(do.call("glm",list(Y ~ X -1,family=quasi(link=power(lambda[i]),variance="mu"),maxit=maxit)))[c("deviance","dispersion","iter")])
	x2<-as.data.frame(summary(do.call("glm",list(Y ~ X -1,family=quasi(link=power(lambda[i]),variance="mu^2"),maxit=maxit)))[c("deviance","dispersion","iter")])
	cat("Lambda ",lambda[i],"(iteracoes): inversa(",x1[["iter"]],") log(",x2[["iter"]],")\n",sep="")
	devi[i,1]<-x1[["deviance"]]/x1[["dispersion"]]
	devi[i,2]<-x2[["deviance"]]/x2[["dispersion"]]
}

na1<-sum(1*is.na(devi[,1]))
na2<-sum(1*is.na(devi[,2]))
if(na1 > 0) {
	lambda1<-lambda[-which.na(devi[,1])]
	devi1<-devi[-which.na(devi[,1]),1]
	t1<-"(retas interrompidas indicam problemas de convergência)"
} else {
	lambda1<-lambda
	devi1<-devi[,1]
	t1<-""
}
if(na2 > 0) {
	lambda2<-lambda[-which.na(devi[,2])]
	devi2<-devi[-which.na(devi[,2]),2]
	t2<-"(retas interrompidas indicam problemas de convergência)"
} else {
	lambda2<-lambda
	devi2<-devi[,2]
	t2<-""
}
if(sep == T) {
	par(mfrow=c(1,2))
	plot(lambda1,devi1,xlab="Potência na Média da Função de Ligação", ylab="Quase-Desvio Escalonado",main="Função de Variância: média",type="l",pch=16,sub=t1)
	plot(lambda2,devi2,xlab="Potência na Média da Função de Ligação", ylab="Quase-Desvio Escalonado",main="Função de Variância: média ^ 2",type="l",pch=16,sub=t2)
	par(mfrow=c(1,1))
} else {
	if(sum(1*is.na(devi)) > 0) {
		t0<-"(retas interrompidas indicam problemas de convergência)"
	} else {
		t0<-""
	}
	#par(mfrow=c(1,1))
	plot(rbind(lambda1,lambda2),rbind(devi1,devi2),xlab="Potência na Média da Função de Ligação", ylab="Quase-Desvio Escalonado",type="n",pch=16,sub=t0)
	lines(lambda1,devi1,type="o",pch=4,col=3,lty=1,lwd=3)
	lines(lambda2,devi2,type="o",pch=5,col=4,lty=1,lwd=3)
	if(is.null(version$language)) {
		key(x=min(rbind(lambda1,lambda2)),y=max(rbind(devi1,devi2)),text=list(c("média","média ^ 2")),lines=list(type=c("o","o")),lty=c(1,1),pch=c(4,5),col=c(3,4),lwd=c(3,3),divide=2,border=T,title="Função de Variância:",cex.title=1.1)
	} else {
		legend(x=min(rbind(lambda1,lambda2)),y=max(rbind(devi1,devi2)),legend=c("média","média ^ 2"),density=0,col=c(3,4),lty=c(1,1),lwd=c(3,3),pch=c(4,5))
	}
}
cat("\n")
dimnames(devi)<-list(lambda,c("mu","mu^2"))
devi
}