A análise de dados binários é usualmente feita através da regressão logística, mas esse modelo possui limitações. Modificar a função de ligação da regressão permite maior flexibilidade na modelagem e diversas propostas já foram feitas nessa área.
Fenômenos que assumem apenas dois estados são chamados de binários e seu estudo sempre foi de grande interesse. Seja prever o resultado de uma eleição ou estudar a prevalência do uso de drogas ou remédios em uma certa população, é comum a presença de dados dicotômicos e o problema de estimar a proporção de ocorrência de um dos eventos faz-se importante. Além disso, pode ser interessante avaliar quanto algumas variáveis influenciam na proporção desse evento. No entanto, técnicas de regressão linear não podem ser utilizadas devido a violação de algumas suposições.
Diversas soluções para esse problema foram propostas e a regressão logística, um caso particular dos modelos lineares generalizados (Nelder e Wedderburn, 1972), é uma das mais utilizadas. Nesses, uma função de ligação é aplicada na combinação linear dos preditores de modo que a resposta seja um valor apropriado para o problema. Para dados binários, é natural utilizar funções quantil (inversa da função de distribuição) como é o caso da regressão logística que utiliza a distribuição logística como o nome sugere. Outras distribuições como a normal padrão (probito) e valor extremo (complementar log-log) também são opções comuns para a função de ligação (Agresti, 2002).
Diversas outras funções de ligações foram propostas na literatura para estender as mencionadas. Prentice (1975) apresentou uma função de ligação baseada no logaritmo da distribuição F-Snedecor. Aranda-Ordaz (1981) propôs uma função assimétrica que tem os modelos logístico e complementar log-log como casos particulares. Stukel (1988) introduziu uma generalização do modelo logístico adicionando dois parâmetros à função de ligação logística que controlam o comportamento das caudas da distribuição. Caron et al. (2009) e Caron (2010) sugerem um modelo baseado na distribuição Weibull que possui o modelo complementar log-log como caso especial e aproxima os modelos logístico e probito.
Apesar das qualidades destes modelos, até onde se sabe, não existe função para a estimação desses nos pacotes estatísticos a não ser pelos três primeiros casos mencionados. Todas as outras funções compartilham a qualidade de serem paramétricas, ou seja, há elementos extras a serem estimados além dos coeficientes de regressão. Prentice (1976) aponta dificuldades em estimar seu modelo com os dois parâmetros livres utilizando máxima verossimilhança e sugere fixar um deles e calcular o outro.
Banco de dados sobre mortalidade de besouros (Bliss, 1935). O objetivo original desse estudo era obter um inseticida eficaz contra besouros. Para isso, 481 animais foram expostos a diferentes concentrações de dissulfeto de carbono (CS2) durante cinco horas e contou-se o número de insetos mortos.
Esse conjunto de dados é conhecido por não ser bem ajustado por modelos simétricos, em particular logito e probito; por conta disso, é amplamente citado em trabalhos que buscam alternativas a esses modelos.
library(VGAM)
data("flourbeetle")
show(flourbeetle)
## logdose CS2mgL exposed killed
## 1 1.690728 49.06 59 6
## 2 1.724194 52.99 60 13
## 3 1.755189 56.91 62 18
## 4 1.784189 60.84 56 28
## 5 1.811307 64.76 63 52
## 6 1.836894 68.69 59 53
## 7 1.860996 72.61 62 61
## 8 1.883888 76.54 60 60
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
grid()
library(rsq)
m <- lapply(c("logit", "probit", "cloglog", "cauchit"), function(type)
glm(cbind(killed,exposed-killed) ~ logdose, data = flourbeetle, family = binomial(link = type)))
names(m) <- c("logit", "probit", "cloglog", "cauchit")
m
## $logit
##
## Call: glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type),
## data = flourbeetle)
##
## Coefficients:
## (Intercept) logdose
## -60.72 34.27
##
## Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
## Null Deviance: 284.2
## Residual Deviance: 11.22 AIC: 41.42
##
## $probit
##
## Call: glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type),
## data = flourbeetle)
##
## Coefficients:
## (Intercept) logdose
## -34.94 19.73
##
## Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
## Null Deviance: 284.2
## Residual Deviance: 10.11 AIC: 40.3
##
## $cloglog
##
## Call: glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type),
## data = flourbeetle)
##
## Coefficients:
## (Intercept) logdose
## -39.57 22.04
##
## Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
## Null Deviance: 284.2
## Residual Deviance: 3.44 AIC: 33.64
##
## $cauchit
##
## Call: glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = type),
## data = flourbeetle)
##
## Coefficients:
## (Intercept) logdose
## -77.32 43.52
##
## Degrees of Freedom: 7 Total (i.e. Null); 6 Residual
## Null Deviance: 284.2
## Residual Deviance: 20.15 AIC: 50.35
Observemos os resultados de cada ajuste.
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = 2)
lines(fitted(m[[2]]) ~ logdose, data = flourbeetle, col = 3)
lines(fitted(m[[3]]) ~ logdose, data = flourbeetle, col = 4)
lines(fitted(m[[4]]) ~ logdose, data = flourbeetle, col = 5)
grid()
Existem muitos casos em que obtemos uma especificação incorreta da função de ligação. O motivo é simples: temos que escolher a função do ligação antes de termos informações suficientes sobre a escolha da ligação. Assim, gostaríamos de descrever uma maneira de melhorar a qualidade do ajuste dos GLMs, reduzindo o desvio. Descobrirá que uma maneira elegante de melhorar os modelos é permitir funções de ligação provenientes das famílias de ligação paramétricas especificadas em Czado (2007). Esta classe paramétrica vantajosa de transformações de ligações foi desenvolvida por Czado (1992). As funções gerais de transformação de potência \(h(\cdot)\) são os elementos-chave para modificar as caudas de um gráfico.
A seguir, vamos apresentar as funções de ligação transformada de potência \(h(\cdot)\) (Czado, 2007) para valores específicos de \(\psi\). Portanto, temos que passar à função um valor inicial \(\eta_0\in\mathbb{R}\) e \(\psi_1\) ou \(\psi_2\in \mathbb{R}\) para uma modificação de cauda única ou os valores do vetor \(\psi\in\mathbb{R}^2\) para uma modificação de ambas as caudas.
Para uma modificação da cauda esquerda, todo ponto \(<\eta_0\) será modificado, enquanto para uma modificação da cauda direita todo ponto é \(\geq\eta_0\). Uma modificação de ambas as caudas modifica ambas as caudas, ou seja, todos os pontos \(<\eta_0\) e \(\geq \eta_0\) são modificados.
Ao definirmos os parâmetros como 1, não obtemos nenhuma modificação, ou seja, uma linha reta. Se definirmos um parâmetro como 1 na modificação de ambas as caudas, obtemos uma única modificação da cauda, por exemplo, \(\psi= (1,\psi_2)\) modifica a cauda esquerda. Numa modificação de cauda direita, o parâmetro \(\psi_1 <1\) diminui a inclinação, enquanto a configuração \(\psi_1> 1\) a aumenta. Na modificação da cauda esquerda ocorre o contrário para \(\psi_2\).
Definição ((Modificação da cauda direita): \[\begin{equation*} h_{\eta_0}(\eta,\psi=\psi_1)=\left\{ \begin{array}{cc} \eta_0+\ln\big(\eta-\eta_0+1\big) & \mbox{ se } \eta\geq \eta_0 \mbox{ e } \psi_1=0 \\[1.2em] \eta_0+\dfrac{(\eta-\eta_0+1)^{\psi_1}-1}{\psi_1} & \mbox{ se } \eta\geq \eta_0 \mbox{ e } \psi_1\neq 0\\[1.2em] \eta & \mbox{ caso contrário, isto é, } \eta<\eta_0\cdot \end{array}\right. \end{equation*}\]
Definição ((Modificação da cauda direita): \[\begin{equation*} h_{\eta_0}(\eta,\psi=\psi_2)=\left\{ \begin{array}{cc} \eta & \mbox{ se } \eta\geq \eta_0\\[1.2em] \eta_0-\ln\big(-\eta+\eta_0+1\big) & \mbox{ se } \eta< \eta_0 \mbox{ e } \psi_2=0 \\[1.2em] \eta_0-\dfrac{(-\eta+\eta_0+1)^{\psi_2}-1}{\psi_2} & \mbox{ caso contrário } \eta< \eta_0 \mbox{ e } \psi_2\neq 0\cdot \end{array}\right. \end{equation*}\]
Definição (Modificação em ambas as caudas): \[\begin{equation*} h_{\eta_0}(\eta,\psi=(\psi_1,\psi_2))=\left\{ \begin{array}{cc} \eta_0+\ln\big(\eta-\eta_0+1\big) & \mbox{ se } \eta\geq \eta_0 \mbox{ e } \psi_1=0 \\[1.2em] \eta_0+\dfrac{(\eta-\eta_0+1)^{\psi_1}-1}{\psi_1} & \mbox{ se } \eta\geq \eta_0 \mbox{ e } \psi_1\neq 0\\[1.2em] \eta_0-\ln\big(-\eta+\eta_0+1\big) & \mbox{ se } \eta< \eta_0 \mbox{ e } \psi_2=0 \\[1.2em] \eta_0-\dfrac{(-\eta+\eta_0+1)^{\psi_2}-1}{\psi_2} & \mbox{ caso contrário } \eta<\eta_0 \mbox{ e } \psi_2\neq 0\cdot \end{array}\right. \end{equation*}\]
psi1GAUSS<-function(psi1 = 1, eta0 = 0)
{linkfun <- function(mu) {hpsi1INV(psi1, mu, eta0)}
linkinv <- function(eta){hpsi1(psi1, eta, eta0)}
mu.eta <- function(eta) {hpsi1DERIV(psi1, eta, eta0)}
valideta<-function(eta) {h <- 1:length(eta)
for (i in 1:length(eta) ) { if (is.finite(linkinv(eta[i]))) {h[i] <- TRUE}
else {h[i] <- FALSE}
}
h
}
link <- paste("psi1GAUSS(", psi1, " , " , eta0, ")")
structure(list(linkfun = linkfun,
linkinv = linkinv,
mu.eta = mu.eta,
valideta = valideta,
name = link),
class = "link-glm")
}
psi2LOGIT<-function(psi2 = 1, eta0 = 0)
{linkfun <- function(mu) hpsi2INV(psi2, log(mu/(1 - mu)), eta0)
linkinv <- function(eta) 1 - (1/(1 + exp(hpsi2(psi2, eta, eta0))))
mu.eta <- function(eta){
(hpsi2DERIV(psi2, eta, eta0)*(1 - (1/(1 + exp(hpsi2(psi2, eta, eta0))))))/
(1 + exp(hpsi2(psi2, eta, eta0)))
}
valideta<-function(eta){h <- 1:length(eta)
for (i in 1:length(eta) ) {
if (is.finite(linkinv(eta[i]))) {h[i] <- TRUE}
else {h[i] <- FALSE}
}
h
}
link <- paste("psi2LOGIT(", psi2, " , " , eta0, ")")
structure(list(linkfun = linkfun,
linkinv = linkinv,
mu.eta = mu.eta,
valideta = valideta,
name = link),
class = "link-glm")
}
hpsi1 <- function(psi1 = stop("Argument 'psi1' is missing"),
eta = stop("Argument 'eta' is missing"), eta0 = 0)
{h <- 1:length(eta)
h[eta < eta0] <- eta[eta < eta0]
if (any(psi1 > -1e-14 && psi1 < 1e-14)) {
h[eta >= eta0] <- eta0 + log(eta[eta >= eta0] - eta0 + 1)
}
else {
h[eta >= eta0] <- ((eta[eta >= eta0] - eta0 + 1)^psi1 - 1)/psi1
h[eta >= eta0] <- h[eta >= eta0] + eta0
}
h
}
hpsi1INV<-function(psi1, y, eta0 = 0)
{h <- 1:length(y)
if (any(psi1 > -1e-14 && psi1 < 1e-14)){
h[y >= eta0] <- eta0 - 1 + exp(y[y >= eta0] - eta0)
}
else {
h[y >= eta0] <- eta0 - 1 + ((1 + (psi1 * (-eta0 + y[y >= eta0])))^(1/psi1))
}
h[y < eta0] <- y[y < eta0]
h
}
hpsi1DERIV<-function(psi1, eta, eta0 = 0)
{h <- 1:length(eta)
if (any(psi1 > -1e-14 && psi1 < 1e-14)){
h[eta >= eta0] <- 1/(eta[eta >= eta0] - eta0 + 1)
}
else {
h[eta >= eta0] <- (1 - eta0 + eta[eta >= eta0])^(psi1 - 1)
}
h[eta < eta0] <- 1
h
}
hpsi1DERIV1<-function(psi1, eta, eta0 = 0)
{h <- eta
temp <- eta
if (any(psi1 > -1e-14 && psi1 < 1e-14)){
h[eta >= eta0] <- ((log(eta[eta >= eta0] - eta0 + 1))^2)/2
}
else {
temp[eta >= eta0] <- (1 - eta0 + eta[eta >= eta0])^psi1
h[eta >= eta0] <- ((temp[eta >= eta0] *
log(eta[eta >= eta0] - eta0 + 1) * psi1) - (temp[eta >= eta0] - 1)) /(psi1^2)
}
h[eta < eta0] <- 0
h
}
hpsi2 <- function (psi2 = stop("Argument 'psi2' is missing"),
eta = stop("Argument 'eta' is missing"), eta0 = 0)
{h <- 1:length(eta)
h[eta >= eta0] <- eta[eta >= eta0]
if (any(psi2 > -1e-14 && psi2 < 1e-14)) {
h[eta < eta0] <- eta0 - log(-(eta[eta < eta0]) + eta0 + 1)
}
else {
h[eta < eta0] <- -((- eta[eta < eta0] + eta0 + 1)^psi2 - 1)/psi2
h[eta < eta0] <- h[eta < eta0] + eta0
}
h
}
hpsi2INV<-function(psi2, y, eta0 = 0)
{h <- 1:length(y)
h[y >= eta0] <- y[y >= eta0]
if (any(psi2 > -1e-14 && psi2 < 1e-14)){
h[y < eta0] <- eta0 + 1 - exp(eta0 - y[y < eta0])
}
else {
h[y < eta0] <- eta0 + 1 - ((1 - (psi2 * (y[y < eta0] - eta0)))^(1/psi2))
}
h
}
hpsi2DERIV<-function(psi2, eta, eta0 = 0)
{h <- 1:length(eta)
h[eta >= eta0] <- 1
if (any(psi2 > -1e-14 && psi2 < 1e-14)){
h[eta < eta0] <- 1/(eta0 + 1 - eta[eta < eta0])
}
else{
h[eta < eta0] <- (eta0 + 1 - eta[eta < eta0])^(psi2 - 1)
}
h
}
hpsi2DERIV2<-function(psi2, eta, eta0 = 0)
{h <- eta
temp <- eta
h[eta >= eta0] <- 0
if (any(psi2 > -1e-14 && psi2 < 1e-14)){
h[eta < eta0] <- - ((log(1 + eta0 - eta[eta < eta0]))^2)/2
}
else {
temp[eta < eta0] <- (1 + eta0 - eta[eta < eta0])^psi2
h[eta < eta0] <- - ((temp[eta < eta0] *
log( - eta[eta < eta0] + eta0 + 1) * psi2) - (temp[eta < eta0] - 1))/(psi2^2)
}
h
}
hpsi12 <- function (psi1 = stop("Argument 'psi1' is missing"),
psi2 = stop("Argument 'psi2' is missing"),
eta = stop("Argument 'eta' is missing"), eta0 = 0)
{h <- 1:length(eta)
if (any(psi1 > -1e-14 && psi1 < 1e-14)) { h[eta >= eta0] <- eta0 + log(eta[eta >= eta0] - eta0 + 1)}
else {
h[eta >= eta0] <- ((eta[eta >= eta0] - eta0 + 1)^psi1 - 1)/psi1
h[eta >= eta0] <- h[eta >= eta0] + eta0
}
if (any(psi2 > -1e-14 && psi2 < 1e-14)) { h[eta < eta0] <- eta0 - log(-(eta[eta < eta0]) + eta0 + 1)}
else {
h[eta < eta0] <- -((- eta[eta < eta0] + eta0 + 1)^psi2 - 1)/psi2
h[eta < eta0] <- h[eta < eta0] + eta0
}
h
}
hpsi12INV<-function(psi1, psi2, y, eta0 = 0)
{h <- 1:length(y)
if (any(psi1 > -1e-14 && psi1 < 1e-14)){
h[y >= eta0] <- eta0 - 1 + exp(y[y >= eta0] - eta0)
}
else {
h[y >= eta0] <- eta0 - 1 + ((1 + (psi1 * (y[y >= eta0] - eta0)))^(1/psi1))}
if (any(psi2 > -1e-14 && psi2 < 1e-14)){ h[y < eta0] <- eta0 + 1 - exp(eta0 - y[y < eta0])}
else {h[y < eta0] <- eta0 + 1 - ((1 - (psi2 * (y[y < eta0] - eta0)))^(1/psi2))}
h
}
hpsi12DERIV<-function(psi1, psi2, eta, eta0 = 0)
{h <- 1:length(eta)
if (any(psi1 > -1e-14 && psi1 < 1e-14)){
h[eta >= eta0] <- 1/(eta[eta >= eta0]- eta0 + 1)
}
else {
h[eta >= eta0] <- (1 - eta0 + eta[eta >= eta0])^(psi1 - 1)
}
if (any(psi2 > -1e-14 && psi2 < 1e-14)){ h[eta < eta0] <- 1/(eta0 + 1 - eta[eta < eta0])
}
else{ h[eta < eta0] <- (1 + eta0 - eta[eta < eta0])^(psi2 - 1)
}
h
}
hpsi12DERIV1<-function(psi1, psi2, eta, eta0 = 0)
{h <- eta
temp <- eta
if (any(psi1 > -1e-14 && psi1 < 1e-14)){ h[eta >= eta0] <- ((log(eta[eta >= eta0] - eta0 + 1))^2)/2
}
else {
temp[eta >= eta0] <- (1 - eta0 + eta[eta >= eta0])^psi1
h[eta >= eta0] <- ((temp[eta >= eta0] *
log(eta[eta >= eta0] - eta0 + 1) * psi1) - (temp[eta >= eta0] - 1))/(psi1^2)
}
if (any(psi2 > -1e-14 && psi2 < 1e-14)){ h[eta < eta0] <- 0 }
else { h[eta < eta0] <- 0 }
h
}
hpsi12DERIV2<-function(psi1, psi2, eta, eta0)
{h <- eta
temp <- eta
if (any(psi1 > -1e-14 && psi1 < 1e-14)){ h[eta >= eta0] <- 0 }
else { h[eta >= eta0] <- 0 }
if (any(psi2 > -1e-14 && psi2 < 1e-14)){
h[eta < eta0] <- - ((log(1 + eta0 - eta[eta < eta0]))^2)/2
}
else {
temp[eta < eta0] <- (1 + eta0 - eta[eta < eta0])^psi2
h[eta < eta0] <- - ((temp[eta < eta0] *
log( - eta[eta < eta0] + eta0 + 1) * psi2) - (temp[eta < eta0] - 1))/(psi2^2)
}
h
}
ajuste <- glm(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
family = binomial(link=psi2LOGIT(psi2=0.1,eta0=0)))
summary(ajuste)
##
## Call:
## glm(formula = cbind(killed, exposed - killed) ~ logdose, family = binomial(link = psi2LOGIT(psi2 = 0.1,
## eta0 = 0)), data = flourbeetle)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0253 -0.2851 0.1766 0.5411 0.9081
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -88.363 10.283 -8.593 <2e-16 ***
## logdose 49.546 5.719 8.664 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.2024 on 7 degrees of freedom
## Residual deviance: 3.1136 on 6 degrees of freedom
## AIC: 33.312
##
## Number of Fisher Scoring iterations: 4
rsq(ajuste)
## [1] 0.4942133
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Logito paramétrico"), col = c("black","red"), fill = c(1,2))
grid()
Aproximação à função de ligação paramétrica Stukel (Stukel (1988)).
ajuste1 = update(m[[1]], formula = . ~ . + ifelse(predict(m[[1]])>0, 0.5*predict(m[[1]])^2,0) +
ifelse(predict(m[[1]])<0, -0.5*predict(m[[1]])^2,0), family = binomial(link="logit"))
summary(ajuste1)
##
## Call:
## glm(formula = cbind(killed, exposed - killed) ~ logdose + ifelse(predict(m[[1]]) >
## 0, 0.5 * predict(m[[1]])^2, 0) + ifelse(predict(m[[1]]) <
## 0, -0.5 * predict(m[[1]])^2, 0), family = binomial(link = "logit"),
## data = flourbeetle)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.26801 0.67253 -0.44440 -0.40939 1.02043 -0.85234 0.03192 0.61935
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -53.7809 16.6477
## logdose 30.1860 9.3888
## ifelse(predict(m[[1]]) > 0, 0.5 * predict(m[[1]])^2, 0) 0.3597 0.2666
## ifelse(predict(m[[1]]) < 0, -0.5 * predict(m[[1]])^2, 0) -0.1765 0.2534
## z value Pr(>|z|)
## (Intercept) -3.231 0.00124 **
## logdose 3.215 0.00130 **
## ifelse(predict(m[[1]]) > 0, 0.5 * predict(m[[1]])^2, 0) 1.349 0.17734
## ifelse(predict(m[[1]]) < 0, -0.5 * predict(m[[1]])^2, 0) -0.696 0.48613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.2024 on 7 degrees of freedom
## Residual deviance: 3.0416 on 4 degrees of freedom
## AIC: 37.24
##
## Number of Fisher Scoring iterations: 4
rsq(ajuste1)
## [1] 0.4922435
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste1) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Ligação Stukel"), col = c("black","red"), fill = c(1,2))
grid()