Ligações paramétricas

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.

Dados

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()

Modelos

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()

Modelos lineares generalizados com ligação paramétrica

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.

Transformações de potência gerias \(h(\cdot)\)

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()