A Seção 1.1 discute como estimar e fazer inferências sobre uma única probabilidade de sucesso \(\pi\). A Seção 1.2 generaliza essa discussão para a situação de duas probabilidades de sucesso que agora dependem de um nível de grupo. Agora completamos a generalização para uma situação em que há muitas possibilidades diferentes de sucesso para estimar e realizar inferências. Além disso, quantificamos como uma variável explicativa com muitos níveis possíveis, talvez contínuo em vez de categórico, afeta a probabilidade de sucesso. Essas generalizações são feitas através do uso de modelos de regressão binária, também chamados de modelos de regressão binomial.


2.1 Modelos de regressão linear


Antes de introduzirmos os modelos de regressão binária, precisamos revisar os modelos de regressão linear normal. Seja \(Y_i\) a variável de resposta para observações \(i = 1,\cdots,n\). Além disso, considere a situação com \(p\) variáveis explicativas \(x_{i1},\cdots,x_{ip}\) que são medidos na \(i\)-ésima observação. Relacionamos as variáveis explicativas com a variável resposta através de um modelo linear: \[ Y_i = \beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}+\epsilon_i, \] onde \(\epsilon_i\) para \(i = 1,\cdots,n\) são independentes e cada um tem uma distribuição de probabilidade normal com média 0 e variância \(\sigma^2\).

Isso leva a cada \(Y_i\) para \(i = 1,\cdots,n\) sendo independente com distribuições normais. Os \(\beta_0\),\(\cdots\),\(\beta_p\) são parâmetros de regressão que quantificam as relações entre as variáveis explicativas e \(Y_i\). Por exemplo, se \(\beta_1 = 0\), isso indica que não existe uma relação linear entre a primeira variável explicativa e a variável de resposta dadas as outras variáveis do modelo. Alternativamente, se \(\beta_1 > 0\), existe uma relação positiva, um aumento na variável explicativa leva a um aumento na variável resposta e se \(\beta_1 < 0\), existe uma relação negativa, um aumento na variável explicativa leva a uma diminuição na a variável de resposta.

Tomando a esperança de \(Y_i\), também podemos escrever o modelo como \[ \mbox{E}(Y_i)=\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}, \]

onde assumimos que as variáveis explicativas são condicionadas ou simplesmente constantes. Essas equações do modelo permitem que diferentes valores possíveis de \(Y_i\) ou \(\mbox{E}(Y_i)\) sejam funções de um conjunto de variáveis explicativas. Assim, se tivermos \(x_{i1},\cdots, x_{ip}\) e de alguma forma conhecer todos os valores dos parâmetros \(\beta_0,\cdots,\beta_p\), poderíamos encontrar o que \(Y_i\) seria em média. Claro, não sabemos \(\beta_0,\cdots,\beta_p\) em aplicações reais, mas podemos estimar esses parâmetros usando estimação de mínimos quadrados ou outros procedimentos mais complicados, se desejado.

No contexto de respostas binárias, a quantidade que queremos estimar é a probabilidade de sucesso, \(\pi\). Sejam \(Y_i\) variáveis de resposta binárias independentes para observações \(i = 1,\cdots,n\), onde um valor de 1 denota um sucesso e um valor de 0 denota uma falha. Uma distribuição normal obviamente não é mais um modelo apropriado para \(Y_i\). Em vez disso, semelhante à Seção 1.1, uma distribuição de Bernoulli descreve \(Y_i\) muito bem, mas agora permitimos que o parâmetro de probabilidade de sucesso \(\pi_i\) seja diferente para cada observação. Assim, \[ P(Y_i = y_i) = \pi_i^{y_i} (1-\pi_i)^{1-y_i}, \qquad \mbox{para} \qquad y_i = 0 \quad \mbox{ou} \quad 1 \] é a função de probabilidade (PMF) para \(Y_i\).

Para encontrar o estimador de máxima verossimilhança de \(\pi_i\), a função de verossimilhança é \[ L(\pi_1,\cdots,\pi_n|y_1,\cdots,y_n) = P(Y_1=y_1)\times \cdots \times P(Y_n=y_n) = \prod_{i=1}^n P(Y_i=y_i)=\prod_{i=1}^n \pi_i^{y_i}(1-\pi_i)^{1-y_i}\cdot \]

Observe que existem \(n\) parâmetros diferentes para as \(n\) observações diferentes. Isso significa que o modelo está saturado: há tantos parâmetros quanto observações. As estimativas de parâmetros neste caso são as mesmas que os dados, \(\widehat{\pi}_i= 0\) ou 1, então não há ganho trabalhando com este modelo.

Suponha, em vez disso, que temos variáveis explicativas \(x_{i1},\cdots, x_{ip}\) como anteriormente, e desejamos relacionar os \(\pi_i\)’s a essas variáveis. Podemos propor uma função matemática, digamos \(\pi_i = (x_{i1},\cdots,x_{ip})\), para descrever a relação. Por exemplo, em uma configuração simples, esta função pode designar dois valores possíveis de \(\pi(x_{i1})\) dependendo de um único valor binário para \(x_{i1}\), \(x_{i1}\) poderia denotar o número do grupo para a observação \(i\) como fizemos na Seção 1.2. Em seguida, substituiríamos \(\pi(x_{i1},\cdots,x_{ip})\) na equação acima para estimar seus parâmetros. Isso nos forneceria um modelo que pode ser usado para estimar uma probabilidade de sucesso em função de quaisquer valores possíveis para um conjunto de variáveis explicativas. A próxima seção mostra que uma escolha apropriada da função \(\pi(x_{i1},\cdots, x_{ip})\) leva a uma gama completa de procedimentos úteis e convenientes de estimação e inferência.


2.2 Modelos de regressão logística


O tipo mais simples de função para \(\pi(x_{i1},\cdots, x_{ip})\) é o modelo linear \(\pi_i = \beta_0+ \beta_1 x_{i1}+\cdots + \beta_p x_{ip}\). No entanto, este modelo pode levar a valores de \(\pi_i\) menores que 0 ou maiores que 1, dependendo dos valores das variáveis explicativas e dos parâmetros de regressão. Obviamente, isso é bastante indesejável ao estimar uma probabilidade. Felizmente, estão disponíveis muitas expressões não lineares que forçam \(\pi_i\) a ficar entre 0 e 1. A expressão mais comumente usada é o modelo de regressão logística: \[ \pi_i=\dfrac{\exp\left( \beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)}{1+\exp\left(\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)}\cdot \]

Observe que \(\exp\left( \beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)\) é sempre positivo e o numerador da equação acima é menor que o denominador, o que leva a \(0 < \pi_i < 1\).

Um modelo de regressão logística também pode ser escrito como \[ \log\left(\dfrac{\pi}{1-\pi_i}\right) = \beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}, \]

através de algumas manipulações algébricas. O lado esquerdo da equação acima é o logaritmo natural para as chances de sucesso, que exploraremos ao interpretar o modelo na Seção 2.2.3. Essa transformação de \(\pi_i\) é frequentemente chamada de transformação logit, assim chamada em analogia ao modelo “probit”, discutida na Seção 2.3 e é simplesmente denotada como \(logit(\pi_i)\) ao escrever o lado esquerdo do modelo. O lado direito da equação acima é uma combinação linear dos parâmetros de regressão com as variáveis explicativas, e isso é frequentemente chamado de preditor linear. Estas equações são declarações equivalentes do modelo de regressão logística.

Frequentemente escreveremos os modelos sem os subscritos \(i\) como \[ \pi=\dfrac{\exp\left( \beta_0+\beta_1 x_{1}+\cdots+\beta_p x_{p}\right)}{1+\exp\left(\beta_0+\beta_1 x_{1}+\cdots+\beta_p x_{p}\right)} \] e \[ logit(\pi)= \beta_0+\beta_1 x_{1}+\cdots+\beta_p x_{p}, \] quando não estamos nos referindo a observações particulares. Ficará claro com base no contexto quando a probabilidade de sucesso for uma função de variáveis explicativas e não como \(\pi\) foi usado na Seção 1.1.1, mesmo que o símbolo \(\pi\) não as mencione explicitamente.


Exemplo 2.1: Gráfico do modelo de regressão logística.


O objetivo deste exemplo é examinar a forma do modelo de regressão logística quando há uma única variável explicativa \(x_1\). Considere o modelo \[ \pi = \dfrac{\exp( \beta_0 + \beta_1 x1)}{1 + \exp(\beta_0 + \beta_1 x_1)}, \]

que é equivalentemente expresso como \[ logit(\pi ) = \beta_0 + \beta_1 x_1\cdot \]

Suponha que \(\beta_0 = 1\) e \(\beta_1 = 0.5\). A figura abaixo mostra este modelo à esquerda. O gráfico à direita é o mesmo, mas com \(\beta_1 = -0.5\).

par( mfrow = c(1 ,2))
beta0 <- 1
beta1 <- 0.5
curve ( expr = exp( beta0 + beta1 *x) / (1+ exp( beta0 + beta1 *x)), xlim = c(-15, 15) , 
        col = " black ", main = expression (pi == frac (e ^{1+0.5* x[1]} , 1+e ^{1+0.5* x [1]}) ), 
        xlab = expression (x [1]) , ylab = expression (pi))
grid()
beta1 <- -0.5
curve ( expr = exp( beta0 + beta1 *x) / (1+ exp( beta0 + beta1 *x)), xlim = c(-15, 15) , 
        col = " black ", main = expression (pi == frac (e ^{1-0.5* x[1]} , 1+e ^{1-0.5* x [1]}) ), 
        xlab = expression (x [1]) , ylab = expression (pi))
grid()


Podemos fazer as seguintes generalizações examinando o modelo e esses gráficos:


Usamos duas novas funções aqui. A função par() define parâmetros gráficos que controlam várias opções de plotagem. Em nosso código, nós o usamos para particionar a janela gráfica em 1 linha e 2 colunas usando o argumento mfrow, que significa “make frame by row”. Existem muitos outros usos da função par(), e nós encorajamos os leitores a investigá-los na ajuda da função.

A segunda nova função é a função curve() que é usada para plotar o modelo. Esta é uma função muito útil para plotar funções matemáticas que variam em uma variável. Em nosso exemplo, o argumento expr contém o modelo de regressão logística onde a letra x deve ser usada como o nome da variável plotada no eixo x. Por padrão, a função matemática é avaliada em 101 valores do eixo x igualmente espaçados dentro do intervalo especificado por xlim. Esses 101 pontos resultantes são então unidos por linhas retas. Também dentro de curve(), usamos a função expression() com os rótulos de título e eixo para incluir letras e frações gregas.


2.2.1 Estimação de parâmetros


A estimação por máxima verossimilhança é usada para estimar os parâmetros de regressão \(\beta_0,\cdots,\beta_p\) do modelo de regressão logística. Substituindo nosso modelo por \(\pi_i\) e tomando o logaritmo natural para obter a função de log-verossimilhança: \[ \begin{array}{rcl} \log L(\beta_0,\cdots,\beta_p|y_1,\cdots,y_n) & = & \displaystyle \log\left(\prod_{i=1}^n \pi_i^{y_i}(1-\pi_i)^{1-y_i}\right) \, = \, \sum_{i=1}^n y_i \log(\pi_i)+(1-y_i)\log(1-\pi_i) \\ & = & \displaystyle \sum_{i=1}^n y_i \log\left( \dfrac{\exp\left( \beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)}{1+\exp\left(\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)}\right) \\ & & \qquad \qquad \qquad + (1-y_i)\log\left( \dfrac{1}{1+\exp\left(\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right)}\right) \\ & = & \displaystyle \sum_{i=1}^n y_i \left(\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip}\right) \\ & & \qquad \qquad \qquad - \log\left(1+\exp(\beta_0+\beta_1 x_{i1}+\cdots+\beta_p x_{ip})\right)\cdot \end{array} \]

Tomamos as derivadas da expressão acima em relação a \(\beta_0,\cdots,\beta_p\), iguale-os a 0 e resolva-os simultaneamente para obter as estimativas dos parâmetros \(\widehat{\beta}_0,\cdots,\widehat{\beta}_p\). Quando as estimativas dos parâmetros são incluídas no modelo, podemos obter a probabilidade estimada de sucesso como \[ \widehat{\pi}=\dfrac{\exp\left(\widehat{\beta}_0+\widehat{\beta}_1 x_{i1}+\cdots+\widehat{\beta}_p x_{ip}\right)}{1+\exp\left(\widehat{\beta}_0+\widehat{\beta}_1 x_{i1}+\cdots+\widehat{\beta}_p x_{ip}\right)}\cdot \]

Infelizmente, existem apenas alguns casos simples em que essas estimativas de parâmetros têm soluções de forma fechada; ou seja, geralmente não podemos escrever as estimativas dos parâmetros em termos dos dados observados como poderíamos para a estimativa de probabilidade única \(\widehat{\pi}\) na Seção 1.1.2.

Em vez disso, usamos procedimentos numéricos iterativos para encontrar sucessivamente estimativas dos parâmetros de regressão que aumentam a função de log-verossimilhança. Quando as estimativas mudam de forma insignificante para iterações sucessivas, isso sugere que atingimos o pico da função de log-verossimilhança e dizemos que elas convergiram. Se as estimativas continuarem a mudar visivelmente até um número máximo de iterações selecionado, o procedimento numérico iterativo não convergiu e essas estimativas de parâmetros finais não devem ser usadas. Discutiremos convergência e não convergência com mais detalhes na Seção 2.2.7.

Dentro do R e da maioria dos pacotes de software estatístico, o algoritmo dos mínimos quadrados reponderados iterativamente (IRLS) é o procedimento numérico iterativo usado para encontrar as estimativas dos parâmetros. Este procedimento utiliza o critério dos mínimos quadrados ponderados, que é comumente utilizado para modelos de regressão linear normal quando há variância não constante. O algoritmo IRLS alterna entre atualizar os pesos e as estimativas dos parâmetros de forma iterativa até que a convergência seja alcançada. A função glm(), “glm” significa “modelo linear generalizado”, dentro de R implementa este procedimento de estimativa de parâmetros e mostramos como usar esta função no próximo exemplo. Mais adiante nesta seção, mostramos um método alternativo para maximizar a função de verossimilhança usando a função optim().


Exemplo 2.2: Placekicking.


Conforme discutido no Capítulo 1, os pontos podem ser marcados no futebol americano por um placekicker chutando uma bola através de uma área-alvo no final do campo. Um sucesso ocorre quando a bola é chutada por cima do travessão e entre as duas traves dos postes da baliza.

A equipe do placekicker recebe 1 ou 3 pontos, onde um ponto após touchdown (PAT) recebe 1 ponto e uma cesta de campo recebe 3 pontos. Um placekick que não for bem sucedido recebe 0 pontos. Bilder and Loughin (1998) usam regressão logística para estimar a probabilidade de sucesso de um placekick (PAT ou field goal) na National Football League (NFL).

Eles examinam uma série de variáveis explicativas, incluindo:

Existem 1.425 observações de placekick coletadas durante a temporada de 1995 da NFL que estão disponíveis no conjunto de dados. Abaixo está como os dados são lidos em R:

placekick <- read.csv(file = "http://leg.ufpr.br/~lucambio/CE073/20222S/Placekick.csv")
head ( placekick )
##   week distance change  elap30 PAT type field wind good
## 1    1       21      1 24.7167   0    1     1    0    1
## 2    1       21      0 15.8500   0    1     1    0    1
## 3    1       20      0  0.4500   1    1     1    0    1
## 4    1       28      0 13.5500   0    1     1    0    1
## 5    1       20      0 21.8667   1    0     0    0    1
## 6    1       25      0 17.6833   0    0     0    0    1


Os dados são salvos em um formato delimitado por vírgulas (.csv), então usamos a função read.table() para lê-los. O argumento file especifica onde o arquivo está localizado. O valor do argumento header = TRUE instrui R a usar os nomes na primeira linha do arquivo como os nomes das variáveis. O valor do argumento sep = “,” indica que o arquivo é delimitado por vírgulas.

Depois de ler o conjunto de dados em R, ele é salvo em um quadro de dados chamado placekick e as primeiras seis observações são impressas usando a função head(). Por exemplo, o primeiro chute no conjunto de dados é da primeira semana da temporada e foi tentado em 21 jardas com 24,7 minutos restantes na metade. Além disso, o primeiro chute foi um field goal bem-sucedido tentado sob condições de mudança de liderança, na grama, em um estádio abobadado e sob condições sem vento. Observe que não consideramos o placekicker como uma variável de interesse, a justificativa para esta decisão é dada na p.22 de Bilder and Loughin, 1998.

Por enquanto, usaremos apenas a distância para estimar a probabilidade de um placekick bem-sucedido. Formalmente, \(Y\) é a variável de resposta com um valor de 1 para um sucesso e 0 para uma falha, e a variável explicativa \(x_1\) que denota a distância em jardas para o placekick. Usaremos os dados observados para estimar os parâmetros no modelo \[ logit(\pi) = \beta_0 + \beta_1 x_1; \]

que também pode ser escrito com “distância” substituindo \(x_1\) para ajudar a tornar a declaração do modelo mais descritiva.

mod.fit <- glm( formula = good ~ distance , family = binomial ( link = logit ), data = placekick )
mod.fit
## 
## Call:  glm(formula = good ~ distance, family = binomial(link = logit), 
##     data = placekick)
## 
## Coefficients:
## (Intercept)     distance  
##       5.812       -0.115  
## 
## Degrees of Freedom: 1424 Total (i.e. Null);  1423 Residual
## Null Deviance:       1013 
## Residual Deviance: 775.7     AIC: 779.7


Os resultados de glm() são salvos em um objeto que chamamos de mod.fit, que é uma versão abreviada de “model fit”, podemos escolher um nome diferente. Os argumentos dentro de glm() são:

Ao imprimir o objeto mod.fit executando mod.fit no prompt de comandos, vemos que o modelo de regressão logística estimado é \(logit(\widehat{\pi}) = 5.8121 - 0.1150\)distance. Como há uma estimativa de parâmetro negativa correspondente à distância, a probabilidade estimada de sucesso diminui à medida que a distância aumenta.

Informações adicionais são impressas na saída, mas na verdade há muito mais itens, chamados de componentes, que são armazenados no objeto mod.fit. Através do uso da função names(), obtemos a seguinte lista desses componentes:

names(mod.fit)
##  [1] "coefficients"      "residuals"         "fitted.values"    
##  [4] "effects"           "R"                 "rank"             
##  [7] "qr"                "family"            "linear.predictors"
## [10] "deviance"          "aic"               "null.deviance"    
## [13] "iter"              "weights"           "prior.weights"    
## [16] "df.residual"       "df.null"           "y"                
## [19] "converged"         "boundary"          "model"            
## [22] "call"              "formula"           "terms"            
## [25] "data"              "offset"            "control"          
## [28] "method"            "contrasts"         "xlevels"
mod.fit$coefficients
## (Intercept)    distance 
##   5.8120798  -0.1150267


O objeto mod.fit tem formato de lista. Como vimos anteriormente, as listas permitem que vários objetos sejam vinculados em um lugar comum. Para acessar os componentes de uma lista, usamos o formato \(\mbox{<object>}\)$\(\mbox{<component>}\) sem os \(<\) \(>\).

O exemplo na saída mostra que mod.fit$coefficients contém o valores de \(\widehat{\beta}_0\) e \(\widehat{\beta}_1\).

Os resultados de muitos dos outros componentes são óbvios por seus nomes. Por exemplo, mod.fit$fitted.values são as probabilidades estimadas para cada observação, \(\widehat{\pi}_i\) para \(i = 1,\cdots,n\) e mod.fit$residuals contém os resíduos do ajuste do modelo, \(y_i-\widehat{\pi}_i\) para \(i = 1,\cdots,n\). Outros componentes não são tão óbvios e discutiremos muitos deles à medida que avançarmos neste capítulo. Observe que as descrições de todos os componentes estão disponíveis na ajuda do glm().

Um resumo do que está dentro do mod.fit é obtido da função summary():

summary (mod.fit )
## 
## Call:
## glm(formula = good ~ distance, family = binomial(link = logit), 
##     data = placekick)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7441   0.2425   0.2425   0.3801   1.6092  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.812080   0.326277   17.81   <2e-16 ***
## distance    -0.115027   0.008339  -13.79   <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: 1013.43  on 1424  degrees of freedom
## Residual deviance:  775.75  on 1423  degrees of freedom
## AIC: 779.75
## 
## Number of Fisher Scoring iterations: 6


A saída acima exibe muitas informações sobre o modelo que descreveremos ao longo deste capítulo. Por enquanto, observe que os valores de \(\widehat{\beta}_0\) e \(\widehat{\beta}_1\) são exibidos na tabela “Coeficientes” sob o cabeçalho “Estimate”. Além disso, foram necessárias 5 iterações para obter essas estimativas, conforme fornecido pela última linha da saída; “Fisher scoring” é equivalente a IRLS para modelos de regressão logística.

Informações adicionais sobre como o R funciona são úteis para entender os resultados de summary(). A classe do objeto mod.fit é glm. Associadas a esta classe estão várias funções de método que ajudam a resumir as informações dentro dos objetos ou que completam cálculos úteis adicionais. Para ver uma lista dessas funções de método, a função methods() é usada:

class (mod.fit )
## [1] "glm" "lm"
methods ( class = glm)
##  [1] add1           anova          coerce         confint        cooks.distance
##  [6] deviance       drop1          effects        extractAIC     family        
## [11] formula        influence      initialize     logLik         model.frame   
## [16] nobs           predict        print          residuals      rstandard     
## [21] rstudent       show           slotsFromS3    summary        vcov          
## [26] weights       
## see '?methods' for accessing help and source code


Usaremos muitas dessas funções de método ao longo deste capítulo. Por enquanto, observe a função summary.glm(). Quando a função genérica summary() é executada, R primeiro encontra a classe do objeto e, em seguida, verifica se existe uma função de método correspondente a essa classe. No caso aqui, ele procura por summary.glm(). Como essa função existe, ela é usada para concluir os cálculos principais.

Se mais de uma variável explicativa for incluída no modelo, os nomes das variáveis podem ser separados por símbolos “+” no argumento da fórmula. Por exemplo, suponha que incluímos a variável de change além de distance no modelo:

mod.fit2 <- glm( formula = good ~ change + distance , 
                  family = binomial ( link = logit ), data = placekick )
mod.fit2$coefficients
## (Intercept)      change    distance 
##   5.8931814  -0.4477832  -0.1128888


O modelo de regressão logística estimado é \[ logit(\widehat{\pi}) = 5.8932 -0.4478\, \mbox{change} - 0.1129\,\mbox{distance}\cdot \]


A matriz de variância-covariância estimada, doravante denominada “matriz de covariâncias” para simplificar, para um vetor de estimativas de parâmetros \(\widehat{\beta}=(\widehat{\beta}_0,\widehat{\beta}_1,\cdots,\widehat{\beta}_p)^\top\), \(\widehat{\mbox{Var}}(\widehat{\beta})\), tem a forma usual \[ \begin{pmatrix} \widehat{\mbox{Var}}(\widehat{\beta}_0) & \widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_1) & \cdots & \widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_p) \\ \widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_1) & \widehat{\mbox{Var}}(\widehat{\beta}_1) & \cdots & \widehat{\mbox{Cov}}(\widehat{\beta}_1,\widehat{\beta}_p) \\ \widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_p) & \widehat{\mbox{Cov}}(\widehat{\beta}_1,\widehat{\beta}_p) & \cdots & \widehat{\mbox{Var}}(\widehat{\beta}_p) \end{pmatrix} \] e pode ser encontrada em R usando a função vcov(). Observe que vcov.glm() está listado como uma função de método na saída methods(class = glm) fornecida no último exemplo. Esta é a função chamada por vcov() quando aplicada a um objeto da classe glm.

Esta matriz é encontrada invertendo uma matriz das segundas derivadas parciais, frequentemente referida como uma “matriz Hessiana”, da função de log-verossimilhança avaliada nas estimativas dos parâmetros e multiplicando esta matriz resultante por -1. No final, a matriz de covariância resultante simplifica para \((X^\top V X)^{-1}\) onde \[ X=\begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix} \] e \(V=\mbox{diag}\big( \widehat{\pi}_i (1-\widehat{\pi}_i)\big)\) uma matriz \(n\times n\) com \(\widehat{\pi}_i (1-\widehat{\pi}_i)\) na diagonal e 0’s em outro lugar. Esta matriz de covariâncias tem a mesma forma que a resultante da estimação dos mínimos quadrados ponderados em um modelo linear. Os elementos individuais da matriz não têm uma forma simples, por isso não os apresentamos aqui.


Exemplo 2.3: Placekicking.


O objetivo deste exemplo é mostrar como obter a matriz de covariância estimada para os estimadores de parâmetros de regressão. Usamos o objeto mod.fit do modelo com distância como única variável explicativa.

Abaixo está a tabela de coeficientes novamente da saída summary(mod.fit):

round ( summary (mod.fit)$coefficients ,4)
##             Estimate Std. Error  z value Pr(>|z|)
## (Intercept)   5.8121     0.3263  17.8133        0
## distance     -0.1150     0.0083 -13.7937        0


Limitamos a saída exibida usando o fato de que summary() cria uma lista com coefficients como um componente. A coluna Std. Error fornece os erros padrão para os estimadores dos parâmetros de regressão, \(\widehat{\mbox{Var}}(\widehat{\beta}_0)^{1/2}\) na linha de (Intercept) e \(\widehat{\mbox{Var}}(\widehat{\beta}_1)^{1/2}\) na linha distance.

A função vcov() produz a matriz de covariância estimada:

vcov (mod.fit )
##             (Intercept)      distance
## (Intercept)  0.10645675 -2.606250e-03
## distance    -0.00260625  6.953996e-05


Podemos extrair a variância estimada para \(\widehat{\beta}_1\) especificando o elemento (2,2) da matriz. Assim, \(\widehat{\mbox{Var}}(\widehat{\beta}_1)= 0.0000695\), que é o quadrado de 0.0083 dado na tabela de coeficientes.

Para leitores interessados nos cálculos reais da matriz, fornecemos um exemplo de como a matriz de covariância é calculada usando \((X^\top VX)^{-1}\):

pi.hat <- mod.fit$fitted.values
V <- diag (pi.hat *(1 - pi.hat))
X <- cbind (1, placekick$distance )
solve (t(X) %*% V %*% X)
##             [,1]          [,2]
## [1,]  0.10645678 -2.606250e-03
## [2,] -0.00260625  6.953997e-05


Criamos a matriz diagonal \(V=\mbox{diag}\big( \widehat{\pi}_i (1-\widehat{\pi}_i)\big)\) usando a função diag(), onde primeiro extraímos \(\widehat{\pi}_i\) de mod.fit. Em seguida, criamos a matriz \(X\) usando a função cbind() que coloca uma coluna de 1’s antes de uma coluna dos valores da variável explicativa. Continuando, a função t() encontra a transposta da matriz \(X\), a sintaxe %*% instrui R a realizar a multiplicação de matrizes e a função solve() encontra a inversa da matriz. A matriz resultante é praticamente a mesma produzida por vcov(), onde as diferenças são devidas ao erro de arredondamento.


Embora a função glm() seja usada sempre que estimarmos um modelo de regressão logística é instrutivo poder visualizar a função de log-verossimilhança e ajustar o modelo de maneira mais geral. No próximo exemplo, ilustramos como a função log-verossimilhança pode ser programada diretamente em R e a função de otimização geral optim() pode ser usada para encontrar as estimativas dos parâmetros. Essas técnicas de programação podem ser utilizadas para problemas mais complicados, onde funções como glm() podem não estar disponíveis.


Exemplo 2.4: Placekicking.


Focamos novamente em ajustar o modelo de regressão logística apenas com a variável explicativa distance. Abaixo, criamos uma função logL() para calcular a função de log-verossimilhança para quaisquer valores de parâmetros fornecidos, valores de variáveis explicativas para \(x_1\) e respostas binárias para \(Y\):

logL <- function (beta , x, Y) {
  pi <- exp( beta [1] + beta [2]* x)/(1 + exp( beta [1] + beta [2]* x)) 
  sum (Y*log(pi) + (1-Y)*log (1- pi))
}
logL ( beta = mod.fit$coefficients , x = placekick$distance , Y = placekick$good )
## [1] -387.8725
logLik (mod.fit )
## 'log Lik.' -387.8725 (df=2)


Dentro da função, encontramos \(\pi_i\) através de nossa expressão de modelo e, em seguida, encontramos

\[ \log \left(L(\beta_0,\beta_1|y_1,\cdots,y_n) \right) = \sum_{i=1}^n y_i\log(\pi_i) +(1-y_i)\log(1-\pi_i) \]

usando as funções sum() e log(), enquanto tira proveito de como R executa adição e multiplicação em vetores de maneira elementar. Quando avaliamos a função usando as estimativas de parâmetros e os dados observados, obtemos \[ \log \left(L(\beta_0,\beta_1|y_1,\cdots,y_n) \right) = -387.87\cdot \]

Como verificação, obtemos o mesmo valor da função genérica logLik(), que extrai o valor máximo da função log-likelihood dos objetos mod.fit, a função logLik.glm() é usada aqui.

Para maximizar a função de log-verossimilhança e encontrar os MLEs correspondentes, usamos a função optim(). Esta é uma função de otimização muito geral que minimiza uma função R em relação a um vetor de parâmetros. Como queremos maximizar a função de log-verossimilhança, usamos o valor do argumento control = list(fnscale = -1) dentro de optim(), que essencialmente instrui R a minimizar o negativo de logL(). Segue abaixo nosso código:

# Find starting values for parameter estimates
reg.mod <- lm( formula = good ~ distance , data = placekick )
reg.mod$coefficients
## (Intercept)    distance 
##  1.25202444 -0.01330212
mod.fit.optim <- optim (par = reg.mod$coefficients , fn = logL , 
                        hessian = TRUE , x = placekick$distance , Y = placekick$good ,
                        control = list ( fnscale = -1) , method = "BFGS")
names (mod.fit.optim )
## [1] "par"         "value"       "counts"      "convergence" "message"    
## [6] "hessian"
mod.fit.optim$par
## (Intercept)    distance 
##   5.8112544  -0.1150046
mod.fit.optim$value
## [1] -387.8725
mod.fit.optim$convergence - solve (mod.fit.optim$hessian )
##              (Intercept)      distance
## (Intercept)  0.106482867 -2.607258e-03
## distance    -0.002607258  6.957463e-05


Os procedimentos de otimização dentro de optim() precisam de valores iniciais para os parâmetros de regressão. Simplesmente ajustamos um modelo de regressão linear aos dados usando a função lm() e tomamos as estimativas dos parâmetros correspondentes como esses valores iniciais.

Dentro da chamada para optim(), especificamos as estimativas de parâmetros iniciais usando o argumento par e especificamos a função a ser maximizada usando o argumento fn. Observe que o primeiro argumento na função nomeada em fn deve corresponder às estimativas dos parâmetros iniciais; é por isso que beta foi dado como o primeiro argumento em logL(). O valor TRUE para o argumento hessiano instrui R a obter uma estimativa numérica da matriz hessiana para os parâmetros. Em outras palavras, ele estima uma matriz de segundas derivadas parciais para a função log-verossimilhança, que pode então ser invertida para obter a matriz de covariância estimada para as estimativas de parâmetros. Finalmente, especificamos os valores correspondentes de x e Y para a função logL().

O objeto produzido por optim() é uma lista que salvamos como mod.fit.optim. Os componentes dentro da lista incluem as estimativas de parâmetro, $par, o valor máximo da função de log-verossimilhança, $value e a matriz de covariância estimada, $hessian; posteriormente usamos solve() para encontrar a matriz de covariância estimada como o inverso da matriz hessiana.

Todos esses valores são praticamente os mesmos produzidos por glm(), onde pequenas diferenças são devidas a diferentes critérios de convergência para os procedimentos numéricos iterativos. O valor 0 para mod.fit.optim$convergence significa que a convergência foi alcançada. Para esta implementação de optim(), usamos o methods = “BFGS”, semelhante a um procedimento de Newton-Raphson, mas outros procedimentos estão disponíveis dentro da função.

A função de log-verossimilhança pode ser plotada usando a função contour() do pacote graphics para plotar os contornos. Podemos ver pelo gráfico que os valores que maximizam esta função são 5.81 para \(\beta_0\) e -0.115 para \(\beta_1\).

# Evaluate the log-likelihood function at a lot of different values for beta0 and beta1
beta0.values<-seq(from = -5, to = 18, by = 0.1)
beta1.values<-seq(from = -0.65, to = 0.25, by = 0.01)
count<-1
save.logL<-numeric(length(beta0.values)*length(beta1.values))
for (beta0 in beta0.values) {
  for (beta1 in beta1.values) {
    save.logL[count]<-logL(beta = c(beta0, beta1), x = placekick$distance, Y = placekick$good)
    count<-count+1
  }
}
max(save.logL)
## [1] -388.0476
save.logL2<-matrix(save.logL, nrow = length(beta0.values), ncol = length(beta1.values), byrow = T)
save.logL[2]  # Verify correct order by checking one value
## [1] -27207.92
save.logL2[1:2,1:2]
##           [,1]      [,2]
## [1,] -27534.45 -27207.92
## [2,] -27408.25 -27081.72
# Contour plot
par(pty = "s")
contour(x = beta0.values, y = beta1.values, z = save.logL2,  xlab = 
expression(beta[0]), ylab = expression(beta[1]), 
                      levels = -c(10000, 7500, 5000, 2500, 1000, 750, 500, 450, 400))
abline(h = mod.fit$coefficients[2], lty = "dashed", lwd = 2, col= "red")
abline(v = mod.fit$coefficients[1], lty = "dashed", lwd = 2, col= "red")
beta0.1.values<-expand.grid(beta1.values, beta0.values)
save500 <- data.frame(beta0.1.values[save.logL < -499 & save.logL > -501,], 
                      logL = save.logL[save.logL < -499 & save.logL > -501]) 
points(x = save500$Var2, y = save500$Var1, pch = 20, col = "darkblue")
grid()


Uma forma alternativa da variável de resposta está disponível quando há vários ensaios de Bernoulli registrados em algumas ou todas as combinações das variáveis explicativas. Em particular, se houver \(J < n\) combinações únicas de variáveis explicativas, podemos agregar os resultados das tentativas individuais como na Seção 1.2, de modo que observemos \(w_j\) sucessos em \(n_j\) tentativas, cada uma com \(x_{j1},\cdots,x_{jp}\) valores da variável explicativa para \(j = 1,\cdots,J\). Em outras palavras, agora temos \(J\) observações de uma distribuição binomial com probabilidades de sucesso \(\pi_j\), \(j = 1,\cdots,J\).

Alternativamente, os dados podem surgir naturalmente na forma de \(w_j\) sucessos de \(n_j\) tentativas. Por exemplo, em um experimento para determinar um nível de dose letal, uma determinada quantidade de pesticida \(x_{j1}\) poderia ser aplicada a um recinto com um número conhecido de formigas \(n_j\), onde o número de formigas exterminadas \(w_j\) seria registrado. Ao aplicar diferentes quantidades de pesticidas a outros recintos, temos uma configuração de resposta binomial para os \(J\) grupos de tratamento.

O mesmo modelo de regressão logística pode ser usado para modelar a probabilidade de sucesso para a configuração de resposta binomial. A função de log-verossimilhança é \[ \log\big( L(\beta_0,\cdots,\beta_p|w_1,\cdots,w_J)\big) = \sum_{j=1}^J \log\binom{n_j}{w_j}+w_j\log(\pi_j)+(n_j-w_j)\log(1-\pi_j)\cdot \]

Como \(\binom{n_j}{w_j}\) é uma constante em relação aos parâmetros, esse termo não altera os valores dos parâmetros estimados que resultam na maximização da log-verossimilhança. Portanto, os MLEs são os mesmos como se os dados consistissem em \(n = \sum_{j = 1}^J n_j\) respostas binárias. Ilustramos isso no próximo exemplo.


Exemplo 2.5: Placekicking.


O objetivo deste exemplo é mostrar como reformar os dados de placekicking para um formato de resposta binomial e mostrar que as estimativas dos parâmetros do modelo são as mesmas para esse formato. Focamos apenas na variável explicativa distance para este exemplo.

A função aggregate() é usada para encontrar o número de sucessos e o número de observações para cada distância:

w <- aggregate ( good ~ distance , data = placekick , FUN = sum )
n <- aggregate ( good ~ distance , data = placekick , FUN = length )
w.n <- data.frame ( distance = w$distance , success = w$good , 
                    trials = n$good , proportion = round ( w$good /n$good ,4) )
head (w.n)
##   distance success trials proportion
## 1       18       2      3     0.6667
## 2       19       7      7     1.0000
## 3       20     776    789     0.9835
## 4       21      19     20     0.9500
## 5       22      12     14     0.8571
## 6       23      26     27     0.9630


O argumento de formula dentro de aggregate() é usado da mesma forma que dentro de glm(). O argumento FUN especifica como a variável de resposta good será resumida, onde a função sum() adiciona as respostas 0 e 1 para cada distância e a função length() conta o número de observações para cada distância. O resultado final está no arquivo de dados w.n. Por exemplo, há 2 sucessos em 3 tentativas a uma distância de 18 jardas, o que resulta em uma proporção observada de sucessos de \(2/3\approx 0.6667\). Observe que a razão para o grande número de observações a 20 jardas é porque a maioria dos PATs são tentados a partir dessa distância.

O modelo de regressão logística é estimado pela função glm() com duas alterações em nosso código anterior:

mod.fit.bin <- glm( formula = success / trials ~ distance , 
                    weights = trials , family = binomial ( link = logit ), data = w.n)
summary (mod.fit.bin)
## 
## Call:
## glm(formula = success/trials ~ distance, family = binomial(link = logit), 
##     data = w.n, weights = trials)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0373  -0.6449  -0.1424   0.5004   2.2758  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.812080   0.326277   17.81   <2e-16 ***
## distance    -0.115027   0.008339  -13.79   <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: 282.181  on 42  degrees of freedom
## Residual deviance:  44.499  on 41  degrees of freedom
## AIC: 148.46
## 
## Number of Fisher Scoring iterations: 5


Primeiro, especificamos a resposta no argumento formula como success/trials (sucesso/ensaios). Observe que isso encontra a proporção observada de sucesso para cada distância. Em segundo lugar, usamos o argumento weights para especificar o número de observações por distância.

O modelo de regressão logística estimado é \[\begin{equation*} logit(\widehat{\pi}) = 5.8121 - 0.1150 \, {\bf \mbox{distance}}, \end{equation*}\] que é essencialmente o mesmo que anteriormente, diferenças muito pequenas podem ocorrer nas estimativas de parâmetros devido ao procedimento numérico iterativo.


2.2.2 Testes de hipóteses para parâmetros de regressão


Testes de hipóteses podem ser usados para avaliar a importância de variáveis explicativas em um modelo. Por exemplo, um teste de \(H_0 : \beta_r = 0\) vs. \(H_1 : \beta_r \neq 0\) avalia o \(r\)-ésimo termo da variável explicativa no modelo \[ logit(\pi ) = \beta_0+\beta_1 x_1 +\cdots +\beta_{r} x_{r} + +\cdots + \beta_p x_p\cdot \] Se \(\beta_r = 0\), isso significa que o termo correspondente é excluído do modelo. Se \(\beta_r \neq 0\), significa que o termo correspondente está incluído no modelo.

De forma equivalente, podemos enunciar as hipóteses como \[ H_0 : \, logit(\pi) = \beta_0+\beta_1 x_1 +\cdots +\beta_{r-1} x_{r-1} +\beta_{r+1} x_{r+1} +\cdots + \beta_p x_p \\ H_1 : \, logit(\pi) = \beta_0+\beta_1 x_1 +\cdots +\beta_{r} x_{r} + +\cdots + \beta_p x_p \cdot \] o que pode ser útil para enfatizar a comparação de dois modelos diferentes. Em geral, as variáveis explicativas no modelo de hipótese nula, também conhecido como modelo reduzido, devem estar todas no modelo de hipótese alternativa, também conhecido como modelo completo.

Usaremos um teste Wald ou um LRT para realizar esses testes. Agora nos concentramos nas especificidades da configuração de regressão logística.


Teste Wald


A estatística Wald \[ Z_0 = \dfrac{\widehat{\beta}_r}{\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_r)}} \]

é usado para testar \(H_0 : \beta_r = 0\) vs. \(H_1 : \beta_r\neq 0\). Se a hipótese nula for verdadeira, \(Z_0\) tem uma distribuição normal padrão aproximada para uma amostra grande e rejeitamos a hipótese nula se \(Z_0\) tiver um valor observado incomum para esta distribuição. Definimos incomum por \(|Z_0| > Z_{1-\alpha/2}\). O \(p\)-valor é \(2P(Z > |Z_0|)\) onde \(Z\) tem uma distribuição normal padrão.

Mais de um parâmetro pode ser testado para igualdade a 0 na hipótese nula. No entanto, renunciamos à discussão sobre isso porque os procedimentos de inferência de Wald aqui geralmente encontram os mesmos tipos de problemas discutidos no Capítulo 1. No contexto de testes de hipóteses, isso significa que a taxa de erro do tipo I declarada para um teste de hipóteses, ou seja \(\alpha\), não é a mesma que a taxa de erro real do tipo I.

O teste da razão de verossimilhanças (LRT) tipicamente tem um desempenho melhor do que o teste de Wald, então nos concentramos neste procedimento a seguir.


Teste da razão de verossimilhanças (LRT)


A estatística LRT pode ser escrita informalmente como \[ \Lambda =\dfrac{\mbox{Máximo da função de verossimilhança sob } H_0}{\mbox{Máximo da função de verossimilhança sob } H_0 \mbox{ ou } H_1}\cdot \]

No contexto de testar a igualdade dos parâmetros de regressão a 0, o denominador na equação acima (\(\Lambda\)) é a função de verossimilhança avaliada nas MLEs para o modelo contendo todos os \(p+1\) parâmetros de regressão. O numerador desta equação é a função de verossimilhança avaliada nas MLEs para o modelo que exclui aquelas variáveis cujos parâmetros de regressão são definidos como 0 na hipótese nula.

Por exemplo, para testar \(H_0 : \beta_r = 0\) vs. \(H_1 : \beta_r\neq 0\), \(\beta_r\) seria mantida igual a 0, o que significa que essa variável explicativa correspondente seria excluída do modelo de hipótese nula. É claro que as estimativas para os parâmetros restantes não precisam ser as mesmas do modelo de hipóteses alternativas devido às diferenças entre os dois modelos.

Se \(q\) parâmetros de regressão forem definidos como 0 na hipótese nula e se a hipótese nula for verdadeira, a estatística \(-2 \log(\Lambda )\) tem uma distribuição aproximada de \(\chi^2_q\) para uma amostra grande e rejeitamos a hipótese nula se \(-2\log(\Lambda )\) tem um valor observado incomumente grande para essa distribuição. Por exemplo, se \(\alpha= 0.05\) para o teste de \(H_0 : \beta_r = 0\) vs. \(H_1 : \beta_r\neq 0\), a região de rejeição é > 3.84, onde \(\chi^2_{1,0.95} = 3.84\) é o quantil 0.95 a partir de uma distribuição qui-quadrado. O \(p\)-valor é \(P(A > -2\log(\Lambda ))\) onde \(A\) tem uma distribuição \(\chi^2_1\) quando \(q = 1\).

Na maioria das vezes, usaremos as funções genéricas anova() do pacote stats e Anova() do pacote car para realizar esses testes, funções comumente usadas por outros incluem a função drop1() do pacote stats e a função lmtest() do pacote lmtest.

Essas funções testam hipóteses de uma estrutura semelhante àquelas normalmente vistas em uma análise de variância (ANOVA), mas podem realizar LRTs em vez dos típicos testes \(F\) ANOVA. As duas funções são baseadas em testes de comparação de modelos de diferentes tipos. Quando usado com um objeto de ajuste de modelo resultante de glm(), a função anova() calcula testes tipo I (sequenciais), enquanto Anova() calcula testes tipo II (parciais).

A diferença entre os tipos está nos modelos de hipótese nula usados para cada termo testado. Com os testes do tipo II, cada modelo de hipótese nula consiste em todas as outras variáveis listadas no lado direito do argumento formula, ignorando quaisquer interações de ordem superior que contenham o termo, consulte a Seção 2.2.5 para obter mais informações sobre interações. Para testes do tipo 1, o modelo de hipótese nula contém apenas as variáveis listadas no argumento formula antes do termo testado. Geralmente, os testes do tipo 2 são preferidos, a menos que haja algum motivo específico para considerar uma sequência de testes, como em modelos polinomiais. Ver Milliken and Johnson (2004) para mais detalhes.

A estatística LRT transformada \(-2\log(\Lambda )\) tem uma forma simplificada. Suponha que a probabilidade estimada de sucesso nos modelos sob a hipótese nula e alternativa sejam denotados como \(\widehat{\pi}^{(0)}_i\) e \(\widehat{\pi}^{(1)}_i\), respectivamente. Da mesma forma, definimos um vetor das estimativas dos parâmetros da regressão como \(\widehat{\beta}^{(0)}\) e \(\widehat{\beta}^{(1)}\). Temos então \[ \begin{array}{rcl} -2\log(\Lambda) & = & -2\log\left( \dfrac{L\big(\widehat{\beta}^{(0)}|y_1,\cdots,y_n\big)}{L\big(\widehat{\beta}^{(1)}|y_1,\cdots,y_n\big)}\right) \\ & = & -2\left(\log\Big(L\big(\widehat{\beta}^{(0)}|y_1,\cdots,y_n\big)\Big)-\log\Big(L\big(\widehat{\beta}^{(1)}|y_1,\cdots,y_n\big)\Big) \right) \\ & = & \displaystyle -2\sum_{i=1}^n y_i \log\big(\widehat{\pi}^{(0)}_i\big)+(1-y_i)\log\big(1-\widehat{\pi}^{(0)}_i\big)-y_i\log\big(\widehat{\pi}^{(1)}_i\big)-(1-y_i)\log\big(\widehat{\pi}^{(1)}_i\big)\cdot \end{array} \] Essa formulação da estatística pode ser facilmente codificada em R. Mostramos a seguir e como usar as funções anova() e Anova().


Exemplo 2.6: Placekicking.


Começamos este exemplo usando o modelo que contém tanto change quanto distance, \[ logit(\pi)=\beta_0+\beta_1 \, \mbox{change} +\beta_2 \, \mbox{distance}, \]

que foi originalmente estimado na Seção 2.2.1 e seus resultados salvos no objeto mod.fit2. Abaixo está o código e saída para testes Wald.

summary (mod.fit2 )
## 
## Call:
## glm(formula = good ~ change + distance, family = binomial(link = logit), 
##     data = placekick)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7061   0.2282   0.2282   0.3750   1.5649  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.893181   0.333184  17.687   <2e-16 ***
## change      -0.447783   0.193673  -2.312   0.0208 *  
## distance    -0.112889   0.008444 -13.370   <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: 1013.4  on 1424  degrees of freedom
## Residual deviance:  770.5  on 1422  degrees of freedom
## AIC: 776.5
## 
## Number of Fisher Scoring iterations: 6


As estatísticas do teste Wald e os \(p\)-valores para cada parâmetro de regressão estão localizados em coefficients na tabela de summary(mod.fit2).

Por exemplo, para testar a importância da variável explicativa change, usamos \(H_0 : \beta_1 = 0\) vs. \(H_1 : \beta_1 \neq 0\) resultando em \(Z_0 = -0.4478/0.1936 = -2.3124\), coluna z value e um \(p\)-valor de \(2P(Z > |-2.3124|) = 0.0208\), coluna Pr(>|Z|). Usando \(\alpha = 0.05\), rejeitaríamos a hipótese nula; no entanto, não o rejeitaríamos usando \(\alpha = 0.01\). Assim, dizemos que há evidências marginais de que change é importante para incluir no modelo, uma vez que distance está no modelo. Observe que é importante incluir a parte “distance está no modelo” porque as hipóteses nula e alternativa são equivalentemente declaradas como: \[ H_0 : logit(\pi) =\beta_0+\beta_1 \, \mbox{distance}\\ H_1: logit(\pi) +\beta_0+\beta_1 \, \mbox{change} +\beta_2 \, \mbox{distance}\cdot \]

Com relação ao teste do parâmetro de distance no modelo que contém change, o \(p\)-valor é dado como <2e-16 na saída, o que significa que o \(p\)-valor é menor que 2\(\times 10^{-16}\). Assim, há fortes evidências da importância de distance, uma vez que change está no modelo. No programa correspondente a este exemplo, mostramos o código R necessário para extrair as informações do mod.fit2 para realizar os testes de Wald sem a ajuda do summary(mod.fit2).

Abaixo está o código e saída de Anova() para os LRTs:

library( package = car)
Anova(mod.fit2 , test = "LR")
## Analysis of Deviance Table (Type II tests)
## 
## Response: good
##          LR Chisq Df Pr(>Chisq)    
## change      5.246  1      0.022 *  
## distance  218.650  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


A função Anova() produz LRTs usando o valor do argumento test = “LR”, que é o padrão. Para o teste de change com \(H_0 : \beta_1 = 0\) vs. \(H_1 : \beta_1\neq 0\), obtemos \(-2\log(\Lambda ) = 5.246\) com um \(p\)-valor de \(P(A > 5.246) = 0.0220\) e chegamos à mesma conclusão do teste de Wald anterior. O \(p\)-valor para o teste de distance é dado como < 2e-16 na saída, o que novamente indica que há fortes evidências de que distance é importante, dado que o modelo inclui change.

A função anova() com o valor do argumento test = “Chisq” também realiza LRTs, mas de forma sequencial com base na ordenação das variáveis explicativas no ajuste do modelo:

anova (mod.fit2 , test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: good
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      1424    1013.43              
## change    1   24.277      1423     989.15 8.343e-07 ***
## distance  1  218.650      1422     770.50 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


A linha de change da saída da função fornece um teste de \[ H_0 : logit(\pi) = \beta_0 \\ H_1 : logit(\pi )= \beta_0+ \beta_1 \, \mbox{change}, \] que resulta em um \(p\)-valor de 8.34\(\times 10^{-7}\). A linha de distance da saída da função fornece um teste de \[ H_0 : logit(\pi) = \beta_0 + \beta_1 \,\mbox{change} \\ H_1 : logit(\pi )= \beta_0 + \beta_1 \, \mbox{change} + \beta_2 \, \mbox{change}, \] o que resulta num \(p\)-valor inferior a 2\(\times 10^{-16}\). As hipóteses e cálculos para o teste de distance são exatamente os mesmos da função Anova().

Outra maneira de usar a função anova() é ajustar os modelos de hipótese nula e alternativa e, em seguida, usar os objetos de ajuste de modelo correspondentes como valores de argumento dentro da função. Isso será especialmente útil mais tarde, quando testarmos hipóteses nulas que não correspondem a simples deleções de um termo do modelo completo.

Para testar a variável change dado que distance está no modelo, podemos usar o seguinte código:

mod.fit.H0 <- glm( formula = good ~ distance , family = binomial ( link = logit ), data = placekick )
anova (mod.fit.H0 , mod.fit2 , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: good ~ distance
## Model 2: good ~ change + distance
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1      1423     775.75                       
## 2      1422     770.50  1   5.2455    0.022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Obtemos os mesmos resultados de antes com \[ -2\log(\Lambda ) = 5.246 \] e um \(p\)-valor de \(P(A > 5.246) = 0.0220\).

A equação acima pode ser programada diretamente em R para realizar LRTs. Embora isso geralmente não seja necessário para aplicações de regressão logística, geralmente é útil passar pelos cálculos dessa maneira para entender melhor quais funções anova() e Anova() estão computando. Além disso, haverá momentos mais adiante em que funções como essas não estarão disponíveis. Demonstramos o LRT para o teste de \(H_0 : logit(\pi ) = \beta_0\) vs. \(H_1 : logit(\pi ) = \beta_0 + \beta_1 \, \mbox{change}\):

mod.fit.Ho <- glm( formula = good ~ 1, family = binomial ( link = logit ), data = placekick )
mod.fit.H1 <- glm( formula = good ~ change , family = binomial ( link = logit ), data = placekick )
anova (mod.fit.Ho , mod.fit.H1 , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: good ~ 1
## Model 2: good ~ change
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1424    1013.43                          
## 2      1423     989.15  1   24.277 8.343e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pi.hat.Ho <- mod.fit.Ho$fitted.values
pi.hat.H1 <- mod.fit.H1$fitted.values
y <- placekick$good
stat <- -2*sum(y*log(pi.hat.Ho/pi.hat.H1) + (1-y)*log ((1 - pi.hat.Ho)/(1 - pi.hat.H1)))
pvalue <- 1- pchisq (q = stat , df = 1)
data.frame (stat , pvalue )
##       stat       pvalue
## 1 24.27703 8.342813e-07


Ajustamos os modelos de hipótese nula e alternativa usando glm(). Observe que usamos formula = good ~ 1 para estimar um modelo com apenas o intercepto. Por outro lado, se quiséssemos excluir o parâmetro de intercepto da estimativa para um modelo, incluiríamos um “-1” como parte de um valor de argumento formula junto com quaisquer nomes de variáveis explicativas.

A função anova() fornece um valor de \(-2\log(\lambda ) = 24.277\) junto com um \(p\)-valor igual a 8.343\(\times 10^{-7}\) como visto antes. Depois de criar objetos para \(\widehat{\pi}^{(0)}_i\), \(\widehat{\pi}^{(1)}_i\) e \(y_i\), codificamos a equação de verossimilhança usando o fato de que R realiza operações matemáticas em vetores elemento a elemento. Lembre-se de que é útil executar segmentos de código para determinar o que o código faz. Por exemplo, pode-se primeiro executar y*log(pi.hat.H0/pi.hat.H1) separadamente para ver que um vetor de valores zero e diferente de zero é retornado. Quando combinado com (1-y)*log((1-pi.hat.H0)/(1-pi.hat.H1)) e então através do uso de -2*sum(), obtemos a equação de verossimilhança.

Para encontrar o \(p\)-valor, usamos o pchisq() e subtraí-lo de 1 porque queremos a área à direita do quantil, mas a função encontra a área à sua esquerda. A estatística de teste resultante e o \(p\)-valor são os mesmos encontrados anteriormente.


Na saída do exemplo anterior, a palavra deviance e sua abreviação “dev” apareceram várias vezes. Essa palavra também aparece na saída summary(). O desvio (deviance) refere-se à quantidade que um determinado modelo se desvia de outro modelo medido pela estatística LRT transformada \(-2 \log(\Lambda )\). Por exemplo, o valor \(-2\log(\Lambda ) = 5.246\) usado para testar change, dado que distance está no modelo é uma medida de quanto a probabilidade estimada de sucesso para o modelo excluindo change se desvia daquelas para o modelo incluindo change.

O desvio residual é uma forma particular da estatística \(-2\log(\Lambda)\) que mede o quanto as probabilidades estimadas de um modelo de interesse se desviam das proporções observadas de sucessos, simplesmente 0/1 ou 1/1 para dados de resposta de Bernoulli. Essas proporções observadas são equivalentes a estimar um modelo de regressão logística onde cada observação é representada por um parâmetro, digamos \(logit(\pi) = \gamma_i\) para \(i = 1,\cdots,n\). Esse modelo é frequentemente chamado de modelo saturado, porque o número de parâmetros é igual ao número de observações, de modo que nenhum parâmetro adicional pode ser estimado.

O desvio nulo denota o quanto as probabilidades estimadas a partir do modelo \(logit(\pi) = \beta_0\) para \(i = 1,\cdots,n\) desviam-se da proporção observada de sucessos. Como \(logit(\pi) = \beta_0\) contém apenas o termo de intercepto e, portanto, apenas um parâmetro, cada \(i\) é estimado como sendo o mesmo valor para esse modelo específico.

As estatísticas de desvio residual são frequentemente calculadas como uma etapa intermediária para realizar um LRT para comparar dois modelos. Por exemplo, considere as hipóteses de \[ H_0 : logit(\widehat{\pi}^{(0)}) = \beta_0 + \beta_1 \, x_1 \\ H_1 : logit(\widehat{\pi}^{(1)}) = \beta_0 + \beta_1 \, x_1 + \beta_2 \, x_2 + \beta_3 \, x_3, \]

onde usamos o sobrescrito para novamente para ajudar a diferenciar entre os dois modelos e suponha que estamos trabalhando com dados de resposta de Bernoulli.

O desvio residual para o modelo \(logit(\pi^{(0)})=\beta_0+\beta_1 x_1\) testando \(H_0 : logit(\pi^{(0)})=\beta_0+\beta_1 x_1\) vs. \(H_1 : \mbox{Modelo saturado}\) é \[ -2\log(\Lambda) = -2\sum_{i=1}^n y_i\log\left(\dfrac{\pi^{(0)}_i}{y_i}\right)+(1-y_i)\log\left(\dfrac{1-\pi^{(0)}_i}{1-y_i}\right) \cdot \]

Valores de \(y_i = 0\) ou 1 nos denominadores das frações não causam problemas devido ao fato de \(y_i\) ser um coeficiente na função logarítmica, ou seja, 0 multiplicado por outra quantidade é 0.

O desvio residual para o modelo \[ logit(\pi^{(1)}) = \beta_0 + \beta_1 x_1 +\beta_2 x_2 + \beta_3 x_3 \] testando \(H_0 : logit(\pi^{(1)}) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3\) vs. \(H_1 : \mbox{Modelo saturado}\) é obtido usando a quantidade: \[ -2\log(\Lambda) = -2\sum_{i=1}^n y_i\log\left(\dfrac{\pi^{(1)}_i}{y_i}\right)+(1-y_i)\log\left(\dfrac{1-\pi^{(1)}_i}{1-y_i}\right) \cdot \]

O desvio residual ou deviance assume uma forma ligeiramente diferente quando temos respostas binomiais \(W_j\) para \(j = 1,\cdots,J\). A estatística é \[ -2\log(\Lambda) = -2\sum_{j=1}^J w_j\log\left(\dfrac{\widehat{\pi}_j}{w_j/n_j}\right)+(n_j-y_j)\log\left(\dfrac{1-\widehat{\pi}_j}{1-w_j/n_j}\right), \cdot \]

onde \(\widehat{\pi}_j\) é a estimativa da probabilidade de sucesso em cada \(j = 1,\cdots,J\). Observe que o modelo saturado nesta configuração estima a probabilidade de sucesso como \(w_j/n_j\). Desvios residuais para modelos de hipóteses nulas e alternativas ainda podem ser usados em testes de hipóteses da mesma maneira que no caso de resposta de Bernoulli anterior.


Exemplo 2.7: Placekicking.


Para realizar o teste de \[ H_0 : logit(\pi) = \beta_0 + \beta_1 \,\mbox{distance} \\ H_1 : logit(\pi) = \beta_0 + \beta_1 \, \mbox{change} + \beta_2 \, \mbox{distance}, \] usamos a função anova() para obter \(-2\log(\Lambda) = 5.246\).

Alternativamente, podemos usar os desvios residuais dos dois modelos:

mod.fit.Ho <- glm( formula = good ~ distance , 
                   family = binomial ( link = logit ), data = placekick )
df <- mod.fit.H0$df.residual - mod.fit2$df.residual
stat <- mod.fit.H0$deviance - mod.fit2$deviance
pvalue <- 1 - pchisq (q = stat , df = df)
data.frame (H0.resid.dev = mod.fit.H0$deviance , 
            H1.resid.dev = mod.fit2$deviance , df = df , stat = round (stat ,4), 
            pvalue = round (pvalue ,4))
##   H0.resid.dev H1.resid.dev df   stat pvalue
## 1      775.745     770.4995  1 5.2455  0.022


Usamos os componentes deviance e df.residual dos objetos de ajuste do modelo para obter os desvios residuais e os graus de liberdade correspondentes, respectivamente. O valor \(-2\log(\Lambda )\) e o \(p\)-valor são os mesmos de antes.


Nesta seção, ilustramos algumas maneiras diferentes de realizar um LRT. Geralmente, as funções Anova() e anova() são as funções mais fáceis de usar e recomendamos seu uso. Incluímos outras maneiras de realizar os testes para ajudar a demonstrar métodos mais gerais e reforçar como os cálculos são realmente feitos. Esses métodos mais gerais também serão úteis para outros tipos de modelos onde uma função do tipo glm() não existe para ajuste de modelo e nenhuma função Anova() ou anova() existe para teste.


2.2.3 Razões de chances (odds ratios)


Em um modelo linear onde \(\mbox{E}(Y) =\beta_0 + \beta_1 x_1 +\cdots + \beta_p x_p\), o parâmetro de regressão \(\beta_r\) é interpretado como a mudança na resposta média para cada aumento de 1 unidade em \(x_r\), mantendo as outras variáveis no modelo constantes. Em um modelo de regressão logística, a interpretação dos parâmetros de regressão precisa levar em conta o fato de que eles estão relacionados à probabilidade de sucesso através de \(logit(\pi ) =\beta_0 + \beta_1 x_1 +\cdots + \beta_p x_p\). Mantendo as outras variáveis constantes, um aumento de 1 unidade em \(x_r\) faz com que \(logit(\pi )\) mude por \(\beta_r\). Essa mudança constante nas probabilidades logarítmicas de um sucesso leva a uma maneira conveniente de usar as probabilidades na interpretação.

Para ver como isso funciona, considere um modelo de regressão logística com apenas uma variável explicativa \(x\). As chances de um sucesso em um determinado valor de \(x\) são \[ Odds_x = \exp(\beta_0+\beta_1 x)\cdot \]

Se \(x\) for aumentado em \(c > 0\) unidades, as chances de sucesso se tornam \[ Odds_{x+c} = \exp\big(\beta_0+\beta_1 (x+c)\big)\cdot \] Para determinar o quanto as chances de sucesso mudaram por esse aumento de \(c\) unidades, encontramos a razão de chances: \[ OR = \dfrac{Odds_{x+c}}{Odds_x} = \dfrac{\exp\big(\beta_0+\beta_1 (x+c)\big)}{\exp\big(\beta_0+\beta_1 x)\big)} = \exp(c \,\beta_1)\cdot \]

Curiosamente, o valor original da variável explicativa \(x\) é cancelado na simplificação; apenas a quantidade de aumento \(c\) e o coeficiente \(\beta_1\) importam. A interpretação padrão desta razão de chances é \[ \mbox{As chances de sucesso mudam em } \exp(c\, \beta_1) \mbox{ vezes para cada } c \mbox{ unidades aumentadas em } x\cdot \]

Também é comum dizer “aumentar” ao invés de “mudar” quando \(\exp(c\, \beta_1) > 1\), e “diminuir” quando \(\exp(c\, \beta_1) < 1\). Quando \(x\) é binário com codificação de 0 e 1, \(c\) é sempre 1. A interpretação da razão de chances ou odds ratio então compara as chances de sucesso no nível codificado “\(x\) = 1” com as chances em “\(x\) = 0” da mesma forma que na Seção 1.2.5.

As estimativas dos parâmetros do modelo podem ser substituídas por seus parâmetros correspondentes em \(OR\) para estimar a razão de chances. A razão de chances estimada torna-se \[ \widehat{OR} = \exp(c \, \widehat{\beta}_1), \]

e sua interpretação é que as chances estimadas de um sucesso mudam por \(\exp(c\, \widehat{\beta}_1)\) vezes para cada \(c\) unidades de aumento em \(x\).

Como a razão de chances estimada é uma estatística, ela varia de amostra para amostra. Portanto, precisamos encontrar um intervalo de confiança para \(OR\) a fim de fazer inferências com um determinado nível de confiança. Os procedimentos baseados em verossimilhanças fornecem a base para os intervalos de confiança discutidos a seguir.

Para encontrar um intervalo de confiança de Wald para \(OR\), primeiro encontramos um intervalo de confiança para \(c\, \beta_1\) e, em seguida, usamos a função exponencial com os pontos finais do intervalo. Para este fim, extraímos \(\widehat{\mbox{Var}}(\widehat{\beta}_1)\) da matriz de covariâncias estimada para os estimadores dos parâmetros e formamos o \((1-\alpha)100\)% intervalo de confiança de Wald para \(c\, \beta_1\) como \[ c\,\widehat{\beta}_1 \pm c Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_1)}, \] onde ustilizamos \(\widehat{\mbox{Var}}(c \,\widehat{\beta}_1) = c^2\widehat{\mbox{Var}}(\widehat{\beta}_1)\).

O \((1-\alpha)100\)% intervalo de confiança de Wald para \(OR\) torna-se \[ \exp\left( c\,\widehat{\beta}_1 \pm c Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_1)}\right)\cdot \]

A interpretação padrão do intervalo de confiança é \[ \mbox{Com } (1-\alpha) \mbox{ 100% de confiança, as chances de sucesso mudam por um valor } \\ \mbox{entre <limite inferior> e <limite superior> vezes para cada aumento de } c \mbox{ unidade em } x, \]

onde os valores numéricos apropriados são inseridos dentro de < >.

O intervalo de confiança de Wald geralmente tem um nível de confiança verdadeiro próximo ao intervalo de confiança declarado somente quando há amostras grandes. Quando o tamanho da amostra não é grande, os intervalos de confiança LR perfilada geralmente apresentam melhor desempenho. Para esta configuração, encontramos o conjunto de valores de \(\beta_1\) tais que \[ -2\log\left( \dfrac{L(\widehat{\beta}_0,\beta_1 \, | \, y_1,\cdots,y_n)}{L(\widehat{\beta}_0,\widehat{\beta}_1 \, | \, y_1,\cdots,y_n)}\right) < \chi^2_{1,1,-\alpha} \]

seja satisfeito, onde \(\widehat{\beta_0}\) é o MLE de \(\beta_0\) seguindo a especificação de um valor de \(\beta_1\).

Na maioria das configurações, não há soluções de forma fechada para os limites inferior e superior, portanto, são necessários procedimentos numéricos iterativos para encontrá-los. Uma vez que os limites do intervalo de confiança para \(\beta_1\) são encontrados, digamos, “lower” e “upper”, usamos a função exponencial e levamos em consideração um valor de \(c\) para encontrar o intervalo de confiança LR perfilada \((1-\alpha)100\)% para \(OR\): \[ \exp(c\times lower) <OR < \exp(c\times upper)\cdot \]

Algumas notas adicionais são necessárias sobre as razões de chances antes de prosseguir para um exemplo:


Exemplo 2.8: Placekicking.


Para este exemplo, utilizamos o objeto mod.fit do modelo de regressão logística que utiliza apenas distance como variável explicativa. A razão de chances estimada para distance é calculada como:

exp(mod.fit$coefficients [2])
##  distance 
## 0.8913424
exp (-10* mod.fit$coefficients [2])
## distance 
## 3.159035


Vemos que \(\exp(\widehat{\beta}_1) = 0.8913\) com \(c = 1\). Como um incremento de 1 jarda é bastante pequeno, as tentativas de field goal variam aproximadamente de 20 a 60 jardas, focamos na mudança nas chances de sucesso para um incremento de 10 jardas. Além disso, como as chances estimadas de sucesso são menores para um aumento na distância \(\exp(c\, \widehat{\beta}_1) < 1\) para \(c > 0\), focamos em \(c = -10\) para nossa interpretação primária.

Descobrimos que as chances estimadas de sucesso mudam em 3.16 vezes para cada diminuição de 10 jardas na distância do chute. Se um treinador de futebol pudesse escolher entre tentar um chute de 50 jardas ou tentar ganhar mais 10 jardas e tentar um chute de 40 jardas ou qualquer outra diminuição de 10 jardas, o fato de que as chances estimadas de sucesso são 3.16 vezes maior para a distância mais curta pode influenciar a decisão do treinador. Portanto, a preferência seria por placekicks mais curtos. Embora os fãs de futebol não se surpreendam que chutes mais curtos tenham uma chance maior de sucesso, a quantidade pela qual as chances de sucesso mudam não é amplamente conhecida.

A função confint() fornece uma maneira conveniente de calcular intervalos de confiança associados a parâmetros de regressão. Esta é uma função genérica que, por padrão, encontra intervalos LR perfilada para nossa configuração. Usamos esta função abaixo para calcular o intervalo correspondente à distância do placekick:

beta.ci <- confint ( object = mod.fit , parm = "distance", level = 0.95)
beta.ci
##       2.5 %      97.5 % 
## -0.13181435 -0.09907103
rev(exp ( -10* beta.ci)) # OR C.I. for c = -10
##   97.5 %    2.5 % 
## 2.693147 3.736478
# Remove labels with as. numeric ()
as.numeric (rev (exp ( -10* beta.ci)))
## [1] 2.693147 3.736478


O intervalo de confiança LR perfilada de 95% para o parâmetro de distância é \[ -0.1318 < \beta_1 < -0.0991\cdot \] Usando \(c = -10\), o intervalo LR perfilada de 95% para a razão de chances é \(2.69 < OR < 3.74\) onde \(OR = \exp(-10\, \beta_1)\). Com 95% de confiança, as chances de sucesso mudam de 2.69 a 3.74 vezes para cada diminuição de 10 jardas na distância do placekick.

Como o intervalo está inteiramente acima de 1, há evidências suficientes de que uma diminuição de 10 jardas na distância aumenta as chances de um chute bem-sucedido. Observe que usamos as.numeric() no código anterior para evitar que rótulos desnecessários sejam impressos.

Para calcular um intervalo Wald, precisamos usar a função específica confint.default():

beta.ci <- confint.default ( object = mod.fit , parm = "distance", level = 0.95)
beta.ci
##               2.5 %     97.5 %
## distance -0.1313709 -0.0986824
rev (1/ exp( beta.ci *10) ) # Invert OR C.I. for c = 10
## [1] 2.682701 3.719946


O intervalo de confiança de Wald de 95% para o parâmetro de distance é \[ -0.1314 < \beta_1 < -0.0987\cdot \] Usando \(c = -10\), obtemos o intervalo de Wald de 95% \(2.68 < OR < 3.72\). Esse intervalo é semelhante ao intervalo LR perfilada devido ao grande tamanho da amostra.

Para ver como os cálculos para o intervalo Wald são realizados, o código abaixo mostra a equação correspondente codificada diretamente em R:

beta.ci <- mod.fit$coefficients [2] + qnorm (p = c (0.025 , 0.975) ) * sqrt ( vcov ( mod.fit)[2 ,2])
beta.ci
## [1] -0.1313709 -0.0986824
rev (1/ exp( beta.ci *10) )
## [1] 2.682701 3.719946


A função vcov() calcula a matriz de covariância estimada para as estimativas dos parâmetros usando as informações dentro do mod.fit. Ao especificar vcov(mod.fit)[2,2], extraímos \(\widehat{\mbox{Var}}(\widehat{\beta}_1)\) da matriz. A sintaxe mod.fit$coeficientes[2] extrai \(\widehat{\beta}_1\) do vetor de estimativas de parâmetros. Juntando esses elementos, calculamos o intervalo de confiança para \(\beta_1\) e, em seguida, o intervalo de confiança desejado para a razão de chances.

Para ver como os cálculos para o intervalo LR perfilada são realizados, incluímos dois conjuntos de código no programa correspondente que calcula o intervalo sem a função confint(). Ambos dependem de reescrever a equação correspondente como

\[ -2\left( \log\Big(L(\widetilde{\beta}_0,\beta_1 \, | \, y) \Big)-\log\Big(L(\widehat{\beta}_0,\widehat{\beta}_1 \, | \, y) \Big)\right) - \chi^2_{1,0.95} < 0\cdot \]

O primeiro método calcula sequencialmente o lado esquerdo da equação para um intervalo de valores possíveis de \(\beta_1\). Os valores menores e maiores de \(\beta_1\) que satisfazem a desigualdade correspondem aos limites de intervalo inferior e superior para \(\beta_1\). O segundo método substitui o sinal de menor que na equação por um sinal de igualdade e resolve iterativamente para valores de \(\beta_1\) usando a função uniroot(). Depois que os limites inferior e superior são encontrados por qualquer método, a função exp() é usada para encontrar os limites de intervalo para a razão de chances.


2.2.4 Probabilidade de sucesso


Uma vez que um modelo de regressão logística é estimado, muitas vezes é interessante estimar a probabilidade de sucesso para um conjunto de valores de variáveis explicativas. Isso pode ser feito simplesmente substituindo as estimativas dos parâmetros no modelo: \[ \widehat{\pi}=\dfrac{\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x_1+\cdots +\widehat{\beta}_p x_p\right)}{1+\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x_1+\cdots +\widehat{\beta}_p x_p\right)}\cdot \]

Como \(\widehat{\pi}\) é uma estatística, ela varia de amostra para amostra. Portanto, precisamos encontrar um intervalo de confiança para fazer inferências com um determinado nível de confiança. Os intervalos Wald e LR perfilados serão discutidos nesta seção.

Para ajudar a explicar como o intervalo de Wald é calculado, considere novamente um modelo de regressão logística com apenas uma variável explicativa. A probabilidade estimada de sucesso é então \[ \widehat{\pi} =\dfrac{\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x\right)}{1+\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x\right)}\cdot \] A distribuição normal é uma melhor aproximação para a distribuição de \(\widehat{\beta}_0\) e \(\widehat{\beta}_1\) do que para \(\widehat{\pi}\), então procedemos de maneira semelhante a encontrar o intervalo de Wald para a razão de chances.

Primeiro encontramos um intervalo de confiança para o preditor linear, \(logit(\pi ) =\beta_0 + \beta_1 x\), e então transformamos as extremidades desse intervalo em um intervalo usando a transformação \(\exp(\cdot)/\big(1 + \exp(\cdot )\big)\). O \((1-\alpha)100\)% intervalo de confiança de Wald para \(\beta_0 +\beta_1 x\) é \[ \widehat{\beta}_0+\widehat{\beta}_1 x \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x\big)}\cdot \]

A variância estimada é encontrada por \[ \widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x\big)=\widehat{\mbox{Var}}\big(\widehat{\beta}_0\big)+ x^2 \widehat{\mbox{Var}}\big(\widehat{\beta}_1 \big)+2x \widehat{\mbox{Cov}}\big(\widehat{\beta}_0,\widehat{\beta}_1 \big), \] onde cada variância e o termo de covariância estão disponíveis a partir da matriz de covariância estimada das estimativas de parâmetros. Esta é outra aplicação do seguinte resultado: \[ \mbox{Var}(a U+bV ) = a^2\mbox{Var}(U)+b^2\mbox{Var}(V )+2ab\mbox{Cov}(U,V ), \] onde \(U\) e \(V\) são variáveis aleatórias e \(a\) e \(b\) são constantes. Usando os limites de intervalo para \(\beta_0+\beta_1 x\), o intervalo de confiança \((1-\alpha)100\)% Wald para \(\pi\) é \[ \dfrac{\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x\big)}\right)}{1+\exp\left(\widehat{\beta}_0+\widehat{\beta}_1 x \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x\big)} \right)}\cdot \] Observe que o limite inferior (superior) para \(\pi\) usa o sinal de menos (mais) na parte \(\pm\) da equação acima.

Quando existem \(p\) variáveis explicativas no modelo, o intervalo de Wald para \(\pi\) é encontrado da mesma maneira. O intervalo é \[ \dfrac{\exp\left( \widehat{\beta}_0+\widehat{\beta}_1 x_1+ \cdots+\widehat{\beta}_p x_p \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x_+\cdots+\widehat{\beta}_p x_p\big)}\right)}{1+\exp\left(\widehat{\beta}_0+\widehat{\beta}_1 x_1 +\cdots+\widehat{\beta}_p x_p \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x_1+\cdots+\widehat{\beta}_p x_p\big)} \right)}, \] onde a expressão de variância é encontrada de maneira semelhante ao caso de variável única: \[ \widehat{\mbox{Var}}\big(\widehat{\beta}_0+\widehat{\beta}_1 x_1+\cdots+\widehat{\beta}_p x_p\big)=\sum_{i=0}^p x_i^2 \widehat{\mbox{Var}}\big(\widehat{\beta}_i\big) +2\sum_{i=1}^{p-1}\sum_{j=i+1}^p x_i x_j \widehat{\mbox{Cov}}\big(\widehat{\beta}_i,\widehat{\beta}_j\big), \] com \(x_0=1\). Embora essa expressão de variancia possa ser longa, ela é calculada automaticamente pela função predict() conforme ilustrado no próximo exemplo.

Os intervalos LR perfilada também são calculados primeiro para \(logit(\pi )\) e depois transformados em intervalos para \(\pi\). A mesma metodologia da Seção 2.2.3 é usada para calcular o primeiro intervalo, mas agora o cálculo é mais difícil devido ao maior número de parâmetros envolvidos. Por exemplo, em um modelo com uma variável explicativa, \(logit(\pi ) =\beta_0 +\beta_1 x\) é uma combinação linear de \(\beta_0\) e \(\beta_1\). O numerador de \(-2\log(\Lambda )\) envolve maximizar a função de verossimilhança com uma restrição para esta combinação linear, que é mais difícil do que restringir apenas um parâmetro.

O problema se torna ainda mais complicado quando há várias variáveis explicativas e, em alguns casos, os procedimentos numéricos iterativos podem demorar muito para serem executados ou até mesmo não convergir. Supondo que a computação seja bem-sucedida e um intervalo para \(\beta_0 +\beta_1 x\) seja encontrado, realizamos a transformação \(\exp(\cdot )/\big(1 + \exp(\cdot )\big)\) para obter um intervalo para \(\pi\).

Usamos o pacote mcprofile, um pacote de contribuição de usuários que não está na instalação padrão do R, para calcular os intervalos de verossimilhança perfilada para \(\pi\) e para a maioria das outras combinações lineares de parâmetros que serão discutidas aqui. Usamos este pacote porque é o pacote mais geral disponível para calcular esses tipos de intervalos. No entanto, versões anteriores produziram resultados questionáveis às vezes. A versão 1.0-1 do pacote, a versão mais atual no momento, foi usada para todos os cálculos nestas notas e esses resultados questionáveis não ocorrem mais. Ainda assim, sugerimos o seguinte curso de ação ao calcular os intervalos de confiança LR perfilada com mcprofile:

Nossas declarações aqui sobre o mcprofile não pretendem alarmar os leitores sobre cálculos realizados por pacotes de contribuição de usuários em geral. Em vez disso, há um grande número de pacotes R contribuídos por usuários, e muitos deles podem realizar cálculos que nenhum outro software pode; no entanto, esses pacotes de contribuição precisam ser usados de maneira criteriosa.

Se possível, as comparações iniciais devem ser feitas entre os cálculos realizados por um pacote de contribuição daqueles computados de alguma outra forma confiável para garantir que eles estejam de acordo. Além disso, como o R é de código aberto, o código para todas as funções está disponível para os usuários examinarem, se necessário. Implementações linha por linha de código dentro de uma função, embora possivelmente árduas, podem fornecer as garantias necessárias de que o código funciona conforme desejado. Nenhum outro software estatístico amplamente utilizado oferece essa oportunidade de verificação.


Exemplo 2.9: Placekicking.


Para este exemplo, usamos novamente o objeto mod.fit do modelo de regressão logística que usa apenas distance como variável explicativa. Abaixo estão três maneiras pelas quais a probabilidade estimada de sucesso pode ser calculada para uma distância de 20 jardas:

linear.pred <- mod.fit$coefficients [1] + mod.fit$coefficients [2] * 20 
linear.pred
## (Intercept) 
##    3.511547
as.numeric (exp ( linear.pred )/(1 + exp( linear.pred )))
## [1] 0.9710145
predict.data <- data.frame ( distance = 20)
predict ( object = mod.fit , newdata = predict.data , type = "link")
##        1 
## 3.511547
predict ( object = mod.fit , newdata = predict.data , type = "response")
##         1 
## 0.9710145
head ( placekick$distance == 20)
## [1] FALSE FALSE  TRUE FALSE  TRUE FALSE
mod.fit$fitted.values [3] # 3rd obs. has distance = 20
##         3 
## 0.9710145


A primeira maneira calcula diretamente o preditor linear como \[ \widehat{\beta}_0 + \widehat{\beta}_1 x = 5.8121 - 0.1150\times 20 = 3.5112 \] resultando em \[ \widehat{\pi} = \dfrac{\exp(3.5112)}{1 + \exp(3.5112)} = 0.9710\cdot \]

Observe que o título “(Intercept)” é um rótulo que sobrou de mod.fit$coeficientes[1], que pode ser removido usando as.numeric().

A segunda maneira de calcular \(\widehat{\pi}\) é usar a função predict() e recomendamos usar essa função na maioria das vezes para esses tipos de cálculos. A função é uma função genérica que realmente usa a função predict.glm() para realizar cálculos. Para usar o predict(), um arquivo de dados deve conter os valores da variável explicativa nos quais as estimativas de \(\pi\) são desejadas. Este arquivo de dados é incluído então no argumento newdata de predict(). Além disso, o argumento ,strong>object especifica onde as informações de ajuste do modelo de glm() estão localizadas, e o valor do argumento type = “response” instrui R a estimar \(\pi\). Alternativamente, o valor do argumento type = “link” instrui R a estimar \(\beta_0 + \beta_1 x\).

Se nenhum arquivo de dados for fornecido no argumento newdata de predict(), \(\pi\) será estimado para cada observação no arquivo de dados. O resultado é exatamente o mesmo dado por mod.fit$fitted.values. Assim, uma terceira maneira de estimar \(\pi\) é observar os valores de mod.fit$fitted.values correspondentes à distância de 20 jardas. É claro que as estimativas de \(\pi\) podem ser de interesse para valores de \(x\) que não estão no conjunto de dados, mas para nosso caso atual a terceira observação corresponde a \(x = 20\).

Para calcular um intervalo de confiança Wald para \(\pi\), a maneira mais fácil é usar a função predict() para calcular \(\widehat{\beta}_0+\widehat{\beta}_1 x\) e \(\widehat{\mbox{Var}}(\widehat{\beta}_0+\widehat{\beta}_1 x)\) primeiro. Em seguida, calculamos o intervalo de confiança para usar o código apropriado para as equações correspondentes:

alpha <- 0.05
linear.pred <- predict ( object = mod.fit , 
                         newdata = predict.data , type = "link", se = TRUE )
linear.pred
## $fit
##        1 
## 3.511547 
## 
## $se.fit
## [1] 0.1732707
## 
## $residual.scale
## [1] 1
pi.hat <- exp( linear.pred$fit ) / (1 + exp ( linear.pred$fit ))
CI.lin.pred <- linear.pred$fit + qnorm (p = c( alpha /2, 1- alpha /2) ) * linear.pred$se
CI.pi <- exp(CI.lin.pred ) /(1+ exp (CI.lin.pred ))
data.frame ( predict.data , pi.hat , lower = CI.pi [1] , upper = CI.pi [2])
##   distance    pi.hat     lower     upper
## 1       20 0.9710145 0.9597647 0.9791871


Usamos o valor do argumento se = TRUE dentro de predict() para encontrar \(\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_0+\widehat{\beta}_1 x)}\), o “erro padrão” de \(\widehat{\beta}_0+\widehat{\beta}_1 x\). Observe que o uso desse argumento adicional faz com que o objeto resultante seja uma lista de vetores componentes, enquanto um único vetor foi produzido quando os erros padrão não foram solicitados.

O intervalo de confiança de Wald de 95% para é \(0.9598 < \pi < 0.9792\); assim, a probabilidade de sucesso para o placekick é bastante alta a uma distância de 20 jardas. No programa R correspondente a este exemplo, também demonstramos como usar o predict() com mais de uma distância e como realizar os cálculos codificando a equação diretamente no R.

Para encontrar um intervalo LR perfilada para uma distância de 20 jardas, usamos o código abaixo:

library ( package = mcprofile )
K <- matrix ( data = c(1, 20) , nrow = 1, ncol = 2)
K
##      [,1] [,2]
## [1,]    1   20
# Calculate -2 log( Lambda )
linear.combo <- mcprofile ( object = mod.fit , CM = K)
# CI for beta_0 + beta_1 * x
ci.logit.profile <- confint ( object = linear.combo , level = 0.95)
ci.logit.profile
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   single-step 
## 
##    Estimate lower upper
## C1     3.51  3.19  3.87
names (ci.logit.profile )
## [1] "estimate"    "confint"     "CM"          "quant"       "alternative"
## [6] "level"       "adjust"
exp(ci.logit.profile$confint )/(1 + exp (ci.logit.profile$confint ))
##       lower    upper
## 1 0.9603165 0.979504


Após a chamada inicial para a função library() do pacote mcprofile, criamos uma matriz \(K\) que contém coeficientes para a combinação linear de interesse. Neste caso, queremos primeiro encontrar um intervalo de confiança para \(1\times \beta_0+20\times \beta_1\). A função mcprofile() calcula \(-2\log(\Lambda)\) para um grande número de valores possíveis da combinação linear, onde o argumento CM é curto para “matriz de contraste”. A função confint() encontra o intervalo LR perfilada de 95% como \(3.19 < \beta_0+20 \, \beta_1 < 3.87\). Usamos a transformação \(\exp(\cdot )/\big(1+\exp(\cdot )\big)\) para encontrar o intervalo para \(\pi\) como \(0.9603 < \pi< 0.9795\), que é muito semelhante ao intervalo de Wald devido ao grande tamanho da amostra.

Incluímos exemplos adicionais de mcprofile em nosso programa correspondente que fazem o seguinte:


Para um modelo linear com apenas uma variável explicativa \(x\), é muito comum examinar um gráfico de dispersão dos dados com o modelo de regressão linear estimado plotado sobre ele. O objetivo desse tipo de gráfico é obter uma avaliação visual das tendências nos dados e avaliar o ajuste do modelo. Como a variável de resposta na regressão logística é binária, construir um gráfico como esse não é muito informativo porque todos os pontos plotados estariam em \(y\) = 0 ou 1 no eixo \(y\). No entanto, podemos traçar a proporção observada de sucessos em cada \(x\) único para obter uma compreensão geral de quão bem o modelo se ajusta aos dados. Espera-se que o modelo estimado passe “próximo” a essas proporções, desde que cada uma seja baseada em um número suficientemente grande de observações. Uma chave para a avaliação com esses tipos de gráficos, então, é conhecer o número de observações que existem em cada \(x\), para que se possa realmente distinguir entre um ajuste ruim e bom. No caso extremo de um \(x\) contínuo, o número de observações seria 1 em cada \(x\) único e esse gráfico geralmente não seria muito útil.


Exemplo 2.10: Placekicking.


Em um exemplo anterior da Seção 2.2.1, encontramos o número de sucessos do número total de tentativas para cada distância. Usando o arquivo de dados w.n resultante, usamos as funções plot() e curve() aqui para plotar a proporção observada de sucessos em cada distância e, em seguida, sobrepor o modelo de regressão logística estimado:

head (w.n)
##   distance success trials proportion
## 1       18       2      3     0.6667
## 2       19       7      7     1.0000
## 3       20     776    789     0.9835
## 4       21      19     20     0.9500
## 5       22      12     14     0.8571
## 6       23      26     27     0.9630
plot (x = w$distance , y = w$good /n$good , xlab = "Distance ( yards )", 
      ylab = " Estimated probability ", panel.first = grid (col = "gray", lty = "dotted"))
curve ( expr = predict ( object = mod.fit , newdata = data.frame ( distance = x), 
                         type = "response") , col = "red", add = TRUE , xlim = c(18 , 66))
grid()


A função curve() avalia a função predict() em 101 distâncias dentro dos limites do eixo \(x\) especificando newdata = data.frame(distance = x) para o argumento expr. O gráfico resultante é dado na figura acima.

Também gostaríamos de traçar intervalos de confiança para \(\pi\) de Wald de 95% para cada distância possível. Essas bandas são adicionadas ao gráfico usando mais duas chamadas para curve():

ci.pi <- function ( newdata , mod.fit.obj , alpha ){
  linear.pred <- predict ( object = mod.fit.obj , newdata = newdata , type = "link", se = TRUE )
  CI.lin.pred.lower <- linear.pred$fit - qnorm (p = 1- alpha /2)* linear.pred$se 
  CI.lin.pred.upper <- linear.pred$fit + qnorm (p = 1- alpha /2)* linear.pred$se
  CI.pi.lower <- exp(CI.lin.pred.lower ) / (1 + exp (CI.lin.pred.lower )) 
  CI.pi.upper <- exp(CI.lin.pred.upper ) / (1 + exp (CI.lin.pred.upper ))
  list ( lower = CI.pi.lower , upper = CI.pi.upper )
}
# Test case
ci.pi( newdata = data.frame ( distance = 20) , mod.fit.obj = mod.fit , alpha = 0.05)
## $lower
##         1 
## 0.9597647 
## 
## $upper
##         1 
## 0.9791871


Como o argumento expr de curve() requer um nome de função como seu valor, escrevemos uma nova função chamada ci.pi() que realiza nossos cálculos de limite de confiança. Essa nova função contém essencialmente o mesmo código usado anteriormente para calcular um intervalo de confiança Wald para \(\pi\), e armazena os limites inferior e superior como componentes separados, para que cada um possa ser plotado com uma chamada separada para curve().

A função ci.pi() pode ser usada de forma mais geral para calcular intervalos de confiança de Wald para sempre que forem necessários. No final do código, usamos a função legend() para esclarecer o que cada linha plotada representa. Os argumentos \(x\) e \(y\) fornecidos em legend() instruem o R a colocar a legenda no local no qual o usuário indica.

# Plot C.I. bands
plot (x = w$distance , y = w$good /n$good , xlab = "Distance ( yards )", pch =19 ,
      ylab = " Estimated probability ", panel.first = grid (col = "gray", lty = "dotted"))
curve ( expr = predict ( object = mod.fit , newdata = data.frame ( distance = x), 
                         type = "response") , col = "red", add = TRUE , xlim = c(18 , 66))
curve ( expr = ci.pi( newdata = data.frame ( distance = x), mod.fit.obj = mod.fit, 
        alpha = 0.05)$lower, col = "blue", lty = "dotdash", add = TRUE , xlim = c(18 , 66))
curve ( expr = ci.pi( newdata = data.frame ( distance = x), mod.fit.obj = mod.fit, 
        alpha = 0.05)$upper, col = "blue", lty = "dotdash", add = TRUE , xlim = c(18 , 66))
# Legend
legend ( 20, 0.4 , legend = c("Logistic regression model", "95% individual C.I .") , 
         lty = c("solid", "dotdash") , col = c("red", "blue") , bty = "n")
grid()


Para avaliar informalmente quão bem nosso modelo se ajusta aos dados, pode-se ficar tentado a comparar quão próximas as proporções apresentadas na figura acima estão do modelo estimado. Há uma série de proporções distantes do modelo estimado; por exemplo, examine a proporção mostrada a uma distância de 18 jardas. No entanto, existem apenas 3 observações a 18 jardas e aconteceu que 1 das 3 é perdida, apesar da alta probabilidade de sucesso. Com resultados binários e um pequeno número de observações para um determinado valor de \(x\), isso pode acontecer e não indica necessariamente um ajuste ruim para o modelo.

Observe que a curva do modelo estimado na figura geralmente passa pelo meio das proporções plotadas. No entanto, existem várias proporções que estão longe do modelo estimado, por exemplo, em 18, 56 e 59 jardas. Essas são distâncias nas quais muito poucas tentativas foram registradas – três em 18 e uma em 56 e 59 – e, portanto, as proporções da amostra podem variar consideravelmente. Não se deve esperar que a curva do modelo estimado chegue tão perto desses pontos quanto deveria dos pontos baseados em um grande número de tentativas.

Uma maneira melhor de abordar um gráfico como a figura acima é incorporar informações sobre o número de tentativas \(n_j\) em cada ponto de plotagem. Isso pode ser feito com um gráfico de bolhas, onde cada ponto de plotagem é proporcional em tamanho ao número de tentativas. A função symbols() produz esses gráficos em R:

symbols (x = w$distance , y = w$good /n$good , circles = sqrt ( n$good ), inches = 0.5 , pch = 19,
         xlab = " Distance ( yards )", ylab = " Estimated probability ", 
         panel.first = grid (col = "gray", lty = "dotted" ))
curve ( expr = predict ( object = mod.fit , newdata = data.frame ( distance = x), 
                         type = "response") , col = "red", add = TRUE , xlim = c(18 , 66))
curve ( expr = ci.pi( newdata = data.frame ( distance = x), mod.fit.obj = mod.fit, 
        alpha = 0.05)$lower, col = "blue", lty = "dotdash", add = TRUE , xlim = c(18 , 66))
curve ( expr = ci.pi( newdata = data.frame ( distance = x), mod.fit.obj = mod.fit, 
        alpha = 0.05)$upper, col = "blue", lty = "dotdash", add = TRUE , xlim = c(18 , 66))
# Legend
legend ( 20, 0.4 , legend = c("Logistic regression model ", "95% individual C.I .") , 
         lty = c("solid", "dotdash") , col = c("red", "blue") , bty = "n")
grid()


Os argumentos \(x\) e \(y\) dentro de symbols() funcionam da mesma forma que em plot(). O argumento circles controla o tamanho do ponto de plotagem onde o raio máximo é especificado em polegadas pelo argumento inches. Para este gráfico, usamos o número de observações em cada distância única para corresponder ao tamanho do ponto de plotagem. A função sqrt() é usada apenas para diminuir as diferenças absolutas entre os valores em n$good, e isso não precisa ser usado em geral.

Existem 789 observações a uma distância de 20 jardas. O próximo maior número de observações é 30, que está a uma distância de 32 jardas. Devido à grande diferença absoluta no número de observações, isso faz com que o ponto de plotagem a 20 jardas seja muito grande e os outros pontos de plotagem sejam muito pequenos. Para diminuir essa diferença absoluta, usamos uma transformação em n$good.

A figura acima fornece o gráfico de bolhas final com estimativas e faixas de intervalo de confiança de Wald incluídas. Podemos ver que, de fato, os placekicks de 18 jardas não devem causar preocupação devido ao pequeno tamanho do ponto de plotagem. Existem pontos semelhantes nas maiores distâncias, onde muitos destes incluem apenas uma observação. Como esperado, o Placekicks de 20 jardas têm o maior ponto de plotagem e parece estar centralizado próximo ao modelo estimado.

Existem alguns pontos de plotagem grandes, como 32 jardas, \(\widehat{\pi}= 0.89\), proporção observada de 23/30 = 0.77 e 51 jardas \(\widehat{\pi} = 0.49\), proporção observada de 11/15 = 0.73, para o qual o modelo pode não se encaixar bem. Como avaliar mais formalmente essas observações e outras será um assunto importante do Capítulo 5 quando examinamos medidas diagnósticas de modelo.

Fornecemos código adicional no programa correspondente que fornece o mesmo gráfico mostrado na figura acima, mas com bandas de intervalo de confiança LR perfilada em vez das bandas de intervalo de confiança de Wald. Esses dois tipos de bandas são indistinguíveis nessas parcelas devido ao grande tamanho da amostra. Uma pequena vantagem de usar bandas de Wald aqui é que elas são calculadas muito mais rapidamente do que aquelas de intervalos LR perfilada.

# profile LR ratio version of the plot
# Easiest way - add the bands to the previous plot
#  Notice the bands are practically the same as those from the Wald intervals
distances<-18:66
K<-cbind(1, distances)
class(K)  # matrix
## [1] "matrix" "array"
head(K)
##        distances
## [1,] 1        18
## [2,] 1        19
## [3,] 1        20
## [4,] 1        21
## [5,] 1        22
## [6,] 1        23
linear.combo<-mcprofile(object = mod.fit, CM = K)  # Calculate -2log(Lambda)
# CI for beta_0 + beta_1 * x
ci.logit.profile<-confint(object = linear.combo, level = 0.95, adjust = "none")  
#head(ci.logit.profile)
profile.lr.int<-exp(ci.logit.profile$confint)/(1 + exp(ci.logit.profile$confint))
plot (x = w$distance , y = w$good /n$good , xlab = "Distance ( yards )", pch =19 ,
      ylab = " Estimated probability ", panel.first = grid (col = "gray", lty = "dotted"))
# Add bands
lines(x = distances, y = profile.lr.int[,1], col = "green")
lines(x = distances, y = profile.lr.int[,2], col = "green")

# Another way to do the sample type of plot using a new ci.pi()-like function
ci.pi2<-function(x, mod.fit.obj, alpha){
  K<-cbind(1, x)
  linear.combo<-mcprofile(object = mod.fit, CM = K)
  ci.logit.profile<-confint(object = linear.combo, level = 1 - alpha, adjust = "none")
  profile.lr.int<-exp(ci.logit.profile$confint)/(1 + exp(ci.logit.profile$confint))
  list(lower = profile.lr.int[,1], upper = profile.lr.int[,2])
}
# Test case
ci.pi2(x = 20, mod.fit.obj = mod.fit, alpha = 0.05)
## $lower
## [1] 0.9603165
## 
## $upper
## [1] 0.979504
symbols(x = w$distance, y = w$good/n$good, circles = sqrt(n$good), inches = 0.5, 
        xlab = "Distance (yards)", ylab = "Estimated probability", 
        panel.first = grid(col = "gray", lty = "dotted"))
# Put estimated logistic regression model on the plot
curve(expr = predict(object = mod.fit, newdata = data.frame(distance = x), 
                     type = "response"), col = "red", add = TRUE, xlim = c(18, 66))
# Plot C.I. bands
curve(expr = ci.pi2(x = x, mod.fit.obj = mod.fit, alpha = 0.05)$lower, 
      col = "blue", lty = "dotdash", add = TRUE, xlim = c(18, 66))
curve(expr = ci.pi2(x = x, mod.fit.obj = mod.fit, alpha = 0.05)$upper, 
      col = "blue", lty = "dotdash", add = TRUE, xlim = c(18, 66))
# Legend
legend( 20, 0.4, legend = c("Logistic regression model", "95% individual C.I."), 
        lty = c("solid", "dotdash"), col = c("red", "blue"), bty = "n")



2.2.5 Interações e transformações para variáveis explicativas


As interações e transformações ampliam a variedade de formas que podem ser propostas para se relacionar \(\pi\) com as variáveis explicativas. Os termos são criados e adicionados ao preditor linear no modelo de regressão logística exatamente da mesma maneira que são usados na regressão linear. Os termos mais comuns e interpretáveis a serem adicionados são interações em pares (two-way) e termos quadráticos. Embora outros tipos de termos também possam ser adicionados, por exemplo, interações de três termos (three-way) e transformações cúbicas, deixamos de explorar essas alternativas como exercícios.

Interações entre variáveis explicativas são necessárias quando o efeito de uma variável explicativa na probabilidade de sucesso depende do valor de uma segunda variável explicativa. Existem algumas maneiras de incluir essas interações em um argumento de fórmula para glm(). Para descrevê-los, suponha que existam duas variáveis explicativas denominadas x1 e x2 em um arquivo de dados representando \(x_1\) e \(x_2\), e o objetivo é estimar o modelo \[ logit(\pi ) = \beta_0+ \beta_1 x_1+\beta_2 x_2+ \beta_3 x_1 x_2\cdot \]

Observe que adotamos a convenção de que as próprias variáveis explicativas, seus efeitos principais, são sempre incluídas no modelo quando a interação é incluída. Abaixo estão maneiras equivalentes de codificar o modelo: \[ \begin{array}{rccccl} 1. & \mbox{formula} & = & \mbox{y} & \sim & \mbox{x1} + \mbox{x2} + \mbox{x1:x2} \\ 2. & \mbox{formula} & = & \mbox{y} & \sim & \mbox{x1*x2} \\ 3. & \mbox{formula} & = & \mbox{y} & \sim & \mbox{(x1+x2)^2} \end{array} \]

No primeiro argumento formula, o símbolo de dois pontos denota uma interação entre x1 e x2. No segundo argumento formula, o asterisco cria automaticamente todos os principais efeitos e interações entre as variáveis conectadas. Para o terceiro argumento formula, o \(\mbox{^2}\) cria todas as combinações de até duas variáveis dentre as variáveis listadas entre parênteses. Se houver uma terceira variável x3, então \(\mbox{(x1 + x2 + x3)^2}\) é equivalente a \(\mbox{x1*x2}\) + \(\mbox{x1*x3}\) + \(\mbox{x2*x3}\). Além disso, \(\mbox{(x1 + x2 + x3)^3}\) é equivalente a \(\mbox{x1*x2*x3}\).

Polinômios quadráticos e de ordem superior são necessários quando a relação entre uma variável explicativa e \(logit(\pi )\) não é linear. Para incluir uma transformação quadrática no argumento formula, a função identidade I() precisa ser aplicada à variável explicativa transformada. I() é uma função de “identidade” que instrui R a interpretar seu argumento como faria normalmente. Esta função é usada para fornecer valores de argumento para outras funções que podem ter uma interpretação especial para um formato de símbolo ou objeto.

Por exemplo, o modelo \(logit(\pi ) = \beta_0 + \beta_1 x_1 + \beta_2 x^2_1\) é codificado como

formula = y ~ x1 + I(x1^2)


A função I() é necessária porque o símbolo ^ tem um significado especial em argumentos formula como acabamos de ver no parágrafo anterior.

A inclusão de termos de interação e/ou transformações faz com que as razões de chances dependam do valor numérico de uma variável explicativa. Por exemplo, considere o modelo \(logit(\pi ) = \beta_0 + \beta_1 x_1 +\beta_2 x_2 + \beta_3 x_1x_2\). A razão de chances associada à mudança de \(x_2\) por \(c\) unidades enquanto mantendo \(x_1\) constante é \[ \begin{array}{rcl} OR & = & \dfrac{Odds_{x_2+c}}{Odds_{x_2}} = \dfrac{\exp\big(\beta_0+\beta_1 x_1+\beta_2 (x_2+c) +\beta_3 x_1 (x_2+c)\big)}{\exp\big(\beta_0+\beta_1 x_1+\beta_2 x_2 +\beta_3 x_1 x_2\big)} \\ & = & \exp\big(c\, \beta_2 + c\, \beta_3 x_1\big) \, = \, \exp\big( c \,(\beta_2+\beta_3 x_1)\big)\cdot \end{array} \]

Assim, o aumento ou diminuição das chances de sucesso para uma mudança de unidade \(c\) em \(x_2\) depende do nível de \(x_1\), que decorre da definição de uma interação. Como outro exemplo, considere o modelo \(logit(\pi ) =\beta_0 +\beta_1 x_1 + \beta_2 x^2_1\). A razão de chances para uma mudança de unidade \(c\) em \(x_1\) é \[ \begin{array}{rcl} OR & = & \dfrac{Odds_{x_1+c}}{Odds_{x_1}} = \dfrac{\exp\big(\beta_0+\beta_1 (x_1+c)+\beta_2 (x_2+c)^2\big)}{\exp\big(\beta_0+\beta_1 x_1+\beta_2 x_1^2\big)} \\ & = & \exp\big(c\, \beta_1 + 2c\, x_1\beta_2 +c^2\beta_2 \big) \, = \, \exp\big( c\, \beta_1+c\, \beta_2(2x_1+c)\big), \end{array} \] que depende do valor de \(x_1\).

Os intervalos de confiança para essas razões de chances podem ser mais complicados de calcular, porque geralmente são baseados em combinações lineares de parâmetros de regressão, e a função genérica confint() não pode ser usada para a tarefa. Para o exemplo de interação, onde \(OR=\exp\big( c\, (\beta_2+\beta_3 x_1)\big)\), o \((1-\alpha)100\)% intervalo de confiança Wald é \[ \exp\left( c\, (\widehat{\beta}_2+\widehat{\beta}_3 x_1) \pm c\, Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_2+\widehat{\beta}_3 x_1)}\right) \] com \[ \widehat{\mbox{Var}}(\widehat{\beta}_2+\widehat{\beta}_3 x_1)=\widehat{\mbox{Var}}(\widehat{\beta}_2)+x_1^2 \, \widehat{\mbox{Var}}(\widehat{\beta}_3)+2x_1 \, \widehat{\mbox{Cov}}(\widehat{\beta}_2,\widehat{\beta}_3)\cdot \]

Outros intervalos de Wald podem ser calculados de maneira semelhante. Geralmente, a única parte desafiadora é encontrar a variância para a combinação linear dos estimadores de parâmetros, então fornecemos código para fazer isso no próximo exemplo.

Os intervalos LR perfilados podem ser calculados para essas razões de chances usando o pacote mcprofile. Semelhante ao intervalo de Wald, a parte mais desafiadora é especificar a combinação linear correta dos estimadores de parâmetros. A função wald() deste pacote também fornece uma maneira um tanto automatizada de calcular os intervalos de confiança de Wald. Discutiremos os detalhes do cálculo no próximo exemplo.


Exemplo 2.11: Placekicking.


Geralmente, pode-se conjecturar que quanto mais tempo uma bola de futebol está no ar, mais ela é suscetível ao efeito do vento. Além disso, chutes mais longos geralmente têm mais tempo no ar do que chutes mais curtos. Assim, seria interessante verificar se existe uma interação entre distance e wind.

A variável explicativa wind no conjunto de dados é uma variável binária para placekicks testados em condições de vento (1) versus condições sem vento (0), onde as condições de vento são definidas como um vento mais forte do que 15 milhas por hora no início em uma pista ao ar livre estádio.

Embora as informações de direção do vento estivessem disponíveis quando coletamos esses dados, a orientação do estádio não estava, tornando impossível contabilizar a direção do vento em relação aos placekicks. Mesmo que a orientação estivesse disponível, poderia não ter sido útil devido aos efeitos do vento em turbilhão que podem ocorrer dentro de um estádio. Além disso, observe que a velocidade do vento foi dicotomizada em condições ventosas versus não ventosas devido à velocidade do vento desconhecida dentro dos estádios abobadados, apenas um estádio abobadado registrou essa informação. Devido aos sistemas de ventilação desses estádios, não é correto supor que a velocidade do vento seja 0.

Em seguida, estimamos um modelo incluindo este termo de interação junto com os efeitos principais correspondentes:

mod.fit.H0 <- glm( formula = good ~ distance + wind , 
                   family = binomial ( link = logit ), data = placekick )
mod.fit.H1 <- glm( formula = good ~ distance + wind + distance:wind , 
                   family = binomial ( link = logit ), data = placekick )
summary (mod.fit.H1)
## 
## Call:
## glm(formula = good ~ distance + wind + distance:wind, family = binomial(link = logit), 
##     data = placekick)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7291   0.2465   0.2465   0.3791   1.8647  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    5.684181   0.335962  16.919   <2e-16 ***
## distance      -0.110253   0.008603 -12.816   <2e-16 ***
## wind           2.469975   1.662144   1.486   0.1373    
## distance:wind -0.083735   0.043301  -1.934   0.0531 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1013.43  on 1424  degrees of freedom
## Residual deviance:  767.42  on 1421  degrees of freedom
## AIC: 775.42
## 
## Number of Fisher Scoring iterations: 6
# LRT - Anova ( mod.fit.Ha) would also work here too.
anova (mod.fit.H0 , mod.fit.H1 , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: good ~ distance + wind
## Model 2: good ~ distance + wind + distance:wind
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1      1422     772.53                       
## 2      1421     767.42  1   5.1097  0.02379 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


O modelo de regressão logística estimado incluindo os principais efeitos e a interação é \[ logit(\widehat{\pi}) = 5.6842 - 0.1103 \, \mbox{distance} + 2.4700 \, \mbox{wind} -0.0837 \, \mbox{distance}\times \mbox{wind}\cdot \]

Para testar \(H_0 : \beta_3 = 0\) vs. \(H_1 : \beta_3 \neq 0\), o teste de Wald dá um \(p\)-valor de 0.0531 e o LRT dá um \(p\)-valor de 0.0238. Ambos os testes sugerem que há evidências marginais para apoiar uma interação wind e distance. O coeficiente negativo no efeito principal da distância indica que quando wind = 0, a probabilidade de sucesso diminui com o aumento da distância. O coeficiente negativo da interação indica que esse efeito é exacerbado em condições de vento.

A figura abaixo apresenta gráficos dos modelos de regressão logística estimados para os casos wind = 0 e 1. O gráfico da esquerda usa o modelo de regressão logística estimado sem a interação, trong>mod.fit.H0 contém as estimativas do modelo e o gráfico da direita usa o modelo de regressão logística estimado com a interação, mod.fit.H1 contém as estimativas do modelo.

O efeito da interação é evidente no gráfico da direita pela redução da probabilidade estimada de sucesso para os placekicks mais longos quando está ventando. Para criar esses gráficos, usamos a função curve() de maneira semelhante à que fizemos na Seção 2.2.4.

beta.hat <- mod.fit.H1$coefficients [2:4]
c <- 1
distance <- seq( from = 20, to = 60, by = 10)
# par( mfrow = c(1,2) )


OR.wind <- exp(c*( beta.hat [2] + beta.hat [3]* distance ))
cov.mat <- vcov ( mod.fit.H1) [2:4 ,2:4]
# Var(beta - hat_2 + distance *beta - hat_3 )
var.log.OR <- cov.mat [2 ,2] + distance^2* cov.mat [3 ,3] + 2* distance *cov.mat [2 ,3]
ci.log.OR.low = c*( beta.hat [2] + beta.hat [3]* distance ) - c* qnorm (p = 0.975) * sqrt (var.log.OR)
ci.log.OR.up = c*( beta.hat [2] + beta.hat [3]* distance ) + c* qnorm (p = 0.975) * sqrt (var.log.OR)
# Plot
par(mfrow = c(1,2))
curve(expr = predict(object = mod.fit.H0, newdata = data.frame(distance = x, wind = 0), 
      type = "response"), col = "red", lty = "solid", xlim = c(20,60), ylim = c(0,1), 
      ylab = "Estimated probability", main = "Without interaction", xlab = "Distance", 
      panel.first = grid(col = "gray", lty = "dotted"), cex.main = 0.9, lwd = 1)
curve(expr = predict(object = mod.fit.H0, newdata = data.frame(distance = x, wind = 1), 
      type = "response"), col = "blue", lty = "dotdash", lwd = 1, add = TRUE) 
legend(x = 20, y = 0.4, legend = c("Wind = 0", "Wind = 1"), lty = c("solid", "dotdash"), 
      col = c("red", "blue"), lwd = c(1,1), bty = "n")
curve(expr = predict(object = mod.fit.H1, newdata = data.frame(distance = x, wind = 0), 
      type = "response"), col = "red", lty = "solid", xlim = c(20,60), ylim = c(0,1), 
      ylab = "Estimated probability", main = "With interaction", xlab = "Distance", 
      panel.first = grid(col = "gray", lty = "dotted"), cex.main = 0.9, lwd = 1)
curve(expr = predict(object = mod.fit.H1, newdata = data.frame(distance = x, wind = 1), 
      type = "response"), col = "blue", lty = "dotdash", lwd = 1, add = TRUE)
legend(x = 20, y = 0.4, legend = c("Wind = 0", "Wind = 1"), lty = c("solid", "dotdash"), 
      col = c("red", "blue"), lwd = c(1,1), bty = "n")


Por causa da interação no modelo, precisamos interpretar o efeito do wind em níveis específicos de distance. Da mesma forma, precisamos interpretar o efeito da distance em níveis específicos de wind. Abaixo está o código R correspondente e a saída para obter estimativas de odds ratio e intervalos de confiança de Wald de 95%:

beta.hat <- mod.fit.H1$coefficients [2:4]
c <- 1
distance <- seq( from = 20, to = 60, by = 10)
OR.wind <- exp(c*( beta.hat [2] + beta.hat [3]* distance ))
cov.mat <- vcov ( mod.fit.H1) [2:4 ,2:4]
# Var(beta - hat_2 + distance *beta - hat_3 )
var.log.OR <- cov.mat [2 ,2] + distance^2* cov.mat [3 ,3] + 2* distance *cov.mat [2 ,3]
ci.log.OR.low <- c*( beta.hat [2] + beta.hat [3]* distance ) - c* qnorm (p = 0.975) * sqrt (var.log.OR)
ci.log.OR.up <- c*( beta.hat [2] + beta.hat [3]* distance ) + c* qnorm (p = 0.975) * sqrt (var.log.OR)
round ( data.frame ( distance = distance , OR.hat = 1/ OR.wind , 
                     OR.low = 1/ exp (ci.log.OR.up), OR.up = 1/ exp(ci.log.OR.low)) ,2)
##   distance OR.hat OR.low OR.up
## 1       20   0.45   0.09  2.34
## 2       30   1.04   0.40  2.71
## 3       40   2.41   1.14  5.08
## 4       50   5.57   1.54 20.06
## 5       60  12.86   1.67 99.13


Começamos extraindo \(\widehat{\beta}_1\), \(\widehat{\beta}_2\) e \(\widehat{\beta}_3\) de mod.fit$coeficientes. Embora isso não seja necessário, fazemos isso aqui para que o vetor resultante tenha os mesmos índices que os subscritos para as estimativas de parâmetros, ou seja, beta.hat[1] é \(\widehat{\beta}_1\).

As razões de chances estimadas estão dentro do objeto OR.wind, onde usamos \(c = 1\) porque wind é uma variável binária e examinamos distâncias de 20, 30, 40, 50 e 60 jardas. Usamos as expressões anteriores para calcular \(\widehat{\mbox{Var}}(\widehat{\beta}_2 + \widehat{\beta}_3 x)\) ou equivalentemente \(\widehat{\mbox{Var}}(\log(\widehat{OR}))\), que é como escolhemos o nome para este objeto e o intervalo de confiança para \(OR\).

Por exemplo, o intervalo de confiança de 95% para a razão de chances do wind a uma distância de 30 jardas é 0.40 < \(OR\) < 2.71. Como 1 está dentro do intervalo, não há evidências suficientes para indicar que condições de vento tenham efeito sobre o sucesso ou fracasso de um placekick. Por outro lado, o intervalo de confiança de 95% para a razão de chances do wind a uma distância de 50 jardas é 1.54 < \(OR\) < 20.06. Como 1 está fora do intervalo, há evidências suficientes de que ,strong<wind tem um efeito no sucesso ou fracasso de um chute a essa distância. Como originalmente conjecturado, quanto maior a distância, mais suscetível ao vento.

Os intervalos LR perfilada para as razões de chances são obtidos usando mcprofile() e, em seguida, confint():

K <- matrix ( data = c(0, 0, 1, 20, 
                       0, 0, 1, 30, 
                       0, 0, 1, 40, 
                       0, 0, 1, 50,
                       0, 0, 1, 60), 
              nrow = 5, ncol = 4, byrow = TRUE )
# K <- cbind (0, 0, 1, distance ) # A quicker way to form K
linear.combo <- mcprofile ( object = mod.fit.H1 , CM = K)
ci.log.OR <- confint ( object = linear.combo , level = 0.95 , adjust = "none")
data.frame ( distance , OR.hat = round (1/ exp(ci.log.OR$estimate ), 2) , 
             OR.low = round (1/ exp(ci.log.OR$confint$upper ), 2) , 
             OR.up = round (1/ exp(ci.log.OR$confint$lower ), 2))
##    distance Estimate OR.low  OR.up
## C1       20     0.45   0.06   1.83
## C2       30     1.04   0.34   2.42
## C3       40     2.41   1.14   5.15
## C4       50     5.57   1.68  22.79
## C5       60    12.86   2.01 130.85


Enquanto os intervalos LR perfilada são um pouco diferentes dos intervalos de Wald, as mesmas conclusões são alcançadas. Observe que os intervalos de confiança de Wald também podem ser calculados com a função wald():

save.wald <- wald ( object = linear.combo )
ci.log.OR.wald <- confint ( object = save.wald , level = 0.95 , adjust = "none")


Exponenciando os limites dos intervalos em ci.log.OR.wald e depois invertendo os intervalos, produz os mesmos intervalos de Wald calculados anteriormente.

As razões de chances estimadas e os intervalos de confiança correspondentes também são encontrados para distance, onde o wind é 0 ou 1. Para intervalos de Wald de 95%, as chances de um chute bem sucedido mudam em uma quantidade entre 2.54 a 3.56 vezes para uma diminuição de 10 jardas em distance quando não há vento, wind = 0, e por uma quantidade entre 3.03 a 15.98 para uma diminuição de 10 jardas em distance quando há vento, wind = 1. Para intervalos LR perfilada de 95%, os intervalos são de 2.55 a 3.58 (sem vento) e 3.40 a 18.79 (vento) para uma diminuição de 10 jardas em distance.

O código R correspondente está disponível no programa para este exemplo.

# A little quicker way to form K
distance<-seq(from = 20, to = 60, by = 10)
K<-cbind(0, 0, 1, distance)
K
##            distance
## [1,] 0 0 1       20
## [2,] 0 0 1       30
## [3,] 0 0 1       40
## [4,] 0 0 1       50
## [5,] 0 0 1       60
linear.combo<-mcprofile(object = mod.fit.H1, CM = K)
ci.log.OR<-confint(object = linear.combo, level = 0.95, adjust = "none")
# ci.log.OR
data.frame(distance, OR.hat = round(1/exp(ci.log.OR$estimate), 2),
    OR.low = round(1/exp(ci.log.OR$confint$upper), 2),
    OR.up = round(1/exp(ci.log.OR$confint$lower), 2))
##    distance Estimate OR.low  OR.up
## C1       20     0.45   0.06   1.83
## C2       30     1.04   0.34   2.42
## C3       40     2.41   1.14   5.15
## C4       50     5.57   1.68  22.79
## C5       60    12.86   2.01 130.85
# profile LR
linear.combo<-mcprofile(object = mod.fit.H1, CM = K)
    ci.log.OR<-confint(object = linear.combo, level = 0.95, adjust = "none")
    exp(ci.log.OR)
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##    Estimate   lower  upper
## C1   2.2150 0.54627 16.061
## C2   0.9588 0.41399  2.930
## C3   0.4150 0.19430  0.876
## C4   0.1796 0.04387  0.594
## C5   0.0778 0.00764  0.497



2.2.6 Variáveis explicativas categóricas


As variáveis explicativas categóricas são representadas em um modelo de regressão logística da mesma forma que em um modelo de regressão linear. Se uma variável tiver \(q\) níveis categóricos, as \(q-1\) variáveis indicadoras podem ser usadas para representá-la em um modelo. Por exemplo, vimos anteriormente que change no conjunto de dados placekicking era uma variável categórica com dois níveis. Foi representado em um modelo de regressão logística com uma variável indicadora: 0 para placekicks de troca sem chumbo e 1 para placekicks de troca de chumbo.

R trata as variáveis categóricas como um tipo de objeto de fator. Uma exceção pode ocorrer se a variável explicativa for codificada numericamente, como change. Discutiremos como lidar com essas situações mais adiante nesta seção.

Por padrão, R ordena os níveis usando uma ordenação numérica e, em seguida, alfabética, em que as letras minúsculas são ordenadas antes das maiúsculas. Para ver a ordenação de qualquer fator, a função levels() pode ser usada.

Essa ordenação é importante porque a maioria das funções em R a usa para construir variáveis indicadoras com o método de construção “definir primeiro nível como 0”. A construção “definir último nível como 0” é usada às vezes por outros softwares estatísticos, como o SAS. Existem outras formas de construção; consulte a Seção 6.3.5 para obter um exemplo.

Por exemplo, a tabela abaixo mostra a codificação R padrão para variáveis indicadoras de uma variável explicativa categórica com o \(q = 4\) níveis codificados como A, B, C e D. O primeiro nível A é o nível base onde todas as variáveis indicadoras são definidas como 0. Os níveis restantes são atribuídos a uma variável indicadora. Neste caso, o nível B é representado por \(x_1\) onde \(x_1 = 1\) para um nível observado de B e \(x_1 = 0\) caso contrário. Os níveis C e D são definidos de maneira semelhante para \(x_2\) e \(x_3\), respectivamente. A função relevel() em R pode ser usada para definir um novo nível base, se desejado, e isso será discutido mais detalhadamente no próximo exemplo. \[ \mbox{Variáveis indicadoras para uma variável explicativa categórica de 4 níveis.} \\ \begin{array}{cccc} \hline \mbox{Nível} & x_1 & x_2 & x_3 \\\hline \mbox{A} & 0 & 0 & 0 \\ \mbox{B} & 1 & 0 & 0 \\ \mbox{C} & 0 & 1 & 0 \\ \mbox{D} & 0 & 0 & 1 \\\hline \end{array} \]

Uma variável explicativa categórica é representada em um modelo de regressão logística usando todas as suas variáveis indicadoras. Por exemplo, um modelo de regressão logística com a variável explicativa categórica de 4 níveis é escrito como \[ logit(\pi) = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3\cdot \]

Uma representação um pouco menos formal substitui os nomes dos níveis para cada uma de suas variáveis indicadoras correspondentes; ou seja,

\[ logit(\pi) = \beta_0 + \beta_1\, \mbox{B} + \beta_2 \, \mbox{C} + \beta_3 \, \mbox{D}\cdot \]

Para testar a importância de uma variável explicativa categórica, todos os seus parâmetros de regressão correspondentes devem ser iguais a 0 na hipótese nula. Por exemplo, testaríamos todas as três variáveis indicadoras na equação acima simultaneamente com \[ H_0 : \beta_1 = \beta_2 = \beta_3 = 0 \qquad \mbox{vs.} \qquad H_1 : \mbox{Nem todas iguais a 0}; \] onde um LRT ou outro procedimento de teste apropriado é usado.


Exemplo 2.12: Controle do vírus da murcha do tomateiro.


Os vírus de plantas são frequentemente transmitidos por insetos. Isso ocorre quando os insetos se alimentam de plantas já infectadas com um vírus e posteriormente se tornam portadores do vírus. Quando esses insetos se alimentam de outras plantas, eles podem transmitir esse vírus para essas novas plantas.

Para entender melhor um vírus em particular, o Tomato Spotted Wilt Virus, e como controlar os tripes que o espalham, pesquisadores da Kansas State University realizaram um experimento em várias estufas. Dados cortesia dos Drs. James Nechols e David Margolies, Departamento de Entomologia, Kansas State University.

Cem plantas de tomate não infectadas foram colocar em cada estufa. Dentro de cada estufa, um dos dois métodos foi usado para introduzir o vírus nas plantas limpas:

Para examinar formas de controlar a propagação do vírus às plantas, os pesquisadores usaram um dos três métodos:

Entre as plantas que estavam originalmente limpas, o número de sintomas de infecção foi registrado após 8 semanas para cada estufa. Abaixo está uma parte dos dados em que cada linha do quadro ou arquivo de dados tomato representado uma estufa:

tomato <- read.csv( file = "http://leg.ufpr.br/~lucambio/CE073/20222S/TomatoVirus.csv")
head ( tomato )
##   Infest Control Plants Virus8
## 1      1       C    100     21
## 2      2       C    100     10
## 3      1       B    100     19
## 4      1       N    100     40
## 5      2       C    100     30
## 6      2       B    100     30


Da linha 1 do conjunto de dados, vemos que 21 de 100 plantas originalmente não infectadas em uma estufa apresentaram sintomas após 8 semanas e essas plantas tiveram o método de infestação #1 aplicado enquanto tentavam controlar os tripes com uma aplicação química.

As variáveis explicativas Control e Infest são de natureza categórica. R automaticamente trata a variável Control como um fator, com um nome de classe factor por causa de seus valores serem caracteres:

tomato$Control = as.factor( tomato$Control)
class ( tomato$Control )
## [1] "factor"
levels ( tomato$Control )
## [1] "B" "C" "N"
contrasts ( tomato$Control )
##   C N
## B 0 0
## C 1 0
## N 0 1


A função contrasts() mostra como R representaria Control com variáveis indicadoras em um modelo, semelhante ao que foi mostrado na tabela acima. Assim, o nível B de Control é o nível base e existem duas variáveis indicadoras que representam os níveis C e N de Control. Observe que Infest é codificado numericamente. Como a variável possui apenas dois níveis, ela já é representada por dois níveis numéricos 1 e 2, em vez de 0 e 1, portanto não há diferença prática entre usá-la como numérica ou como fator.

No entanto, se tivesse mais de dois níveis, teríamos que redefinir sua classe para factor usando a função factor(). Um uso adicional de factor() é alterar a ordem dos níveis, semelhante ao relevel(), mas com mais flexibilidade e alterar os nomes dos níveis. Exemplos são fornecidos no programa correspondente a este exemplo.

Abaixo está o código correspondente:

class ( tomato$Infest )
## [1] "integer"
levels ( tomato$Infest )
## NULL
class ( factor ( tomato$Infest ))
## [1] "factor"
levels ( factor ( tomato$Infest ))
## [1] "1" "2"
contrasts ( factor ( tomato$Infest ))
##   2
## 1 0
## 2 1
tomato$Infest <- factor ( tomato$Infest )
class ( tomato$Infest )
## [1] "factor"


A função factor() pode ser usada em toda a análise de maneira semelhante à acima. Por exemplo, podemos usá-lo no argumento de formula de glm() ao especificar a variável. Alternativamente, simplesmente substituímos a versão original do Infest no arquivo de dados tomat0 por uma nova versão do fator. Novamente, isso normalmente não seria necessário para uma variável explicativa de dois níveis, mas fazemos isso aqui para fins de demonstração.

Abaixo está a saída da estimativa do modelo de regressão logística com os fatores de Infest e control:

mod.fit <- glm( formula = Virus8 / Plants ~ Infest + Control , 
                family = binomial ( link = logit ), data = tomato , weights = Plants )
summary (mod.fit )
## 
## Call:
## glm(formula = Virus8/Plants ~ Infest + Control, family = binomial(link = logit), 
##     data = tomato, weights = Plants)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.288  -2.425  -1.467   1.828   8.379  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6652     0.1018  -6.533 6.45e-11 ***
## Infest2       0.2196     0.1091   2.013   0.0441 *  
## ControlC     -0.7933     0.1319  -6.014 1.81e-09 ***
## ControlN      0.5152     0.1313   3.923 8.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 278.69  on 15  degrees of freedom
## Residual deviance: 183.27  on 12  degrees of freedom
## AIC: 266.77
## 
## Number of Fisher Scoring iterations: 4


Como a variável de resposta é dada em uma forma binomial, usamos o argumento weights junto com a formulação sucesso/ensaios no argumento de formula, consulte a Seção 2.2.1. O modelo de regressão logística estimado é \[ logit(\widehat{\pi}) = -0.6652 + 0.2196\, \mbox{Infest2} -0.7933 \, \mbox{ControlC} + 0.5152 \, \mbox{ControlN}, \] onde usamos a notação do R de Infest2 para a variável indicadora que representa o nível 2 de Infest e ControlC e ControlN para as variáveis indicadoras que representam os níveis C e N de Control, respectivamente.

Com base no parâmetro positivo estimado para Infest2, estima-se que a probabilidade de apresentar sintomas seja maior em estufas onde o método de infestação #2 é usado. Isso é um pouco esperado porque os tripes já são portadores de vírus com esse método. Além disso, com base nos parâmetros estimados para ControlC e ControlN, a probabilidade estimada de apresentar sintomas é mais baixa para o método de controle químico e mais alta para quando nenhum método de controle é usado. Este modelo e interpretação assumem que não há interação entre os métodos de infestação e controle. A seguir, examinaremos como incluir essas interações e avaliaremos sua importância.



Interações


Para representar uma interação entre uma variável explicativa categórica e uma variável explicativa contínua em um modelo, todos os produtos aos pares dos termos variáveis são necessários. Por exemplo, considere a equação \[ logit(\pi) = \beta_0+\beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 \]

com um termo adicional de \(\beta_4 v\) adicionado a ela, onde \(v\) denota genericamente uma variável explicativa contínua.

A interação entre a variável explicativa categórica de 4 níveis e \(v\) é representada pela inclusão de três produtos pareados entre \(v\) e \(x_1\), \(x_2\) e \(x_3\) no modelo com parâmetros de regressão como coeficientes.

Para representar uma interação pareada entre duas variáveis explicativas categóricas, todos os produtos pareados entre os dois conjuntos de variáveis indicadoras são incluídos no modelo com parâmetros de regressão apropriados como coeficientes. Por exemplo, suponha que haja um fator \(X\) de 3 níveis e um fator \(Z\) de 3 níveis, representados por variáveis indicadoras apropriadas \(x_1\) e \(x_2\) para \(X\) e \(z_1\) e \(z_2\) para \(Z\). O modelo de regressão logística com a interação entre \(X\) e \(Z\) é \[ logit(\pi)=\beta_0 +\beta_1 x_1 + \beta_2 x_2 + \beta_3 z_1 + \beta_4 z_2 + \beta_5 x_1 z_1 + \beta_6 x_1 z_2 + \beta_7 x_2 z_1 + \beta_8 x_2 z_2\cdot \]

Para testar a interação, a hipótese nula é que todos os parâmetros de regressão para termos de interação são 0, ou seja, \[ H_0 : \beta_5 = \beta_6 = \beta_7 = \beta_8 = 0, \]

em nosso exemplo e a alternativa é que pelo menos um desses parâmetros não seja 0. Exemplos são dados a seguir.


Exemplo 2.13: Controle do vírus da murcha do tomateiro.


Na aplicação real, é importante considerar a possibilidade de interação entre os métodos de infestação e controle. O código abaixo mostra como incluir essa interação e avaliar sua importância:

mod.fit.inter <- glm( formula = Virus8 / Plants ~ Infest + Control + Infest : Control, 
                      family = binomial ( link = logit ), data = tomato , weights = Plants )
summary (mod.fit.inter )
## 
## Call:
## glm(formula = Virus8/Plants ~ Infest + Control + Infest:Control, 
##     family = binomial(link = logit), data = tomato, weights = Plants)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.466  -2.705  -1.267   2.811   6.791  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -1.0460     0.1316  -7.947 1.92e-15 ***
## Infest2            0.9258     0.1752   5.283 1.27e-07 ***
## ControlC          -0.1623     0.1901  -0.854    0.393    
## ControlN           1.1260     0.1933   5.826 5.68e-09 ***
## Infest2:ControlC  -1.2114     0.2679  -4.521 6.15e-06 ***
## Infest2:ControlN  -1.1662     0.2662  -4.381 1.18e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 278.69  on 15  degrees of freedom
## Residual deviance: 155.05  on 10  degrees of freedom
## AIC: 242.55
## 
## Number of Fisher Scoring iterations: 4
library ( package = car)
Anova (mod.fit.inter )
## Analysis of Deviance Table (Type II tests)
## 
## Response: Virus8/Plants
##                LR Chisq Df Pr(>Chisq)    
## Infest            4.060  1     0.0439 *  
## Control          91.584  2  < 2.2e-16 ***
## Infest:Control   28.224  2  7.434e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


O modelo de regressão logística estimado é \[ logit(\widehat{\pi}) = -1.9718 + 0.9258 \, \mbox{Infest2} - 0.1623 \, \mbox{ControlC} + 1.1260 \, \mbox{ControlN} \\ \qquad \qquad \qquad \qquad \qquad \qquad - 1.2114 \, \mbox{Infest2} \times \mbox{ControlC} - 1.1662 \, \mbox{Infest2}\times \mbox{ControlN}\cdot \]

Um LRT para avaliar a importância do termo de interação testa os parâmetros de regressão correspondentes a Infest2\(\times\)ControlC e Infest2\(\times\)ControlN, digamos \(\beta_4\) e \(\beta_5\). As hipóteses são \(H_0 : \beta_4 = \beta_5 = 0\) vs. \(H_1 : \beta_4\neq 0\) e/ou \(\beta_5\neq 0\). De forma equivalente, também poderíamos escrever as hipóteses em termos de comparações de modelos: \[ \begin{array}{rcl} H_0 : logit(\pi) & = & \beta_0 + \beta_1 \, \mbox{Infest2} + \beta_2 \, \mbox{ControlC} + \beta_3 \, \mbox{ControlN} \\ H_1 : logit(\pi) & = & \beta_0 + \beta_1 \, \mbox{Infest2} + \beta_2 \, \mbox{ControlC} + \beta_3 \, \mbox{ControlN} +\\ & & \qquad \qquad \beta_4 \, \mbox{Infest2} \times \mbox{ControlC} + \beta_5 \, \mbox{Infest2} \times \mbox{ControlN}\cdot \end{array} \]

A estatística de teste é \(-2\log(\Lambda) = 28.224\), e o \(p\)-valor é \(7.4\times 10^{-7}\) usando uma aproximação de \(\chi_2^2\). Assim, há fortes evidências de uma interação entre os métodos de infestação e controle.


Odds ratios


Considere novamente o modelo dado na equação \[ logit(\pi) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3\cdot \] A chance de sucesso no nível B é \(\exp(\beta_0 + \beta_1)\) porque \(x_1 = 1\), \(x_2 = 0\) e \(x_3 = 0\) e as chances de sucesso no nível A são \(\exp(\beta_0)\) porque \(x_1 = 0\), \(x_2 = 0\) e \(x_3 = 0\). A razão de chances resultante comparando o nível B com o A é \[ \dfrac{Odds_{x_1=1, \, x_2=0, \, x_3=0}}{Odds_{x_1=0, \, x_2=0, \, x_3=0}} = \dfrac{\exp(\beta_0 + \beta_1)}{\exp(\beta_0)} = \exp(\beta_1)\cdot \]

De maneira semelhante, pode-se mostrar que a razão de chances comparando C com A é \(\exp(\beta_2)\) e a razão de chances comparando D com A é \(\exp(\beta_3)\). Podemos usar a mesma técnica para comparar níveis não básicos. Por exemplo, a razão de chances comparando o nível B ao nível C é \[ \dfrac{Odds_{x_1=1, \, x_2=0, \, x_3=0}}{Odds_{x_1=0, \, x_2=1, \, x_3=0}} = \dfrac{\exp(\beta_0 + \beta_1)}{\exp(\beta_0+\beta_2)} = \exp(\beta_1-\beta_2)\cdot \]

As estimativas e os intervalos de confiança são formados da mesma forma discutida nas Seções 2.2.3 e 2.2.5. Por exemplo, a razão de chances estimada comparando o nível \(B\) com \(C\) é \(\exp(\widehat{\beta}_1-\widehat{\beta}_2)\), e o \((1-\alpha)100\)% intervalo de confiança de Wald é \[ \exp\left( \widehat{\beta}_1-\widehat{\beta}_2 \pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\beta}_1-\widehat{\beta}_2\big)}\right) \]

onde \[ \widehat{\mbox{Var}}\big(\widehat{\beta}_1-\widehat{\beta}_2\big)=\widehat{\mbox{Var}}(\widehat{\beta}_1)+\widehat{\mbox{Var}}(\widehat{\beta}_2)-2\widehat{\mbox{Cov}}(\widehat{\beta}_1,\widehat{\beta}_2)\cdot \]

Quando há duas variáveis explicativas categóricas \(X\) e \(Z\), a comparação de dois níveis para uma das variáveis é novamente realizada pela formação das razões das chances. Quando há interação entre as variáveis, essas razões de chance dependem do nível da outra variável explicativa. Considere novamente a equação \[ logit(\pi) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 z_1 + \beta_4 z_2 + \beta_5 x_1 z_1 + \beta_6 x_1 z_2 + \beta_7 x_2 z_1 + \beta_8 x_2 z_2, \] que modela as probabilidades logarítmicas de sucesso em relação a duas variáveis explicativas categóricas de 3 níveis, digamos, com os níveis A, B e C e sua interação.

Uma análise típica com foco em \(X\) é calcular as razões de chances para comparar B com A, C com A e B com C. Devido à interação, cada uma dessas razões de chances precisa ser calculada condicionalmente a um nível de \(Z\). Assim, há são 9 razões de chances diferentes entre os níveis de \(X\). Essas razões de chances são exibidas na tabela abaixo.

\[ \mbox{Razões de chances comparando níveis de X condicionais em um nível de Z.} \\ \begin{array}{ccc} \hline{Z} & \mbox{X} & \mbox{OR} \\\hline \mbox{A} & \mbox{B com A} & \exp(\beta_1) \\ & \mbox{C com A} & \exp(\beta_2) \\ & \mbox{B com C} & \exp(\beta_1-\beta_2) \\\hline \mbox{B} & \mbox{B com A} & \exp(\beta_1+\beta_5) \\ & \mbox{C com A} & \exp(\beta_2+\beta_7) \\ & \mbox{B com C} & \exp(\beta_1-\beta_2+\beta_5-\beta_7) \\\hline \mbox{C} & \mbox{B com A} & \exp(\beta_1+\beta_6) \\ & \mbox{C com A} & \exp(\beta_2+\beta_8) \\ & \mbox{B com C} & \exp(\beta_1-\beta_2+\beta_6-\beta_8) \\\hline \end{array} \]

Por exemplo, a comparação dos níveis B e C de \(X\) no nível B de \(Z\) é encontrada por \[ \dfrac{Odds_{x_1=1, \, x_2=0, \, z_1=1, \, z_2=0}}{Odds_{x_1=0, \, x_2=1, \, z_1=1, \, z_2=0}}=\dfrac{\exp(\beta_0+\beta_1+\beta_3+\beta_5)}{\exp(\beta_0+\beta_2+\beta_3+\beta_7)}=\exp(\beta_1-\beta_2+\beta_5-\beta_7)\cdot \]

As razões de chances para comparar níveis de \(Z\) em cada nível de \(X\) podem ser encontradas da mesma maneira, se necessário.

Observe que se não houver interação, então os parâmetros que definem a interação entre \(X\) e \(Z\), \(\beta_5\), \(\beta_6\), \(\beta_7\) e \(\beta_8\) são todos 0. É então evidente na tabela acima que as razões de chances comparando os níveis de \(X\) são as mesmas em todos os níveis de \(Z\) e, portanto, precisam ser computadas apenas uma vez.


Exemplo 2.14: Controle do vírus da murcha do tomateiro.


Como a interação entre os métodos de infestação e controle é significativa, normalmente nos concentraríamos apenas no modelo que inclui a interação. No entanto, é instrutivo primeiro ver como os cálculos são realizados usando o modelo sem a interação. Mais tarde, realizaremos as investigações mais complicadas com o modelo, incluindo a interação.

Para entender o efeito do Control em plantas que apresentam sintomas de infecção, calculamos as razões de chances estimadas para todas as comparações possíveis entre os níveis:

exp(mod.fit$coefficients [3:4])
## ControlC ControlN 
## 0.452342 1.674025
# Control N vs. Control C
exp(mod.fit$coefficients [4] - mod.fit$coefficients [3])
## ControlN 
## 3.700795


Por exemplo, a razão de chances estimada comparando nenhum controle a um controle biológico, o nível base é \[ \exp(0.5152) = 1.67\cdot \]

Assim, as chances estimadas de plantas apresentarem sintomas são 1.67 vezes maiores para não usar métodos de controle do que usar um controle biológico, onde o método de infestação é mantido constante. Invertemos essa razão de chances porque preferimos fazer uma afirmação sobre os efeitos do controle biológico em relação a nenhum controle, e não vice-versa. Alternativamente, poderíamos ter definido o nível base para Control = “N” antes de estimar o modelo, o que levaria a esses tipos de comparações mais rapidamente.

Assim, as chances estimadas de plantas apresentarem sintomas são 1/1.67 = 0.5973 vezes maiores para usar um método de controle biológico do que usar métodos sem controle, onde o método de infestação é mantido constante. Estima-se que o uso de ácaros (controle biológico) reduz em aproximadamente 40% as chances de uma planta apresentar sintomas.

Abaixo está o código para encontrar os intervalos de confiança correspondentes:

# Wald interval
exp( confint.default ( object = mod.fit , parm = c("ControlC", "ControlN") , level = 0.95) )
##              2.5 %    97.5 %
## ControlC 0.3492898 0.5857982
## ControlN 1.2941188 2.1654579
# Profile likelihood ratio interval
exp( confint ( object = mod.fit , parm = c("ControlC", "ControlN") , level = 0.95) )
##              2.5 %    97.5 %
## ControlC 0.3486313 0.5848759
## ControlN 1.2945744 2.1666574
# Wald interval for Control N vs. Control C
beta.hat <- mod.fit$coefficients [ -1] # Matches up beta indices with [i] to help avoid mistakes
exp( beta.hat [3] - beta.hat [2])
## ControlN 
## 3.700795
cov.mat <- vcov ( mod.fit) [2:4 ,2:4]
var.N.C <- cov.mat [3 ,3] + cov.mat [2 ,2] - 2* cov.mat [3 ,2]
CI.betas <- beta.hat [3] - beta.hat [2] + qnorm (p = c (0.025 , 0.975) )* sqrt (var.N.C)
exp(CI.betas )
## [1] 2.800422 4.890649


Por exemplo, o intervalo de confiança de Wald de 95% comparando o nível “N” ao nível “B” é de 1.29 a 2.17. Assim, com 95% de confiança, as chances de plantas apresentarem sintomas são entre 1.29 e 2.17 vezes maiores quando não se usam métodos de controle em vez de usar um controle biológico, mantendo o método de infestação constante.

Alternativamente, também poderíamos dizer com 95% de confiança, as chances de as plantas apresentarem sintomas são entre 0.46 e 0.77 vezes maior ao usar um método de controle biológico em vez de usar nenhum método de controle (mantendo o método de infestação constante). Assim, o uso de ácaros (controle biológico) reduz as chances de uma planta apresentar sintomas em aproximadamente 23% a 54%.

Para comparar nenhum controle com o controle químico, precisamos codificar \[ \exp\Big(\widehat{\beta}_3-\widehat{\beta}_2 \pm Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_3-\widehat{\beta}_2)}\Big), \] onde \[ \widehat{\mbox{Var}}(\widehat{\beta}_3-\widehat{\beta}_2) = \widehat{\mbox{Var}}(\widehat{\beta}_3)+\widehat{\mbox{Var}}(\widehat{\beta}_2)-2\widehat{\mbox{Cov}}(\widehat{\beta}_3,\widehat{\beta}_2) \] no R.

Uma maneira mais simples de encontrar o intervalo de confiança é usar a função relevel() e alterar o nível base para “C”. Após reestimar o modelo, podemos usar as funções confint.default() e confint() para obter o intervalo de confiança correspondente. Este processo é demonstrado por código no programa correspondente a este exemplo. Também mostramos exemplos no programa de uso do pacote mcprofile para encontrar o intervalo de confiança sem relevel().

Agora considere o modelo de regressão logística estimado que inclui a interação. Para entender o efeito do Control na resposta, precisaremos calcular as razões de chances entre pares de níveis de Control, mantendo o nível de Infest2 fixo em 0 (Infest = 1) ou 1 (Infest = 2).

A razão de chances comparando o nível “N” ao nível “B” com Infest2 = 0 é \(\exp(\beta_3)\): \[ \dfrac{Odds_{\mbox{ControlC=0},\, \mbox{ControlN=1}, \, \mbox{Infest2=0}}}{Odds_{\mbox{ControlC=0},\, \mbox{ControlN=0}, \, \mbox{Infest2=0}}} = \dfrac{\exp(\beta_0+\beta_3)}{\exp(\beta_0)}=\exp(\beta_3)\cdot \]

A razão de chances comparando o nível “N” ao nível “B” com Infest2 = 1 é \[ \dfrac{Odds_{\mbox{ControlC=0},\, \mbox{ControlN=1}, \, \mbox{Infest2=1}}}{Odds_{\mbox{ControlC=0},\, \mbox{ControlN=0}, \, \mbox{Infest2=1}}} = \dfrac{\exp(\beta_0+\beta_1+\beta_3+\beta_5)}{\exp(\beta_0+\beta_1)}=\exp(\beta_3+\beta_5)\cdot \]

Outras razões de chances podem ser calculadas de maneira semelhante. Abaixo estão todas as razões de chances estimadas para Control mantendo o Infest2 constante:

beta.hat <- mod.fit.inter$coefficients [ -1]
N.B.Infest2.0 <- exp( beta.hat [3])
N.B.Infest2.1 <- exp( beta.hat [3] + beta.hat [5])
C.B.Infest2.0 <- exp( beta.hat [2])
C.B.Infest2.1 <- exp( beta.hat [2] + beta.hat [4])
N.C.Infest2.0 <- exp( beta.hat [3] - beta.hat [2])
N.C.Infest2.1 <- exp( beta.hat [3] - beta.hat [2] + beta.hat [5] - beta.hat [4])
comparison <- c("N vs. B", "N vs. B", "C vs. B", "C vs. B", "N vs. C", "N vs. C")
data.frame ( Infest2 = c(0, 1, 0, 1, 0, 1) , Control = comparison , 
             OR.hat = round (c(N.B.Infest2.0, N.B.Infest2.1, C.B.Infest2.0, 
                               C.B.Infest2.1, N.C.Infest2.0, N.C.Infest2.1) , 2))
##   Infest2 Control OR.hat
## 1       0 N vs. B   3.08
## 2       1 N vs. B   0.96
## 3       0 C vs. B   0.85
## 4       1 C vs. B   0.25
## 5       0 N vs. C   3.63
## 6       1 N vs. C   3.79


Por exemplo, a razão de chances estimada comparando o nível “N” ao nível “B” com Infest2 = 1 é \[ \exp(1.1260-1.1662) = 0.96\cdot \] Assim, as chances estimadas de plantas apresentarem sintomas são 0.96 vezes maiores para não usar nenhum controle do que usar um controle biológico quando tripes infectados são liberados em casa de vegetação.

Também podemos ver porque a interação entre Infest e Control foi significativa. Os odds ratios “N” vs. “B” e “C” vs. “B” diferem bastante nos dois níveis de Infest. Podemos usar métodos semelhantes aos anteriores para encontrar intervalos de confiança. Como há muitos intervalos para calcular envolvendo muitas variâncias e covariâncias, é mais fácil usar o pacote mcprofile e sua codificação de formulação matricial:

K <- matrix ( data = c(0, 0, 0, 1, 0, 0, 
                       0, 0, 0, 1, 0, 1,
                       0, 0, 1, 0, 0, 0,
                       0, 0, 1, 0, 1, 0,
                       0, 0, -1, 1, 0, 0,
                       0, 0, -1, 1, -1, 1) , 
              nrow = 6, ncol = 6, byrow = TRUE )
linear.combo <- mcprofile ( object = mod.fit.inter , CM = K)
ci.log.OR <- confint ( object = linear.combo , level = 0.95 , adjust = "none")
data.frame ( Infest2 = c(0, 1, 0, 1, 0, 1) , comparison , 
             OR = round (exp (ci.log.OR$estimate ) ,2) , 
             OR.CI = round (exp (ci.log.OR$confint ) ,2))
##    Infest2 comparison Estimate OR.CI.lower OR.CI.upper
## C1       0    N vs. B     3.08        2.12        4.52
## C2       1    N vs. B     0.96        0.67        1.37
## C3       0    C vs. B     0.85        0.58        1.23
## C4       1    C vs. B     0.25        0.17        0.36
## C5       0    N vs. C     3.63        2.47        5.36
## C6       1    N vs. C     3.79        2.54        5.71
save.wald <- wald ( object = linear.combo )
ci.logit.wald <- confint ( object = save.wald , level = 0.95 , adjust = "none")
data.frame ( Infest2 = c(0, 1, 0, 1, 0, 1) , comparison , 
             OR = round (exp (ci.log.OR$estimate ) ,2) , 
             lower = round (exp (ci.logit.wald$confint [ ,1]) ,2) , 
             upper = round (exp (ci.logit.wald$confint [ ,2]) ,2))
##    Infest2 comparison Estimate lower upper
## C1       0    N vs. B     3.08  2.11  4.50
## C2       1    N vs. B     0.96  0.67  1.38
## C3       0    C vs. B     0.85  0.59  1.23
## C4       1    C vs. B     0.25  0.17  0.37
## C5       0    N vs. C     3.63  2.46  5.34
## C6       1    N vs. C     3.79  2.53  5.68


As colunas de K são ordenadas correspondendo aos 6 parâmetros de regressão estimados pelo modelo. Cada linha identifica a combinação linear de parâmetros necessários para uma determinada razão de chances. Por exemplo, a linha 2 corresponde à estimativa da razão de chances comparando o nível “N” ao nível “B” com Infest2 = 1, que é \(\exp(\beta_3 +\beta_5)\), \(\beta_3\) e \(\beta_5\) são o quarto e sexto parâmetros do modelo, respectivamente.

O intervalo LR perfilada de 95% comparando o nível N ao nível B com Infest2 = 1 é (0.67, 1.37). Assim, com 95% de confiança, as chances de plantas apresentarem sintomas são entre 0.67 e 1.37 vezes maiores para não usar controle do que usar um controle biológico quando tripes infectados são liberados em casa de vegetação.

Como 1 está dentro do intervalo, não há evidências suficientes para concluir que um controle biológico é eficaz nesse cenário. Interpretações semelhantes podem ser construídas para as outras razões de chance. Observe que o intervalo LR perfilada para comparar o nível N ao nível B com Infest2 = 0 é (2.12, 4.52). Como o intervalo é superior a 1, há evidências suficientes para concluir que o controle biológico reduz as chances de plantas apresentarem sintomas ao intercalar plantas infectadas com tripes não infectados.

Quando há muitas razões de chances diferentes calculadas, deve-se considerar o controle do nível geral de confiança familiar. Por exemplo, o ajuste de Bonferroni pode ser aplicado a \(g = 6\) intervalos de confiança possíveis para Control usando um nível de confiança individual 1 -0.05/6 = 0.9917. Isso garantiria que o nível geral de confiança familiar permanecesse \(\ge\) 95%. Usando o argumento adjust = “bonferroni” e ainda usando level = 0.95 em confint.mcprofile() faz esse ajuste automaticamente.

Uma melhor forma de controle da taxa de erro está disponível através do argumento adjust = “single-step”, que é baseado em lógica comparável à estatística de intervalo estudantil de Tukey para comparações de pares em ANOVA, conforme descrito em Westfall and Young (1993). Os resultados de ambos os métodos são fornecidos abaixo para os intervalos LR perfilada:

# Bonferroni
ci.log.bon <- confint ( object = linear.combo , level = 0.95 , adjust = "bonferroni")
# Single Step
ci.log.ss <- confint ( object = linear.combo , level = 0.95 , adjust = "single-step")
data.frame ( Infest2 = c(0, 1, 0, 1, 0, 1) , comparison , 
             bon = round (exp (ci.log.bon$confint ) ,2) , 
             ss = round (exp (ci.log.ss$confint ) ,2))
##   Infest2 comparison bon.lower bon.upper ss.lower ss.upper
## 1       0    N vs. B      1.86      5.16     1.87     5.12
## 2       1    N vs. B      0.59      1.56     0.60     1.55
## 3       0    C vs. B      0.51      1.40     0.52     1.39
## 4       1    C vs. B      0.15      0.41     0.15     0.41
## 5       0    N vs. C      2.17      6.14     2.18     6.09
## 6       1    N vs. C      2.22      6.59     2.24     6.54


O Bonferroni e os intervalos ajustados de um passo são muito semelhantes, embora os intervalos de um passo sejam mais apertados em ambas as extremidades do que o Bonferroni. Ambos os conjuntos de intervalos são mais amplos do que suas contrapartes não ajustadas fornecidas na caixa de saída anterior, refletindo seus maiores níveis de confiança individual. A vantagem desses intervalos ajustados é que, com qualquer um dos conjuntos, podemos afirmar que estamos 95% confiantes de que o processo cobriu todas as razões de chances, em vez de ter apenas 95% de confiança separadamente em cada uma.

Embora tenhamos aplicado esses métodos de controle de nível de confiança aos intervalos de confiança do LR perfilada aqui, a mesma sintaxe de ajuste pode ser aplicada aos intervalos Wald produzidos com este pacote.



2.2.7 Convergência dos estimadores dos parâmetros


A função glm() usa procedimentos computacionais que iteram até que a convergência seja alcançada ou até que o número máximo de iterações seja alcançado. O critério utilizado para determinar se há convergência é a mudança nos valores sucessivos da deviance residual, ao invés da mudança nas estimativas sucessivas dos parâmetros de regressão. Isso cria um critério único e compacto que equilibra equitativamente a convergência de todas as estimativas de parâmetros simultaneamente.

Se consideramos \(G^{(k)}\) denotar o desvio residual na iteração \(k\), então a convergência ocorre quando \[ \dfrac{|G^{(k)}-G^{(k-1)}|}{0.1+|G^{(k)}|} < \epsilon, \] onde \(\epsilon\) é algum número pequeno especificado maior que 0.

O numerador \(|G^{(k)}-G^{(k-1)}|\) dá uma medida geral de quanto \(\widehat{\pi}_1,\cdots,\widehat{\pi}_n\) mudam, assim, quanto \(\widehat{\beta}_0,\cdots,\widehat{\beta}_p\) mudam da iteração anterior para a iteração atual. O denominador \(0.1+|G^{(k)}|\) converte isso aproximadamente em uma mudança proporcional.

A função glm() fornece algumas maneiras de controlar como a convergência é determinada. Primeiro, seu argumento epsilon permite que o usuário declare o valor de \(\epsilon\). O valor padrão é epsilon = 10^(-8). Segundo, o argumento maxit indica o número máximo de iterações permitidas para o procedimento numérico, onde o padrão é maxit = 25. Terceiro, o valor do argumento trace = TRUE pode ser usado para ver os valores reais de \(G^{(k)}\) para cada iteração. O padrão é trace = FALSE.


Exemplo 2.15: Placekicking.


Embora os valores padrão funcionem bem para o modelo de regressão logística usando a distância como a única variável explicativa, exploramos as mudanças que ocorrem quando alteramos os valores dos argumentos para epsilon, maxit e trace:

mod.fit <- glm( formula = good ~ distance , 
                family = binomial ( link = logit ), data = placekick , 
                trace = TRUE , epsilon = 0.0001 , maxit = 50)
## Deviance = 836.7715 Iterations - 1
## Deviance = 781.1072 Iterations - 2
## Deviance = 775.8357 Iterations - 3
## Deviance = 775.7451 Iterations - 4
## Deviance = 775.745 Iterations - 5
mod.fit$control
## $epsilon
## [1] 1e-04
## 
## $maxit
## [1] 50
## 
## $trace
## [1] TRUE
mod.fit$converged
## [1] TRUE


O valor do argumento trace = TRUE resulta na impressão dos desvios residuais junto com o número de iteração. O valor do critério de convergência para iteração \(k = 5\) é \[ \dfrac{|G^{(k)}-G^{(k-1)}|}{0.1+|G^{(k)}|} = \dfrac{|775.745 - 775.7451|}{0.1+|775.745|} = 1.3\times 10^{-7}, \]

que é menor que o indicado 0.0001, então o procedimento numérico iterativo foi interrompido. Para a iteração \(k = 4\), o valor do critério de convergência é 0.00012, que é maior que 0.0001, por isso o procedimento continua. Os componentes control e converged do objeto mod.fit também são impressos na saída mostrando os valores de controle de convergência e que a convergência foi obtida.

Incluímos outro exemplo no programa correspondente onde maxit = 3. Como a convergência não é alcançada em tão poucas iterações, R imprime “Warning message: glm.fit: algorithm did not converge” após a iteração \(k = 3\). Enquanto a convergência não é alcançada, R ainda permite acesso às estimativas de parâmetros obtidas na última iteração. Como esperado, essas estimativas de parâmetros são diferentes daquelas quando a convergência foi obtida e não devem ser utilizadas.


Quando a convergência não ocorre, a primeira solução possível é tentar um número maior de iterações. Por exemplo, se as 25 iterações padrão não forem suficientes, tente 50 iterações. Se isso não funcionar, pode haver algum problema fundamental com os dados ou o modelo que está interferindo nos procedimentos numéricos iterativos.

O problema mais comum ocorre quando uma variável explicativa separa perfeitamente os dados entre os valores de \(y_i = 0\) e 1; isso é muitas vezes referido como separação completa. Além de uma mensagem de aviso de convergência, glm() também pode relatar “glm.fit: fitted probabilities numerically 0 or 1 occurred”; no entanto, esta afirmação por si só nem sempre é indicativa de problemas de ajuste de modelo. Ilustramos o problema de separação completa no próximo exemplo.


Exemplo 2.16: Separação completa.


Começamos criando um conjunto de dados simples com uma variável explicativa \(x_1\) tal que \(y = 0\) quando \(x_1\leq 5\) e \(y = 1\) quando \(x_1\geq 6\). Assim, \(x_1\) separa perfeitamente os dois valores possíveis de \(y\).

set1<-data.frame(x1 = c(1,2,3,4,5,6,7,8,9,10), y = c(0,0,0,0,0, 1,1,1,1,1))
mod.fit1<-glm(formula = y ~ x1, data = set1, family = binomial(link = logit), trace = FALSE)
summary(mod.fit1)
## 
## Call:
## glm(formula = y ~ x1, family = binomial(link = logit), data = set1, 
##     trace = FALSE)
## 
## Deviance Residuals: 
##        Min          1Q      Median          3Q         Max  
## -1.983e-05  -2.110e-08   0.000e+00   2.110e-08   1.983e-05  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -245.8   337834.2  -0.001    0.999
## x1              44.7    61172.1   0.001    0.999
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.3863e+01  on 9  degrees of freedom
## Residual deviance: 7.8648e-10  on 8  degrees of freedom
## AIC: 4
## 
## Number of Fisher Scoring iterations: 25
mod.fit1$coefficients
## (Intercept)          x1 
##  -245.84732    44.69951
plot(x = set1$x1, y = set1$y, main = "Plot for set1", ylab = "Estimated probability", 
     xlab = expression(x[1]), panel.first = grid(col = "gray", lty = "dotted"))
curve(expr = predict(object = mod.fit1, newdata = data.frame(x1 = x), type = "response"),
      col = "red", add = TRUE, lwd = 2, n = 1000)
# Firth model
library(package = logistf)
mod.fit.firth1<-logistf(formula = y ~ x1, data = set1)
# Controls printing in R Console window - the width argument can also be used
options(digits = 4)  
summary(mod.fit.firth1)
## logistf(formula = y ~ x1, data = set1)
## 
## Model fitted by Penalized ML
## Coefficients:
##                coef se(coef) lower 0.95 upper 0.95 Chisq        p method
## (Intercept) -5.3386   3.0309   -30.7818    -0.8313 6.374 0.011579      2
## x1           0.9706   0.5268     0.1962     5.5642 7.759 0.005345      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=7.759 on 1 df, p=0.005345, n=10
## Wald test = 3.395 on 1 df, p = 0.06538
names(mod.fit.firth1)  # Check what is in the object
##  [1] "coefficients"      "alpha"             "terms"            
##  [4] "var"               "df"                "loglik"           
##  [7] "iter"              "n"                 "y"                
## [10] "formula"           "call"              "conv"             
## [13] "firth"             "linear.predictors" "predict"          
## [16] "hat.diag"          "method"            "conflev"          
## [19] "ci.lower"          "ci.upper"          "prob"             
## [22] "method.ci"         "pl.iter"           "betahist"         
## [25] "pl.conv"           "flic"              "control"          
## [28] "modcontrol"        "model"
data.frame(x = set1$x1, pi.hat = mod.fit.firth1$predict)
##     x  pi.hat
## 1   1 0.01252
## 2   2 0.03238
## 3   3 0.08116
## 4   4 0.18908
## 5   5 0.38100
## 6   6 0.61900
## 7   7 0.81092
## 8   8 0.91884
## 9   9 0.96762
## 10 10 0.98748
options(digits = 7)  # Go back to default
# Put model on plot
beta.hat<-mod.fit.firth1$coefficients
curve(expr = exp(beta.hat[1]+beta.hat[2]*x)/(1 + exp(beta.hat[1]+beta.hat[2]*x)), 
      col = "blue", add = TRUE, n = 1000)
legend(x = 0.5, y = 0.8, legend = c("glm()", "logistf()"), lty = c(1,1), 
       lwd = c(2,1), col = c("red", "blue"), bty = "n")

# Check class and availability of method functions
class(mod.fit.firth1)
## [1] "logistf"
methods(class = logistf)
##  [1] add1       anova      backward   confint    drop1      extractAIC
##  [7] family     flac       flic       forward    nobs       predict   
## [13] print      profile    summary    terms      vcov      
## see '?methods' for accessing help and source code
# LRT with penalized likelihood
# logistftest(formula = y ~ x1, test = ~ x1 - 1, values = 0, data = set1)  # Old syntax
logistftest(object = mod.fit.firth1, test = ~ x1 - 1, values = 0)  # New syntax
## logistftest(object = mod.fit.firth1, test = ~x1 - 1, values = 0)
## Model fitted by Penalized ML 
## 
## Factors fixed as follows:
## (Intercept)          x1 
##          NA           0 
## 
## Likelihoods:
## Restricted model       Full model       difference 
##        -4.960074        -1.080698         3.879376 
## 
## Likelihood ratio test=7.758753 on 1 df, p=0.005345286


R indica que a convergência não ocorreu e que o modelo está tentando criar algumas estimativas de que são praticamente iguais a 0 ou 1. O gráfico na figura acima mostra o modelo estimado na última iteração. Como há uma separação entre os valores de \(y = 0\) e 1, a probabilidade estimada de sucesso é apropriadamente igual a 0 até \(x_1\) = 5 e 1 para \(x_1 = 6\) e além. A inclinação da linha entre \(x = 5\) e 6 continuará a aumentar à medida que as iterações continuarem. Se o argumento maxit foi definido para um valor mais alto que o padrão de 25, glm() indicará que a convergência foi alcançada na iteração 26 com \(\widehat{\beta}_0 = -256.85\) e \(\widehat{\beta}_1 = 46.70\).

No entanto, essas estimativas não devem ser usadas, pois ainda serão diferentes para um número maior de iterações. Por exemplo, \(\widehat{\beta}_0 = -338.02\) e \(\widehat{\beta}_1 = 60.55\) para 34 iterações em que um valor menor para epsilon permitia um número maior de iterações. As estimativas continuariam a mudar para ainda mais iterações, onde a precisão do software leva a um eventual número máximo de iterações que podem ser implementadas, 34 para este conjunto de dados em R.

Como ilustração, considere o que acontece se invertermos os valores de \(y\) em \(x_1\) = 5 e 6. A convergência ocorre em 6 iterações, e a figura abaixo mostra um gráfico do modelo estimado. A inclinação estimada é muito menor do que antes, levando a maioria dos \(\widehat{\pi}_i\) a não ser 0 ou 1.

# Example #2
set2<-data.frame(x1 = c(1,2,3,4,6,5,7,8,9,10), y = c(0,0,0,0,0, 1,1,1,1,1))
mod.fit2<-glm(formula = y ~ x1, data = set2, family = binomial(link = logit), trace = TRUE)
## Deviance = 5.958246 Iterations - 1
## Deviance = 5.187599 Iterations - 2
## Deviance = 5.029584 Iterations - 3
## Deviance = 5.018104 Iterations - 4
## Deviance = 5.018017 Iterations - 5
## Deviance = 5.018017 Iterations - 6
summary(mod.fit2)
## 
## Call:
## glm(formula = y ~ x1, family = binomial(link = logit), data = set2, 
##     trace = TRUE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4633  -0.2426   0.0000   0.2426   1.4633  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -7.159      4.759  -1.504    0.133
## x1             1.302      0.840   1.550    0.121
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13.863  on 9  degrees of freedom
## Residual deviance:  5.018  on 8  degrees of freedom
## AIC: 9.018
## 
## Number of Fisher Scoring iterations: 6
mod.fit1$coefficients
## (Intercept)          x1 
##  -245.84732    44.69951
plot(x = set2$x1, y = set2$y, main = "Plot for set2", ylab = "Estimated probability",
     panel.first = grid(col = "gray", lty = "dotted"), xlab = expression(x[1]))
curve(expr = predict(object = mod.fit2, newdata = data.frame(x1 = x), type = "response"),
      col = "red", add = TRUE, lwd = 2, n = 1000)
mod.fit.firth2<-logistf(formula = y ~ x1, data = set2)
summary(mod.fit.firth2)
## logistf(formula = y ~ x1, data = set2)
## 
## Model fitted by Penalized ML
## Coefficients:
##                   coef  se(coef)  lower 0.95 upper 0.95    Chisq          p
## (Intercept) -3.8589313 2.2069938 -12.3193578 -0.3685249 4.962374 0.02590463
## x1           0.7016239 0.3749585   0.1141795  2.1879602 6.127290 0.01331108
##             method
## (Intercept)      2
## x1               2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=6.12729 on 1 df, p=0.01331108, n=10
## Wald test = 3.501404 on 1 df, p = 0.06131681
beta.hat<-mod.fit.firth2$coefficients
curve(expr = exp(beta.hat[1]+beta.hat[2]*x)/(1 + exp(beta.hat[1]+beta.hat[2]*x)), 
      col = "blue", add = TRUE, n = 1000)
legend(x = 0.5, y = 0.8, legend = c("glm()", "logistf()"), lty = c(1,1), 
       lwd = c(2,1), col = c("red", "blue"), bty = "n")


Heinze and Schemper, 2002 esboçam uma série de opções possíveis para o que fazer quando ocorre a separação completa. As opções mais desejáveis são usar

A primeira abordagem usa a distribuição exata de cada estimador dos parâmetro de regressão condicional a certas características dos dados. Fornecemos detalhes deste procedimento e outros métodos exatos na Seção 6.2. Em segundo lugar, como a função de probabilidade aumenta sem limites durante a estimação, a função pode ser modificada, às vezes chamada de “penalizada” para potencialmente evitar que esse problema aconteça. Detalhamos essa abordagem a seguir.

Firth (1993) propõe uma modificação na função de verossimilhança para diminuir o viés das estimativas de parâmetros em tamanhos de amostra fixos. O viés mede a diferença entre o valor esperado de uma estatística e o parâmetro correspondente que ela está estimando. Idealmente, gostaríamos que o viés fosse 0. Os estimadores de máxima verossimilhança são aproximadamente imparciais, o viés é 0 para amostras grandes, mas pode haver viés significativo presente com amostras pequenas.

Como a separação completa leva a um caso extremo de viés devido a infinitas estimativas de parâmetros, a proposta de Firth é útil aqui. Para ajudar a explicar esse procedimento, escrevemos a primeira derivada em relação a \(\beta_r\) da função de log-verossimilhança como \[ \sum_{i=1}^n (y_i-\pi_i) x_{ir}, \] para \(r = 0,\cdots,p\) e \(x_{i0} = 1\).

Essa equação é frequentemente chamada de função escore. Conforme discutido na Seção 2.2.1, normalmente definiríamos a função escore igual a 0 e resolveríamos os parâmetros de regressão para produzir suas estimativas. Firth sugere adicionar uma penalidade à função escore para produzir \[ \sum_{i=1}^n \big(y_i-\pi_i +h_i(0.5-\pi_i)\big) x_{ir}, \] onde \(h_i\) é o \(i\)-ésimo elemento diagonal da matriz hat \[ H = V^{1/2}X(X^\top VX)^{-1}X^\top V^{1/2}\cdot \] Observe que \(h_i\) ajuda a medir a influência da \(i\)-ésima observação nas estimativas dos parâmetros. Discutimos a influência e a matriz chapéu (hat) com mais detalhes na Seção 5.2.3.

Novamente definimos essa função escore modificada para 0 e resolvemos os parâmetros de regressão para produzir suas estimativas. Esta penalidade tem o efeito de mudar \(y_i\); \(i = 1,\cdots, n\), de 0 ou 1 a valores ligeiramente mais próximos de 0.5. Portanto, evita que qualquer \(\widehat{\beta}_r\) se afaste muito de 0, de modo que os procedimentos numéricos iterativos sempre convergem para valores finitos (Heinze and Schemper, 2002). Os procedimentos de inferência de Wald e LR são usados da mesma forma que a estimação de máxima verossimilhança comum. Heinze and Schemper (2002) mostram que a convergência para estimativas de parâmetros finitos sempre ocorre usando esta função escore modificada.


Exemplo 2.17: Dados simulados.


Abaixo está nosso código usado para estimar o mesmo modelo do último exemplo, mas agora usando a função escore modificada:

library ( package = logistf )
mod.fit.firth1 <- logistf ( formula = y ~ x1 , data = set1 )
options ( digits = 4) # Controls printing in R Console window
summary (mod.fit.firth1 )
## logistf(formula = y ~ x1, data = set1)
## 
## Model fitted by Penalized ML
## Coefficients:
##                coef se(coef) lower 0.95 upper 0.95 Chisq        p method
## (Intercept) -5.3386   3.0309   -30.7818    -0.8313 6.374 0.011579      2
## x1           0.9706   0.5268     0.1962     5.5642 7.759 0.005345      2
## 
## Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
## 
## Likelihood ratio test=7.759 on 1 df, p=0.005345, n=10
## Wald test = 3.395 on 1 df, p = 0.06538
options ( digits = 7) # Default


A função logistf() do pacote com o mesmo nome estima o modelo de regressão logística usando o método de Firth. Na saída summary(), limitamos o número de dígitos exibidos usando o argumento digits da função options(). Caso contrário, a saída é bastante ampla. Observe que também poderíamos ter usado o argumento width de options() para controlar a largura do texto na janela do R Console.

As estimativas dos parâmetros convergem e o modelo estimado é \(logit(\widehat{\pi}) = -5.3386 + 0.9706 x_1\). Juntamente com as estimativas dos parâmetros, a saída summary() fornece informações como: \(\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_1)}=0.5268\), um intervalo de verossimilhança perfilada de 95% de 0.19 < \(\beta_1\) < 5.56 e uma estatística LRT para \(H_0 : \beta_1 = 0\) vs. \(H_1 : \beta_1\neq 0\) com \(p\)-valor = 0.0053, onde todos os cálculos utilizam a verossimilhança modificada em suas formulações.

Objetos resultantes de logistf() possuem uma classe também chamada logistf. A maioria das funções usadas com objetos de classe glm, como predict() e anova(), não estão mais disponíveis para uso. No entanto, vários itens úteis, como as estimativas dos valores \(x\) observados, estão disponíveis no objeto resultante, use names(mod.fit.firth1) para ver esses itens. Além disso, a função logistftest() dentro do pacote executa LRTs da mesma maneira que anova(), mas agora com a função de verossimilhança penalizada.


Concluímos esta seção com três notas adicionais sobre a separação completa:


2.2.8 Simulações Monte Carlo


Através do uso da simulação de Monte Carlo, as propriedades do modelo de regressão logística e os procedimentos de inferência associados podem ser investigados. Por exemplo, nesta seção examinaremos se uma distribuição normal é uma boa aproximação para a distribuição de probabilidade das estimativas dos parâmetros de regressão. Isso nos ajudará a determinar se os intervalos de confiança de Wald são apropriados.


Exemplo 2.18: Dados simulados.


Antes de podermos investigar a distribuição de nossos estimadores, precisamos ser capazes de simular respostas em um cenário de regressão logística. Consideramos arbitrariamente o modelo de regressão logística \[ logit(\pi ) = \beta_0 + \beta_1 x_1, \] onde simulamos \(x_1\) a partir de uma distribuição de probabilidade uniforme que tem valor mínimo e máximo de 0 e 1, respectivamente.

Escolhemos valores de \(\beta_0\) e \(\beta_1\) para nosso “modelo verdadeiro” de modo que \(\pi= 0.01\) quando \(x\) = 0 e \(\pi= 0.99\) quando \(x = 1\). Abaixo está nosso código para encontrar os valores de \(\beta_0\) e \(\beta_1\) que satisfazem esses requisitos.

# x=0: logit (pi) = beta_0 + beta_1 *0 = beta_0
pi.x0 <- 0.01
beta0 <- log(pi.x0 /(1 - pi.x0))
beta0
## [1] -4.59512
# x=1: logit (pi) = beta_0 + beta_1 *1 = beta_0 + beta_1
pi.x1 <- 0.99
beta1 <- log(pi.x1 /(1 - pi.x1)) - beta0
beta1
## [1] 9.19024
# Check
exp( beta0 + beta1 *c(0 ,1)) / (1 + exp( beta0 + beta1 *c(0 ,1)))
## [1] 0.01 0.99
curve ( expr = exp( beta0 + beta1 *x) / (1 + exp( beta0 + beta1 *x)), 
        xlim = c(0 ,1) , ylab = expression (pi), n = 1000 , lwd = 3, 
        xlab = expression (x [1]) , main = "True model")
grid()


Para simular dados do modelo, utilizamos muito da mesma metodologia do Capítulo 1, onde a função rbinom() foi usada para simular observações de uma distribuição de Bernoulli. A principal diferença agora é que permitimos que o parâmetro de probabilidade de sucesso varie em função de uma variável explicativa.

Primeiro, amostramos \(n = 500\) valores da distribuição uniforme usando a função runif() para encontrar valores para \(x\). Esses valores \(x\) são inseridos em nosso modelo para encontrar \(\pi\). Em seguida, usamos a função rbinom() para simular uma resposta binária \(y\) para cada \(x\). Abaixo está nosso código correspondente:

set.seed (8238)
x1 <- runif (n = 500 , min = 0, max = 1)
pi <- exp ( beta0 + beta1 *x1) / (1 + exp ( beta0 + beta1 *x1))
set.seed (1829)
y <- rbinom (n = length (x1), size = 1, prob = pi)
head ( data.frame (y, x1 , pi))
##   y        x1         pi
## 1 1 0.4447197 0.37565323
## 2 1 0.9523950 0.98459621
## 3 1 0.5019858 0.50456242
## 4 0 0.1425656 0.03609259
## 5 0 0.2042915 0.06194091
## 6 0 0.2137520 0.06718933
mod.fit <- glm( formula = y ~ x1 , family = binomial ( link = logit ))
mod.fit$coefficients
## (Intercept)          x1 
##     -5.1434     10.7140
beta.hat0 <- mod.fit$coefficients [1]
beta.hat1 <- mod.fit$coefficients [2]
exp( beta0 + beta1 *c(0 ,1)) / (1 + exp( beta0 + beta1 *c(0 ,1)))
## [1] 0.01 0.99
curve ( expr = exp( beta0 + beta1 *x) / (1 + exp( beta0 + beta1 *x)), 
        xlim = c(0 ,1) , ylab = expression (pi), n = 1000 , lwd = 3, 
        xlab = expression (x [1]) )
curve ( expr = exp( beta.hat0 + beta.hat1*x) / (1 + exp( beta.hat0 + beta.hat1 *x)), 
        xlim = c(0 ,1) , ylab = expression (pi), add = TRUE , col = "red", n = 1000)
legend (x = 0, y = 1, legend = c(" True ", " Estimated ") , lty = c(1 ,1) , 
        col = c("black", "red") , lwd = c(3, 1) , bty = "n")
grid()


Por exemplo, o primeiro valor simulado para \(x_1\) é 0.4447, o que resulta em \(\pi= 0.3757\). A resposta simulada para este \(x_1\) é um sucesso, \(y = 1\). O modelo de regressão logística é ajustado aos 500 pares \((x, y)\), resultando na estimativa, \[ logit(\widehat{\pi} ) = -5.14 + 10.71 \,x_1\cdot \]

Como esperado, isso é um pouco semelhante ao verdadeiro modelo \(logit(\pi ) = -4.60 + 9.19 \, x_1\). Observe que não usamos um argumento data em glm() porque \(x\) e \(y\) são definidos fora de um arquivo de dados.


As distribuições de probabilidade de \(\widehat{\beta}_0\) e \(\widehat{\beta}_1\) podem ser encontradas exatamente gerando aleatoriamente infinitos conjuntos de dados e, em seguida, estimando os parâmetros de regressão separadamente em cada conjunto. As estimativas podem ser resumidas em um histograma ou alguma outra forma que nos permita examinar sua distribuição. Como isso é impraticável, podemos aproximar as distribuições mostrais usando um grande número de conjuntos de dados, onde “grande” normalmente depende do tempo computacional necessário para concluir a tarefa e da precisão desejada com a qual a distribuição amostral deve ser estimada.

Em seguida, realizamos esse processo para nosso modelo de regressão logística. Usando o processo de simulação de dados do último exemplo, o próximo exemplo discute como repeti-lo. Isso nos permite investigar a distribuição de nossos estimadores.


Exemplo 2.19:
Distribuição de probabilidade de um parâmetro de regressão.


Repetimos o processo de simulação do conjunto de dados 10.000 vezes e estimamos o modelo para cada conjunto de dados simulado de tamanho \(n = 500\). Geralmente, um número menor de réplicas para o processo de simulação seria suficiente; entretanto, um valor maior é escolhido aqui para obter boas estimativas de alguns quantis de distribuição extremos.

Observe que a maior parte do nosso código está embutido em uma função que criamos, chamada sim.all(), para que possamos repetir o processo de simulação para diferentes tamanhos de amostra.

Gostaríamos de destacar a estrutura da simulação, pois a mesma estrutura é usada frequentemente para simulações em R. Os dados são armazenados em uma matriz, onde cada coluna da matriz é um conjunto de dados amostrais de \(n\) observações. Nossa matriz de dados, y.mat, armazena as respostas binárias em 500 linhas e 10.000 colunas. Uma função de análise é aplicada a cada coluna da matriz em uma etapa usando a função apply().

A sintaxe para apply() é: apply(X = , MARGIN = , FUN = , …). O argumento \(X\) é usado para a matriz de dados, o valor do argumento MARGIN especifica se a função de análise deve ser aplicada às linhas (=1) ou colunas (=2) da matriz e FUN indica a função de análise a ser aplicada. Alternativamente, poderíamos ter usado uma função for() para fazer um loop sobre cada conjunto de dados, mas apply() é mais conveniente.

# Simulate more than one data set and calculate true confidence level 
sim.all<-function(sample.size, number.data.sets, x, seed.numb, color = TRUE, 
                  linetype = "solid") {
  # Allows for black-and-white plots
  if(color == TRUE) {line.color = "red"} else {line.color = "black"}  
  # Set seed number to reproduce results and simulate responses
  set.seed(seed.numb)
  y<-rbinom(n = length(x)*number.data.sets, size = 1, prob = pi)
  # Check out simulated data
  all.x<-rep(x, times = number.data.sets)
  # Restructure y for apply() function
  y.mat<-matrix(data = y, nrow = length(x), ncol = number.data.sets)
  # y.mat[1:3,1:5]  # Check
  # Function to estimate model and obtain pi^'s - could put this outside of sim.all() too
  calc2<-function(y, x) {
    mod.fit<-glm(formula = y ~ x, family = binomial(link = logit))
    den.logLik<-logLik(mod.fit)
    mod.fit2<-glm(formula = y ~ offset(beta1*x), family = binomial(link = logit))  
    num.logLik<-logLik(mod.fit2)
    profileLR.beta<-confint(mod.fit, parm = "x", level = 0.95)  # profile LR interval
    wald.beta<-confint.default(mod.fit, parm = "x", level = 0.95)  # Wald interval
    cov.mat<-vcov(mod.fit)
    c(as.numeric(profileLR.beta),  wald.beta, mod.fit$coefficients, mod.fit$converged,
      cov.mat[1,1], cov.mat[2,2], den.logLik, num.logLik)
  }
  # round(calc2(y = y.mat[,1], x = x),4)  # Check
  # Apply the function to every column of the y.mat matrix.
  save.results<-apply(X = y.mat, MARGIN = 2, FUN = calc2, x = x)
  #Calculate the true confidence level and investigate problems
  profile.ck<-ifelse(test = beta1 > save.results[1,], yes = ifelse(test = beta1 < save.results[2,], 
                                                                   yes = 1, no = 0), no = 0)
  cat("Profile est. true conf. level", sum(profile.ck, na.rm = TRUE)/sum(!is.na(profile.ck)), "\n")  
  # Replaced number.data.sets with sum(!is.na(profile.ck)) due to the data sets where NA is produced 
  # for a confidence interval limit.
  # Number of data sets where profile LR interval did not work
  cat("Profile problems = ", sum(is.na(profile.ck)), "\n")  
  # sum(is.na(save.results[2,]))  # Problem with upper limit?
  # sum(is.na(save.results[1,]))  # Problem with lower limit?
  wald.ck<-ifelse(test = beta1 > save.results[3,], 
                  yes = ifelse(test = beta1 < save.results[4,], yes = 1, no = 0), no = 0)
  cat("Wald est. true conf. level", sum(wald.ck)/number.data.sets, "\n")
  # Plots of the first 100 intervals
  par(mfrow = c(1,2))
  plot(x = save.results[1,1:100], y = 1:100, type = "p", pch = "(", 
       xlim = c(min(save.results[c(1,3),1:100], na.rm = TRUE), 
                max(save.results[c(2,4),1:100], na.rm = TRUE)),
       xlab = expression(beta[1]), ylab = "Data set number", 
       main = "Profile likelihood ratio", cex = 0.75) 
  points(x = save.results[2,1:100], y = 1:100, pch = ")", cex = 0.75)
  segments(x0 = save.results[1,1:100], x1 =  save.results[2,1:100], 
           y0 = 1:100, y1 = 1:100, lty = "dotted")
  abline(v = beta1, lty = "solid", col = line.color)
  cat("First 100 for profile = ", sum(profile.ck[1:100])/100, "\n")
  grid()
  plot(x = save.results[3,1:100], y = 1:100, type = "p", pch = "(", 
       xlim = c(min(save.results[c(1,3),1:100], na.rm = TRUE), 
                max(save.results[c(2,4),1:100], na.rm = TRUE)),
       xlab = expression(beta[1]), ylab = "Data set number", main = "Wald", cex = 0.75) 
  points(x = save.results[4,1:100], y = 1:100, pch = ")", cex = 0.75)
  segments(x0 = save.results[3,1:100], x1 =  save.results[4,1:100], 
           y0 = 1:100, y1 = 1:100, lty = "dotted")
  abline(v = beta1, lty = "solid", col = line.color)
  cat("First 100 for wald = ", sum(wald.ck[1:100])/100, "\n")
  grid()
  # Histogram of parameter estimates with normal distribution overlay
  par(mfrow = c(1,2))
  # beta^_0
  hist(x = save.results[5,], xlab = expression(hat(beta)[0]), 
       main = paste("n = ", sample.size), freq = FALSE)
  curve(expr = dnorm(x = x, mean = mean(save.results[5,]), 
                     sd = sqrt(var(save.results[5,]))), col = line.color, add = TRUE, lwd = 2)
  grid()
  # beta^_1
  hist(x = save.results[6,], xlab = expression(hat(beta)[1]), 
       main = paste("n = ", sample.size), freq = FALSE)
  curve(expr = dnorm(x = x, mean = mean(save.results[6,]), 
                     sd = sqrt(var(save.results[6,]))), col = line.color, add = TRUE, lwd = 2)
  grid()
 # Z = (beta^_1 - beta_1)/sqrt(var(beta^_1))
  par(mfrow = c(1,1))
  z<-(save.results[6,]-beta1)/sqrt(save.results[9,])
  hist(x = z,  freq = FALSE, 
    main = paste("Figura A: Histograma e distribuição normal padrão, n = ", sample.size),
    xlab = expression((hat(beta)[1] - beta[1]) / sqrt(widehat(Var)(hat(beta)[1])) ),
    xlim = c(-4,4))
  box()
  curve(expr = dnorm(x = x, mean = 0, sd = 1), col = line.color, add = TRUE, lwd = 2)
  cat("z and standard normal quantiles: \n")
  print(quantile(x = z, probs = c(0.005, 0.05, 0.025, 0.95, 0.975, 0.995)))
  print(qnorm(p = c(0.005, 0.05, 0.025, 0.95, 0.975, 0.995)))
  grid()
  # -2log(Lambda)
  tran.stat<--2*(save.results[11,] - save.results[10,])
  hist(x = tran.stat, xlab = expression(-2*log(L(tilde(beta)[0], beta[1])/L(hat(beta)[0], 
    hat(beta)[1]))), main = paste("Figura B: n = ", sample.size), freq = FALSE, nclass = 50) 
  box()
  # For n = 500, helpful to have more classes than default
  curve(expr = dchisq(x = x, df=1), col = line.color, add = TRUE, lwd = 2, 
        n = 1000, xlim = c(0.01,16))  # Starting at 0.01 eliminates the line being plotted at 0
  cat("-2log(Lambda) and chi^2_1 quantiles: \n")
  print(quantile(x = tran.stat, probs = c(0.9, 0.95, 0.995)))
  print(qchisq(p = c(0.9, 0.95, 0.995), df = 1))
  grid()
  # Plot of beta.hat's
  par(mfrow = c(1,2))
  beta.lab<-c("Intercept", "Slope")
  plot.data<-data.frame(name = rep(x = beta.lab, each = number.data.sets), 
                        beta.hat = c(save.results[5,], save.results[6,]))
  grid()
  par(mfrow = c(1,2))
  boxplot(formula = beta.hat ~ name, data = plot.data, col = "lightblue", 
          main = paste("Box plot, n = ", sample.size), ylab = expression(hat(beta)), xlab = " ")
  grid()
  stripchart(x = beta.hat ~ name, data = plot.data, method = "jitter", vertical = TRUE, pch = 1,
    main = paste("Dot plot, n = ", sample.size), ylab = expression(hat(beta)), xlab = "")
  grid()
  # paste("n = ", sample.size)
  # Plot of first 100 models
  save.results100<-save.results[,1:100]
  par(mfrow = c(1,1))
  curve(expr = exp(beta0 + beta1*x) / (1 + exp(beta0 + beta1*x)), xlim = c(0,1), ylab = expression(pi), 
    lwd = 1, main = "", type = "n", xlab = expression(x[1]))  
  grid()
  # Estimated models
  for (i in 1:100) {
    beta.hat0<-save.results[5,i]
    beta.hat1<-save.results[6,i]
    curve(expr = exp(beta.hat0 + beta.hat1*x) / (1 + exp(beta.hat0 + beta.hat1*x)), 
          xlim = c(0,1), ylab = expression(pi), add = TRUE, col = line.color, 
          lty = linetype)  # Use solid line for book
  }
  # Plot true model on top of every other models plotted - this makes sure it is not covered
  curve(expr = exp(beta0 + beta1*x) / (1 + exp(beta0 + beta1*x)), xlim = c(0,1), 
        ylab = expression(pi), lwd = 5, n = 1000, add = TRUE)
  mtext("Figura C: Modelo verdadeiro (preto, linha grossa) e os primeiros 100 modelos estimados", 
        side =3)
  grid()
  # Check convergence
  cat("Convergence information: ", sum(save.results[7,]), "\n")
  save.results
}
sample.size<-500
set.seed(8238)
x<-runif(n = sample.size, min = 0, max = 1)
pi<-exp(beta0 + beta1*x) / (1 + exp(beta0 + beta1*x))
# Used for colored version of book
# save.it<-sim.all(sample.size = sample.size, number.data.sets = 10000, x = x, seed.numb = 1829)
# Used for black-and-white version of book
# save.it<-sim.all(sample.size = sample.size, number.data.sets = 10000, x = x, seed.numb = 1829,
#   color = FALSE, linetype = "dotted")
# Smaller number of data sets for testing purposes
save.it<-sim.all(sample.size = sample.size, number.data.sets = 100, x = x, seed.numb = 1829,
    color = FALSE, linetype = "dotted")
## Profile est. true conf. level 0.96 
## Profile problems =  0 
## Wald est. true conf. level 0.97
## First 100 for profile =  0.96

## First 100 for wald =  0.97

## z and standard normal quantiles: 
##      0.5%        5%      2.5%       95%     97.5%     99.5% 
## -2.452325 -1.430523 -1.743850  1.637402  1.771129  1.862260 
## [1] -2.575829 -1.644854 -1.959964  1.644854  1.959964  2.575829

## -2log(Lambda) and chi^2_1 quantiles: 
##      90%      95%    99.5% 
## 2.851862 3.513461 5.384444 
## [1] 2.705543 3.841459 7.879439

## Convergence information:  100
#For those intervals that achieve the correct true confidence level, we would expect
# their estimated true confidence levels to be within
number.data.sets<-1000
0.95 + qnorm(p = c(0.005, 0.995)) * sqrt(0.95*(1-0.95)/number.data.sets)
## [1] 0.9322473 0.9677527


A Figura C o modelo verdadeiro com os modelos estimados dos primeiros 100 conjuntos de dados. A variabilidade entre os modelos estimados é evidente a partir do gráfico, mas eles mostram claramente a mesma tendência geral do modelo real do qual foram amostrados. Observe que esta imagem é altamente dependente do tamanho da amostra. Por exemplo, os leitores podem executar novamente a simulação com um tamanho de amostra menor e ver como a variabilidade aumenta!

Essa variabilidade é exibida na Figura A. A figura fornece um histograma dos 10.000 valores diferentes de \[ z_0 =\dfrac{\widehat{\beta}_1-\beta_1}{\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_1)}}, \] onde incluímos uma sobreposição de distribuição normal padrão. A tabela abaixo fornece uma comparação quantílica para os valores de \(z_0\) simulados e para uma distribuição normal padrão. No geral, a aproximação de distribuição normal padrão não é tão ruim. Há alguma assimetria à esquerda na distribuição simulada, e isso leva a pequenas diferenças nos quantis extremos.

\[ \mbox{Quantis para a distribuição normal padrão e para o simulado distribuição de } z_0.\\ \begin{array}{ccccccc}\hline & 0.005 & 0.025 & 0.05 & 0.90 & 0.95 & 0.995\\\hline \mbox{Normal padrão} & -2.576 & -1.960 & -1.645 & 1.645 & 1.960 & 2.576 \\ z_0 \, \mbox{simulada} & -2.691 & -2.004 & -1.666 & 1.615 & 1.897 & 2.382 \\\hline \end{array} \]

Saber que a distribuição de probabilidade das estimativas dos parâmetros é bem aproximada por uma distribuição normal é importante porque precisamos de normalidade para um intervalo de confiança de Wald para \(\beta_1\). Nesse caso, onde a aproximação normal parece razoavelmente boa, o nível de confiança verdadeiro estimado do intervalo de Wald é 0.9520 quando o nível indicado for 0.9.

Uma aproximação ruim levaria a um intervalo com um nível de confiança verdadeiro bem diferente do nível declarado. Os detalhes do cálculo são fornecidos no programa correspondente.

A Figura B exibe um histograma dos 10.000 valores diferentes de \[ -2\log(\Lambda)=-2\log\left(\dfrac{L(\widehat{\beta}_0,\beta_1 \, | \, y_1,\cdots,y_{500})}{L(\widehat{\beta}_0,\widehat{\beta}_1 \, | \, y_1,\cdots,y_{500})}\right), \] onde incluímos uma sobreposição da distribuição \(\chi_1^2\). A tabela abaixo fornece uma comparação dos quantis para os valores simulados de \(-2\log(\Lambda )\) e para uma distribuição \(\chi_1^2\). Semelhante aos resultados de \(z_0\), vemos que a aproximação da distribuição \(\chi_1^2\) não é tão ruim, mas há mais assimetria à direita na distribuição simulada do que em uma distribuição \(\chi^2_1\). Isso leva a pequenas diferenças nos quantis extremos. O nível de confiança verdadeiro estimado para o intervalo LR do perfil é 0.9482 quando o nível declarado é 0.95.

\[ \mbox{Quantis para a distribuição } \chi_1^2 \mbox{ e para o simulado distribuição de } -2\log(\Lambda).\\ \begin{array}{ccccccc}\hline & 0.90 & 0.95 & 0.995\\\hline \chi^2_1 & 2.706 & 3.841 & 7.879\\ -2\log(\Lambda) \, \mbox{simulada} & 2.750 & 3.896 & 8.408 \\\hline \end{array} \]

Escrevemos o programa para que os leitores possam explorar outros tamanhos de amostra. Por exemplo, as distribuições normal e \(\chi^2\) não funcionam tão bem para tamanhos menores e funcionam melhor para tamanhos de amostra maiores.


Existem muitas outras investigações que podem ser realizadas para modelos de regressão logística e procedimentos associados. Por exemplo, podemos investigar a consistência das estimativas dos parâmetros. Os estimadores são consistentes se continuarem a se aproximar de seus parâmetros correspondentes à medida que o tamanho da amostra cresce. Também podemos investigar o desempenho de procedimentos de inferência para outros parâmetros. Por exemplo, podemos examinar os intervalos Wald e LR perfilada para um determinado valor de \(x\).


2.3 Modelos lineares generalizados


Os modelos de regressão logística se enquadram em uma família de modelos chamados modelos lineares generalizados (GLMs). Cada modelo linear generalizado tem três partes diferentes:

Observe que “linear” em modelos lineares generalizados vem do uso de uma combinação linear dos parâmetros de regressão com as variáveis explicativas no componente sistemático. Componentes sistemáticos alternativos podem envolver formas funcionais mais complexas, como \(x^\beta\). Estes seriam então chamados de modelos não lineares generalizados.


Funções de ligação usadas com modelos de regressão binária


Outros modelos lineares generalizados também são usados para modelar respostas binárias. Esses modelos de regressão binária têm os mesmos componentes aleatórios e sistemáticos do modelo de regressão logística, mas suas funções de ligação são diferentes do logit. O aspecto mais importante da função ligação nesses casos é que sua inversa deve garantir que \(\mbox{E}(Y)\) esteja entre 0 e 1. Por exemplo, vimos na Seção 2.2 que \[ \pi=\dfrac{\exp(\beta_0 + \beta_1 x_1 + \cdots \beta_p x_p)}{1+\exp(\beta_0 + \beta_1 x_1 + \cdots \beta_p x_p)}, \] está sempre entre 0 e 1.

Essa garantia é obtida usando uma função de distribuição acumulada (CDF) como o inverso da função de ligação. De fato, o CDF de uma distribuição de probabilidade logística é usado para regressão logística, o que resulta no nome do modelo.

Para revisar CDFs, suponha que \(X\) seja uma variável aleatória contínua com uma função de densidade de probabilidade (PDF) \(f(x)\). O CDF \(F(x)\) fornece a área sob o PDF à esquerda de \(x\). Mais formalmente, o CDF de \(X\) é \[ F(x) = P(X\leq x) = \int_{-\infty}^x f(u)\mbox{d}u, \] \(f(u)\) em vez de \(f(x)\) é usado como o integrando para evitar confusão entre o que está sendo integrado e os limites da integração.

Como todas as probabilidades estão entre 0 e 1, \(0\leq F(x)\leq 1\) para \(-\infty < x < +\infty\). Já usamos CDFs muitas vezes em R através do uso de funções como pnorm() e pchisq(). Por exemplo, pnorm(q = 1.96) pode ser expresso como \(F(1.96)\) onde \(F(\cdot )\) é a CDF de uma distribuição normal padrão. Equivalentemente, \(Z_{0.975} = 1.96\).


Exemplo 2.20: Distribuição de probabilidade logística.


Uma variável aleatória \(X\) com uma distribuição de probabilidade logística tem um PDF de \[ f(x)=\dfrac{\exp\big(-(x-\mu)/\sigma\big)}{\sigma \Big(1+\exp\big(-(x-\mu)/\sigma\big)\Big)^2}, \] para \(-\infty < x < +\infty\) e parâmetros para \(-\infty < \mu < +\infty\) e para \(\sigma >0\). A CDF é \[ F(x) = \int_{-\infty}^x \dfrac{\exp\big(-(u-\mu)/\sigma\big)}{\sigma \Big(1+\exp\big(-(u-\mu)/\sigma\big)\Big)^2}\mbox{d}u = \dfrac{1}{1+\exp\big(-(x-\mu)/\sigma\big)}=\dfrac{\exp\big((x-\mu)/\sigma\big)}{ 1+\exp\big((x-\mu)/\sigma\big)}\cdot \]

A figura abaixo exibe gráficos de \(f(x)\) e \(F(x)\) onde \(\mu=-2\) e \(\sigma= 2\). O código usado para produzir os gráficos está a seguir:

mu <- -2
sigma <- 2
# Examples for f( -2) and F( -2)
dlogis (x = -2, location = mu , scale = sigma )
## [1] 0.125
plogis (q = -2, location = mu , scale = sigma )
## [1] 0.5
par( mfrow = c(1 ,2))
curve ( expr = dlogis (x = x, location = mu , scale = sigma ), ylab = "f(x)", 
        xlab = "x", xlim = c(-15 , 15) , main = "PDF", col = "black", n = 1000)
grid()
abline (h = 0)
curve ( expr = plogis (q = x, location = mu , scale = sigma ), ylab = "F(x)", 
        xlab = "x", xlim = c(-15 , 15) , main = "CDF", col = "black", n = 1000)
grid()


As funções dlogis() e plogis() avaliam a função de densidade (Probablity Density Function, PDF) e a função de distribuição (Cumulative Density Function, CDF) para uma distribuição logística, respectivamente. Vemos que a distribuição é simétrica e centrada em \(\mu\). Pode-se mostrar que \(\mbox{E}(X) =\mu\) e \(\mbox{Var}(X) = \sigma^2\pi^2/3\).

Em comparação com uma distribuição normal com a mesma média e variância, as distribuições são centradas no mesmo local, mas a distribuição logística é mais espalhada, caudas mais grossas. Observar que \[ F(x) = \dfrac{\exp\big((x-\mu)/\sigma\big)}{1 + \exp\big((x-\mu)/\sigma\big)} = \dfrac{\exp(\beta_0+\beta_1 x)}{1 + \exp(\beta_0+\beta_1 x)} = \pi \] quando \(\beta_0 = -\mu/\sigma\) e \(\beta_1 = 1/\sigma\). Por exemplo, com \(\mu= -2\) e \(\sigma= 2\), obtemos \(\beta_0 = 1\) e \(\beta_1 = 0.5\).

Portanto, vemos que a regressão logística usa a CDF de uma distribuição logística para relacionar a probabilidade de sucesso ao preditor linear. Se invertermos essa relação para resolver o preditor linear, obtemos \[ \beta_0 + \beta_1 x = \log\left( \dfrac{F(x)}{1-F(x)}\right)= F^{-1}(\pi )\cdot \]

Podemos denotar a expressão na equação acima por \(F^{-1}(\pi )\), porque é a CDF inversa para uma distribuição logística. Portanto, a função de ligação de um modelo de regressão logística é o CDF inverso da distribuição de probabilidade logística.


Como o CDF sempre produz um valor entre 0 e 1, outros CDFs inversos são usados como funções de ligação para modelos de regressão binária. Embora a função de ligação logit seja a mais usada para regressão binária, existem duas outras que são comuns:


Exemplo 2.21: Comparando três modelos de regressão binária.


Suponha que o preditor linear tenha apenas uma variável explicativa da forma \(\beta_0+ \beta_1 x\). A figura abaixo fornece gráficos dos modelos de regressão logística, probit e log-log complementar.

Escolhemos diferentes valores de \(\beta_0\) e \(\beta_1\) para cada modelo de modo que a média seja 2 e a variância 3 para variáveis aleatórias com os CDFs correspondentes.

# All distributions - holding mean and variance constant
rm(pi)  # I used pi to mean something else than pi = 3.14 above
M<-2  # Mean
V<-3  # ariance
# x-axis limits
xlim.val<-c(qnorm(0.001, mean = M, sd = sqrt(V)),  qnorm(0.999, mean = M, sd = sqrt(V))) 
rd<-2  # Number of digits for rounding.
# Logistic
# Note: If X ~ Logistic(mu, sigma) (mu is location parameter and sigma is scale parameter), then
#       the mean (M) is M = mu and the variance (V) is V = sigma^2 * pi^2 / 3. This leads to
#       sigma = sqrt(V) * sqrt(3) / pi. (X-mu)/sigma = beta0 + beta1*X which leads to the beta0 and
#       beta1 as given below
beta0<--M*pi/(sqrt(3) * sqrt(V))  # -mu/sigma
beta1<-pi/(sqrt(3) * sqrt(V))     # 1/sigma
curve(expr = plogis(beta0 + beta1*x), xlim = xlim.val, col = "red", lwd = 2,
    lty = 1, ylim = c(0,1), xlab =  "x", ylab = expression(pi), 
    panel.first = grid(col = "gray", lty = "dotted"))
logistic.beta<-data.frame(name = "logistic", beta0 = round(beta0,rd), beta1 = round(beta1,rd))
# Normal
#  (X-mu)/sigma = beta0 + beta1*X which leads to the beta0 and beta1 as given below
beta0<--M/sqrt(V)  # -mu/sigma
beta1<-1/sqrt(V)   # 1/sigma
curve(expr = pnorm(beta0 + beta1*x), add = TRUE, col = "blue", lwd = 2, lty = "dashed")
  normal.beta<-data.frame(name = "probit", beta0 = round(beta0,rd), beta1 = round(beta1,rd))
# Complementary log-log
# Note: If X ~ Gumbel(mu, sigma) (mu is location parameter and sigma is scale parameter), then
#       Y = 1-X has a CDF of F(y) = 1 - exp(-(exp(y+mu-1)/sigma)). We can set (y+mu-1)/sigma equal to
#       beta0 + beta1*x to find the beta0 and beta1 as shown below. Also, the mean (M) of Y is
#       M = 1 - mu - euler*sigma, and the variance (V) of Y is V = pi^2 * sigma^2 / 6. Solving for
#       mu and sigma, we obtain the parameter values as shown below that are written in terms of M and V.
euler<--digamma(1)  #Euler's constant
sigma<-sqrt(V) * sqrt(6) / pi
mu<-1 - M - euler*sigma
beta0<-(mu-1)/sigma
beta1<-1/sigma
curve(expr = 1 - exp(-exp(beta0 + beta1*x)), add = TRUE, col = "green", lwd = 2, lty = "dotdash")
  cloglog.beta<-data.frame(name = "cloglog", beta0 = round(beta0,rd), beta1 = round(beta1,rd))
# Compare beta0 and beta1 among CDFs
rbind.data.frame(logistic.beta, normal.beta, cloglog.beta)
##       name beta0 beta1
## 1 logistic -2.09  1.05
## 2   probit -1.15  0.58
## 3  cloglog -2.06  0.74
# Use in legend
logistic<-substitute(paste("Logistic: ", beta[0] == beta0, ", ", beta[1] == beta1),
    list(beta0 = as.numeric(logistic.beta[2]), beta1 = as.numeric(logistic.beta[3])))
probit<-substitute(paste("Probit: ", beta[0] == beta0, ", ", beta[1] == beta1),
    list(beta0 = as.numeric(normal.beta[2]), beta1 = as.numeric(normal.beta[3])))
cloglog<-substitute(paste("Comp. log-log: ", beta[0] == beta0, ", ", beta[1] == beta1),
    list(beta0 = as.numeric(cloglog.beta[2]), beta1 = as.numeric(cloglog.beta[3])))
legend(x = -3.5, y = 1, legend = c(as.expression(logistic), as.expression(probit), 
                                    as.expression(cloglog)), lty = c(1,2,4), lwd = c(2,2,2),
       col = c("red", "blue", "green"), bty = "n", cex = 1, seg.len = 4)


Ambos os modelos de regressão logística e probit são simétricos em torno de \(\pi= 0.5\); ou seja, acima \(\pi= 0.5\) é uma imagem espelhada de abaixo \(\pi= 0.5\). O modelo log-log complementar não possui essa mesma simetria, pois uma distribuição de Gumbel é uma distribuição assimétrica. A principal diferença entre os modelos logístico e probit é que o modelo logístico sobe um pouco mais lentamente para \(\pi= 1\) e cai um pouco mais lentamente para \(\pi= 0\). Essa característica se deve a distribuição logística ter mais probabilidade em suas caudas do que uma distribuição normal.



Estimação e inferência para modelos de regressão binária


Os modelos probit e log-log complementares são estimados da mesma forma que o modelo de regressão logística. A diferença é que \(\pi_i\) é representado na função de log-verossimilhança pela correspondente especificação do modelo probit ou log-log complementar. Procedimentos numéricos iterativos são novamente necessários para encontrar as estimativas dos parâmetros. Uma vez encontradas as estimativas dos parâmetros, os mesmos procedimentos de inferência utilizados para o modelo de regressão logística estão disponíveis para os modelos probit e log-log complementares. Estes incluem os testes de Wald e LRTs da Seção 2.2.1 e os intervalos Wald e LR perfilada das Seções 2.2.3 e 2.2.4.

Por exemplo, podemos usar a função anova() para LRTs envolvendo parâmetros de regressão e a função predict() para intervalos de confiança de Wald para \(\pi\). Uma das diferenças mais importantes entre os modelos de regressão logística, probit e log-log complementar surge ao calcular as razões de chance.

Na Seção 2.2.3, mostramos que a razão de chances para \(x\) no modelo \(logit(\pi ) =\beta_0+ \beta_1 x\) é \(OR = e^{c \, \beta_1}\) para um aumento de \(c\)-unidade em \(x\). Um aspecto muito importante dessa razão de chances é que ela é a mesma para qualquer valor de \(x\). Infelizmente, isso não ocorre para modelos probit e log-log complementares. Por exemplo, considere o modelo \(probit(\pi ) =\beta_0 +\beta_1 x\). As chances de sucesso são então \[ Odds_{x} =\dfrac{\Phi(\beta_0+\beta_1 x)}{1-\Phi(\beta_0+\beta_1 x)}, \] em um determinado valor da variável explicativa \(x\). Se \(x\) for aumentado em \(c > 0\) unidades, as chances de sucesso se tornam \[ Odds_{x+c} =\dfrac{\Phi(\beta_0+\beta_1 x +\beta_1 \, c)}{1-\Phi(\beta_0+\beta_1 x+ \beta_1 \, c)}\cdot \]

Quando a razão dessas duas chances é tomada, o resultado não tem uma forma fechada e depende do valor de \(x\). As razões de chances do modelo log-log complementar também dependem de \(x\). Essa é uma das principais razões pelas quais os modelos de regressão logística são mais usados entre os modelos de regressão binária.


Exemplo 2.22: Placekicking.


O objetivo deste exemplo é comparar os modelos de regressão logística, probit e log-log complementares estimados em relação ao conjunto de dados placekicking. Usamos distance como a única variável explicativa no modelo. Nosso código para estimar esses modelos é mostrado abaixo:

mod.fit.logit <- glm( formula = good ~ distance , 
                      family = binomial ( link = logit ), data = placekick )
round ( summary (mod.fit.logit ) $coefficients , 4)
##             Estimate Std. Error  z value Pr(>|z|)
## (Intercept)   5.8121     0.3263  17.8133        0
## distance     -0.1150     0.0083 -13.7937        0
mod.fit.probit <- glm( formula = good ~ distance , 
                       family = binomial ( link = probit ), data = placekick )
round ( summary (mod.fit.probit ) $coefficients , 4)
##             Estimate Std. Error  z value Pr(>|z|)
## (Intercept)   3.2071     0.1570  20.4219        0
## distance     -0.0628     0.0043 -14.5421        0
mod.fit.cloglog <- glm( formula = good ~ distance , 
                        family = binomial ( link = cloglog ), data = placekick )
round ( summary (mod.fit.cloglog ) $coefficients , 4)
##             Estimate Std. Error  z value Pr(>|z|)
## (Intercept)   2.3802     0.1184  20.1011        0
## distance     -0.0522     0.0037 -14.0763        0


Para estimar os modelos probit e log-log complementares, alteramos o valor do argumento link em glm() para o nome da função de ligação correspondente. Os modelos estimados são:

Cada modelo resulta em uma estimativa negativa do efeito que a distância tem na probabilidade de sucesso e o teste de Wald de \(H_0 : \beta_1 = 0\) vs. \(H_1 : \beta_1\neq 0\) em cada caso resulta em \(p\)-valores muito pequenos indicando forte evidência da importância da distância.

As estimativas da probabilidade de sucesso podem ser encontradas usando as estimativas de parâmetros de um modelo e o CDF correspondente. Por exemplo, usando o modelo probit e uma distância de 20 jardas produz:

lin.pred <- as.numeric (mod.fit.probit$coefficients [1] + mod.fit.probit$coefficients [2]*20)
lin.pred
## [1] 1.951195
pnorm (q = lin.pred )
## [1] 0.9744831


Assim, \[ \widehat{\pi} = \Phi(3.2071 - 0.0628\times 20) = \Phi(1.9512) = 0.9745 \] quando distance = 20. Essas estimativas de probabilidade podem ser encontradas mais facilmente usando a função predict():

predict.data <- data.frame ( distance = c(20 , 35, 50))
logistic.pi <- predict ( object = mod.fit.logit , newdata = predict.data , type = "response")
probit.pi <- predict ( object = mod.fit.probit , newdata = predict.data , type = "response")
cloglog.pi <- predict ( object = mod.fit.cloglog , newdata = predict.data , type = "response")
round ( data.frame ( predict.data , logistic.pi , probit.pi , cloglog.pi) ,4)
##   distance logistic.pi probit.pi cloglog.pi
## 1       20      0.9710    0.9745     0.9777
## 2       35      0.8565    0.8436     0.8239
## 3       50      0.5152    0.5268     0.5477


No geral, vemos que as diferenças na probabilidade de sucesso entre os modelos são pequenas – especialmente para os modelos logísticos e probit. Isso pode ser surpreendente porque as estimativas de parâmetros para os modelos estimados são bastante diferentes. Por exemplo, \(\widehat{\beta}_1\) para o modelo log-log complementar é menos da metade do valor de \(\widehat{\beta}_1\) para o modelo logístico.

No entanto, esta é uma maneira inadequada de comparar os modelos devido às diferenças na forma como os parâmetros de regressão são usados por suas funções de ligação inversa para formar probabilidades.

As semelhanças entre os modelos estimados são vistas mais adiante na figura abaixo. Nesta figura, fornecemos um gráfico de bolhas com o ponto de plotagem proporcional ao número de observações em cada distância. Adicionados ao gráfico estão os três modelos estimados. Exceto pelas distâncias muito grandes, há pouca diferença entre esses modelos.

# Bubble plot - Much of the same code from before is included here
w<-aggregate( good ~ distance, data = placekick, FUN = sum)
n<-aggregate( good ~ distance, data = placekick, FUN = length)
w.n<-data.frame(distance = w$distance, success = w$good, trials = n$good, 
                proportion = round(w$good/n$good,4))
# Plot of the observed proportions with logistic regression model
symbols(x = w$distance, y = w$good/n$good, circles = sqrt(n$good), 
        inches = 0.5, xlab="Distance (yards)", ylab="Estimated probability",
panel.first = grid(col = "gray", lty = "dotted"))
# Estimated logistic regression model
curve(expr = predict(object = mod.fit.logit, newdata = data.frame(distance = x), 
                     type = "response"), col = "red", lwd = 2, add = TRUE, lty = 1, xlim = c(18,66))
# Estimated probit model
curve(expr = predict(object = mod.fit.probit, newdata = data.frame(distance = x), 
                     type = "response"), col = "blue", lwd = 2, add = TRUE, lty = 2, xlim = c(18,66))
# Estimated complementary log-log model
curve(expr = predict(object = mod.fit.cloglog, newdata = data.frame(distance = x), 
                     type = "response"), col = "green", lwd = 2, add = TRUE, lty = 4, xlim = c(18,66))
# Legend
legend(x = 18, y = 0.42, legend = c("Logistic", "Probit", "Complementary log-log"), 
       lty = c(1, 2, 4), lwd = c(2,2,2), bty = "n", col=c("red", "blue", "green"), cex = 1)


Terminamos este exemplo examinando as razões de chances estimadas para os três modelos:

pi.hat <- data.frame ( predict.data , logistic.pi , probit.pi , cloglog.pi)
odds.x20 <- pi.hat [1, 2:4]/(1 - pi.hat [1, 2:4])
odds.x35 <- pi.hat [2, 2:4]/(1 - pi.hat [2, 2:4])
odds.x50 <- pi.hat [3, 2:4]/(1 - pi.hat [3, 2:4])
OR.20.35 <- odds.x20/ odds.x35
OR.35.50 <- odds.x35/ odds.x50
data.frame (OR = c ("20 vs. 35" , "35 vs. 50") , round ( rbind (OR.20.35 , OR.35.50) , 2))
##          OR logistic.pi probit.pi cloglog.pi
## 1 20 vs. 35        5.61      7.08       9.36
## 2 35 vs. 50        5.61      4.84       3.86


O código calcula as chances estimadas de sucesso para uma distância de 20, 35 e 50 jardas. Formando razões de chances consecutivas, obtemos as razões de chances estimadas para incrementos de 15 jardas de distância. Como esperado, as razões de chance são exatamente as mesmas para o modelo de regressão logística, mas são diferentes para os modelos probit e log-log complementar.



Componentes aleatórios usados com modelos lineares generalizados


A distribuição de Bernoulli é o componente aleatório padrão para modelos de regressão de respostas binárias. Existem muitas outras distribuições de probabilidade que podem ser usadas, dependendo do tipo de resposta. Por exemplo, o modelo de regressão linear padrão \[ \mbox{E}(Y)=\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \] usa a distribuição normal para \(Y\) juntamente com uma função de ligação de identidade, ou seja, \(\mbox{E}(Y)\) é igual ao preditor linear.

Este tipo de modelo é apropriado com uma variável de resposta contínua que possui uma distribuição normal aproximada condicionada às variáveis explicativas. Além disso, no Capítulo 4, examinaremos modelos para respostas de contagem que surgem de outros mecanismos além da contagem de sucessos em repetidas tentativas de Bernoulli. Nesse cenário, um componente aleatório de Poisson é frequentemente utilizado, onde sua média está relacionada às variáveis explicativas por meio de uma ligação logaritmo \[ \log\big( \mbox{E}(Y)\big) = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p\cdot \]

Esta função de ligação garante que \(\mbox{E}(Y) > 0\) para qualquer combinação de parâmetros de regressão e valores de variáveis explicativas.


2.3.1 Modelos lineares generalizados com lingação paramétrica (em elaboração)


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.


Exemplo 2.23: Mortalidade de besouros.


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

Vamos utilizar diferentes modelos diferentes nas funções de ligação.

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

Percebemos que algumas curvas subestimam algumas das probabilidades observadas de mortalidade e outras superestimam as probabilidades. Por exemplo, na logdose 1.724194 todas as curvas estimadas subestimam a probabilidade de mortalidade observada, sendo as que pior estimam essa probabilidade as obtidas com as ligações probit, logit e cauchit nessa ordem. Logo percebemos que a probabilidade estimada na logdose 1.784189 todas as curvas estimadas superestimam essa probabilidade observada.


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.

Uma melhoria nas funções de ligação na regressão logística seria incorporar parámetros a estas funções com o intuito claro de diminuir os erros de estimação.


Funções de ligação paramétricas


Modelos lineares generalizados (GLM’s) permitem o tratamento de problemas de regressão nos quais a resposta pode ser distribuída de forma não normal. Mais especificamente, a distribuição da resposta pode ser qualquer distribuição numa família exponencial. Isso inclui gama, Poisson, binomial, normal e respostas gaussianas inversas.

Além disso, uma função de ligação conectando a resposta média ao preditor linear deve ser escolhida. GLM’s com ligações canônicas, como a ligação logit na regressão binomial, garantem o máximo de informação e uma interpretação simples dos parâmetros de regressão porque neste caso obtemos um modelo linear para o parâmetro natural da família exponencial subjacente.

Por exemplo, a ligação logit fornece uma representação simples das probabilidades que auxilia na interpretação dos resultados. Além disso, a concavidade da função de verossimilhança garante a unicidade dos estimadores de máxima verossimilhança (MLE).

No entanto, as ligações canônicas nem sempre fornecem o melhor ajuste disponível para um determinado conjunto de dados. Neste caso, a ligação pode ser mal especificada, o que pode produzir um viés substancial para as estimativas dos parâmetros de regressão, bem como para as estimativas da resposta média, por exemplo, ver Czado e Santner (1992) no caso da resposta binomial.

A abordagem mais comum para se proteger contra tal especificação incorreta é incorporar a função de ligação em uma ampla classe paramétrica de ligações \(\mathcal{F}=\{F(\cdot,\psi), \psi\in\Psi\}\), que inclui as ligações canônicas como um caso especial quando \(\psi = \psi_*\) digamos.

Muitas dessas classes paramétricas de funções de ligação para dados de regressão binária foram propostas na literatura. Por exemplo, Van Montford e Otten (1976), Copenhaver e Mielke (1977), Aranda-Ordaz (1981), Guerrero e Johnson (1982), Morgan (1983) e Whittmore (1983) propuseram famílias de um parâmetro enquanto Prentice (1976), Pregibon (1980), Stukel (1988) e Czado (1992) consideraram duas famílias de parâmetros. O caso geral de funções de ligação em GLMs foi estudado por Pregibon (1980) e Czado (1992). Czado (1995) desenvolveu critérios sobre como escolher tal família.

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 modelos lineares generalizados (GLM’s), reduzindo o desvio.

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\).

Todas as funções de ligação paramétricas que apresentamo a seguir estão programadas no pacote glmx, o qual depende do pacote gld.

library(glmx); library(gld)

Exemplificamos a continuação da utiliazação de diferents funções de ligação paramétricas aplicadas aos dados do Exemplo 2.23.

Pregibon

A família de funções de liagação proposta por Pregibon (1980) depende de dois parámetros e é definida como \[ g(\mu;\alpha,\beta) = \dfrac{\mu^{\alpha-\beta}-1}{\alpha-\beta}-\dfrac{(1-\mu)^{\alpha+\beta}-1}{\alpha+\beta}\cdot \]

Quando \(\alpha=\beta=0\) obtemos a ligação logit. Para \(\beta=0\) esta função é simétrica e \(\beta\) controla a assimetria, o peso das caudas é controlado por \(\alpha\). A implementação usa a distribuição lambda generalizada.

## Pregibon model
pregibin <- function(shape) binomial(link = pregibon(shape[1], shape[2]))
ajuste1 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = pregibin, xstart = c(0, 0), xlink = "logit")
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","Pregibon"), col = c("black","red"), fill = c(1,2))
grid()

Gosset

A família Gosset é dada pelo inverso da função de distribuição acumulada da distribuição \(t\)-Student. Os graus de liberdade \(\nu\) controlam o peso das caudas e está restrito a valores positivos. Para \(\nu=1\) se obtém a ligação Cauchy (cauchit) e para \(\nu\to\infty\) esta função de ligação converga para a probit.

## Gosset model
gossbin <- function(nu) binomial(link = gosset(nu))
ajuste2 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = gossbin, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste2) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Gosset"), col = c("black","red"), fill = c(1,2))
grid()

Aranda-Ordaz

A função Aranda-Ordaz (1981) simétrica é definida como \[ g(\mu;\phi) = \dfrac{2}{\phi}\dfrac{\mu^\phi-(1-\mu)^\phi}{\mu^\phi+(1-\mu)^\phi} \] e está implementada na função ao1. A função Aranda-Ordaz (1981) assimétrica é dada por \[ g(\mu;\phi) = \log\left( \dfrac{(1-\mu)^{-\phi}-1}{\phi} \right)\cdot \] Ambas contêm a função logit para \(\phi=0\) e \(\phi=1\) respectivamente, onde esta última também inclui o complementar log-log (cloglog) para \(\phi=0\).

## Aranda-Ordaz model asymmetric
aranda.ordaz <- function(phi) binomial(link = ao2(phi))
ajuste3 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = aranda.ordaz, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste3) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Aranda-Ordaz"), col = c("black","red"), fill = c(1,2))
grid()

Guerrero-Johnson

A seguinte função de ligação foi proposta em Guerrero, V. and Johnson, R. (1982). “Use of the Box-Cox Transformation with Binary Response Models.” Biometrika, 69, 309–314. A família Guerrero-Johnson é definida como \[ g(\mu;\phi) = \frac{1}{\phi}\left( \Big(\dfrac{\mu}{1-\mu}\Big)^\phi-1\right), \] sendo simétrica e contendo a função logística quando \(\phi=0\).

## Guerrero-Johnson model
guerrero.johnson = function(phi) binomial(link = gj(phi))
ajuste4 <- glmx(formula = cbind(killed,exposed-killed) ~ logdose, data = flourbeetle,
           family = guerrero.johnson, xstart = 0, xlink = "logit")
plot(killed/exposed ~ logdose, data = flourbeetle, pch = 19)
lines(fitted(ajuste3) ~ logdose, data = flourbeetle, col = "red")
lines(fitted(m[[1]]) ~ logdose, data = flourbeetle, col = "black")
legend("topleft", legend = c("Logito","Guerrero-Johnson"), col = c("black","red"), fill = c(1,2))
grid()

Aproximação à função de ligação paramétrica Stukel (Stukel (1988)).

A família de funções de ligação Stukel, proposta por Thérèse A. Stukel em 1988 também generaliza a ligação logito, mas preservando a propriedade de que se \(\eta=0\) então \(\mu=0.5\). Define-se como \[ \log\left(\dfrac{\mu}{1-\mu}\right) = h_\alpha(\eta), \] onde \(\alpha=(\alpha_1,\alpha_2)\) é o vetor de parâmetros que definem a forma da função de ligação. Desta forma, o modelo proposto é \[ \mu_\alpha(\eta)=\dfrac{\exp\big(h_\alpha(\eta)\big)}{1+\exp\big(h_\alpha(\eta)\big)}\cdot \]

A função \(h_\alpha(\eta)\) é definida como \[ h_\alpha(\eta)=\left\{ \begin{array}{ccl} \alpha_1^{-1}\Big(\exp(\alpha_1|\eta|)-1\Big) & \mbox{se} & \alpha_1>0 \\ \eta & \mbox{se} & \alpha_1=0 \\ -\alpha_1^{-1}\log\big(1-\alpha_1|\eta|\big) & \mbox{se} & \alpha_1<0 \end{array}\right. \] se \(\eta\geq 0\) ou, equivalentemente, se \(\mu\geq 0.5\) e \[ h_\alpha(\eta)=\left\{ \begin{array}{ccl} -\alpha_2^{-1}\Big(\exp(\alpha_2|\eta|)-1\Big) & \mbox{se} & \alpha_2>0 \\ \eta & \mbox{se} & \alpha_2=0 \\ \alpha_2^{-1}\log\big(1-\alpha_2|\eta|\big) & \mbox{se} & \alpha_2<0 \end{array}\right. \] se \(\eta\leq 0\) ou, equivalentemente, se \(\mu\leq 0.5\).

Quando \(\alpha_1 = \alpha_2 = 0\) obtemos o modelo logístico (logit). Caso contrário cada parâmetro governa o comportamento da curva de maneira diferente. Se \(\alpha_1 = \alpha_2\) a curva obtida é simétrica, se \(\alpha_1\neq \alpha_2\) a função de ligação é assimetrica. O modelo probito é obtido quando \(\alpha_1 = \alpha_2 ≈ 0.165\), quando \(\alpha_1 = −0.037\) e \(\alpha_2 = 0.62\) se obtém a chamada função de ligação log-log e quando \(\alpha_1 = 0.62\) e \(\alpha_2 = −0.037\) obtem-se a função de ligação complemetar log-log. No caso de \(\alpha_1 = \alpha_2 ≈ −0.077\) a função de ligação chama-se Laplace.

Podemos utilizar também uma aproximação à função de ligação paramétrica Stukel, também proposta por Stukel (1988). Suponhamos ajustamos um modelo com a função de ligação \(g_1(\mu)=g(\mu;\alpha_1,\delta_1)\) quando a função de ligação correta, mas desconhecida, é \(g_0(\mu)=g(\mu;\alpha_1,\delta_1)\).

Utilizando a expansão em série de Taylor de primeira órdem, no ponto \(g_1(\mu)\), temos a seguinte relação aproximada \[ g_0(\mu) \approx g_1(\mu)+(\alpha_0-\alpha_1)D_\alpha(g_1)+(\delta_0-\delta_1) D_\delta(g_1), \] onde \[ D_\alpha(g_1)=\left. \dfrac{\partial}{\partial\alpha} g(\mu;\alpha,\delta)\right|_{\alpha=\alpha_1,\delta=\delta_1} \] e similarmente para \(D_\delta(g_1)\).

Observemos que agora podemos aproximar a função de ligação correta \(g_0(\mu) = X\beta\) por \[ g_1(\mu)=X\beta+Z\gamma, \] onde \(Z=\big(D_\alpha(g_1),D_\delta(g_1)\big)\) e \(\gamma=-(\alpha_0-\alpha_1,\delta_0-\delta_1)\).

A forma de aproximar as estimativas dos parâmetros do modelo ampliado é escrevendo como preditor linear \[ \eta = \beta_0+\beta_1 \mbox{logdose}+\dfrac{1}{2}\alpha_1\widehat{\eta}^2\mbox{I}_{\widehat{\eta}>0}-\dfrac{1}{2}\alpha_2\widehat{\eta}^2\mbox{I}_{\widehat{\eta}<0}, \] para o caso do conjunto de dados considerado.

Neste modelo ampliado as estimativas dos coeficientes \(\alpha_1\) e \(\alpha_2\) representam aproximações de primeira ordem das estimativas destes mesmos parâmetros se utilizados diretamente na definição do modelo de regressão. Este modelo permite além de encontrar aproximações dos parâmetros que definem a nova função de ligação permite também testar se devemos ou não utilizar funções de ligação paramétricas

ajuste5 = 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(ajuste5)
## 
## 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
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()

Percebemos claramente a melhoria na qualidade na estimação das probabilidades de mortalidade.

Selecionando a função de ligação

Um problema comdiversas destas funções de ligação é a falta de um critério para escolher qual à mais adequada. Estamos dizendo que na maioria destes ajustes o AIC, AICc e BIC não estão disponíveis nos objeto de ajuste, nem está claro se estes critérios devem ser utilizados nestas situações. Devemos lembrar que estamos adicionando parámetros aos modelos, os quais devem ser considerados no cálculo desses critérios.

Como alternativa propomos utilizar diversas medidas de precisão, sendo estas:

Um erro de estimação é a diferença entre um valor observado e seu estimado. Pode ser escrito como \[ e_i = y_i-\widehat{y}_i \] sendo \(y_i\) o valor observado da resposta e \(\widehat{y}_i\) sua estimativa pelo modelo.

Observe que os erros de estimação são diferentes dos resíduos de duas maneiras. Podemos medir a precisão da estimação utilizando os erros dependentes de escala.

Os erros de estimação estão na mesma escala dos dados. As medidas de precisão baseadas apenas em \(e_i\) são, portanto, dependentes da escala e não podem ser usadas para fazer comparações entre modelos que envolvem unidades diferentes.

As duas medidas dependentes de escala mais comumente usadas são baseadas em erros absolutos ou erros quadrados:

Erro médio absoluto: MAE

\[ MAE = \dfrac{1}{n}\sum_{i=1}^n |y_i-\widehat{y}_i| \]

Erro absoluto médio: RMSE

\[ RMSE = \sqrt{\dfrac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2} \]

Ao comparar modelos o MAE é popular porque é fácil de entender e calcular. Um método de precisão que minimize o MAE levará a precisões medianas, enquanto a minimização do RMSE levará a previsões médias. Consequentemente, o RMSE também é amplamente utilizado, apesar de ser mais difícil de interpretar.

MAE = function(obs,adj) mean(abs((obs-adj)))
RMSE = function(obs,adj) sqrt(mean((obs-adj)^2))
obs = flourbeetle$killed/flourbeetle$exposed
m1.erro = MAE(obs,m[[1]]$fitted.values)
m2.erro = MAE(obs,m[[2]]$fitted.values)
m3.erro = MAE(obs,m[[3]]$fitted.values)
m4.erro = MAE(obs,m[[4]]$fitted.values)
erro1 = MAE(obs,ajuste1$fitted.values)
erro2 = MAE(obs,ajuste2$fitted.values)
erro3 = MAE(obs,ajuste3$fitted.values)
erro4 = MAE(obs,ajuste4$fitted.values)
erro5 = MAE(obs,ajuste5$fitted.values)
m5.erro = RMSE(obs,m[[1]]$fitted.values)
m6.erro = RMSE(obs,m[[2]]$fitted.values)
m7.erro = RMSE(obs,m[[3]]$fitted.values)
m8.erro = RMSE(obs,m[[4]]$fitted.values)
erro6 = RMSE(obs,ajuste1$fitted.values)
erro7 = RMSE(obs,ajuste2$fitted.values)
erro8 = RMSE(obs,ajuste3$fitted.values)
erro9 = RMSE(obs,ajuste4$fitted.values)
erro0 = RMSE(obs,ajuste5$fitted.values)
## adicionando espaço extra à margem direita do gráfico dentro do quadro
par(mar=c(4, 6, 1, 0) + 0.1, cex.axis=0.6)
## Plotando o primeiro conjunto de dados e desenhando seu eixo
plot(1:9,cbind(m1.erro,m2.erro,m3.erro,m4.erro,erro1,erro2,erro3,erro4,erro5), xlab= "Código do modelo",
     pch=16, axes=FALSE, ylim=c(0.02,0.06), ylab="", type="b", col="black", main="")
axis(2, ylim=c(0.02,0.06), col="black", las=1)  ## las=1 faz etiquetas horizontais
mtext("MAE", side=2, line=2.5)
box()
par(new=TRUE)
## Plotando o segundo gráfico e colocando a escala do eixo à direita
plot(1:9,cbind(m5.erro,m6.erro,m7.erro,m8.erro,erro6,erro7,erro8,erro9,erro0), xlab= "Código do modelo",
     pch=15, ylab="", ylim=c(0.02,0.06), axes=FALSE, type="b", col="red")
axis(2, ylim=c(0.02,0.06), lwd=2, line=3.5, col="red", col.axis = "red")
axis(1, at = 1:9, labels = c("logit", "probit", "cloglog", "cauchit","Pregibon","Gosset",
                             "Aranda-Ordaz","Guerrero \n Johnson","Stukel"), tick = FALSE)
mtext("RMSE", side=2, line=5.3, col="red")
grid()

Concluímos então que, utilizando ambas as medida de precisão, o modelo no qual minimizamos o erro de estimação é àquele no qual utilizamos a função de ligação Stukel, ou seja, o modelo no objeto R ajuste5. Observemos que neste caso é onde foi necessário utilizar mais parámetros.

Decidimos utilizar somente algumas das funções de ligação disponíveis no pacote de funções glmx


2.3.2 Modelos lineares generalizados com ligação não paramétrica (em elaboração)



2.3.3 \(R^2\) para modelos lineares generalizados (em elaboração)


Para modelos de regressão linear, o coeficiente de determinação, também conhecido como \(R^2\), é bem definido para medir a proporção de variância na variável dependente explicada pelos preditores incluídos no modelo, seguindo a lei da variância total. Embora seja direto estender tal medida para modelos lineares mistos, a heterocedasticidade natural de modelos lineares generalizados desafia sua extensão. Seguindo uma medida baseada em função de variância recentemente proposta para modelos lineares generalizados, define-se o coeficiente de determinação apropriado para modelos mistos lineares generalizados, medindo a proporção de variância na variável dependente modelada por efeitos fixos ou aleatórios ou ambos. Como a medida original definida para modelos lineares generalizados, esta nova medida precisa apenas conhecer as funções de média e variância, portanto aplicáveis a quase-modelos mais gerais. É consistente com a medida clássica de incerteza usando variância e reduz à definição clássica do coeficiente de determinação quando modelos de regressão linear mista são considerados.


Introdução


Para um par de variáveis aleatórias \(X\) e \(Y\) , a lei da variância total afirma que \[ \mbox{Var}(Y) = \mbox{Var}\big(\mbox{E}(Y\, | \, X)\big) + \mbox{E}\big(\mbox{Var}(Y\,|\, X)\big), \] decompondo a variância total de \(Y\) em duas partes: a primeira parte \(\mbox{Var}\big(\mbox{E}(Y\, | \, X)\big)\) para a variação em \(Y\) explicada por \(X\) e a segunda parte \(\mbox{E}\big(\mbox{Var}(Y\,|\, X)\big)\) para a variação em \(Y\) não explicada por \(X\). Com a razão \[\begin{equation} \tag{1} \dfrac{\mbox{Var}\big(\mbox{E}(Y\, | \, X)\big)}{\mbox{Var}(Y)}=1-\dfrac{\mbox{E}\big(\mbox{Var}(Y\,|\, X)\big)}{\mbox{Var}(Y)}, \end{equation}\] medindo a proporção da variação em \(Y\) explicada por \(X\), a lei da variância total fornece a base teórica para definir o coeficiente de determinação para modelos lineares.

Quando um modelo linear misto (McCulloch et al., 2008) é considerado para a variável de resposta observada \(Y_{ij}\) do \(j\)-ésimo indivíduo dentro do \(i\)-ésimo agrupamento, geralmente o modelamos com efeitos fixos e aleatórios, para \(j = 1,\cdots, n_i\) dentro de cada \(i = 1,\cdots, m\). Para simplificar, escrevemos o modelo linear misto correspondente como, \[\begin{equation} \tag{2} Y_{ij}=\eta^F_{ij}+\eta^R_{ij}+\epsilon_{ij}, \qquad \epsilon\sim N(0,\sigma^2), \end{equation}\] onde \(\eta^F_{ij}\) e \(\eta^R_{ij}\), respectivamente, resumem todos os efeitos fixos e aleatórios na variável de resposta com \(\eta^R_{ij} \, | \, \tau^2_{ij} \sim N (0, \tau^2_{ij})\).

Segundo a equação (\(\ref{1}\)), a proporção de variação em \(Y_{ij}\) modelada pelos efeitos fixos pode ser definido como \[\begin{equation} \tag{3} \rho^2_F = 1-\dfrac{\mbox{E}\big(\mbox{Var}(Y_{ij} \,|\, \eta^F_{ij},\tau_{ij})\big)}{\mbox{Var}(Y_{ij})} = 1- \dfrac{\mbox{E}\Big(\big(Y_{ij}-\mbox{E}(Y_{ij} \,|\, \eta^F_{ij},\tau^2_{ij})\big)^2 \Big)}{\mbox{E}\big((Y_{j}-\mbox{E}(Y_{ij}))^2\big)}\cdot \end{equation}\]

Com \(\eta^F_{ij}\) estimado por \(\widehat{\eta}^F_{ij}\) e \(\mbox{E}(Y_{ij})\) estimado por \(\overline{Y}_{\cdot\cdot}\), temos a seguinte estimativa de \(\rho^2_F\), \[\begin{equation} \tag{4} R^2_F=1-\dfrac{\displaystyle \sum_{i=1}^n \big( Y_{ij}-\widehat{\eta}^F_{ij}\big)^2}{\displaystyle \sum_{i=1}^n \big(Y_{ij}-\overline{Y}_{\cdot\cdot}\big)^2}\cdot \end{equation}\]

Como \(\mbox{E}\Big(\big(Y_{ij}-\mbox{E}(Y_{ij} \,|\, \eta^F_{ij},\tau^2_{ij})\big)^2 \Big) = \mbox{E}\Big(\mbox{E} \big( (Y_{ij}-\eta^F_{ij})^2 \, | \, \tau^2_{ij}\big)\Big) = \mbox{E}(\tau^2_{ij}) + \sigma^2\), também podemos considerar componentes de variância estimados para construir \(R^2_F\) para estimar \(\rho^2_F\), ver, por exemplo, Xu (2003); Nakagawa and Schielzeth (2013); Nakagawa et ai. (2017); Jaeger et ai. (2017). Como nossa revisão da definição de \(R^2\) em modelos mistos lineares é para lançar luz sobre suas extensões para modelos mistos lineares generalizados (McCullagh and Nelder, 1989), não seguiremos esse caminho, pois ele não pode gerenciar a heterogeneidade em modelos mistos lineares generalizados.


2.4 Exercícios


1- Começando com o modelo de regressão logística escrito como \[ \pi = \dfrac{\exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}{1+\exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}, \]

mostre que leva ao modelo escrito como \[ logit(\pi) =\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p, \] e vice versa.

2- Mostre que \(\partial \pi/\partial x_1 = \beta_1 \pi(1-\pi)\) para o modelo \[ \pi = \dfrac{\exp(\beta_0+\beta_1 x_1)}{1+\exp(\beta_0+\beta_1 x_1)}\cdot \] Discuta como a inclinação muda em função de \(x\).

3- Modificando o código adequado, descreva o que acontece com o modelo de regressão logística quando \(\beta_1\) é aumentado ou diminuído no caso de \(\beta_1 < 0\).

4- A falha de um O-ring nos foguetes de reforço do ônibus espacial Challenger levou à sua destruição em 1986. Usando dados de lançamentos anteriores de ônibus espaciais, Dalal et al. (1989) examinaram a probabilidade de uma falha do O-ring em função da temperatura no lançamento e da pressão de combustão. Os dados de seu artigo estão incluídos no arquivo:

challenger = read.csv(file = "http://leg.ufpr.br/~lucambio/CE073/20222S/challenger.csv")


Descrição das variáveis:

Flight: número do voo
Temp: Temperatura (F) no lançamento
Pressure: Pressão de combustão (psi)
O.ring: Número de falhas de O-ring do campo primário
Number: Número total de O-rings de campo primário (seis no total, três cada para os dois foguetes de reforço)

A variável de resposta é O.ring e as variáveis explicativas são Temp e Pressure. Responda o seguinte:

  1. Os autores usam a regressão logística para estimar a probabilidade de um O-ring falhar. Para usar este modelo, os autores precisaram assumir que cada O-ring é independente para cada lançamento. Discuta por que essa suposição é necessária e os possíveis problemas com ela. Observe que uma análise posterior ajudou a aliviar as preocupações dos autores sobre a independência.

  2. Estime o modelo de regressão logística usando as variáveis explicativas de forma linear.

  3. Realize testes da razão de verossimilhanças (LRTs) para julgar a importância das variáveis explicativas no modelo.

  4. Os autores optaram por remover Pressure do modelo com base nos LRTs. Com base em seus resultados, discuta por que você acha que isso foi feito. Há algum problema potencial com a remoção dessa variável?

5- Continuando o Exercício 4, considere o modelo simplificado \(logit(\pi ) =\beta_0 +\beta_1 \, \mbox{Temp}\), onde \(\pi\) é a probabilidade de falha do O-ring. Responda o seguinte:

  1. Estime o modelo.

  2. Construa dois gráficos: (1) \(\pi\) vs. Temp e (2) Número esperado de falhas vs. Temp. Use uma faixa de temperatura de 31°F a 81°F no eixo \(x\), mesmo que a temperatura mínima no conjunto de dados seja 53°F.

  3. Inclua as bandas do intervalo de confiança de Wald de 95% para no gráfico. Por que as bandas são muito mais largas para temperaturas mais baixas do que para temperaturas mais altas?

  4. A temperatura era 31°F no lançamento do Challenger em 1986. Estime a probabilidade de uma falha do O-ring usando esta temperatura e calcule um intervalo de confiança correspondente. Discuta quais suposições precisam ser feitas para aplicar os procedimentos de inferência.

  5. Em vez de usar intervalos Wald ou LR perfilada para a probabilidade de falha, Dalal et al. (1989) usam um bootstrap paramétrico para calcular intervalos. Seu processo foi:

(1) simular um grande número de conjuntos de dados, por exemplo 500 ou 1000 com $n = 23$ para cada, a partir do modelo estimado de $logit(\widehat{\pi}) = \widehat{\beta}_0 + \widehat{\beta}_1 \, \mbox{Temp}$; 

(2) estimar novos modelos para cada conjunto de dados, digamos $logit(\widehat{\pi}^*) = \widehat{\beta}_0^* + \widehat{\beta}_1^* \, \mbox{Temp}$ e 

(3) calcular $\widehat{\pi}^*$ a uma temperatura específica de interesse. Os autores usaram os quantis observados de 0.05 e 0.95 da distribuição simulada de  $\widehat{\pi}^*$ como seus limites de intervalo de confiança de 90%. Usando o bootstrap paramétrico, calcule intervalos de confiança de 90% separadamente em temperaturas de 31&#176;F e 72&#176;F.

Os métodos bootstrap usados aqui correspondem ao que é conhecido como intervalo “percentil”. Existem maneiras melhores de formar intervalos de confiança usando o bootstrap, ou seja, o BCa e os intervalos estudantis. Por favor, veja Davison and Hinkley (1997) para uma discussão e as funções boot() e boot.ci() no pacote boot para implementação.

  1. Determine se um termo quadrático é necessário no modelo para a temperatura.


6- Continuando o Exercício 5, investigue se a probabilidade estimada de falha teria mudado significativamente se um modelo de regressão probit ou log-log complementar fosse usado em vez de um modelo de regressão logística.

7- Berry and Wood (2004) apresentaram dados para determinar se uma estratégia de “congelar o chutador” implementada pela equipe adversária reduziria a probabilidade de sucesso de uma cesta de campo. Dados adicionais coletados para esta investigação estão incluídos no arquivo:

placekick.BW = read.csv ( file = "http://leg.ufpr.br/~lucambio/CE073/20222S/placekick.BW.csv")


Abaixo estão as descrições das variáveis disponíveis neste arquivo:

GameNum: Identifica o ano e o jogo
Kicker: Sobrenome do kicker
Good: Variável de resposta (“Y” = sucesso, “N” = falha)
Distance: Comprimento em jardas do field goal
Weather: Níveis de “Clouds”, “Inside”, “SnowRain” e “Sun”
Wind15: 1 se a velocidade do vento for \(\geq\) 15 milhas por hora e o placekick for ao ar livre, 0 caso contrário.
Temperature: Níveis de “Nice” (40°F < temperatura < 80° ou dentro de uma cúpula), “Cold” (temperatura 40° e ao ar livre) e “Hot” (temperatura \(\geq\) 80° e ao ar livre)
Grass: 1 se chutar em um campo de grama, 0 caso contrário
Pressure: “Y” se a tentativa for nos últimos 3 minutos de um jogo e uma cesta de campo bem-sucedida causar uma mudança de liderança, “N” caso contrário
Ice: 1 se Pressure = 1 e um tempo limite é chamado antes da tentativa, 0 caso contrário

Observe que essas variáveis são semelhantes, mas nem todas são exatamente as mesmas fornecidas para os dados de placekicking descritos na Seção 2.2.1, por exemplo, as informações foram coletadas apenas em metas de campo, portanto, não há variável PAT.

Usando esse novo conjunto de dados, conclua o seguinte:

  1. Ao usar um valor de argumento formula = Good ~ Distance em glm(), como você sabe se R está modelando a probabilidade de sucesso ou falha? Explique.

  2. Estime o modelo da parte (a) e plote-o usando a função curve().

  3. Adicione ao gráfico na parte (b) o modelo \(logit(\widehat{\pi} ) = 5.8121 - 0.1150\, \mbox{distance}\), estimado na Seção 2.2.1. Observe que os modelos são bastante semelhantes. Por que isso é desejável?

8- Continuando o Exercício 6, use as variáveis explicativas Distance, Weather, Wind15, Temperature, Grass, Pressure e Ice como termos lineares em um novo modelo de regressão logística e complete o seguinte:

  1. Estimar o modelo e definir adequadamente as variáveis indicadoras utilizadas nele.

  2. Os autores usam “Sun” como a categoria de nível base para Weather, que não é o nível padrão que R usa. Descreva como “Sun” pode ser especificado como o nível base em R.

  3. Realize os testes da razão de verossimilhanças (LRTs) para todas as variáveis explicativas para avaliar sua importância dentro do modelo. Discuta os resultados.

  4. Estime uma razão de chances apropriada para a distância e calcule o intervalo de confiança correspondente. Interprete a razão de chances.

  5. Estime todas as seis razões de chances possíveis para Weather e calcule os intervalos de confiança correspondentes. Interprete as razões de chance. Um valor de 1 está dentro de algum dos intervalos? Relacione este resultado com o teste da razão de verossimilhanças (LRT) realizado anteriormente para Weather.

  6. Continuando a parte (e), discuta como a taxa geral de erro familiar pode ser controlada para todos os intervalos de confiança. Faça os cálculos necessários.

  7. O objetivo desta parte é estimar a probabilidade de sucesso de uma cesta de campo tentada pelo Dallas Cowboys em 12 de dezembro de 2011, em seu jogo contra o Arizona Cardinals. Com o placar empatado e sete segundos restantes, o placekicker do Dallas, Dan Bailey, tentou um field goal de 49 jardas que teve sucesso. Infelizmente para Bailey, seu treinador, Jason Garrett, pediu um tempo antes da tentativa, então Bailey teve que chutar novamente. Bailey não teve sucesso na segunda vez, então Garrett essencialmente “congelou” seu próprio chutador e os Cardinals venceram o jogo na prorrogação. Uma descrição jogada a jogada do jogo está disponível em http://www.nfl.com/liveupdate/gamecenter/55351/ARI_Gamebook.pdf. Curiosamente, se o Dallas tivesse vencido um jogo adicional durante a temporada regular, o New York Giants não teria chegado aos playoffs e, posteriormente, não teria vencido o Super Bowl.

  1. Discuta a adequação de usar o modelo aqui para estimar a probabilidade de sucesso.

  2. Os valores das variáveis explicativas para este placekick são: Distance = 49 jardas, Wind15 = 0, 5 milhas por hora, Grass = 1, Pressure = “Y” e Ice = 1. O jogo foi disputado no estádio com teto retrátil do Arizona e informações sobre se o teto estava ou não aberto não estão disponíveis. Para esta parte, assuma então Weather = “Sun” e Temperature = “Nice”. Encontre a probabilidade estimada de sucesso para esta cesta de campo.

  3. Calcule o intervalo de confiança correspondente para a probabilidade de sucesso.

  1. Há evidências para concluir que congelar o kicker é uma boa estratégia a seguir?

9- Continuando o Exercício 6, considere o modelo de \[ logit(\pi ) = \beta_0+ \beta_1 \, \mbox{Distance}+\beta_2 \, \mbox{Wind15}+ \beta_3 \, \mbox{Distance} \times \mbox{Wind15}, \] onde \(\pi\) é a probabilidade de sucesso.

Complete o seguinte:

  1. Estime os modelos com e sem a interação \(\mbox{Distance} \times \mbox{Wind15}\).

  2. Construa um gráfico para ver o efeito da interação. A interação parece ser importante?

  3. Usando o modelo com o termo de interação, calcule e interprete os odds ratios para entender melhor o efeito do $ $ no sucesso ou fracasso de um placekick.

  4. Realize um teste da razão de verossimilhanças (LRT) para examinar a importância do termo de interação. Compare os resultados aqui com o que foi encontrado anteriormente na Seção 2.2.5. Sugira possíveis razões para as semelhanças e/ou diferenças.

10- O Exercício 4 do Capítulo 1 examinou a eficácia do Aptima Combo 2 Assay da Gen-Probe usado para detectar clamídia e gonorreia em indivíduos. Para os resultados do teste de clamídia em relação à sensibilidade, responda o seguinte:

  1. Use um modelo de regressão logística com Gender, Specimen e Symtoms_Status como variáveis explicativas para estimar a sensibilidade do ensaio. Use apenas termos lineares no modelo. Descreva o efeito que cada variável explicativa tem na sensibilidade do ensaio.

  2. O BD Probetec Assay é outro ensaio utilizado para detectar clamídia e gonorreia em indivíduos. Infelizmente, sabe-se que o ensaio tem uma sensibilidade muito menor para mulheres quando uma amostra de urina é testada do que quando uma amostra de swab é testada. O mesmo efeito não ocorre para os machos. Descreva por que esse efeito corresponde a uma interação entre gênero e tipo de espécime.

  3. Construa um gráfico da sensibilidade observada versus sexo condicionando o tipo de amostra e o estado sintomático. Abaixo está o código que o ajudará a construir o gráfico, os dados originais estão em um quadro ou arquivo de dados chamado aptima. O que esse gráfico sugere sobre a possibilidade de uma interação do tipo Gender e Specimen?

aptima <- read.table( file = "http://leg.ufpr.br/~lucambio/CE073/20222S/Aptima_combo.csv", 
                      header = TRUE , sep = ",")
c.table <- array ( data = c( aptima [1 ,6] , aptima [1 ,7], 
      aptima [1 ,9] , aptima [1 ,8]) , dim = c(2 ,2) , 
      dimnames = list ( True = c("+" , "-") , Assay = c("+" , "-")))
c.table
##     Assay
## True   +   -
##    + 190   7
##    -  15 464
Se.hat <- c.table [1 ,1]/ sum(c.table [1 ,]) # Sensitivity
Sp.hat <- c.table [2 ,2]/ sum(c.table [2 ,]) # Specificity
aptima$gender.numb <- ifelse ( test = aptima$Gender == "Female", yes = 0, no = 1)
aptima$sensitivity <- aptima$True_positive / Se.hat 
plot (x = aptima$gender.numb , y = aptima$sensitivity , xlab = "Gender", ylab = "Sensitivity", 
      type = "n", xaxt = "n")
axis ( side = 1, at = c(0 ,1) , labels = c("Female", "Male"))
Swab.Symptomatic <- aptima$Specimen == "Swab" & aptima$Symptoms_Status == "Symptomatic"
lines (x = aptima$gender.numb [ Swab.Symptomatic ], aptima$sensitivity [ Swab.Symptomatic ], 
       type = "o", lty = "solid", col = "black")
grid()

  1. Realize um teste de hipótese para determinar se o mesmo tipo de interação que ocorre com o ensaio BD Probetec também ocorre com o ensaio Aptima Combo 2 da Gen-Probe.

  2. O gráfico em (c) sugere que outras possíveis interações podem existir? Em caso afirmativo, realize testes de hipótese para essas interações.

11- Continuando o Exercício 10 usando o modelo com todas as três variáveis explicativas como termos lineares e uma interação Gender e Specimen, complete o seguinte:

  1. Construa intervalos de confiança para sensibilidade em todas as combinações possíveis de níveis de Gender, Specimen e Symptoms_Status. Controle o nível de erro da família em um nível apropriado.

  2. A metodologia discutida na Seção 1.1.2 também pode ser usada para calcular intervalos para cada combinação de Gender, Specimen e Symptoms_Status sem o uso de um modelo de regressão logística. Discuta as diferenças entre esses métodos e aqueles usados para a parte (a) aqui.

  3. Discutir as diferenças de metodologia entre os intervalos calculados para a parte (c) do Exercício 4 no Capítulo 1 e os intervalos construídos para a parte (a) aqui.

12- Continuando o Exercício 10, conclua o seguinte:

  1. Realize o mesmo tipo de análise das partes (a), (c) e (e) para a especificidade ao testar a clamídia.

  2. Realize o mesmo tipo de análise das partes (a), (c) e (e) para a sensibilidade e especificidade ao testar gonorréia.

13- Os dados de arremesso de lance livre de Larry Bird discutidos na Seção 1.2 podem ser examinados em um contexto de regressão logística. O objetivo aqui é estimar a probabilidade de sucesso para a segunda tentativa de lance livre dado o que aconteceu na primeira tentativa de lance livre. Complete o seguinte:

  1. O objeto c.table criado anteriormente contém o formato de tabela de contingência dos dados. Para criar o formato de quadro de dados necessário para glm(), podemos reinserir os dados em um novo quadro de dados com
bird <- data.frame ( First = c("made", "missed") , success = c(251 , 48) , trials = c(285 , 53))
bird
##    First success trials
## 1   made     251    285
## 2 missed      48     53


ou transformá-lo diretamente com 

c.table <- array ( data = c(251 , 48, 34, 5) , dim = c(2 ,2) , 
                    dimnames = list ( First = c(" made ", " missed ") , 
                                      Second = c(" made ", " missed ")))
bird1 <- as.data.frame (as.table (c.table ))
trials <- aggregate ( Freq ~ First , data = bird1 , FUN = sum)
success <- bird1 [ bird1$Second == " made ", ]
bird2 <- data.frame ( First = success$First , success = success$Freq , trials = trials$Freq )
bird2
##      First success trials
## 1    made      251    285
## 2  missed       48     53

O segundo método é mais geral e pode funcionar com tabelas de contingência maiores, como as examinadas nos Capítulos 3 e 4. Ipplemente cada linha de código no segundo método acima e descreva o que ocorre.

  1. Estime um modelo de regressão logística para a probabilidade de sucesso na segunda tentativa, onde First é a variável explicativa.

  2. Estime a razão de chances comparando os dois níveis de First. Calcule os intervalos Wald e LR perfilada para a razão de chances. Compare esses valores calculados com aqueles obtidos no exemplo de Larry Bird da Seção 1.2.5. Se você usar confint() para calcular os intervalos, certifique-se de usar o valor correto para o argumento parm, mesmo nome da variável indicadora fornecida por summary(mod.fit).

  3. Realize um teste de hipótese de \(H_0 \, : \, \beta_1 = 0\) vs. \(H_1: \, : \, \beta_1\neq 0\) usando as estatísticas Wald e da razão de verossimilhanças (LR). Compare esses resultados com aqueles encontrados para o exemplo de Larry Bird na Seção 1.2.3.

  4. Discuta por que ocorrem semelhanças e/ou diferenças entre os cálculos aqui usando a regressão logística e os cálculos correspondentes no Capítulo 1.

14- Semelhante ao Exercício 13, os dados do ensaio clínico da vacina Salk discutidos na Seção 1.2 podem ser examinados em um contexto de regressão logística. Complete o seguinte:

  1. A partir do formato da tabela de contingência dos dados fornecidos como
c.table <- array ( data = c(57 , 142 , 200688 , 201087) , dim = c(2 ,2) , 
                   dimnames = list ( Treatment = c("vacina", "placebo") , 
                                     Result = c("pólio ", " polio free ")))


crie um quadro de dados em um formato semelhante ao bird2 no Exercício 13. Nomeie esse novo quadro ou arquivo de dados polio.

  1. Use a função levels() para examinar a ordem dos níveis em polio$treatment. Observe que eles não estão em ordem alfabética devido à ordem fornecida originalmente na função array().

  2. Estime um modelo de regressão logística para a probabilidade de estar livre da pólio (segunda coluna da tabela de contingência) onde o tipo de tratamento é a variável explicativa. Indique especificamente a codificação da variável indicadora para a variável explicativa.

  3. Estime uma razão de chances comparando os dois tratamentos e calcule os intervalos de confiança correspondentes. Compare seus resultados com os encontrados na Seção 1.2.5.

  4. Determine se há uma diferença na probabilidade de estar livre da poliomielite com base no tratamento.

  5. Suponha que variáveis explicativas adicionais, como localização geográfica e status socioeconômico, estejam disponíveis para cada indivíduo e também sejam incluídas no modelo. Discuta por que usar um modelo de regressão logística seria uma maneira mais fácil de avaliar a eficácia do tratamento em vez de um formato de tabela(s) de contingência.

15- O risco relativo e a diferença entre duas probabilidades de sucesso podem ser calculados usando um modelo de regressão logística. Para este problema, considere o modelo de regressão logística de \(logit(\pi ) = \beta_0 + \beta_1 x_1\), onde \(x_1\) é 0 ou 1 e conclua o seguinte:

  1. Encontre a diferença das duas probabilidades de sucesso em termos do modelo.

  2. Encontre o risco relativo em termos do modelo.

  3. Através do método delta, pode-se encontrar a variância das estatísticas estimadas em (a) e (b). Encontre essas variações estimadas em termos de \(\widehat{\mbox{Var}}(\widehat{\beta}_0)\), \(\widehat{\mbox{Var}}(\widehat{\beta}_1)\) e \(\widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_1)\).

16- A função deltaMethod() do pacote car fornece uma maneira conveniente de aplicar o método delta para encontrar variâncias de funções matemáticas. A sintaxe da função é \[ \begin{array}{c} \mbox{deltaMethod ( objeto = <objeto de ajuste do modelo >, g = < função matemática >,} \\ \mbox{parameterNames = <Vetor de nomes de parâmetro >)}\cdot \end{array} \] O argumento parameterNames usa os mesmos nomes de parâmetro fornecidos no argumento g. Esses nomes não precisam corresponder aos dados por glm(), mas precisam estar na mesma ordem do objeto de ajuste do modelo. Aplique deltaMethod() aos dados das vacinas Larry Bird e Salk dos Exercícios 13 e 14, respectivamente, como segue:

  1. Sejam \(\pi_\mbox{First=missed}\) e \(\pi_\mbox{First=made}\) as probabilidade de sucesso na segunda tentativa dado o que aconteceu na primeira tentativa para Larry Bird. O código a seguir estima \[ \pi_\mbox{First=missed}-\pi_\mbox{ First=made} \] e \[ \widehat{\mbox{Var}}^{1/2}\big(\widehat{\pi}_\mbox{First=missed}-\widehat{\pi}_\mbox{First=made}\big)\cdot \] \[ \begin{array}{l} \mbox{library ( package = car)} \\ \mbox{deltaMethod ( object = mod.fit , g = "exp( beta0 + beta1 *1) /(1+ exp ( beta0 + beta1 *1) ) } \\ \qquad \mbox{ - exp( beta0 ) /(1+ exp( beta0 ))", parameterNames = c("beta0", "beta1"))} \end{array} \] onde mod.fit contém os resultados do ajuste do modelo com glm(). Execute este código e use-o para encontrar o intervalo de confiança de Wald correspondente para a diferença das duas probabilidades. Compare sua estimativa e intervalo de confiança com aqueles do exemplo de Larry Bird na Seção 1.2.2.

  2. Sejam \(\pi_\mbox{vaccine}\) e \(\pi_\mbox{placebo}\) as probabilidade de estar livre da poliomielite dado o tipo de tratamento recebido. Use a função deltaMethod() para estimar o risco relativo dado por \((1-\pi_\mbox{placebo})/(1-\pi_\mbox{vaccine})\) e sua variância correspondente. Calcule um intervalo de confiança de Wald para o risco relativo com esses resultados. Compare sua estimativa e intervalo de confiança com aqueles do exemplo da vacina Salk na Seção 1.2.4.

  3. Repita a parte (b), mas agora usando \(\log\big((1-\pi_\mbox{placebo})/(1-\pi_\mbox{vaccine})\big)\).

17- Considere o modelo \(logit(\pi) =\beta_0 + \beta_1 \, \mbox{distance}\) para o conjunto de dados placekicking introduzido na Seção 2.2.1. Supondo que os resultados da estimativa do modelo glm() sejam armazenados em mod.fit, conclua o seguinte:

  1. Mostre que \(\widehat{\mbox{Var}}^{1/2}(\widehat{\pi})\) é calculado como o mesmo valor usando predict() e deltaMethod() como mostrado abaixo: \[ \begin{array}{l} \mbox{predict ( object = mod.fit , newdata = data . frame ( distance = x), } \\ \qquad \quad \mbox{ type = " response ", se. fit = TRUE )}\\ \mbox{library ( package = car)} \\ \mbox{deltaMethod ( object = mod.fit , g = "exp( beta0 + beta1 *x)/(1 + exp ( beta0 + beta1 *x))", } \\ \qquad \quad \mbox{ parameterNames = c(" beta0 ", " beta1 "))} \end{array} \] onde \(x\) é substituído por um valor particular de distance.

  2. Com a ajuda da função deltaMethod(), estime \(\pi_{x+c}-\pi_x\), onde \(\pi_x\) é a probabilidade de sucesso a uma distância de \(x\) jardas e \(c\) é um incremento constante na distância. As estimativas mudam em função de \(x\) para uma constante \(c\)? Calcule os intervalos de confiança de Wald correspondentes para \(\pi_{x+c}-\pi_x\).

  3. Com a ajuda da função deltaMethod(), estime \(\pi_{x+c}/\pi_x\) para alguns valores diferentes de \(x\) e calcule os intervalos de confiança de Wald correspondentes. As estimativas mudam em função de \(x\) para uma constante \(c\)?

18- Considere o modelo \(logit(\pi) =\beta_0 + \beta_1 \, \mbox{distance} + \beta_2 \, \mbox{distance}^2\) para o conjunto de dados placekicking introduzido na Seção 2.2.1. Use as razões de chances para interpretar o efeito que a distância tem sobre o sucesso ou falha de um placekick.

19- Thorburn et al. (2001) examinam a prevalência da hepatite C entre profissionais de saúde na área de Glasgow, na Escócia. Esses profissionais de saúde foram categorizados nos seguintes grupos ocupacionais: (1) propensos à exposição (por exemplo, cirurgiões, dentistas, enfermeiras cirúrgicas), (2) contato com fluidos (por exemplo, enfermeiras não cirúrgicas), (3) equipe de laboratório, (4) contato com o paciente (por exemplo, farmacêuticos, assistentes sociais) e (5) nenhum contato com o paciente (por exemplo, administrativo). Os dados coletados estão disponíveis no arquivo:

Healthcare_worker = read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/healthcare_worker.csv")


Existe evidência de um efeito de grupo ocupacional no status de hepatite? Se houver evidência suficiente de um efeito, use as razões de chances apropriadas para fazer comparações entre os grupos. Se não houver evidência suficiente de um efeito, discuta por que esse pode ser o resultado preferido.

20- Um exemplo na Seção 2.2.5 examinado usando \[ logit(\pi ) = \beta_0 + \beta_1 \, \mbox{distance} + \beta_2 \, \mbox{wind} + \beta_3 \, \mbox{distance}\times \mbox{wind} \] para o conjunto de dados placekicking e os intervalos de confiança de Wald foram calculados para interpretar a razão de chances para o vento. Por causa da interação, intervalos de confiança foram encontrados em algumas distâncias específicas usando o código R que programou diretamente as equações de intervalo. Uma forma diferente de encontrar esses intervalos de confiança em R é por meio do pacote multcomp (Bretz et al., 2011).

  1. Implemente a seguinte continuação do código do exemplo. Verifique se o arquivo de dados contém a razão de chances de vento estimada e o intervalo de confiança Wald de 95% quando a distância é de 20 jardas.

\[ \begin{array}{l} \mbox{library ( package = multcomp )} \\ \mbox{K <- matrix ( data = c(0, 0, 1, 1*20) , nrow = 1, ncol = 4,} \\ \qquad \quad \mbox{byrow = TRUE )} \\ \mbox{linear.combo <- glht ( model = mod.fit.Ha, linfct = K)} \\ \mbox{ci.log.OR <- confint ( object = linear.combo, calpha = qnorm(0.975) )}\\ \mbox{data.frame ( estimate = exp(-ci.log.OR\$confint [1]) , lower = exp (-ci.log.OR\$confint [3]),} \\ \qquad \quad \mbox{upper = exp (-ci.log.OR\$confint [2]) )} \end{array} \]

  1. Descreva o que cada linha do código em (a) calcula dados os valores de argumento especificados. Você pode querer examinar a ajuda das funções junto com as vinhetas online para o pacote em http://cran.r-project.org/web/packages/multcomp/index.html. Observe que os elementos da função c() são colocados em uma ordem correspondente a \(\widehat{\beta}_0\), \(\widehat{\beta}_1\), \(\widehat{\beta}_2\) e \(\widehat{\beta}_3\).

  2. Modifique o código para completar os mesmos cálculos para distâncias de 30, 40, 50 e 60 jardas.

  3. Continue a implementar o código abaixo:

\[ \begin{array}{l} \mbox{K <- matriz (dados = c(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, } \\ \quad \qquad \mbox{0, 0, 0, 1) , nrow = 4, ncol = 4, byrow = TRUE )} \\ \mbox{linear.combo <- glht ( model = mod.fit.inter , linfct = K)} \\ \mbox{summary ( objeto = linear.combo , teste = ajustado ( tipo = "none"))} \end{array} \]

e compare com o summary(mod.fit.Ha).

  1. O pacote multcomp contém vários ajustes integrados para inferência simultânea, incluindo aqueles usados com o mcprofile descrito anteriormente. O padrão para testes e intervalos de confiança é usar o ajuste de etapa única. Portanto, a sintaxe para calcular intervalos de confiança com uma taxa de erro familiar fixa para essas distâncias é apenas confint(object = linear.combo). Calcule esses intervalos e explique que afirmação pode ser feita sobre o nível de confiança familiar para cobrir as verdadeiras razões de chances. Compare os intervalos com os intervalos não ajustados e comente as diferenças.

21- Em vez de encontrar a probabilidade de sucesso em um valor de variável explicativa, muitas vezes é interessante encontrar o valor de uma variável explicativa dada uma probabilidade de sucesso desejada. Isso é chamado de previsão inversa. Uma aplicação da predição inversa envolve encontrar a quantidade de um pesticida ou herbicida necessária para ter uma taxa de morte desejada quando aplicada a pragas ou plantas. O nível de dose letal \(x_\pi\), comumente chamado de “LDz”, onde \(z = 100\pi\) é definido como

\[ x_\pi = \big(logit(\pi)-\beta_0 \big)/\beta_1 \]

para o modelo de regressão logística \(logit(\pi) =\beta_0 + \beta_1 x\). O valor estimado de \(x_\pi\) é \(\widehat{x}_\pi = \big(logit(\pi) - \widehat{\beta}_0\big)/\widehat{\beta}_1\). Complete o seguinte:

  1. Mostre como \(x_\pi\) é derivado resolvendo \(x\) no modelo de regressão logística.

  2. Um intervalo de confiança de \((1-\alpha)100\%\) para \(x_\pi\) é o conjunto de todos os valores possíveis de \(x\) tais que

\[ \dfrac{\big|\widehat{\beta}_0+\widehat{\beta}_1 \, x - logit(\pi)\big|}{\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_0)+x^2 \, \widehat{\mbox{Var}}(\widehat{\beta}_1) + 2x \, \widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_0)}} < Z_{1-\alpha/2} \]

é satisfeito. Descreva como esse intervalo de confiança para \(x_\pi\) é derivado. Observe que geralmente não há solução de forma fechada para os limites do intervalo de confiança, o que leva ao uso de procedimentos numéricos iterativos.

  1. A função dose.p() do pacote MASS e a função deltaMethod() do pacote car podem calcular \(\widehat{x}_\pi\) e \(\widehat{\mbox{Var}}^{1/2}(\widehat{x}_\pi)\), onde a variância estimada é encontrada usando o método delta. Usando essa estatística e sua variância correspondente, mostre como um intervalo de confiança Wald de \((1-\alpha)100\%\) para \(x_\pi\) pode ser formado.

  2. Derive a forma de \(x_\pi\) usando modelos probit e completar log-log.

22- Turner et al. (1992) usa a regressão logística para estimar a taxa na qual o picloram, um herbicida, mata o larkspur, uma erva daninha. Seus dados foram coletados pela aplicação de quatro níveis diferentes de picloram em parcelas separadas e o número de ervas daninhas mortas dentro da parcela foi registrado. Os dados estão no arquivo

dados = read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/picloram.csv")
head(dados)
##   picloram rep kill total
## 1      0.0   1    0    18
## 2      1.1   1   13    24
## 3      2.2   1   28    29
## 4      4.5   1   24    24
## 5      0.0   2    0    29
## 6      1.1   2   11    21

O experimento foi realizado em três execuções separadas, que são indicadas pela variável rep no arquivo de dados. Os autores estimam um modelo de regressão logística sem levar em conta a variável rep e nós o mesmo aqui. Incluímos rep em uma reanálise dos dados no Exercício 2 do Capítulo 6.

Complete o seguinte:

  1. Estime um modelo de regressão logística usando a quantidade de picloram como variável explicativa e o número de ervas daninhas mortas como variável de resposta.

  2. Trace a proporção observada de ervas daninhas mortas e o modelo estimado. Descreva o quão bem o modelo se ajusta aos dados.

  3. Estime o nível de taxa de morte de 0.9 (“LD90”) para o picloram. Adicione linhas ao gráfico em (b) para ilustrar como ele é encontrado, a função segment() pode ser útil para esse propósito.

  4. Para calcular um intervalo de confiança para a dosagem que produz uma taxa de morte de 0.9, segue abaixo um conjunto de códigos R:

\[ \begin{array}{l} \mbox{root.func <- function (x, mod.fit.obj , pi0 , alpha ) \{ } \\ \qquad \quad \mbox{beta.hat <- mod.fit.obj\$coefficients } \\ \qquad \quad \mbox{cov.mat <- vcov (mod.fit.obj) } \\ \qquad \quad \mbox{var.den <- cov.mat [1 ,1] + x^2* cov.mat [2 ,2] + 2*x* cov.mat [1 ,2]} \\ \qquad \quad \mbox{abs ( beta.hat [1] + beta.hat [2]* x - log (pi0 /(1 - pi0))) / } \\ \qquad \quad \quad \mbox{ sqrt ( var.den) - qnorm (1- alpha /2) \} } \\ \mbox{lower <- uniroot (f = root.func , interval = c(min( dados\$picloram ), LD.x), }\\ \quad \quad \mbox{mod.fit.obj = mod.fit , pi0 = 0.9 , alpha = 0.95)} \\ \mbox{lower # lower\$root contains the lower bound } \\ \mbox{upper <- uniroot (f = root.func , interval = c(LD.x, max ( dados\$picloram )), } \\ \qquad \quad \mbox{mod.fit.obj = mod.fit , pi0 = 0.9 , alpha = 0.05)} \\ \mbox{upper # upper\$root contains the upper bound} \end{array} \] Descreva o que cada linha do código R faz e execute-a para encontrar um intervalo de confiança de 95% para \(x_{0.9}\).

  1. Que quantidade de picloram deve ser usada para obter uma taxa de morte de 0.9?

  2. Os dados para este problema consistem em apenas quatro diferentes níveis de dosagem de picloram. Que suposições são necessárias para que o modelo forneça uma boa estimativa de \(x_\pi\)?

23- Repita a análise do Exercício 22 usando os modelos probit e log-log complementar. Os resultados mudam em comparação com o uso da regressão logística?

24- Potter (2005) e Heinze (2006) examinam um estudo de incontinência urinária que examinou três variáveis explicativas e sua relação com a incontinência após a administração de um medicamento. Infelizmente, o autor não fornece uma descrição exata do que a codificação da variável de resposta representa, provavelmente, a resposta \(y = 1\) denota continente e \(y = 0\) denota incontinente e não fornece descrições detalhadas das variáveis explicativas; eles apenas dizem que \(x_1\), \(x_2\) e \(x_3\) são medidas no trato urinário inferior. No entanto, o principal motivo para examinar esses dados aqui é ver um exemplo de onde ocorre a separação completa quando essas variáveis explicativas são incluídas em um modelo de regressão logística. Usando os dados do arquivo

incontinence = read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/incontinence.csv")

complete o seguinte.

  1. Estime os modelos de regressão logística \(logit(\pi ) =\beta_0+ \beta_1 \, x_r\) para \(r\) = 1, 2, 3 e \(\pi= P(Y = 1)\) usando a estimativa de máxima verossimilhança. Mostre graficamente os modelos estimados juntamente com os valores observados. Há algum problema com a separação completa?

  2. Estime o modelo de regressão logística \(logit(\pi ) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3\) usando a estimativa de máxima verossimilhança. Usando as estimativas de parâmetro da última iteração de glm(), estime a probabilidade de sucesso para cada observação. Discuta os problemas que ocorrem.

  3. Estime o mesmo modelo de (b), mas agora usando a função de verossimilhança modificada conforme discutido na Seção 2.2.7. Ocorre algum problema com a estimativa? Compare os resultados aqui com aqueles encontrados em (b).

  4. Usando o modelo ajustado de (c), interprete a relação entre as variáveis explicativas e a resposta.

25- A Seção 2.2.7 discute os critérios de como R determina se e quando ocorre a convergência das estimativas dos parâmetros. Discuta por que o deviance residual é usado para esse propósito. Sugira outras maneiras de determinar a convergência.

26- O que o argumento na.action da função glm() controla e qual é seu padrão? Discuta os possíveis problemas que podem ocorrer com esse valor padrão.

27- Considere o modelo de regressão log-log complementar de \(\log\big(-\log(1-\pi)\big) = \beta_0+ \beta_1 x_1\). Mostre que a razão de chances comparando diferentes níveis de \(x_1\) não pode ser escrita sem que a variável explicativa seja apresentada.

28- A Seção 2.2.1 dá um exemplo de como estimar um modelo de regressão logística criando uma função R para calcular a função de log-verossimilhança e então maximizá-la usando optim(). Seguindo este exemplo, escreva uma função R que calcule a função log-verossimilhança para o modelo de regressão probit. Use esta função com optim() para estimar

\[ probit(\pi ) = \beta_0 + \beta_1 \, \mbox{distance}, \]

com o conjunto de dados placekicking. Compare os erros padrão estimados resultantes do uso de optim() com aqueles obtidos de glm(). Repita este processo com o modelo de regressão log-log complementar correspondente.

29- Para o modelo simples \(logit(\pi)=\beta_0 + \beta_1 x_{i1}\), pode ser mostrado que:

  1. \(\widehat{\mbox{Var}}(\widehat{\beta}_0)=\displaystyle \dfrac{1}{m}\sum_{i=1}^n x_{i1}^2 \widehat{\pi}_i(1-\widehat{\pi}_i)\),

  2. \(\widehat{\mbox{Var}}(\widehat{\beta}_1)=\displaystyle \dfrac{1}{m}\sum_{i=1}^n \widehat{\pi}_i(1-\widehat{\pi}_i)\),

  3. \(\widehat{\mbox{Cov}}(\widehat{\beta}_0,\widehat{\beta}_1)=\displaystyle \dfrac{1}{m}\sum_{i=1}^n x_{i1} \widehat{\pi}_i(1-\widehat{\pi}_i)\), onde \[ m = \left(\sum_{i=1}^n x_{i1}^2 \widehat{\pi}_i(1-\widehat{\pi}_i)\right)\left(\sum_{i=1}^n \widehat{\pi}_i(1-\widehat{\pi}_i) \right)-\left(\sum_{i=1}^n x_{i1} \widehat{\pi}_i(1-\widehat{\pi}_i)\right)^2\cdot \] Complete o seguinte:

  1. Mostre como essas fórmulas são encontradas usando \((X^\top VX)^{-1}\). Observe que a inversa de uma matriz 2\(\times\) 2 \[ \begin{pmatrix} a & b \\ c & d \end{pmatrix} = \dfrac{1}{m}\begin{pmatrix} d & -b \\ -c & a \end{pmatrix}, \] onde \(m=ad-bc\).

  2. Usando essas fórmulas, calcule a matriz de covariância estimada para o modelo \[ logit(\widehat{\pi}) = 5.8121 -0.1150 \, \mbox{distance} \] usado com o exemplo de placekicking. Verifique se esta matriz é a mesma como se vcov() fosse usado para calculá-la.

30- Em vez de usar o intervalo de confiança de Wald para \(\pi\), uma aproximação normal para \(\widehat{\pi}\) pode ser feita diretamente, resultando em um intervalo Wald de \(\widehat{\pi}\pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}(\widehat{\pi})}\), onde \(\widehat{\mbox{Var}}(\widehat{\pi})\) é encontrado usando o método delta. Para amostras grandes e não próximas de 0 ou 1, esse intervalo será muito semelhante à intervalos encontrados anteriormente. No entanto, quando isso não acontece, esse intervalo pode ter problemas. Discuta quais são esses problemas e por que esse intervalo não é desejável para uso na prática geral.

31- Suponha que um modelo de regressão logística com três variáveis explicativas binárias seja usado para estimar a probabilidade de sucesso. Este modelo inclui todos os três termos lineares, as três interações two-way e a interação de three-way. Deduza a razão de chances para examinar as chances de sucesso em dois níveis de uma das variáveis explicativas e forneça sua interpretação. Dê a variância do estimador da razão de chances.