Análise de Dados Categóricos
Capítulo 3- Analisando uma resposta multicategórica
2023-10-25
- 3.1 Distribuição de probabilidade multinomial
- 3.2 Tabelas de contingência \(I\times J\) e procedimentos de inferência
- 3.2.1 Uma distribuição multinomial
- 3.2.2 \(I\) distribuições multinomiais
- 3.2.3 Teste de independência
- 3.3 Modelos de regressão de resposta nominal
- 3.3.1 Odds ratios
- 3.3.2 Tabelas de contingência
- 3.4 Modelos de regressão de resposta ordinal
- 3.4.1 Odds ratio
- 3.4.2 Tabelas de contingência
- 3.4.3 Modelo de probabilidades não proporcionais
- 3.5 Modelos de regressão adicionais
- 3.6 Exercícios
Os capítulos 1 e 2 investigaram uma variável de resposta que tinha duas opções de categorias possíveis. O Capítulo 3 generaliza isso para uma configuração em que o valor da variável de resposta é escolhido a partir de um conjunto fixo de mais de duas opções.
Por exemplo, as opções de resposta podem ter o formato:
Escala Likert de cinco níveis - Discordo totalmente, discordo, neutro, concordo ou concordo totalmente,
Compostos químicos em experimentos de descoberta de medicamentos - Positivo, bloqueador ou nenhum,
Colocação de prateleiras de cereais em uma mercearia - fundo, meio ou topo,
Afiliação a partidos políticos canadenses - Conservador, Novo Democrata, Liberal, Bloc Quebecois ou Verde, e
Graus de carne bovina — Prime, Choice, Select, Standard, Utility e Commercial.
Para esses exemplos, algumas respostas são ordinais, por exemplo, escala Likert e outras não, por exemplo, compostos químicos. Investigaremos respostas multicategorias ordinais e nominais, não ordenadas neste capítulo.
Em cada um dos exemplos acima, uma unidade observada se encaixa exatamente em uma categoria. Por exemplo, um composto químico não pode ser tanto positivo quanto bloqueador. Existem outras situações em que uma unidade pode se encaixar simultaneamente em mais de uma categoria, como nas perguntas da pesquisa “escolha todas as que se aplicam”. Investigaremos esses problemas de “resposta múltipla” separadamente na Seção 6.4.
3.1 Distribuição de probabilidade multinomial
A distribuição de probabilidade multinomial é a extensão da distribuição binomial para situações em que há mais de duas categorias para uma resposta. Seja \(Y\) a variável aleatória resposta categórica com níveis \(j = 1,\cdots,J\), onde cada categoria tem probabilidade \(\pi_j = P(Y = j)\) tal que \(\sum_{j=1}^J \pi_j=1\).
Se houverem \(n\) tentativas idênticas com respostas \(Y_1,\cdots,Y_n\), então podemos definir variáveis aleatórias \(N_j\), \(j = 1,\cdots,J\), tal que \(N_j\) conte o número de tentativas respondendo com a categoria \(j\). Ou seja, \(N_j = \sum_{i=1}^n \pmb{I}(Y_i = j)\); onde \(\pmb{I}(\cdot ) = 1\) quando a condição entre parênteses for verdadeira e = 0 caso contrário.
Seja \(n_1,\cdots, n_J\) denotando a contagem de respostas observadas para a categoria \(j\) com \(\sum_{j=1}^J n_j = n\). A função de probabilidade (PMF) para observar um determinado conjunto de contagens \(n_1,\cdots, n_J\) é \[ P(N_1=n_1,\cdots,N_J=n_j)=\dfrac{n!}{\displaystyle \prod_{j=1}^J n_j!}\prod_{j=1}^J \pi_j^{n_j}, \] conhecida como distribuição de probabilidade multinomial.
Observe que quando \(J = 2\), a distribuição simplifica para a distribuição binomial conforme descrito na Seção 1.1.1, onde \(n_1 = w\), \(n_2 = n-w\), \(\pi_1 =\pi\) e \(\pi_2 = 1-\pi\) na notação dessa seção. Usamos a estimalção por máxima verossimilhança para obter estimativas de \(\pi_1,\cdots,\pi_J\). A função de verossimilhança é simplesmente a equação acima e o MLE para cada \(\pi_j\) é \(\widehat{\pi}_j = n_j/n\), ou seja, a proporção observada para cada categoria.
Exemplo 3.1: Amostra simulada multinomial.
Para imaginar como podem ser os dados de uma distribuição multinomial, simulamos uma amostra de n = 1.000 tentativas de uma distribuição multinomial com \[ \pi_1 = 0.25, \quad \pi_2 = 0.35, \quad \pi_3 = 0.2, \quad \pi_4 = 0.1 \quad \mbox{e} \quad \pi_5 = 0.1\cdot \]
Usamos a função rmultinom() para simular um conjunto (\(n = 1\)) de 1.000 tentativas (tamanho = 1.000) de uma distribuição multinomial com probabilidades listadas no argumento prob:
pi.j <- c(0.25 , 0.35 , 0.2 , 0.1 , 0.1)
set.seed (2195) # Set a seed to be able to reproduce the sample
n.j <- rmultinom (n = 1, size = 1000 , prob = pi.j)
data.frame (n.j, pihat.j = n.j/1000 , pi.j)## n.j pihat.j pi.j
## 1 242 0.242 0.25
## 2 333 0.333 0.35
## 3 188 0.188 0.20
## 4 122 0.122 0.10
## 5 115 0.115 0.10As contagens simuladas são salvas em um objeto chamado \(n.j\). Depois de dividir essas contagens por \(n\), vemos que as proporções observadas são muito semelhantes aos valores reais dos parâmetros. Se quiséssemos simular \(m\) conjuntos de 1000 tentativas da mesma distribuição multinomial, poderíamos alterar o valor do argumento \(n\) para \(m\). Nesta situação, haveria \(m\) colunas de observações produzidas por rmultinom(). Por exemplo, o resultado para \(m = 5\) é
set.seed (9182)
n.j <- rmultinom (n = 5, size = 1000 , prob = pi.j)
n.j## [,1] [,2] [,3] [,4] [,5]
## [1,] 259 259 237 264 247
## [2,] 341 346 374 339 341
## [3,] 200 188 198 191 210
## [4,] 92 106 89 108 107
## [5,] 108 101 102 98 95Como podemos ver, há variabilidade nas contagens de um conjunto de amostras para o próximo. Este mesmo tipo de variabilidade deve ser esperado em amostras reais com respostas multinomiais.
3.2 Tabelas de contingência \(I\times J\) e procedimentos de inferência
A Seção 1.2 discute como as contagens de duas distribuições binomiais podem ser analisadas na forma de uma tabela de contingência \(2\times 2\). Expandimos esta discussão agora para permitir mais de duas linhas e/ou colunas. Assumimos que duas variáveis categóricas, \(X\) e \(Y\), são medidas em cada unidade com níveis de \(i = 1,\cdots,I\) para \(X\) e \(j = 1,\cdots,J\) para \(Y\). Nossa amostra é composta por \(n\) unidades que são classificadas cruzadamente de acordo com seus níveis de \(X\) e \(Y\). As contagens de unidades para as quais a combinação \((X = i; Y = j)\) é observada são denotadas por \(n_{ij}\) e representam observações das variáveis aleatórias \(N_{ij}\), \(i = 1,\cdots,I\), \(j = 1,\cdots,J\).
Nosso objetivo é usar as contagens observadas para fazer inferências sobre a distribuição do \(N_{ij}\) e identificar quaisquer relações entre \(X\) e \(Y\). Para isso, podemos criar uma tabela de contingência resumindo essas contagens, sendo \(X\) a variável de linha e \(Y\) a variável de coluna. A tabela abaixo fornece a tabela geral de contingência. Observe que usamos um subscrito + para indicar uma soma; por exemplo, \(n_{+j} = \sum_{i=1}^I n_{ij}\) é o total da coluna \(j\). Vamos denotar \(n_{++}\) simplesmente por \(n\).
\[ \mbox{Tabela de contingência } I\times J\\ \begin{array}{ccccccc} \hline \mbox{} & & & & Y & & \\ \mbox{} & & 1 & 2 & \cdots & J & \mbox{Total}\\\hline \mbox{} & 1 & n_{11} & n_{12} & \cdots & n_{1J} & n_{1+} \\ \mbox{} & 2 & n_{21} & n_{22} & \cdots & n_{2J} & n_{2+} \\ \mbox{X} & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mbox{} & I & n_{I1} & n_{I2} & \cdots & n_{IJ} & n_{2+} \\\hline & \mbox{Total} & n_{+1} & n_{+2} & \cdots & n_{+J} & n \\\hline \end{array} \]
Na discussão a seguir, examinamos dois dos modelos que são normalmente usados para representar contagens em uma tabela de contingência.
3.2.1 Uma distribuição multinomial
Primeiro, considere a situação em que um tamanho de amostra fixo de \(n\) unidades é amostrado de uma grande população. Definimos \(\pi_{ij} = P(X = i; Y = j)\). Esta é a mesma configuração multinomial da Seção 3.1 com \(n\) tentativas e \(I\times J\) respostas possíveis, exceto que agora cada resposta possível é uma combinação de duas variáveis em vez de apenas uma. Assumimos que cada unidade amostrada tem uma e apenas uma combinação de categorias \(X\) e \(Y\), de modo que \(\sum_{i=1}^I \sum_{j=1}^J \pi_{ij}=1\).
O estimador de máxima verossimilhança para essas probabilidades procede de forma análoga ao que foi feito na Seção 3.1, com pequenas modificações na notação para levar em conta a segunda variável. A função de probabilidades para \(N_{11},\cdots,N_{IJ}\) torna-se \[ P(N_{11}=n_{11},\cdots,N_{IJ}=n_{IJ})=\dfrac{n!}{\displaystyle \prod_{i=1}^I \prod_{j=1}^J n_{ij}!}\prod_{i=1}^I \prod_{j=1}^J \pi_{ij}^{n_{ij}}, \] que também é a função de verossimilhança para uma amostra de tamanho \(n\). O MLE de \(\pi_{ij}\) é a proporção estimada \(\widehat{\pi}_{ij} = n_{ij}/n\).
Distribuições marginais para \(X\) e para \(Y\) também podem ser encontradas. A probabilidade marginal para o nível \(i\) de \(X\) é \(\pi_{i+} = P(X = i)\) para \(i = 1,\cdots,I\). A distribuição marginal de \(X\) é, portanto, multinomial com \(n\) tentativas e probabilidades \(\pi_{1+},\cdots,\pi_{I+}\) e resulta nas contagens marginais \(n_{1+},\cdots,n_{I+}\). Da mesma forma, a probabilidade marginal para a categoria \(j\) de \(Y\) é \(\pi_{+j}=\sum_{i=1}^I \pi_{ij}\) para \(j=1,\cdots,J\).
A distribuição marginal de \(Y\) é multinomial com \(n\) tentativas, probabilidades \(\pi_{+1},\cdots,\pi_{+J}\), resultando nas contagens marginais \(n_{+1},\cdots,n_{+J}\). Observe que \(\sum_{i=1}^I \pi_{i+}=1\) e \(\sum_{j=1}^J \pi_{+j} = 1\). Os estimadores de máxima verossimilhança (MLEs) de \(\pi_{i+}\) e \(\pi_{+j}\) são as proporções de linha e coluna correspondentes, \(\widehat{\pi}_{i+} = n_{i+}/n\) e \(\widehat{\pi}_{+j} = n_{+j}/n\), respectivamente.
Quando o resultado de \(X\) não afeta as probabilidades dos resultados de \(Y\) , dizemos que \(Y\) é independente de \(X\). Como resultado da independência, a probabilidade de qualquer resultado conjunto de fatores \((X = i, Y = j)\) nas probabilidades marginais para \(X = i\) e \(Y = j\): \(\pi_{ij} = \pi_{i+}\pi_{+j}\).
A independência simplifica a estrutura das probabilidades dentro de uma tabela de contingência, reduzindo o número de parâmetros desconhecidos para \((I-1) + (J-1) = I + J - 2\) probabilidades marginais. Observe que as partes “-1” ocorrem devido às restrições \(\sum_{i=1}^I \pi_{i+}=1\) e \(\pi_{+j}=\sum_{j=1}^J \pi_{j+}=1\) nas probabilidades. Sem independência, existem \(IJ - 1\) parâmetros de probabilidade desconhecidos, onde a parte “-1” ocorre devido à restrição \(\sum_{i=1}^I\sum_{j=1}^J \pi_{ij}=1\) nas probabilidades. Assim, há uma redução de \((IJ- 1) - (I + J - 2) = (I - 1)(J- 1)\) parâmetros em comparação com a mesma tabela sem independência. Examinaremos em breve como realizar um teste de hipótese para independência.
3.2.2 \(I\) distribuições multinomiais
Um modelo alternativo é necessário quando amostras de tamanhos \(n_{i+}\), \(i = 1,\cdots,I\) sou deliberadamente tiradas de cada um dos diferentes grupos. Neste caso, as contagens marginais \(n_{i+}\) são fixas por planejamento, então temos uma distribuição multinomial de \(J\) categorias separada em cada um dos \(I\) grupos, onde \(n =\sum_{i=1}^I n_{i+}\). Cada uma dessas distribuições tem seu próprio conjunto de parâmetros de probabilidade.
Definindo \(P(Y = j \, | \, X = i) = \pi_{j\,|\,i}\) como a probabilidade condicional de observar a categoria de resposta \(j\) dado que uma unidade é do grupo \(i\). Observe que \(\sum_{j=1}^J \pi_{j \,|\, i}=1\) para cada \(i = 1,\cdots,I\). A distribuição conjunta condicional de \(N_{i1},\cdots,N_{iJ}\) tem função de probabilidade \[ P(N_{i1}=n_{i1},\cdots,N_{iJ}=n_{iJ} \, | \, N_{i+}=n_{i+}) = \dfrac{n_{i+!}}{\displaystyle \prod_{j=1}^ J n_{ij}!}\prod_{j=1}^J \pi_{j\ |\, i}^{n_{ij}}, \]
para cada \(i = 1,\cdots,I\).
Assumindo que \(I\) amostras diferentes são tomadas independentemente umas das outras, a verossimilhança para os parâmetros é apenas o produto de \(I\) distribuições multinomiais, \[ \prod_{i=1}^I \dfrac{n_{i+!}}{\displaystyle \prod_{j=1}^ J n_{ij}!}\prod_{j=1}^J \pi_{j\ |\, i}^{n_{ij}}\cdot \]
Como resultado, este modelo é muitas vezes referido como o modelo produto de multinomiais. O estimador de máxima verossimilhança (MLE) de \(\pi_{j\ |\, i}\) é \(\widehat{\pi}_{j\ |\, i} = n_{ij}/n_{i+}\).
Observe que as mesmas estimativas resultam da aplicação da definição de probabilidade condicional, \[ P(Y = j \, | \, X = i) = P(X = i, Y = j)/P(X = i), \]
às estimativas de um modelo multinomial: \(\widehat{\pi}_{j\ |\, i} = \widehat{\pi}_{ij}/\widehat{\pi}_{i+} = (n_{ij}/n)/(n_{i+}/n) = n_{ij}/n_{i+}\).
A independência de \(X\) e \(Y\) no contexto de um modelo produto de multinomiais significa que as probabilidades condicionais para cada \(Y\) são iguais nas linhas da tabela. Ou seja, para cada \(j\), \(\pi_{j\,|\,1} = \cdots =\pi_{j\,|\,I} = \pi_{+j}\). Observe que essa condição é matematicamente equivalente à independência conforme definido para o modelo multinomial. Para mostrar essa equivalência, novamente precisamos usar a definição de probabilidade condicional. Como \(X\) não é aleatório no modelo produto de multinomiais, definimos \(\pi_{i+}\) como a proporção fixa da amostra total que é retirada do grupo \(i\). Então \(\pi_{j\,|\,i} =\pi_{ij}/\pi_{i+}\) e \(\pi_{j\,|\,i} =\pi_{+j}\) juntos implicam que \(\pi_{ij} =\pi_{i+}\pi_{+j}\).
Vimos agora que as estimativas de parâmetros dos modelos multinomiais um e produto são as mesmas, as definições de independência nos dois modelos são equivalentes e os dois modelos levam exatamente às mesmas distribuições condicionais para \(Y\) dado \(X = i\). Como consequência, as análises realizadas com base em cada modelo geralmente produzem os mesmos resultados. Portanto, ao desenvolver testes de independência e outras análises em tabelas de contingência, assumimos o modelo para a tabela que for mais conveniente.
Exemplo 3.2: Amostra simulada multinomial.
Simulamos amostras aqui para modelos multinomiais e produto de multinomiais para ajudar os leitores a entender as relações entre eles. Os dados são simulados usando rmultinom() como na Seção 3.1. A principal diferença está em como definimos probabilidades para corresponder a células de uma \(I\times J\) tabela de contingência.
Abaixo está como simulamos uma amostra de tamanho \(n = 1000\) para uma tabela de contingência \(2\times 3\) correspondente ao modelo multinomial onde \(\pi_{11} = 0.2\), \(\pi_{21} = 0.3\), \(\pi_{12} = 0.2\), \(\pi_{22} = 0.1\), \(\pi_{13} = 0.1\) e \(\pi_{23} = 0.1\):
# Probabilities entered by column for array ()
pi.ij <- c(0.2 , 0.3 , 0.2 , 0.1 , 0.1 , 0.1)
pi.table <- array ( data = pi.ij , dim = c(2 ,3) , dimnames = list (X = 1:2 , Y = 1:3) )
pi.table # pi_ij## Y
## X 1 2 3
## 1 0.2 0.2 0.1
## 2 0.3 0.1 0.1set.seed (9812)
save <- rmultinom (n = 1, size = 1000 , prob = pi.ij)
c.table1 <- array ( data = save , dim = c(2 ,3) , dimnames = list (X = 1:2 , Y = 1:3) )
c.table1## Y
## X 1 2 3
## 1 191 206 94
## 2 311 95 103c.table1 /sum(c.table1 )## Y
## X 1 2 3
## 1 0.191 0.206 0.094
## 2 0.311 0.095 0.103Por exemplo, \(n_{11} = 191\) e \(n_{23} = 103\). Vemos que cada \(\widehat{\pi}_{ij}\) é bastante semelhante ao seu \(\pi_{ij}\) correspondente, como seria esperado para uma amostra tão grande. Além disso, observe que os totais marginais das linhas variam, \(n_{1+} = 491\) e \(n_{2+} = 509\) e estão muito próximos do que seria esperado considerando que \(\pi_{1+} = \pi_{2+} = 0.5\).
Para a configuração de \(I\) multinomiais, simulamos novamente uma amostra para uma tabela de contingência \(2\times 3\), portanto, \(I = 2\). Com este modelo, precisamos desenhar amostras de tamanho fixo separadamente para cada linha. Também precisamos reexpressar as probabilidades das células como probabilidades condicionais que somam 1 em cada linha, o que resulta em \(\pi_{1\,|\,1} = 0.4\); \(\pi_{2\,|\,1} = 0.4\), \(\pi_{1\,|\,2} = 0.6\), \(\pi_{2\,|\,2} = 0.2\) e \(\pi_{3\,|\,2} = 0.2\).
Ao contrário do modelo multinomial, os totais das linhas não são especificados pelas probabilidades na tabela, então podemos escolhê-los de acordo com as necessidades do problema.
Selecionamos arbitrariamente \(n_{1+} = 400\) e \(n_{2+} = 600\). Duas chamadas separadas para a função rmultinom() são usadas para gerar amostras independentemente das duas distribuições de probabilidade separadas:
pi.cond <- pi.table / rowSums (pi.table )
pi.cond # pi_j |i## Y
## X 1 2 3
## 1 0.4 0.4 0.2
## 2 0.6 0.2 0.2set.seed (8111)
save1 <- rmultinom (n = 1, size = 400 , prob = pi.cond [1 ,])
save2 <- rmultinom (n = 1, size = 600 , prob = pi.cond [2 ,])
c.table2 <- array ( data = c( save1 [1] , save2 [1] , save1 [2] ,
save2 [2] , save1 [3] , save2 [3]) , dim = c(2 ,3) , dimnames = list (X = 1:2 , Y = 1:3) )
c.table2## Y
## X 1 2 3
## 1 162 159 79
## 2 351 126 123rowSums (c.table2 )## 1 2
## 400 600c.table2 / rowSums (c.table2 ) # Estimate of pi_j |i## Y
## X 1 2 3
## 1 0.405 0.3975 0.1975
## 2 0.585 0.2100 0.2050round (c.table1 / rowSums (c.table1 ) ,4) # From 1 multinomial## Y
## X 1 2 3
## 1 0.389 0.4196 0.1914
## 2 0.611 0.1866 0.2024Vemos novamente que cada \(\widehat{\pi}_{j\,|\,i}\) é bastante semelhante ao seu respectivo \(\pi_{j|,|\,i}\), como seria de esperar.
Se quiséssemos simular os dados sob independência, poderíamos encontrar \(\pi_{ij}\)’s que satisfaçam \(\pi_{ij} =\pi_{i+}\pi_{+j}\) para cada par \(i, j\) e usar a função rmultinom() para simular uma amostra de uma distribuição multinomial de 6 categorias. Equivalentemente, poderíamos encontrar \(\pi_{+j}\) para \(j = 1\), 2 e 3, e simular dados de duas distribuições multinomiais de três categorias usando duas chamadas para a função rmultinom(). Demonstraremos este processo a seguir:
# Simulate n_ij under independence for a 2x3 contingency table
# 1 multinomial under independence
pi.i<-rowSums(pi.table)
pi.j<-colSums(pi.table)
pi.ij.ind<-pi.i%o%pi.j # Quick way to find pi_i+ * pi_+j
pi.ij.ind## 1 2 3
## 1 0.25 0.15 0.1
## 2 0.25 0.15 0.1set.seed(9218)
save.ind<-rmultinom(n = 1, size = 1000, prob = pi.ij.ind)
save.ind## [,1]
## [1,] 251
## [2,] 245
## [3,] 141
## [4,] 153
## [5,] 104
## [6,] 106c.table1.ind<-array(data = save.ind, dim = c(2,3), dimnames = list(X = 1:2, Y = 1:3))
c.table1.ind/sum(c.table1.ind) # pi^_ij is similar to pi^_i+ * pi^_+j## Y
## X 1 2 3
## 1 0.251 0.141 0.104
## 2 0.245 0.153 0.106# Using methods found later in the chapter
chisq.test(x = c.table1.ind, correct = FALSE) # Do not reject independence##
## Pearson's Chi-squared test
##
## data: c.table1.ind
## X-squared = 0.51746, df = 2, p-value = 0.772# I multinomials under independence
set.seed(7718)
save1.ind<-rmultinom(n = 1, size = 400, prob = pi.j)
save2.ind<-rmultinom(n = 1, size = 600, prob = pi.j)
c.table2.ind<-array(data = c(save1.ind[1], save2.ind[1], save1.ind[2], save2.ind[2],
save1.ind[3], save2.ind[3]), dim = c(2,3), dimnames = list(X = 1:2, Y = 1:3))
c.table2.ind## Y
## X 1 2 3
## 1 194 124 82
## 2 299 175 126rowSums(c.table2.ind)## 1 2
## 400 600c.table2.ind/rowSums(c.table2.ind) # pi^_j|1 is similar to pi^_j|2## Y
## X 1 2 3
## 1 0.4850000 0.3100000 0.205
## 2 0.4983333 0.2916667 0.210chisq.test(x = c.table2.ind, correct = FALSE) # Do not reject independence##
## Pearson's Chi-squared test
##
## data: c.table2.ind
## X-squared = 0.38518, df = 2, p-value = 0.82483.2.3 Tetes de independência
Um teste de independência, \[ H_0 \, : \, \pi_{ij}=\pi_{i+}\pi_{j+} \qquad \mbox{para cada} \; i,j \\ H_1 \, : \, \pi_{ij}\neq \pi_{i+}\pi_{j+} \qquad \mbox{para algum} \; i,j \] pode ser realizado usando um teste qui-quadrado de Pearson e um teste da razão de verossimilhanças (LRT).
Esses testes já foram mostrados na Seção 1.2.3 como testes para a igualdade de duas probabilidades de sucesso binomial, o que equivale a um teste de independência em um modelo produto de multinomiais com \(I = J = 2\). A estatística do teste qui-quadrado de Pearson é novamente formado pela soma \[ \dfrac{\mbox{(contagem observada - contagem esperada estimada)}^2}{\mbox{(contagem esperada estimada)}} \]
em todas as células da tabela de contingência.
A contagem observada é \(n_{ij}\). A contagem esperada sob independência é \(n\widehat{\pi}_{i+}\widehat{\pi}_{+j} = n_{i+}n_{+j}/n\) para o modelo multinomial, ou equivalentemente \(n_{i+}\widehat{\pi}_{+j} = n_{i+}n_{+j}/n\) sob o modelo produto de multinomiais. Isso leva à estatística de teste \[ X^2= \sum_{i=1}^I \sum_{j=1}^J \dfrac{\big(n_{ij}-n_{i+}n_{+j}/n\big)^2}{n_{i+}n_{+j}/n}\cdot \]
A estatística \(X^2\) tem uma distribuição aproximada de \(\chi^2_{(I-1)(J-1)}\) em grandes amostras quando a hipótese nula é verdadeira.
Quando a hipótese nula é falsa, esperamos grandes desvios entre as contagens observadas e esperadas em relação ao tamanho de \(n_{i+}n_{+j}/n\), que levam a grandes valores da estatística \(X^2\). Portanto, rejeitamos a hipótese nula de independência entre \(X\) e \(Y\) quando \(X^2>\chi^2_{(I-1)(J-1),1-\alpha}\).
A razão de verossimilhança é formada da maneira usual 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 \]
O cálculo de \(\Lambda\) é baseado em \(\widehat{\pi}_{ij}=\widehat{\pi}_{i+}\widehat{\pi}_{+j} = n_{i+}n_{+j}/n\) no numerador e \(\widehat{\pi}_{ij} = n_{ij}/n\) no denominador para estimar cada \(\pi_{ij}\). Aplicando a transformação usual com \(\Lambda\), temos \[ -2\log(\Lambda) = 2\sum_{i=1}^I \sum_{j=1}^J n_{ij}\left( \dfrac{n_{ij}}{n_{i+}n_{+j}/n}\right), \]
onde assumimos \(0\times \log(0) = 0\) por convenção. Assim como \(X^2\), a estatística LRT transformada tem uma distribuição amostral em amostras grandes \(\chi^2_{(I-1)(J-1)}\) quando \(H_0\) é verdadeiro e usa a mesma regra de rejeição.
Os graus de liberdade usados para o teste de independência são encontrados calculando \[ (\mbox{Número de parâmetros sob } H_1) - (\mbox{Número de parâmetros sob } H_0)\cdot \]
Esta é uma maneira geral de encontrar graus de liberdade para qualquer teste de comparação de modelos. Conforme mostrado na Seção 3.2.1, precisamos estimar os \(I-J-2\) parâmetros quando a hipótese nula de independência é válida e \(IJ-1\) quando não. Assim, os graus de liberdade para o teste de independência são \((IJ - 1) - (I + J - 2) = (I - 1)(J - 1)\).
Tanto o teste LRT quanto o teste qui-quadrado de Pearson geralmente fornecem resultados semelhantes em amostras grandes. No entanto, às vezes, seus valores podem diferir consideravelmente em amostras menores, levando a ambiguidade se seus valores estiverem em lados opostos da região de rejeição. Tem havido uma série de recomendações sobre o que constitui uma amostra “grande o suficiente” para obter uma boa aproximação \(\chi^2\).
Os critérios mais comuns são \(n_{i+}n_{+j}/n > 1\) ou > 5 para todas as células da tabela de contingência. Esses critérios podem não ser atendidos quando há contagens de células muito pequenas em muitas células da tabela. Por exemplo, uma linha para a qual a contagem marginal \(n_{i+}\) não é muito maior do que o número de colunas não pode ter contagens de células esperadas “grandes” em todas as colunas. Nesses casos, uma simulação de Monte Carlo ou os procedimentos de teste exatos descritos na Seção 6.2 podem fornecer uma avaliação visual para determinar se a aproximação à distribuição \(\chi^2_{(I-1)(J-1)}\) é apropriada. Esta é a nossa abordagem preferida sempre que houver alguma dúvida em relação à aproximação \(\chi^2\) e fornecemos um exemplo de sua implementação.
Exemplo 3.3: Crackers enriquecidos com fibra.
A fibra dietética é um composto saudável encontrado em muitos vegetais e grãos. Alimentos altamente processados são pobres em fibras, por isso às vezes é adicionado a esses alimentos para torná-los mais nutritivos. A Data and Story Library (DASL) descreve os resultados de um estudo em que os indivíduos receberam um novo tipo de biscoito enriquecido com fibras.
Os participantes comeram os biscoitos e depois uma refeição. Pouco depois, os participantes foram instruídos a descrever qualquer inchaço que experimentassem. A tabela abaixo fornece os dados resultantes em uma tabela de contingência \(4\times 4\). Havia quatro tipos diferentes de biscoitos com base em sua fonte de fibra e havia quatro níveis de gravidade do inchaço relatados.
\[ \mbox{Severidade do inchaço depois de comer um biscoito enriquecido} \\ \mbox{com fibras; a fonte de dados é DASL.} \\ \mbox{(http://lib.stat.cmu.edu/DASL/Stories/HighFiberDietPlan.html).}\\ \qquad \qquad \qquad \qquad \qquad \qquad \mbox{Gravidade do inchaço} \\ \begin{array}{cccccc}\hline & & \mbox{Nenhum} & \mbox{Baixo} & \mbox{Médio} & \mbox{Alto} \\ \mbox{Fonte de fibra} & \mbox{Nenhum} & 6 & 4 & 2 & 0 \\ & \mbox{Farelo} & 7 & 4 & 1 & 0 \\ & \mbox{Chiclete} & 2 & 2 & 3 & 5 \\ & \mbox{Ambos} & 2 & 5 & 3 & 2 \\\hline \end{array} \]
Embora não especificado na descrição dos dados na DASL, assumimos que cada participante foi atribuído a apenas um grupo de fonte de fibra. Além disso, a descrição dos dados não menciona se o fato de haver 12 participantes para cada fonte de fibra foi por design ou por acaso, portanto, não está claro se um modelo multinomial ou produto de multinomiais pode ser mais apropriado. Felizmente, conforme observado no final da Seção 3.2.2, essa distinção não tem relação com os resultados da análise.
Finalmente, enquanto os dados fornecidos na DASL estavam em uma tabela de contingência \(4\times 4\), podemos ver que a fonte de fibra pode ser dividida em duas variáveis explicativas separadas, uma indicando se o farelo (bran) está presente e a outra indicando se a chiclete (gum) está presente. Exploramos métodos alternativos de análise que levam em conta essa estrutura “fatorial” na Seção 3.3.
Começamos lendo os dados em R e resumindo os dados em uma forma de tabela de contingência usando a função xtabs():
diet <- read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/Fiber.csv")
head(diet)## fiber bloat count
## 1 bran high 0
## 2 gum high 5
## 3 both high 2
## 4 none high 0
## 5 bran medium 1
## 6 gum medium 3diet$fiber <- factor (x = diet$fiber , levels = c("none", "bran", "gum", "both"))
diet$bloat <- factor (x = diet$bloat , levels = c("none", "low", "medium", "high"))
diet.table <- xtabs ( formula = count ~ fiber + bloat , data = diet )
diet.table## bloat
## fiber none low medium high
## none 6 4 2 0
## bran 7 4 1 0
## gum 2 2 3 5
## both 2 5 3 2Usamos a função factor() para alterar a ordem dos níveis de ambas as variáveis para corresponder ao dado na tabela acima. Em vez disso, poderíamos ter inserido a tabela de contingência diretamente no R usando a função array() e organizado o array na ordem desejada.
Gostaríamos de determinar se a gravidade do inchaço está relacionada ao tipo de fibra. Se houver uma relação, isso indicaria ao fabricante do biscoito que alguns tipos de fibra podem causar inchaço mais grave e devem ser evitados. Usamos um teste de independência para fazer essa determinação. As hipóteses para o teste são \[ H_0 \, : \, \pi_{ij} = \pi_{i+}\pi_{+j} \qquad \mbox{para cada} \; i,j \\ vs.\\ H_1 \, : \, \pi_{ij} \neq \pi_{i+}\pi_{+j} \qquad \mbox{para algum} \; i, j, \]
onde \(\pi_{ij}\) é a probabilidade de que uma pessoa selecionada aleatoriamente seja designada para a fonte de fibra \(i\) e experimente o nível de inchaço \(j\).
O R fornece várias funções para testar essas hipóteses, incluindo:
chisq.test(),
assocstats() do pacote vcd e
a função genérica summary() para resumir um objeto de tabela.
Normalmente, usaríamos apenas uma dessas funções, mas mostramos todas as três para fins ilustrativos:
ind.test <- chisq.test (x = diet.table , correct = FALSE )
ind.test##
## Pearson's Chi-squared test
##
## data: diet.table
## X-squared = 16.943, df = 9, p-value = 0.04962library ( package = vcd)
assocstats (x = diet.table )## X^2 df P(> X^2)
## Likelihood Ratio 18.880 9 0.026230
## Pearson 16.943 9 0.049621
##
## Phi-Coefficient : NA
## Contingency Coeff.: 0.511
## Cramer's V : 0.343class ( diet.table )## [1] "xtabs" "table"summary ( diet.table )## Call: xtabs(formula = count ~ fiber + bloat, data = diet)
## Number of cases in table: 48
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 16.943, df = 9, p-value = 0.04962
## Chi-squared approximation may be incorrectqchisq (p = 0.95 , df = 9)## [1] 16.91898Observe que incluímos o valor do argumento correct = FALSE em chisq.test() para evitar que uma correção de continuidade seja aplicada. Consulte a Seção 1.2.3 para saber por que evitamos essas correções. A estatística do teste qui-quadrado de Pearson é \(X^2 = 16.94\), o valor crítico no nível = 0.05 é \(\chi^2_{0.95;9} = 16.92\), e o \(p\)-valor é \(P(A > 16.94 ) = 0.0496\) onde \(A\sim \chi^2_9\).
Como o \(p\)-valor é pequeno, mas não extremamente, diríamos que há evidência moderada contra a independência, portanto, evidência moderada de dependência. O LRT leva a uma conclusão semelhante com \(-2\log(\Lambda ) = 18.88\) e um \(p\)-valor de 0.0262.
Há várias contagens de células muito pequenas na tabela de contingência, portanto, pode haver alguma preocupação sobre a adequação da aproximação de \(\chi^2_9\). As contagens de células esperadas sob independência para esta tabela são:
ind.test$expected## bloat
## fiber none low medium high
## none 4.25 3.75 2.25 1.75
## bran 4.25 3.75 2.25 1.75
## gum 4.25 3.75 2.25 1.75
## both 4.25 3.75 2.25 1.75# Calculate -2log(Lambda) long way
# Because log(0) is undefined, a small adjustment is made to these cells
diet.table.adj<-diet.table + ifelse(test = diet.table == 0, yes = 0.01, no = 0)
# Because 0*log(0.01) = 0, the adjustment does not affect the final outcome
test.stat<- 2*sum( diet.table * log(diet.table.adj/ind.test$expected) )
test.stat## [1] 18.88033qchisq(p = 0.95, df = 9)## [1] 16.91898ind.test$residuals # Pearson residuals easier way## bloat
## fiber none low medium high
## none 0.8488747 0.1290994 -0.1666667 -1.3228757
## bran 1.3339459 0.1290994 -0.8333333 -1.3228757
## gum -1.0914103 -0.9036961 0.5000000 2.4567691
## both -1.0914103 0.6454972 0.5000000 0.1889822abs(ind.test$residuals)>qnorm(p = 0.995)## bloat
## fiber none low medium high
## none FALSE FALSE FALSE FALSE
## bran FALSE FALSE FALSE FALSE
## gum FALSE FALSE FALSE FALSE
## both FALSE FALSE FALSE FALSEc.table<-array(data = c(0, 5, 2, 0,
1, 3, 3, 2,
4, 2, 5, 4,
7, 2, 2, 6),
dim=c(4,4), dimnames = list(fiber = c("bran", "gum", "combo", "control"),
bloat = c("high", "medium", "low", "none")))
c.table## bloat
## fiber high medium low none
## bran 0 1 4 7
## gum 5 3 2 2
## combo 2 3 5 2
## control 0 2 4 6ind.test<-chisq.test(c.table, correct = FALSE)
ind.test##
## Pearson's Chi-squared test
##
## data: c.table
## X-squared = 16.943, df = 9, p-value = 0.04962que satisfazem o critério \(n_{i+}n_{+j}/n > 1\) dado anteriormente, mas não \(n_{i+}n_{+j}/n > 5\).
As funções chisq.test() e summary() até imprimem uma mensagem de aviso sobre a adequação da aproximação \(\chi^2_9\). Este aviso é dado pelas funções sempre que \(n_{i+}n_{+j}/n < 5\) para qualquer célula.
Uma maneira de examinar a aproximação distribucional é através de uma simulação de Monte Carlo. Isso envolve a simulação de um grande número de tabelas de contingência de tamanho \(n\) cujas probabilidades satisfazem a hipótese nula de independência. Idealmente, deve-se seguir o mesmo tipo de método de amostragem do modelo assumido. Para um modelo multinomial, podemos usar uma chamada para rmultinom(), onde definimos cada \(\pi_{ij}\), argumento prob, para \(n_{i+}n_{+j}/n\), a contagem esperada estimada sob a hipótese nula dividida pelo tamanho da amostra.
Para um modelo produto de multinomiais, podemos usar \(I\) chamadas diferentes para rmultinom() definindo \(\pi_{j\,|\,i}\) igual a \(\widehat{\pi}_{+j}\). Optamos por usar o modelo multinomial aqui devido à sua simplicidade e à natureza desconhecida do método de amostragem real.
#Monte Carlo simulation
n<-sum(diet.table) # Sample size
pi.hat<-ind.test$expected/n # Estimated expected pi^_ij under Ho
pi.hat## bloat
## fiber high medium low none
## bran 0.03645833 0.046875 0.078125 0.08854167
## gum 0.03645833 0.046875 0.078125 0.08854167
## combo 0.03645833 0.046875 0.078125 0.08854167
## control 0.03645833 0.046875 0.078125 0.08854167c(pi.hat) # Notice the order here - need to know later for calc.stat()## [1] 0.03645833 0.03645833 0.03645833 0.03645833 0.04687500 0.04687500
## [7] 0.04687500 0.04687500 0.07812500 0.07812500 0.07812500 0.07812500
## [13] 0.08854167 0.08854167 0.08854167 0.08854167sum(pi.hat) # Double check that it adds to 1## [1] 1B<-10000 # Number of simulated data sets
set.seed(1298)
table.count<-rmultinom(n = B, size = n, prob = pi.hat)
# Each column contains the counts for one contingency table
table.count[,1:2] # First two sets of contingency table counts## [,1] [,2]
## [1,] 3 2
## [2,] 2 1
## [3,] 0 1
## [4,] 2 1
## [5,] 1 1
## [6,] 2 3
## [7,] 3 2
## [8,] 2 3
## [9,] 5 3
## [10,] 2 3
## [11,] 4 4
## [12,] 6 2
## [13,] 3 9
## [14,] 3 5
## [15,] 9 5
## [16,] 1 3rowMeans(table.count)/n # Close to c(pi.hat), which shows the data simulated is o.k.## [1] 0.03646458 0.03639167 0.03612917 0.03608958 0.04712083 0.04693542
## [7] 0.04673958 0.04639375 0.07812292 0.07775000 0.07898542 0.07841667
## [13] 0.08834583 0.08860833 0.08909167 0.08841458# Function calculates X^2 for a data set and checks the contingency table size
calc.stat<-function(data) {
c.table<-array(data = data, dim=c(4,4), dimnames = list(fiber =
c("bran", "gum", "combo", "control"), bloat = c("high", "medium", "low", "none")))
ind.test<-chisq.test(c.table, correct = FALSE)
ck.0.row<-sum(rowSums(c.table) == 0) # Number of rows with all 0 entries
ck.0.col<-sum(colSums(c.table) == 0) # Number of columns with all 0 entries
c(ind.test$statistic, ck.0.row, ck.0.col)
}
# Example application of the function (could do this with the observed data too)
calc.stat(data = table.count[,1])## X-squared
## 11.31912 0.00000 0.00000# Calculate X^2* for each simulated data set
save.star<-apply(X = table.count, MARGIN = 2, FUN = calc.stat)
# The warnings are for "Chi-squared approximation may be incorrect"
# We do not want to include contingency tables that have an all 0 row or all 0 column
# because this changes the contingency table size (and degrees of freedom for the chi^2 distribution)
sum(save.star[2,])## [1] 0sum(save.star[3,])## [1] 6# The non 4x4 table
c.table<-array(data = table.count[,691], dim=c(4,4), dimnames =
list(fiber = c("bran", "gum", "combo", "control"), bloat = c("high", "medium", "low", "none")))
# Just obtain the 4x4 tables
X.sq.star<-save.star[1, save.star[2,] == 0 & save.star[3,] == 0]
# Quantiles
p<-c(0.01, 0.05, 0.10)
quantile(x = X.sq.star, probs = 1 - p, type = 1)## 99% 95% 90%
## 21.17181 16.67327 14.62051qchisq(p = 1 - p, df = 9)## [1] 21.66599 16.91898 14.68366A estatística do teste qui-quadrado de Pearson, digamos \(X^{2*}\), é calculada para cada tabela de contingência de contagens simuladas. Esses valores \(X^{2*}\) podem ser resumidos por meio de um gráfico de várias maneiras, como um histograma, um gráfico empírico da função de distribuição e um gráfico QQ-plot, a fim de visualizar uma distribuição estimada para a estatística. Usando uma sobreposição da \(\chi^2_9\) ou examinando o desvio de uma linha reta no gráfico QQ-plot, pode-se julgar se a aproximação de distribuição é apropriada. O mesmo procedimento também pode ser aplicado à estatística do LRT.
# Plots
par(mfrow=c(1,2))
# Histogram
df<-9
hist(x = X.sq.star, main = "Histogram", freq = FALSE, xlab = expression(X^"2*"))
curve(expr = dchisq(x = x, df = df), col = "red", add = TRUE)
box()
grid()
# segments(x0 = ind.test$statistic, y0 = -10, x1 = ind.test$statistic, y1 = 10)
# Compare CDFs
plot.ecdf(x = X.sq.star, verticals = TRUE, do.p = FALSE, main = "CDFs", lwd = 2, col = "black",
xlab = expression(X^"2*"), panel.first = grid(col="gray", lty="dotted"), ylab = "CDF")
curve(expr = pchisq(q = x, df = df), col = "red", add = TRUE, lwd = 1)
legend(x = 15, y = 0.4, legend = c(expression(X^"2*"), expression(chi[9]^2)),
lwd = c(2,1), col = c("black", "red"), bty = "n")
# QQ-Plot
par(mfrow=c(1,1))
B.adj<-length(X.sq.star)
chi.quant<-qchisq(p = seq(from = 1/(B.adj+1), to = 1-1/(B.adj+1), by = 1/(B.adj+1)), df = df)
plot(x = sort(X.sq.star), y = chi.quant, main = "QQ-Plot", xlab = expression(X^"2*"),
ylab = expression(paste(chi[9]^2, "quantile")))
abline(a = 0, b = 1)
grid()
# Bootstrap p-value
mean(X.sq.star > ind.test$statistic) # bootstrap p-value## [1] 0.04592756Para uma amostra de observações \(w_1,w_2,\cdots,w_m\), a função de distribuição empírica em \(w\) é a proporção de observações em ou abaixo de \(w\): \(\widehat{F}(w) = \#\{w_j \leq w\}/m\), onde “#” significa “número de”. Os valores observados ordenados de uma estatística são plotados em relação aos quantis correspondentes de uma distribuição de probabilidade particular de interesse. Por exemplo, a mediana observada, desde que o número de valores observados seja ímpar, seria plotada em relação ao quantil 0.5 da distribuição de probabilidade. Se os pontos plotados seguem próximos a uma linha reta de \(y = x\), onde \(y(x)\) indica o que é plotado no eixo \(y\) (eixo \(x\)), a estatística segue aproximadamente a distribuição de probabilidade.
A simulação de Monte Carlo usando 10.000 conjuntos de dados simulados é implementada com código acima. A figura apresentada mostra o histograma e o gráfico da função de distriuição empírica com a sobreposição da distribuição \(\chi^2_9\). Vemos uma concordância geral em ambos os gráficos entre a distribuição de \(X^{2*}\) e a distribuição \(\chi^2_9\). Os quantis em \(\alpha=\) 0.90, 0.95 e 0.99 são 14.57, 16.61 e 21.04 para \(X^{2*}\) e 14.68, 16.92 e 21.67 para \(\chi^2_9\), respectivamente. Assim, há uma concordância geral entre as duas distribuições.
É interessante notar que a simulação de Monte Carlo realizada aqui é na verdade uma forma de bootstrap. O bootstrap é aplicado aqui como um procedimento não paramétrico que fornece estimativas da função de distribuição para uma estatística. Para este exemplo, podemos obter uma estimativa bootstrap do \(p\)-valor para o teste calculando a proporção de valores \(X^{2*}\) simulados maiores ou iguais a \(X^2 = 16.94\). O resultado é 0.0446, o que está de acordo com o \(p\)-valor do teste qui-quadrado de Pearson. Para uma discussão mais aprofundada do bootstrap, veja Davison and Hinkley (1997).
Quando a independência é rejeitada, existem várias maneiras de examinar a associação entre \(X\) e \(Y\). Na Seção 1.2 com tabelas de contingência \(2\times 2\), poderíamos simplesmente calcular uma razão de chances para resumir o grau de associação; ou seja, a quantidade pela qual a razão de chances é > 1 ou < 1. Agora, com tabelas de contingência \(I\times J\), uma razão de chances não resume a associação. Em vez disso, uma abordagem é calcular as razões de chances para \(2\times 2\) partes de uma tabela de contingência que são de interesse especial. A menos que haja uma orientação firme do contexto do problema, há potencial para criar muitas razões de chances diferentes e, portanto, potencial para taxas de erro infladas para vários testes ou intervalos de confiança.
Como alternativa às razões de chances, várias formas de resíduos podem ser calculadas em cada célula para medir até que ponto a contagem de células observada se desvia do que seria esperado sob independência. O resíduo básico na célula \((i, j)\) tem a forma \(n_{ij}-n_{i+}n_{+j}/n\). Uma versão diferente do resíduo, o resíduo padronizado de Pearson, é \[ \dfrac{n_{ij}-n_{i+}n_{+j}/n}{\sqrt{(1-n_{i+}/n)(1-n_{+j}/n) \, n_{i+}n_{+j}/n}}, \] onde o denominador é \(\sqrt{\widehat{\mbox{Var}}(N_{ij}-N_{i+}N_{+j}/N)}\).
O componente stdres de um objeto resultante de chisq.test() fornece esses resíduos Pearson padronizados. A vantagem de um resíduo de Pearson padronizado é que sua distribuição pode ser aproximada pelo padrão normal. Isso significa que valores além de, digamos, $$2 podem ser considerados “incomuns” e devem ocorrer apenas em cerca de 5% das células. As células onde ocorrem grandes resíduos de Pearson padronizados indicam fontes potenciais da associação entre \(X\) e \(Y\).
Resíduos e odds ratios são um tanto ineficazes em tabelas onde \(I\) e/ou \(J\) são muito maiores que 2, porque não levam facilmente a uma compreensão sistemática da relação entre \(X\) e \(Y\). No entanto, exibir os resíduos nas mesmas linhas e colunas da tabela de contingência às vezes pode revelar padrões de valores positivos ou negativos que sugerem uma causa maior para a associação. Por exemplo, se os resíduos em uma determinada linha mostram uma tendência crescente ou decrescente em níveis crescentes de um \(Y\) ordinal, isso indica que as respostas dessa linha tendem a níveis maiores ou menores de \(Y\) do que em outras linhas onde isso não ocorre. ocorrer.
Alternativamente, se uma célula tiver um resíduo grande e de sinal oposto de outras na mesma linha e coluna, essa célula parece violar a independência. Infelizmente, os padrões de resíduos nem sempre são aparentes. A modelagem de regressão é uma abordagem mais sistemática para avaliar a associação entre as variáveis. As seções 3.3 e 3.4 cobrem uma variedade de estruturas de modelos que nos permitem entender melhor a associação entre \(X\) e \(Y\). Além disso, as propriedades ordinais das categorias de resposta multinomial podem ser levadas em conta mais facilmente com esses modelos e apontam para fontes de associação que podem não ser notadas de outra forma. Razões de chances e resíduos ainda podem ser examinados dentro desses modelos de regressão.
3.3 Modelos de regressão de resposta nominal
Passamos agora ao problema de modelar as probabilidades de uma variável de resposta categórica \(Y\) com categorias de resposta \(j = 1,\cdots,J\) usando variáveis explicativas \(x_1,\cdots,x_p\). Na Seção 1.2.5, definimos probabilidades para uma resposta binária como \(P(\mbox{sucesso})/P(\mbox{fracasso})\). De maneira mais geral, para uma resposta multinomial, podemos definir as probabilidades como uma comparação de qualquer par de categorias de resposta.
Por exemplo, \(\pi_j/\pi_{j'}\) é a probabilidade da categoria \(j\) relativa a \(j'\). Um modelo de regressão popular para respostas multinomiais é desenvolvido então selecionando uma categoria de resposta como o nível base e formando as chances das \(J-1\) categorias restantes contra este nível. Por exemplo, considere que \(j = 1\) represente a categoria de nível básico e forme as chances (odds) \(\pi_j/\pi_1\) para \(j = 2,\cdots,J\). Essas probabilidades são então modeladas em função de variáveis explicativas usando uma forma generalizada de regressão logística.
Observe que se \(J = 2\), temos \(\log(\pi_2/\pi_1) = \log\big(\pi_2/(1-\pi_2)\big)\), que é equivalente a \(\log\big(\pi/(1-\pi)\big)\) como na regressão logística, a categoria de resposta 2 torna-se a categoria de “sucesso”.
Especificamente, um modelo de regressão multinomial, também conhecido como modelo logit de categoria de linha de base relaciona um conjunto de variáveis explicativas a cada probabilidade logarítmica de acordo com \[ \log(\pi_j/\pi_1) = \beta_{j0} + \beta_{j1} x_1 +\cdots + \beta_{jp} x_p, \]
para \(j = 2,\cdots,J\).
Observe que o primeiro subscrito nos parâmetros \(\beta\) corresponde à categoria de resposta, o que permite que os logaritmos das probabilidades (log-odds) de cada resposta se relacionem com as variáveis explicativas de uma maneira diferente. Além disso, observe que, subtraindo os log-odds apropriadas, podemos reescrever a equação acima para comparar qualquer par de categorias de resposta.
Por exemplo, para achar \(\log(\pi_j/\pi_{j'})\) onde \(j'\neq 1\) e \(j'\neq j\), temos \[ \begin{array}{rcl} \log(\pi_j/\pi_{j'}) & = & \log(\pi_j)-\log(\pi_{j'}) \, = \, \log(\pi_j) - \log(\pi_{j'}) +\log(\pi_1)-\log(\pi_1) \\ & = & \log(\pi_j/\pi_1)-\log(\pi_{j'}/\pi_1) \\ & = & (\beta_{j0}-\beta_{j' \,0})+(\beta_{j1}-\beta_{j' \, 1})x_1 + \cdots (\beta_{jp}-\beta_{j'\, p})x_p\cdot \end{array} \] Assim, a escolha do nível de base não é importante e pode ser feita com base na conveniência ou interpretação.
As probabilidades para cada categoria individual também podem ser encontradas em termos do modelo. Podemos reescrever \[ \pi_j=\pi_1 \exp\big( \beta_{j0}+ \beta_{j1}x_1+\cdots+\beta_{jp}x_p \big) \] usando propriedades de logaritmos. Observando que \(\pi_1 +\pi_2 +\cdots +\pi_J = 1\), temos \[ \pi_1+\pi_1 \exp\big( \beta_{20}+ \beta_{21}x_1+\cdots+\beta_{2p}x_p \big)+\cdots+\pi_1 \exp\big( \beta_{J0}+ \beta_{J1}x_1+\cdots+\beta_{Jp}x_p \big) \, = \, 1\cdot \]
Fatorando o \(\pi_1\) comum em cada termo, obtemos uma expressão para \(\pi_1\): \[ \pi_1=\dfrac{1}{\displaystyle 1+\sum_{j=2}^J \exp\big( \beta_{j0}+ \beta_{j1}x_1+\cdots+\beta_{jp}x_p \big)}\cdot \] Isso leva a uma expressão geral para \(\pi_j\): \[ \pi_j=\dfrac{\exp\big( \beta_{j0}+ \beta_{j1}x_1+\cdots+\beta_{jp}x_p \big)}{\displaystyle 1+\sum_{\ell=2}^J \exp\big( \beta_{\ell 0}+ \beta_{\ell 1}x_1+\cdots+\beta_{\ell p}x_p \big)}, \]
para \(j=2,\cdots,J\).
Os parâmetros para o modelo são estimados usando a máxima verossimilhança. Para uma amostra de observações, \(y_i\) que denotam a categoria resposta e as correspondentes variáveis explicativas \(x_{i1},\cdots, x_{ip}\) para \(i = 1,\cdots,m\), a função de verossimilhança é simplesmente o produto de \(m\) distribuições multinomiais com parâmetros de probabilidade dados anteriormente.
Procedimentos numéricos iterativos são então usados para encontrar as estimativas dos parâmetros. A função multinom() do pacote nnet executa os cálculos necessários. Conforme mencionado na ajuda desta função, redimensionar as variáveis explicativas para que fiquem entre 0 e 1 pode ajudar às vezes com a convergência das estimativas dos parâmetros.
A matriz de covariância estimada para as estimativas dos parâmetros é encontrada usando procedimentos de verossimilhança padrão. Os métodos de inferência baseados nas estatísticas de Wald e LR são executados da mesma forma que os procedimentos de verossimilhança nos capítulos anteriores. No entanto, os intervalos de razão de verossimilhança perfilada são difíceis de obter em R, porque nem uma função do método confint() nem o pacote mcprofile foram estendidos para calcular esses tipos de intervalos para objetos resultantes de multinom(). Por esta razão, focamos no uso de intervalos de Wald nesta seção.
Exemplo 3.4: Grãos de trigo (Kernels).
A presença de grãos germinados ou doentes no trigo pode reduzir o valor de toda a safra de um produtor de trigo. É importante identificar esses grãos depois de colhidos, mas antes da venda. Para facilitar esse processo de identificação, sistemas automatizados foram desenvolvidos para separar kernels saudáveis dos demais. Melhorar esses sistemas requer uma melhor compreensão das maneiras mensuráveis pelas quais os grãos saudáveis diferem dos grãos que brotaram prematuramente ou estão infectados com um fungo (“Scab”).
Para isso, Martin et al. (1998) realizaram um estudo examinando inúmeras propriedades físicas de grãos - densidade, dureza, tamanho, peso e teor de umidade - medidos em uma amostra de grãos de trigo de duas classes diferentes de trigo, inverno vermelho duro (hrw) e inverno vermelho suave (srw). A condição de cada kernel também foi classificada como “Healthy” (Saudável), “Sprout” ou “Scab” por inspeção visual humana.
Nos dados fornecidos pelos autores deste artigo, temos medições de 275 grãos de trigo. Segue abaixo uma parte dos dados:
wheat <- read.csv(file = "http://leg.ufpr.br/~lucambio/ADC/wheat.csv", stringsAsFactors = TRUE)
# NOTE 2021-10-19: stringsAsFactors = TRUE was added due to changes in R version 4.0.0
head(wheat, n = 3) # n argument gives the number of rows to print## class density hardness size weight moisture type
## 1 hrw 1.349253 60.32952 2.30274 24.6480 12.01538 Healthy
## 2 hrw 1.287440 56.08972 2.72573 33.2985 12.17396 Healthy
## 3 hrw 1.233985 43.98743 2.51246 31.7580 11.87949 Healthytail(wheat, n = 3)## class density hardness size weight moisture type
## 273 srw 0.8491887 34.06615 1.40665 12.0870 11.92744 Scab
## 274 srw 1.1770230 60.97838 1.05690 9.4800 12.24046 Scab
## 275 srw 1.0305543 -9.57063 2.05691 23.8185 12.64962 ScabComeçamos com gráficos úteis para examinar inicialmente os dados. A figura abaixo mostra um gráfico de coordenadas paralelas construído usando a função parcoord() do pacote MASS. O eixo horizontal dá o nome de cada variável, e acima desses nomes estão os valores correspondentes para cada kernel. Como é difícil identificar números de kernel específicos, fornecemos código no programa correspondente que mostra como destacar um kernel. Por exemplo, pode-se ver que o kernel #269 é o kernel com o menor peso.
Esses valores são dimensionados para que os valores mínimo e máximo de cada variável apareçam na parte inferior e superior, respectivamente, do eixo vertical.
Uma linha é desenhada para cada kernel indicando sua posição nas variáveis. Por exemplo, o kernel que tem o menor peso (linha mais baixa para peso) também tem o segundo menor tamanho (segunda linha mais baixa para tamanho). Usamos diferentes tipos de linha para representar os três tipos de kernel, para que possamos entender algumas das características que diferenciam as condições “Healthy”, “Sprout” ou “Scab”.
# Parallel coordinate plot
library(package = MASS) # Location of parcoord() function
# Reorder variables because class is binary (may distort plot)
# Create indicator variable for class
wheat2<-data.frame(kernel = 1:nrow(wheat), wheat[,2:6],
class.new = ifelse(test = wheat$class == "hrw", yes = 0, no = 1))
head(wheat2)## kernel density hardness size weight moisture class.new
## 1 1 1.349253 60.32952 2.30274 24.6480 12.01538 0
## 2 2 1.287440 56.08972 2.72573 33.2985 12.17396 0
## 3 3 1.233985 43.98743 2.51246 31.7580 11.87949 0
## 4 4 1.336534 53.81704 2.27164 32.7060 12.11407 0
## 5 5 1.259040 44.39327 2.35478 26.0700 12.06487 0
## 6 6 1.300258 48.12066 2.49132 33.2985 12.18577 0# Colors by condition:
wheat.colors<-ifelse(test = wheat$type=="Healthy", yes = "black",
no = ifelse(test = wheat$type=="Sprout", yes = "red", no = "green"))
# Line type by condition:
wheat.lty<-ifelse(test = wheat$type=="Healthy", yes = "solid",
no = ifelse(test = wheat$type=="Sprout", yes = "longdash", no = "dotdash"))
parcoord(x = wheat2, col = wheat.colors, lty = wheat.lty) # Plot
legend(x=6.15, y=0.75, legend = c("Healthy", "Sprout", "Scab"), lty = c("solid", "longdash", "dotdash"),
col=c("black", "red", "green"), cex=0.8, bty="n")
# legend(locator(1), legend = c("Healthy", "Sprout", "Scab"), lty = c("solid", "longdash", "dotdash"),
# col=c("black", "red", "green"), cex=0.8, bty="n") #Original legend using locator(1)Gráfico de coordenadas paralelas dos dados de trigo. A variável do kernel corresponde ao número de observação no conjunto de dados. A variável class.new é 1 para trigo de inverno vermelho macio e 0 para trigo de inverno vermelho duro.
# What observation has the largest weight
wheat[wheat$weight == min(wheat2$weight),] # 269## class density hardness size weight moisture type
## 269 srw 0.9343233 48.66988 0.88496 8.532 11.81367 Scaborder(wheat$size) # 268 has the 2nd largest size## [1] 268 269 140 139 274 27 141 171 68 209 33 29 266 273 206 43 253 179
## [19] 264 265 78 63 174 124 36 234 177 65 70 17 208 190 250 97 239 73
## [37] 143 207 19 142 152 64 175 231 267 101 132 271 127 210 26 251 162 122
## [55] 98 105 31 145 184 218 66 154 157 126 48 61 112 205 56 67 168 197
## [73] 204 96 11 211 72 215 227 170 71 110 203 100 263 15 182 128 262 237
## [91] 49 246 91 176 272 80 138 38 53 172 173 192 275 12 169 34 135 107
## [109] 270 45 111 134 83 90 202 241 167 248 22 149 188 212 255 214 10 144
## [127] 130 252 76 115 224 117 30 46 21 191 258 259 261 13 183 81 121 77
## [145] 4 35 181 245 158 233 247 18 24 164 1 108 44 151 155 50 199 106
## [163] 32 39 5 223 55 228 256 161 41 118 136 129 89 109 133 153 213 243
## [181] 88 92 189 260 103 14 194 232 60 69 74 180 186 219 82 85 28 235
## [199] 238 6 220 198 226 54 3 165 229 57 102 163 193 201 62 113 166 93
## [217] 42 221 249 225 187 79 146 195 120 131 147 148 59 75 99 94 217 95
## [235] 137 40 236 240 196 257 150 200 52 8 51 9 2 16 104 123 47 84
## [253] 86 25 114 37 116 185 20 242 23 230 160 7 222 58 119 254 159 244
## [271] 216 178 87 125 156Os grãos de Scab geralmente têm valores menores de density (densidade), size (tamanho) e weight (peso),
Grãos Healthy (saudáveis) podem ter densidades mais altas,
Há muita sobreposição para grãos Healthy e Sprout (germinados), e
O moisture (teor de umidade) parece estar relacionado à classe de trigo de inverno vermelho duro ou macio.
Para examinar esses dados de forma mais formal, estimamos um modelo de regressão multinomial usando a classe do trigo e as cinco medidas em cada grão como variáveis explicativas:
levels ( wheat$type ) # The 3 response categories## [1] "Healthy" "Scab" "Sprout"library ( package = nnet )
mod.fit <- multinom ( formula = type ~ class + density + hardness + size + weight + moisture , data = wheat )## # weights: 24 (14 variable)
## initial value 302.118379
## iter 10 value 234.991271
## iter 20 value 192.127549
## final value 192.112352
## convergedsummary (mod.fit )## Call:
## multinom(formula = type ~ class + density + hardness + size +
## weight + moisture, data = wheat)
##
## Coefficients:
## (Intercept) classsrw density hardness size weight
## Scab 30.54650 -0.6481277 -21.59715 -0.01590741 1.0691139 -0.2896482
## Sprout 19.16857 -0.2247384 -15.11667 -0.02102047 0.8756135 -0.0473169
## moisture
## Scab 0.10956505
## Sprout -0.04299695
##
## Std. Errors:
## (Intercept) classsrw density hardness size weight
## Scab 4.289865 0.6630948 3.116174 0.010274587 0.7722862 0.06170252
## Sprout 3.767214 0.5009199 2.764306 0.008105748 0.5409317 0.03697493
## moisture
## Scab 0.1548407
## Sprout 0.1127188
##
## Residual Deviance: 384.2247
## AIC: 412.2247A saída de levels(wheat$type) nos mostra que “Healthy” é armazenado como o primeiro nível de type, então multinom() o usa como o nível base. Além disso, observe que uma variável indicadora é criada para class, onde 1 é para “srw” e 0 para “hrw”. Essas escolhas de R seguem a convenção usual de tornar o primeiro nível de um fator o nível base. Da tabela de coeficientes na saída, obtemos o modelo estimado como \[ \log\big(\widehat{\pi}_\mbox{Scab}/\widehat{\pi}_\mbox{Healthy}\big) = 30.55 -0.65 \, \mbox{Scab} -21.60 \, \mbox{density} -0.01 \, \mbox{hardness} +1.07 \, \mbox{size} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad -0.29 \, \mbox{weight} + 0.11 \, \mbox{moisture} \]
e \[ \log\big(\widehat{\pi}_\mbox{Sprout}/\widehat{\pi}_\mbox{Healthy}\big) = 19.17 -0.22 \, \mbox{srw} -15.12 \, \mbox{density} -0.021 \, \mbox{hardness} +0.88 \, \mbox{size} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad -0.047 \, \mbox{weight} - 0.043 \, \mbox{moisture}\cdot \]
Testes de hipóteses da forma \(H_0 \, : \, \beta_{jr} = 0\) vs. \(H_1 \, : \, \beta_{jr}\neq 0\) podem ser realizados por testes de Wald ou LRTs. Por exemplo, um teste de Wald de \(H_0 \, : \, \beta_{21} = 0\) vs. \(H_1 \, : \, \beta_{21}\neq 0\), que corresponde à variável explicativa classe no log-odds para scab vs. healty, resulta em \[ Z_0 = \dfrac{\widehat{\beta}_{21}}{\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_{21})}}=\dfrac{-0.6481}{0.6631}=-0.98 \]
e um \(p\)-valor de \(2P(Z > |-0.98|) = 0.33\). Assim, não há evidências suficientes de que o trigo de inverno vermelho duro e mole tenham efeitos diferentes no scab ou no estado healthy dos grãos, uma vez que as outras variáveis explicativas estão no modelo.
As hipóteses mais interessantes exploram os efeitos de uma determinada variável explicativa em todas as categorias de resposta. Por exemplo, a classe de trigo tem algum efeito sobre as probabilidades de grãos saudáveis, sarna (scab) ou brotados? Para investigar isso, precisamos testar se todos os parâmetros correspondentes à variável explicativa de interesse, digamos \(x_r\), são zero: \(H_0 \, : \, \beta_{jr} = 0\), \(j = 2,\cdots,J\).
No caso da variável explicativa class, usamos \(H_0 \, : \, \beta_{21} = \beta_{31} = 0\) vs. \(H_1 \, : \, \beta_{21}\neq 0\) e/ou \(\beta_{31}\neq 0\). Este teste é facilmente realizado com métodos LR usando a função Anova():
library ( package = car)
Anova (mod.fit )## Analysis of Deviance Table (Type II tests)
##
## Response: type
## LR Chisq Df Pr(>Chisq)
## class 0.964 2 0.6175
## density 90.555 2 < 2.2e-16 ***
## hardness 7.074 2 0.0291 *
## size 3.211 2 0.2008
## weight 28.230 2 7.411e-07 ***
## moisture 1.193 2 0.5506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qchisq(p = 0.95, df = 2)## [1] 5.991465# Wald tests for one beta - usually not done since it tests only one parameter
sum.fit<-summary(mod.fit)
test.stat<-sum.fit$coefficients/sum.fit$standard.errors
p.value<-2*(1-pnorm(q = abs(test.stat), mean = 0, sd = 1))
round(test.stat,2)## (Intercept) classsrw density hardness size weight moisture
## Scab 7.12 -0.98 -6.93 -1.55 1.38 -4.69 0.71
## Sprout 5.09 -0.45 -5.47 -2.59 1.62 -1.28 -0.38round(p.value,2)## (Intercept) classsrw density hardness size weight moisture
## Scab 0 0.33 0 0.12 0.17 0.0 0.48
## Sprout 0 0.65 0 0.01 0.11 0.2 0.70# LRT for class:
mod.fit.H0<-multinom(formula = type ~ density + hardness + size + weight + moisture, data=wheat)## # weights: 21 (12 variable)
## initial value 302.118379
## iter 10 value 237.492991
## iter 20 value 192.609457
## final value 192.594412
## convergedanova(mod.fit.H0, mod.fit) ## Likelihood ratio tests of Multinomial Models
##
## Response: type
## Model Resid. df Resid. Dev
## 1 density + hardness + size + weight + moisture 538 385.1888
## 2 class + density + hardness + size + weight + moisture 536 384.2247
## Test Df LR stat. Pr(Chi)
## 1
## 2 1 vs 2 2 0.9641189 0.6175104# Additional code showing the LRT
G.sq.H0<-mod.fit.H0$deviance
G.sq.H1<-mod.fit$deviance
G.sq<-G.sq.H0-G.sq.H1
p.value<-1-pchisq(q = G.sq, df = 2)
data.frame(G.sq.H0, G.sq.H1, G.sq, p.value, df = 2)## G.sq.H0 G.sq.H1 G.sq p.value df
## 1 385.1888 384.2247 0.9641189 0.6175104 2A estatística de teste transformada é \(-2\log(\Lambda ) = 0.964\) e o \(p\)-valor correspondente é 0.6175 usando uma aproximação da distribuição \(\chi_2^2\). Por causa do grande \(p\)-valor, não há evidências suficientes para indicar que class de trigo é importante, uma vez que as demais variáveis estão no modelo. Testes separados para density, hardness e weight na saída indicam pelo menos evidências marginais de importância para essas variáveis explicativas. Observe que a função anova() também pode ser usada para realizar esses testes.
Podemos usar names(mod.fit) para obter uma lista dos objetos dentro da lista mod.fit. Estes incluem deviance (desvio residual) e convergence (indicador de convergência). O componente contém fitted.values, as probabilidades estimadas de healthy (saudável), scab (sarna) ou sprout (broto) para cada kernel. Alternativamente, a função predict() pode calcular essas probabilidades estimadas usando o valor do argumento type = “probs”:
pi.hat <- predict ( object = mod.fit , newdata = wheat , type = "probs")
head (pi.hat)## Healthy Scab Sprout
## 1 0.8552110 0.046396827 0.09839221
## 2 0.7492553 0.021572158 0.22917255
## 3 0.5172800 0.068979903 0.41374011
## 4 0.8982064 0.006740716 0.09505287
## 5 0.5103245 0.176260796 0.31341473
## 6 0.7924907 0.015304122 0.19220522Por exemplo,a probabilidade estimada de ser healthy (saudável) para o primeiro kernel é
\[ \begin{array}{l} \qquad \widehat{\pi}_\mbox{healthy} \, = \, \\ \displaystyle \dfrac{1}{1+\exp(30.55-0.65\times 0+\cdots+0.11\times 12.02)+\exp(19.77-0.22\times 0+\cdots-0.043\times 12.02)} \\ \, = \, 0.8552\cdot \end{array} \]
A função predict() também pode ser usada para encontrar a categoria de resposta mais provável para cada kernel usando o valor do argumento type = “class”, que é o padrão.
Calcular intervalos de confiança para \(\widehat{\pi}_\mbox{Healthy}\), \(\widehat{\pi}_\mbox{Scab}\) e \(\widehat{\pi}_\mbox{Sprout}\) não é tão fácil quanto se poderia esperar. A função predict() não fornece as variâncias estimadas necessárias para calcular um intervalo de confiança Wald para \(\pi_j\). A razão é que essas probabilidades formam uma distribuição de probabilidades. Uma região de confiança conjunta seria necessária nas dimensões \(J-1\), o “-1” é devido a \(\sum_{j=1}^J \pi_j = 1\) para obter adequadamente um “intervalo” para essas probabilidades.
Porque uma região de confiança conjunta pode ser difícil para calcular e interpretar, uma alternativa são intervalos de confiança um de cada vez para cada \(\pi_j\) e use-os como uma aproximação. Isso permite ainda levar em conta a variabilidade do estimador de alguma maneira.
Esses intervalos de confiança de Wald um de cada vez podem ser calculados como \[ \widehat{\pi}_j\pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}(\widehat{\pi}_j)} \] com a ajuda da função deltaMethod(). Abaixo estão os cálculos para a primeira observação do conjunto de dados e \({\pi}_\mbox{Healthy}\):
# Obtain observation values , ";" ends a command
x1 <- 0; x2 <- wheat [1 ,2]; x3 <- wheat [1 ,3]
x4 <- wheat [1 ,4]; x5 <- wheat [1 ,5]; x6 <- wheat [1 ,6]
g.healthy <- "1/(1 + exp(b20 + b21*x1 + b22*x2 + b23*x3 + b24 *x4 + b25*x5 + b26*x6) +
exp (b30 + b31*x1 + b32*x2 + b33 *x3 + b34 *x4 + b35*x5 + b36*x6))"
calc.healthy <- deltaMethod ( object = mod.fit , g = g.healthy,
parameterNames = c("b20", "b21", "b22", "b23", "b24", "b25",
"b26", "b30", "b31", "b32", "b33", "b34", "b35", "b36"))
calc.healthy$Estimate # pi - hat_Healthy## [1] 0.855211calc.healthy$SE # sqrt (Var -hat(pi - hat_Healthy ))## [1] 0.05999678alpha <- 0.05
calc.healthy$Estimate + qnorm (p = c( alpha /2, 1- alpha /2)) * calc.healthy$SE## [1] 0.7376194 0.9728025O objeto g.healthy contém uma cadeia de caracteres para a fórmula usada para calcular \({\pi}_\mbox{Healthy}\), onde usamos bjr para denotar \(\beta_{jr}\) e xr para denotar \(x_r\) para \(j = 2, 3\); \(r = 0,\cdots, 6\). Pode haver maneiras mais eficientes de inserir cadeias de caracteres como essas em R. A função deltaMethod() calcula \(\widehat{\pi}_\mbox{Healthy}\) e usa o método delta para calcular \(\sqrt{\widehat{\mbox{Var}}(\widehat{\pi}_j)}\). O intervalo de confiança de 95% para \({\pi}_\mbox{Healthy}\) é (0.7376, 0.9728). Os intervalos de confiança de 95% para \({\pi}_\mbox{Scab}\) e \({\pi}_\mbox{Sprout}\) são (-0.0067, 0.0995) e (0.0143, 0.1825), respectivamente. Observe que esses tipos de intervalos podem estar abaixo de 0 ou acima de 1.
Além disso, observe que combinações de valores para \({\pi}_\mbox{Healthy}\), \({\pi}_\mbox{Scab}\) e \({\pi}_\mbox{Sprout}\) que todos estejam dentro de seus respectivos intervalos poderia resultar em \(\sum_{j=1}^3 \pi_j\neq 1\), pois esses intervalos não formam uma região de confiança conjunta.
O modelo estimado pode ser apresentado incluindo a probabilidade estimada para cada tipo de trigo no eixo \(y\) e uma das variáveis explicativas no eixo \(x\), mantendo as outras variáveis dentro do modelo constantes. Como o gráfico é condicional a esses valores de variáveis mantidos constantes, mais de um gráfico geralmente deve ser feito em diferentes valores constantes. Em vez disso, para este exemplo, reestimamos um modelo com density como a única variável explicativa: \[ \log\big(\widehat{\pi}_{\mbox{Scab}}/\widehat{\pi}_{\mbox{Healthy}}\big) = 29.38 - 24.56 \, \mbox{density} \] e \[ \log\big(\widehat{\pi}_{\mbox{Sprout}}/\widehat{\pi}_{\mbox{Healthy}}\big) = 19.12 - 15.48 \, \mbox{density}\cdot \]
# Estimate model with density only
mod.fit.nom.density<-multinom(formula = type ~ density, data = wheat)## # weights: 9 (4 variable)
## initial value 302.118379
## iter 10 value 229.769334
## iter 20 value 229.712304
## final value 229.712290
## convergedsummary(mod.fit.nom.density)## Call:
## multinom(formula = type ~ density, data = wheat)
##
## Coefficients:
## (Intercept) density
## Scab 29.37827 -24.56215
## Sprout 19.12165 -15.47633
##
## Std. Errors:
## (Intercept) density
## Scab 3.676892 3.017842
## Sprout 3.337092 2.691429
##
## Residual Deviance: 459.4246
## AIC: 467.4246beta.hat<-coefficients(mod.fit.nom.density)
curve(expr = 1/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) + exp(beta.hat[2,1] + beta.hat[2,2]*x)),
ylab = expression(hat(pi)), xlab = "Density", xlim = c(min(wheat$density), max(wheat$density)),
col = "black", lty = "solid", lwd = 2, n = 1000, type = "n",
panel.first = grid(col = "gray", lty = "dotted"))
# Plot each pi_j
curve(expr = 1/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) + exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "black", lty = "solid", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Healthy"]),
max(wheat$density[wheat$type == "Healthy"]))) # Healthy
curve(expr = exp(beta.hat[1,1] + beta.hat[1,2]*x)/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) +
exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "green", lty = "dotdash", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Scab"]), max(wheat$density[wheat$type == "Scab"]))) # Scab
curve(expr = exp(beta.hat[2,1] + beta.hat[2,2]*x)/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) +
exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "red", lty = "longdash", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Sprout"]),
max(wheat$density[wheat$type == "Sprout"]))) # Sprout
legend(x = 1.4, y = 0.8, legend=c("Healthy", "Sprout", "Scab"), lty=c("solid","longdash","dotdash"),
col=c("black","red","green"), bty="n", lwd = c(2,2,2), seg.len = 4)
# Verify plot is correct
density.values<-seq(from = 0.8, to = 1.6, by = 0.1)
data.frame(density.values, round(predict(object = mod.fit.nom.density,
newdata = data.frame(density = density.values), type = "probs"), 2))## density.values Healthy Scab Sprout
## 1 0.8 0.00 0.95 0.05
## 2 0.9 0.00 0.89 0.11
## 3 1.0 0.01 0.76 0.24
## 4 1.1 0.05 0.54 0.41
## 5 1.2 0.27 0.25 0.48
## 6 1.3 0.69 0.05 0.25
## 7 1.4 0.92 0.01 0.07
## 8 1.5 0.98 0.00 0.02
## 9 1.6 1.00 0.00 0.00O modelo é mostrado na figura, onde as probabilidades estimadas para cada condição de kernel são desenhadas entre o menor e o maior valor de densidade observado para aquela condição. A partir do gráfico, vemos que a probabilidade de sarna (scab) estimada é a maior para os grãos de menor densidade. A probabilidade saudável estimada é a maior para os grãos de alta densidade. Para níveis de densidade no meio, o broto (sprout) tem a maior probabilidade estimada. O gráfico de coordenadas paralelas apresenta achados semelhantes onde os níveis de densidade tendem a seguir a ordenação sarna < broto < saudável.
Transformações e termos de interação podem ser incluídos em multinom() de maneira semelhante a glm(). Fornecemos um exemplo no programa correspondente que mostra como estimar parâmetros em \[ log\big(\pi_j/\pi_{\mbox{Healthy}}\big) = \beta_{j0} + \beta_{j1} \, \mbox{SRW} + \beta_{j2}\, \mbox{density} + \beta_{j3} \, \mbox{density}^2 + \beta_{j4} \, \mbox{SRW}\times \mbox{density}\cdot \]
Após 100 iterações, a convergência não é alcançada. No entanto, quando o número de iterações é aumentado para 500, usando o argumento maxit, a convergência é alcançada.
# Additional investigations
# Quadtratic and interaction terms example
mod.fit.trans1<-multinom(formula = type ~ class + density + I(density^2) +
density:class, data=wheat) # Does not converge (notice says "stopped after 100")## # weights: 18 (10 variable)
## initial value 302.118379
## iter 10 value 231.793812
## iter 20 value 225.610537
## iter 30 value 218.263988
## iter 40 value 217.052033
## iter 50 value 216.875672
## iter 60 value 216.282105
## iter 70 value 216.157837
## iter 80 value 215.874353
## iter 90 value 215.834023
## iter 100 value 215.715350
## final value 215.715350
## stopped after 100 iterationssummary(mod.fit.trans1)## Call:
## multinom(formula = type ~ class + density + I(density^2) + density:class,
## data = wheat)
##
## Coefficients:
## (Intercept) classsrw density I(density^2) classsrw:density
## Scab 92.31329 1.466333 -133.39721 46.78379 -1.056592
## Sprout -11.06620 11.273002 33.83208 -20.06648 -9.289201
##
## Std. Errors:
## (Intercept) classsrw density I(density^2) classsrw:density
## Scab 20.76329 5.805945 31.33423 11.80758 4.682795
## Sprout 37.54842 6.719231 60.98022 24.80113 5.455298
##
## Residual Deviance: 431.4307
## AIC: 451.4307mod.fit.trans2<-multinom(formula = type ~ class + density + I(density^2) + density:class,
data=wheat, maxit = 1000) # Converges## # weights: 18 (10 variable)
## initial value 302.118379
## iter 10 value 231.793812
## iter 20 value 225.610537
## iter 30 value 218.263988
## iter 40 value 217.052033
## iter 50 value 216.875672
## iter 60 value 216.282105
## iter 70 value 216.157837
## iter 80 value 215.874353
## iter 90 value 215.834023
## iter 100 value 215.715350
## iter 110 value 215.545845
## iter 120 value 215.355028
## iter 130 value 215.269680
## iter 140 value 215.182222
## iter 150 value 214.996559
## iter 160 value 214.823422
## iter 170 value 214.771269
## iter 180 value 214.723800
## iter 190 value 214.544479
## iter 200 value 214.394734
## iter 210 value 214.337652
## iter 220 value 214.201358
## iter 230 value 214.171255
## final value 214.168438
## convergedsummary(mod.fit.trans2)## Call:
## multinom(formula = type ~ class + density + I(density^2) + density:class,
## data = wheat, maxit = 1000)
##
## Coefficients:
## (Intercept) classsrw density I(density^2) classsrw:density
## Scab 112.11720 0.3947469 -162.15852 57.04640 -0.1771958
## Sprout 31.72094 9.9762283 -34.49787 7.15602 -8.2035424
##
## Std. Errors:
## (Intercept) classsrw density I(density^2) classsrw:density
## Scab 22.56297 5.875188 33.64421 12.50196 4.729433
## Sprout 36.64202 6.748108 58.77756 23.62749 5.461176
##
## Residual Deviance: 428.3369
## AIC: 448.3369# Make the last level the base level
wheat.original<-wheat # Save original format
wheat$type<-relevel(x = wheat$type, ref = "Sprout")
levels(wheat$type)## [1] "Sprout" "Healthy" "Scab"mod.fit.relevel<-multinom(formula = type ~ class + density + hardness + size + weight + moisture, data=wheat)## # weights: 24 (14 variable)
## initial value 302.118379
## iter 10 value 214.152215
## iter 20 value 192.195956
## final value 192.112352
## convergedsummary(mod.fit.relevel)## Call:
## multinom(formula = type ~ class + density + hardness + size +
## weight + moisture, data = wheat)
##
## Coefficients:
## (Intercept) classsrw density hardness size weight
## Healthy -19.16859 0.2247587 15.116925 0.021020145 -0.8756472 0.04731197
## Scab 11.37795 -0.4234045 -6.480573 0.005113742 0.1935585 -0.24233291
## moisture
## Healthy 0.04299084
## Scab 0.15255994
##
## Std. Errors:
## (Intercept) classsrw density hardness size weight
## Healthy 3.767230 0.5009202 2.764335 0.008105752 0.5409349 0.03697498
## Scab 2.641839 0.6040948 1.852886 0.008597373 0.6774422 0.05581095
## moisture
## Healthy 0.1127190
## Scab 0.1409521
##
## Residual Deviance: 384.2247
## AIC: 412.2247wheat<-wheat.original # Change back to original formatUsando a codificação padrão de “Healthy” como o nível básico da resposta foi ideal neste exemplo para que as comparações pudessem ser feitas mais facilmente com cada uma das duas condições não saudáveis. No entanto, pode haver outras situações em que faça sentido alterar a ordem padrão. A função relevel() pode ser usada como na Seção 2.2.6 para este propósito. Por exemplo, podemos tornar “Sprout” o nível base com wheat$type <- relevel(x = wheat$type, ref = “Sprout”) e, em seguida, reestimar o modelo usando multinom().
A estimativa de modelos de regressão multinomial também pode ser feita usando a função geral de ajuste de modelo vglm() do pacote VGAM (Yee, 2010). Fornecemos código mostrando como usar esta função no programa correspondente para este exemplo. Uma vantagem de usar vglm() é que às vezes ele pode obter uma convergência mais rápida de estimativas de parâmetros. Além disso, a função pode ser usada para ajustar uma ampla variedade de tipos de modelos para dados categóricos. Focamos na função multinom() aqui porque seu pacote, nnet, faz parte da instalação padrão do R.
# Using vglm() from the VGAM (Vector Generalized Additive Models) package
library(package = VGAM)
# This will match mod.fit with multinom()
mod.fit2<-vglm(formula = type ~ class + density + hardness + size + weight + moisture,
family = multinomial(refLevel = 1), data = wheat)
summary(mod.fit2)##
## Call:
## vglm(formula = type ~ class + density + hardness + size + weight +
## moisture, family = multinomial(refLevel = 1), data = wheat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 30.545829 4.289769 7.121 1.07e-12 ***
## (Intercept):2 19.167800 3.767134 5.088 3.62e-07 ***
## classsrw:1 -0.648242 0.663090 -0.978 0.32827
## classsrw:2 -0.224810 0.500918 -0.449 0.65358
## density:1 -21.596968 3.116126 -6.931 4.19e-12 ***
## density:2 -15.116309 2.764258 -5.468 4.54e-08 ***
## hardness:1 -0.015907 0.010275 -1.548 0.12157
## hardness:2 -0.021021 0.008106 -2.593 0.00951 **
## size:1 1.069246 0.772281 1.385 0.16620
## size:2 0.875660 0.540932 1.619 0.10549
## weight:1 -0.289652 0.061702 -4.694 2.67e-06 ***
## weight:2 -0.047315 0.036975 -1.280 0.20066
## moisture:1 0.109590 0.154840 0.708 0.47909
## moisture:2 -0.042978 0.112718 -0.381 0.70299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1])
##
## Residual deviance: 384.2247 on 536 degrees of freedom
##
## Log-likelihood: -192.1124 on 536 degrees of freedom
##
## Number of Fisher scoring iterations: 6
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1', 'weight:1'
##
##
## Reference group is level 1 of the responsenames(attributes(mod.fit2)) # Just list object names within it.## [1] "extra" "family" "iter" "predictors"
## [5] "assign" "callXm2" "contrasts" "df.residual"
## [9] "df.total" "dispersion" "effects" "offset"
## [13] "qr" "R" "rank" "ResSS"
## [17] "smart.prediction" "terms" "Xm2" "Ym2"
## [21] "xlevels" "call" "coefficients" "constraints"
## [25] "control" "criterion" "fitted.values" "misc"
## [29] "model" "na.action" "post" "preplot"
## [33] "prior.weights" "residuals" "weights" "x"
## [37] "y" "class"class(mod.fit2)## [1] "vglm"
## attr(,"package")
## [1] "VGAM"# Provides sqrt( Var^(log(pi^_j / pi^_1)) )
save.pred<-predictvglm(object = mod.fit2, newdata = wheat[1,], type = "link", se.fit = TRUE)
save.pred## $fitted.values
## log(mu[,2]/mu[,1]) log(mu[,3]/mu[,1])
## 1 -2.914012 -2.162308
##
## $se.fit
## log(mu[,2]/mu[,1]) log(mu[,3]/mu[,1])
## 1 0.6392031 0.5011193
##
## $df
## [1] 536
##
## $sigma
## [1] 1# predictvglm(object = mod.fit2, newdata = wheat[1,], type = "response", se.fit = TRUE) # Not allowed
predictvglm(object = mod.fit2, newdata = wheat[1,], type = "response", se.fit = FALSE)## Healthy Scab Sprout
## 1 0.8552002 0.0464011 0.09839873.3.1 Odds ratios
Como as probabilidades log-odds ou logaritmos da probabilidades de chances são modeladas diretamente em um modelo de regressão multinomial, as razões de chances são úteis para interpretar a relação de uma variável explicativa com a resposta. Conforme descrito na Seção 2.2.3 para modelos de regressão logística, as razões de chances para variáveis explicativas numéricas representam a mudança nas chances correspondentes a um aumento de \(c\) unidades em uma variável explicativa específica. A única diferença com os modelos multinomiais é que as probabilidades agora são formadas como uma comparação entre duas das \(J\) categorias de resposta.
Por exemplo, suponha que os termos na equação \[ \log(\pi_j/\pi_1) = \beta_{j0}+\beta_{j1} x_1+\cdots+\beta_{jp} x_p, \] sejam variáveis explicativas distintas, não transformações ou interações. As chances de uma resposta da categoria \(j\) versus uma resposta da categoria 1 são \(\exp(\beta_{j0} + \beta_{j1} x_1 +\cdots +\beta_{jp} x_p)\). Essas probabilidades mudam por \(\exp(c \beta_{jr})\) vezes para cada \(c\) unidade de aumento em \(x_r\), mantendo as outras variáveis constantes.
Da mesma forma, a partir da equação \[ \log(\pi_j/\pi_{j'})=(\beta_{j0}-\beta_{j'0})+(\beta_{j1}-\beta_{j'1}) x_1+\cdots+(\beta_{jp}-\beta_{j'p}) x_p, \]
as chances de uma categoria \(j\) vs. uma resposta de categoria \(j'\), \(j\neq j'\), \(j > 1\) e \(j' > 1\), mudam por \(\exp\big(c(\beta_{jr}-\beta_{j'r})\big)\) vezes para cada \(c\) unidades de aumento em \(x_r\) mantendo as outras variáveis no modelo constantes.
As estimativas de máxima verossimilhança das razões de chances são obtidas substituindo os parâmetros de regressão por suas estimativas correspondentes. Os métodos de inferência baseados em Wald e LR para odds ratios são usados da mesma forma que discutimos nos capítulos anteriores.
Exemplo 3.5: Grãos de trigo (Kernels).
As chances neste problema são construídas como \(P(Y = j)/P(Y = 1)\), \(j\) = 2, 3, onde categoria 1 = saudável (healthy); 2 = sarna (scab) e 3 = broto (sprout). As razões de chances estimadas para sarna vs. saudável ou broto vs. saudável são calculadas para cada variável explicativa como \(\widehat{OR} = \exp(c \widehat{\beta}_{jr})\) para \(j\) = 2, 3 e \(r = 1,\cdots, 6\). Para escolher os valores apropriados de \(c\), definimos \(c\) igual a um desvio padrão para cada variável explicativa contínua.
Cada um desses desvios padrão é encontrado usando a função apply(), que aplica a função sd() a colunas específicas (MARGIN = 2) do data frame do trigo. Para class, escolhemos \(c = 1\) porque é uma variável binária para trigo de inverno vermelho macio ou duro. Abaixo estão os desvios padrão estimados, seguidos pelas razões de chances:
sd.wheat <- apply (X = wheat [,-c(1 ,7 ,8)], MARGIN = 2, FUN = sd)
c.value <- c(1, sd.wheat ) # class = 1 is first value
round (c.value ,2)## density hardness size weight moisture
## 1.00 0.13 27.36 0.49 7.92 2.03# beta.hat_jr for r=1,...,6 and j = 2,3, where beta.hat_jr = coefficients (mod.fit)[j -1,r+1]
beta.hat2 <- coefficients (mod.fit) [1 ,2:7]
beta.hat3 <- coefficients (mod.fit) [2 ,2:7]
# Odds ratios for j = 2 vs. j = 1 ( scab vs. healthy )
round (exp(c.value * beta.hat2 ), 2)## density hardness size weight moisture
## 0.52 0.06 0.65 1.69 0.10 1.25round (1/ exp(c.value * beta.hat2 ), 2)## density hardness size weight moisture
## 1.91 17.04 1.55 0.59 9.90 0.80# Odds ratios for j = 3 vs. j = 1 ( sprout vs. healthy )
round (exp(c.value * beta.hat3 ), 2)## density hardness size weight moisture
## 0.80 0.14 0.56 1.54 0.69 0.92round (1/ exp(c.value * beta.hat3 ), 2)## density hardness size weight moisture
## 1.25 7.28 1.78 0.65 1.45 1.09Usamos a função coefficients() para obter os valores dos parâmetros de regressão estimados porque não há componente de coefficients em objetos produzidos por multinom(). As interpretações das razões de chances estimadas incluem:
1- As chances estimadas de scab vs. healthy mudam em 0.06 vezes para um aumento de 0.13 na densidade mantendo as outras variáveis constantes. De forma equivalente, podemos dizer que as chances estimadas de scab vs. healthy mudam em 17.04 vezes para uma diminuição de 0.13 na densidade mantendo as outras variáveis constantes. As chances estimadas de sprout vs. healthy mudam em 7.28 vezes para uma diminuição de 0.13 na densidade mantendo as outras variáveis constantes.
2- As chances estimadas de scab vs. healthy mudam em 9.90 vezes para uma diminuição de 7.92 no peso mantendo as outras variáveis constantes. As chances estimadas de sprout vs. healthy mudam em 1.45 vezes para uma diminuição de 7.92 no peso mantendo as outras variáveis constantes.Assim como na função glm(), existem várias funções de método associadas a objetos resultantes de multinom():
class (mod.fit )## [1] "multinom" "nnet"methods ( class = multinom )## [1] add1 anova Anova brief coef confint
## [7] Confint deltaMethod drop1 extractAIC logLik model.frame
## [13] predict print S summary vcov
## see '?methods' for accessing help and source codesqrt ( vcov (mod.fit) [2 ,2]) # sqrt (Var -hat(beta - hat_21 ))## [1] 0.6630948Observe que anteriormente usamos coefficients(), o equivalente a coef.multinom() para extrair os parâmetros de regressão estimados.
Os intervalos de confiança podem ser construídos para os parâmetros de regressão usando confint(). Infelizmente, a função do método confint.multinom() que é chamada pelo confint() genérico não produz intervalos LR perfilada. Os intervalos de Wald para as razões de chances são dados abaixo:
conf.beta <- confint ( object = mod.fit , level = 0.95)
conf.beta # Results are in a 3D array## , , Scab
##
## 2.5 % 97.5 %
## (Intercept) 22.13851497 38.954475222
## classsrw -1.94776958 0.651514098
## density -27.70474380 -15.489565975
## hardness -0.03604523 0.004230411
## size -0.44453927 2.582767006
## weight -0.41058295 -0.168713512
## moisture -0.19391723 0.413047326
##
## , , Sprout
##
## 2.5 % 97.5 %
## (Intercept) 11.78496433 26.552173165
## classsrw -1.20652328 0.757046542
## density -20.53461137 -9.698731394
## hardness -0.03690744 -0.005133494
## size -0.18459306 1.935820104
## weight -0.11978643 0.025152642
## moisture -0.26392179 0.177927888ci.OR2 <- exp(c.value * conf.beta [2:7 ,1:2 ,1])
ci.OR3 <- exp(c.value * conf.beta [2:7 ,1:2 ,2])
round ( data.frame ( low = ci.OR2 [,1], up = ci.OR2 [ ,2]) , 2)## low up
## classsrw 0.14 1.92
## density 0.03 0.13
## hardness 0.37 1.12
## size 0.80 3.55
## weight 0.04 0.26
## moisture 0.67 2.32round ( data.frame ( low = 1/ ci.OR2 [,2], up = 1/ ci.OR2 [ ,1]) , 2)[c(2 ,5) ,]## low up
## density 7.64 38.00
## weight 3.80 25.79round ( data.frame ( low = ci.OR3 [,1], up = ci.OR3 [ ,2]) , 2)## low up
## classsrw 0.30 2.13
## density 0.07 0.28
## hardness 0.36 0.87
## size 0.91 2.59
## weight 0.39 1.22
## moisture 0.58 1.44round ( data.frame (low = 1/ ci.OR3 [,2], up = 1/ ci.OR3 [ ,1]) , 2)[c(2 ,3) ,]## low up
## density 3.57 14.82
## hardness 1.15 2.74A função confint() armazena os limites de confiança em uma matriz tridimensional usando o formato [linha, coluna, tabela]. A linha \(r + 1\) de cada tabela contém os limites inferior e superior da variável explicativa \(r\) nas colunas 1 e 2, respectivamente. Por exemplo, para acessar os intervalos scab vs. healthy para density (\(r = 2\)), precisamos de elementos [3, 1:2, 1].
Para encontrar os intervalos para as razões de chances, multiplicamos cada ponto final do intervalo pela constante apropriada e, em seguida, usamos a função exp(). Observe que o segmento de código c.value*conf.beta[2:7,1:2,1] multiplica cada linha de conf.beta[2:7,1:2,1] por seu elemento correspondente em c.value.
As razões de chances de densidade podem ser interpretadas da seguinte
forma:
Com 95% de confiança, as chances de
scab em vez de um kernel healthy mudam
de 7.64 a 38.00 vezes quando a density é diminuída em
0.13 mantendo as outras variáveis constantes. Além disso, com 95% de
confiança, as chances de sprout em vez de um kernel
healthy mudam em 3.57 e 14.82 vezes quando a densidade
é diminuída em 0.13 mantendo as outras variáveis constantes.
Para a comparação de sarna vs. saudável, apenas os intervalos de confiança da razão de chances de densidade e peso não incluem 1. Para a comparação de broto versus saudável, apenas os intervalos de confiança de razão de chances de densidade e dureza não incluem 1. As comparações de razões de chances podem ser feitas para sarna vs broto também. Este é o assunto do Exercício 11.
Para a comparação sarna (scab) vs. saudável (healthy), apenas os intervalos de confiança da razão de chances de densidade (density) e peso (weight) não incluem 1. Para a comparação de broto (sprout) vs. saudável (healthy), apenas os intervalos de confiança de razão de chances de densidade e dureza (hardness) não incluem.
3.3.2 Tabelas de contingência
O modelo de regressão multinomial fornece uma maneira conveniente de realizar o teste da razão de verossimilhanças (LRT) para independência descrito na Seção 3.2.3. Podemos tratar a variável de linha \(X\) como uma variável explicativa categórica (consulte a Seção 2.2.6) construindo \(I-1\) variáveis indicadoras \(x_2,\cdots,x_I\) representando os níveis \(2,\cdots,I\) de \(X\).
Observe que excluímos o rótulo \(x_1\) por conveniência, para que os índices das variáveis indicadoras correspondam aos níveis de \(X\). Usando \(Y\) como variável de resposta com probabilidades de categorias \(\pi_1,\cdots,\pi_J\), o modelo de regressão multinomial para a tabela é \[ \log\big(\pi_j/\pi_1\big) = \beta_{j0} + \beta_{j2} x_2+\cdots+\beta_{jI} x_I, \]
para \(=j=2,\cdots,J\).
Este modelo implica que \(\log\big(\pi_j/\pi_1\big) = \beta_{j0}\) quando \(X = 1\) e \(\log\big(\pi_j/\pi_1\big) = \beta_{j0}+\beta_{j1}\) quando \(X = i\) para \(i = 2,\cdots,I\). Portanto, o log-odds entre duas colunas depende da linha da tabela. Além disso, observe que existem \(I(J-1)\) parâmetros de regressão no total para o modelo. Porque \(\sum_{j=1}^J \pi_j = 1\) para cada \(X\), ou seja, apenas as \(J-1\) probabilidades de resposta precisam ser estimadas em cada linha, o modelo está saturado.
A independência entre \(X\) e \(Y\) remove essa dependência de linha, de modo que o log-odds entre duas colunas seja constante entre as linhas. Isso implica no modelo simplificado \[ \log\big(\pi_j/\pi_1\big) = \beta_{j0}\cdot \]
Assim, testar a independência é equivalente a testar \[ H_0 \, : \, \beta_{j2} = \cdots= \beta_{jI} = 0, \] para cada \(j = 2,\cdots,J\) vs. \[ H_1 \, : \, \mbox{Pelo menos um } \beta_{ji} \neq 0 \mbox{ para alguns } j \mbox{ e } i\cdot \]
Este teste é facilmente realizado usando o teste da razão de verossimilhanças (LRT). Quando a hipótese nula é rejeitada, o próximo passo é investigar as fontes de dependência.
Como na Seção 3.2.3, o cálculo das razões de chances para seções \(2\times 2\) selecionadas da tabela de contingência pode às vezes ajudar a explicar a natureza da associação. Essas razões de chances são facilmente obtidas a partir dos parâmetros do modelo de regressão, conforme descrito na seção anterior. Por exemplo, as chances estimadas para comparar as categorias de resposta \(Y = j\) a \(Y = 1\) são \(\exp\big(\widehat{\beta}_{ji}\big)\) vezes maiores para \(X = i\) do que \(X = 1\). As razões de chances (odds ratios) envolvendo linhas e colunas diferentes da primeira são encontradas como exponenciais de combinações lineares de parâmetros de regressão. Por exemplo, as probabilidades estimadas para comparar as categorias de resposta 2 e 3 são \[ \exp\big( (\widehat{\beta}_{24}-\widehat{\beta}_{34}) - (\widehat{\beta}_{25}-\widehat{\beta}_{35})\big) \] vezes maiores na linha 4 do que na linha 5. Os intervalos de confiança podem ser formados para essas razões de chances, conforme descrito na Seção 3.3.1.
Exemplo 3.6: Crackers enriquecidos com fibra.
Usamos regressão multinomial para repetir os cálculos para o teste da razão de verossimilhanças (LRT) para independência que foram apresentados no exemplo da Seção 3.2.3.
Começamos ajustando o modelo usando multinom() com a gravidade do inchaço como variável de resposta e fonte de fibra como variável explicativa:
library ( package = nnet )
mod.fit.nom <- multinom ( formula = bloat ~ fiber , weights = count , data = diet )## # weights: 20 (12 variable)
## initial value 66.542129
## iter 10 value 54.519963
## iter 20 value 54.197000
## final value 54.195737
## convergedsummary (mod.fit.nom)## Call:
## multinom(formula = bloat ~ fiber, data = diet, weights = count)
##
## Coefficients:
## (Intercept) fiberbran fibergum fiberboth
## low -0.4057626 -0.1538545 0.4055575 1.322135
## medium -1.0980713 -0.8481379 1.5032639 1.503764
## high -12.4401085 -4.1103893 13.3561038 12.440403
##
## Std. Errors:
## (Intercept) fiberbran fibergum fiberboth
## low 0.6455526 0.8997698 1.190217 1.056797
## medium 0.8163281 1.3451836 1.224593 1.224649
## high 205.2385583 1497.8087307 205.240263 205.240994
##
## Residual Deviance: 108.3915
## AIC: 132.3915A função reconhece que fiber é um fator com níveis ordenados como “none” (nenhum), “bran” (farelo), “gum” (goma) e “both” (ambos). Portanto, ele se encaixa no modelo \[ \log\big(\pi_j/\pi_\mbox{none}\big) = \beta_{j0} + \beta_{j1} \, \mbox{bran} + \beta_{j2} \, \mbox{gum} + \beta_{j3} \, \mbox{both}, \] onde bran, gum e both representam variáveis indicadoras para as categorias de fontes de fibra correspondentes e \(j\) = 2, 3, 4 representa as categorias de resposta low (baixa), medium (média) e high (alta), respectivamente. O argumento weights = count é usado em multinom() porque cada linha da dieta representa contagens de tabelas de contingência, armazenadas na variável “count” em vez de observações de testes individuais.
A função summary() relata as estimativas de parâmetros, os erros padrão correspondentes, \(\sqrt{\widehat{\mbox{Var}}(\widehat{\beta}_{ji})}\). Observe que o desvio residual (deviance) é dado como 108.3915 em vez de 0, mesmo que um modelo saturado esteja sendo ajustado. Este é um efeito colateral de como multinom() usa o argumento weights para criar uma função de verossimilhança. Para esta configuração da tabela de contingência, a função de verossimilhança deve ser um produto de 4 funções de probabilidade multinomiais, uma para cada fonte de fibra, onde o número de tentativas é 12 para cada função de probabilidade. Então as probabilidades de categoria de resposta estimadas pelo modelo serão idênticas às proporções amostrais na tabela e isso causa um desvio residual (deviance) de zero.
A função multinom() reformula as contagens em tentativas multinomiais individuais, ou seja, o formato de “dados brutos” conforme examinado na Seção 1.2.1. Isso leva a 48 linhas no quadro de dados, onde cada linha representa os valores de fibra e inchaço de um indivíduo. A função de verossimilhança torna-se o produto de 48 funções de probabilidade multinomiais com uma tentativa para cada função de probabilidade. As proporções amostrais de cada indivíduo para as categorias de resposta de inchaço são 1, se o indivíduo experimentou esse nível de inchaço ou 0 caso contrário. As probabilidades estimadas do modelo não correspondem a essas proporções, portanto, a deviance residual difere de 0. No entanto, apesar dessa formulação diferente da função de verossimilhança por multinom(), não teremos problemas na análise que se segue. Para esses dados e em geral, estimativas de parâmetros, erros padrão e resultados de LRT para variáveis explicativas são os mesmos, não importa qual função de verossimilhança e formato de dados seja usado. Discutiremos esses dois formatos de dados mais adiante no Capítulo 5.
O teste para independência da fonte de fibra e gravidade do inchaço compara \[ H_0 \, : \, \beta_{j2} = \beta_{j3} = \beta_{j4} = 0 \] para cada \(j = 2, 3, 4\) vs. \[ H_1 \, : \, \mbox{Nem todos os } \beta_{ji} = 0 \mbox{ para alguns } j \mbox{ e } i\cdot \] O teste da razão de verossimilhnças (LRT) é realizado por Anova():
library ( package = car)
Anova (mod.fit.nom)## # weights: 8 (3 variable)
## initial value 66.542129
## final value 63.635876
## converged## Analysis of Deviance Table (Type II tests)
##
## Response: bloat
## LR Chisq Df Pr(>Chisq)
## fiber 18.88 9 0.02623 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O LRT resulta em \(-2\log(\lambda) = 18.9\) com um \(p\)-valor de 0.026. Assim, há evidência moderada de dependência entre o tipo de fibra e a gravidade do inchaço. Observe que esses valores correspondem ao que foi encontrado anteriormente usando a função assocstats().
Para examinar as fontes de associação, podemos começar examinando os parâmetros estimados do modelo. Por exemplo, os log-odds estimadas comparando a categoria de gravidade baixa de inchaço com nenhum inchaço são \[ \log\big(\widehat{\pi}_\mbox{low}/\widehat{\pi}_\mbox{none}\big) = -0.41 - 0.15 \, \mbox{bran} + 0:41 \, \mbox{gum} + 1.32 \, \mbox{both}\cdot \]
As estimativas dos parâmetros sugerem que o uso de gum (goma), com ou sem bran (farelo), leva a maiores chances de inchaço baixo em relação à ausência de inchaço do que o uso sem fibra. Por outro lado, parece haver um efeito um pouco menor na direção oposta devido ao uso de farelo em vez de nenhuma fibra. Calculamos as razões de chances estimadas da seguinte forma:
round (exp( coefficients (mod.fit.nom )[ , -1]) , 2)## fiberbran fibergum fiberboth
## low 0.86 1.5 3.75
## medium 0.43 4.5 4.50
## high 0.02 631658.3 252812.49Por exemplo, as chances estimadas de ter baixo inchaço ao invés de nenhum é \(\exp(-0.15) = 0.86\) vezes maior para usar farelo como fonte de fibra do que não usar nenhuma fibra. É claro que as conclusões finais sobre os efeitos da fibra precisam levar em conta a variabilidade dos estimadores de parâmetros, o que faremos em breve.
Observe que a linha de alto inchaço na tabela de razões de chances contém duas estimativas muito grandes. Além disso, observe que a tabela de erros padrão de summary(mod.fit.nom) contém valores muito grandes na linha de alto inchaço. De fato, se alguém diminuir o valor do critério de convergência dado pelo argumento reltol, o padrão é reltol = 10^(-8) para multinom(), pode-se ver que as estimativas de parâmetro estão aumentando em direção ao infinito para as colunas fibergum e fiberboth da linha de alto inchaço apesar da mensagem convergente impressa na saída.
Esses problemas ocorrem devido às contagens 0 presentes nas células correspondentes na tabela de dados. Vimos na Seção 1.2.5 que 0 contagens podem causar dificuldades na estimativa de odds ratios e esses problemas também ocorrem aqui. Uma solução ad-hoc semelhante à que foi usada na Seção 1.2.5 é adicionar 0.5 às contagens de 0 células. Alternativamente, uma pequena constante, como 0.5, pode ser adicionada a todas as células da tabela de contingência. Em ambos os casos, as estimativas dos parâmetros de regressão são tendenciosas para 0 e, portanto, as razões de chances são tendenciosas para 1. No entanto, a redução na variância das estimativas é substancial e não há necessariamente soluções melhores no âmbito do uso da estimação por máxima verossimilhança para o modelo de regressão multinomial.
O código abaixo mostra como 0.5 é adicionado às células 0 usando a função ifelse(). Em seguida, reestimamos o modelo:
diet$count2 <- ifelse ( test = diet$count == 0, yes = diet$count + 0.5 , no = diet$count )
mod.fit.nom2 <- multinom ( formula = bloat ~ fiber , weights = count2 , data = diet )## # weights: 20 (12 variable)
## initial value 67.928424
## iter 10 value 58.549878
## final value 58.394315
## convergedsum.fit <- summary (mod.fit.nom2 )
round (sum.fit$coefficients , 4)## (Intercept) fiberbran fibergum fiberboth
## low -0.4055 -0.1541 0.4055 1.3218
## medium -1.0986 -0.8473 1.5041 1.5040
## high -2.4849 -0.1542 3.4012 2.4849round (sum.fit$standard.errors , 4)## (Intercept) fiberbran fibergum fiberboth
## low 0.6455 0.8997 1.1902 1.0567
## medium 0.8165 1.3452 1.2247 1.2247
## high 1.4719 2.0759 1.6931 1.7795Vemos que os parâmetros estimados e os erros padrão correspondentes são muito semelhantes aos anteriores, exceto que os altos valores de inchaço estão muito mais próximos de 0 devido à constante adicionada. Abaixo estão as razões de chances estimadas e alguns dos intervalos de confiança correspondentes:
round (exp( coefficients (mod.fit.nom2 )[ , -1]) , 2)## fiberbran fibergum fiberboth
## low 0.86 1.5 3.75
## medium 0.43 4.5 4.50
## high 0.86 30.0 12.00conf.beta <- confint ( object = mod.fit.nom2 , level = 0.95)
round (exp( conf.beta [2:4 , ,3]) ,2) # compare high to no bloating## 2.5 % 97.5 %
## fiberbran 0.01 50.13
## fibergum 1.09 828.46
## fiberboth 0.37 392.52Essas razões de chances são as mesmas que encontramos anteriormente, exceto pelo inchaço alto. Os intervalos de confiança para essa categoria de resposta são bastante amplos, questionando sua utilidade. Isso acontece porque as razões de chances e os intervalos de confiança são baseados em contagens tão pequenas que há uma incerteza considerável em relação aos valores da população.
É interessante notar que há uma tendência aparente entre as razões de chances estimadas para goma vs. nenhuma e ambas vs. nenhuma. Vemos que, à medida que a gravidade do inchaço aumenta, o mesmo acontece com as razões de chances estimadas. Isso pode não ser uma coincidência. Observe que a resposta da gravidade do inchaço é na verdade uma variável ordinal. Na Seção 3.4 mostraremos como adaptar nossa modelagem para variáveis de resposta ordinais que permitem uma análise potencialmente mais poderosa, porém parcimoniosa, a ser conduzida. Uma vantagem adicional de modelar a estrutura ordinal é que o impacto das contagens 0 individuais é substancialmente reduzido.
Curiosamente, as probabilidades, chances e razões de chances estimadas produzidas a partir da equação \[ \log\big(\pi_j/\pi_1\big) = \beta_{j0} + \beta_{j2} x_2+\cdots+\beta_{jI} x_I, \]
para \(=j=2,\cdots,J\) são exatamente as mesmas obtidas diretamente da tabela de contingência, desde que não haja contagens 0.
Por exemplo, no último exemplo, a razão de chances estimada de ter baixo inchaço em vez de nenhum é \[ n_{11}n_{22}/(n_{12}n_{21}) = 6\times 4/(4\times 7) = 0.86 \] vezes maior para usar farelo como fonte de fibra do que usando nenhuma fibra em tudo.
Este valor de odds ratio estimado é exatamente o mesmo que \(\exp(-0.15) = 0.86\) calculado usando o modelo de regressão multinomial estimado. Nesse sentido, o modelo multinomial parece ser uma forma complicada de realizar a análise. Isso é verdade para o problema relativamente simples de relacionar uma resposta multinomial a uma única variável explicativa categórica. A força do procedimento de modelagem é que ele pode ser usado em problemas mais complexos onde a simples análise baseada em tabela não pode.
Em particular, o modelo de regressão multinomial pode ser usado para avaliar os efeitos de muitas variáveis arbitrariamente na resposta, sujeito ao tamanho adequado da amostra. As variáveis podem ser numéricas, conforme mostrado anteriormente nesta seção, categóricas, conforme mostrado nesta subseção ou qualquer combinação destes.
Quando variáveis explicativas categóricas adicionais estão disponíveis, podemos examinar os dados em tabelas de contingência de dimensão superior. Por exemplo, considere uma configuração com uma tabela de contingência tridimensional resumindo as contagens para variáveis categóricas \(X\), \(Y\) e \(Z\), com níveis \(i = 1,\cdots,I\), \(j = 1,\cdots,J\) e \(k = 1,\cdots,K\), respectivamente.
Uma análise usando as ferramentas mais simples da Seção 3.2 requer que tabelas de contingência separadas resumindo qualquer par de variáveis, por exemplo, \(X\) e \(Y\) sejam formadas para cada nível da terceira variável, por exemplo, \(Z\). Em alguns casos, existe uma maneira “óbvia” de fazer isso, como quando \(X\) representa uma variável de tratamento, \(Y\) uma variável de resposta e \(Z\) uma variável de bloqueio ou estratificação. Então, modelar a relação entre \(X\) e \(Y\) controlando o nível de \(Z\) é a abordagem natural para a análise.
Um teste de Pearson ou da razão de verossimilhanças (LR) para independência pode ser realizado em cada tabela e, posteriormente, as razões de chances podem ser estimadas assim como outras análises. No entanto, esta abordagem é limitada em seu escopo. Permite avaliações apenas da associação entre \(X\) e \(Y\) condicionalmente ao nível de \(Z\). Não faz nenhuma tentativa para comparar essas associações entre os níveis \(Z\), nem combiná-las em uma única associação que se aplique a todos os níveis \(Z\). Embora existam métodos baseados em tabelas para resolver esta última questão, eles ainda não permitem considerar a associação entre \(Z\) e \(Y\) , o que é de interesse quando \(Z\) é outra variável explicativa em vez de do que uma variável de estratificação.
A estrutura de regressão multinomial e a estrutura do modelo de regressão de Poisson, a ser introduzida na Seção 4.2.4, permitem que essas questões mais gerais sejam respondidas. Por exemplo, um modelo geral que permite a associação de \(X\) e \(Z\) com \(Y\) é \[ \log\big(\pi_j/\pi_1 \big)=\beta_{j0}+\beta_{j2}^X x_2+\cdots+\beta_{jI}^X x_I+\beta_{j2}^Z z_2+\cdots+\beta_{jK}^Z z_K \] \(j=2,\cdots,J\) onde usamos \(x_2,\cdots,x_I\) e \(z_2,\cdots,z_K\) como variáveis indicadoras para os níveis dados em seus índices.
O sobrescrito \(X\) ou \(Z\) é adicionado aos parâmetros \(\beta\) para deixar claro qual associação da variável com \(Y\) é descrita pelo parâmetro. A exclusão das variáveis indicadoras \(X\) e/ou \(Z\) do modelo corresponde à independência entre essa variável explicativa e \(Y\). Uma interação entre \(X\) e \(Z\) também pode ser incluída no modelo para permitir que a associação de cada variável com \(Y\) varie entre os níveis da outra variável explicativa. A adição de variáveis explicativas categóricas adicionais é uma extensão direta deste modelo, assim como a adição de variáveis explicativas numéricas.
Exemplo 3.7: Crackers enriquecidos com fibra.
Conforme mencionado na Seção 3.2, a variável explicativa fibra é, na verdade, construída a partir de combinações de duas variáveis separadas: farelo “sim” ou “não” e goma “sim” ou “não”. Essa estrutura é frequentemente chamada de fatorial \(2\times 2\). Em geral, um fatorial \(a\times b\) consiste em todas as combinações de uma variável com \(a\) níveis e outra com \(b\) níveis. As variáveis que compõem as combinações são chamadas de “fatores” ver, por exemplo, Milliken and Johnson (2004) para detalhes.
É típico nessas estruturas fazer perguntas sobre os efeitos separados de cada variável constituinte nas respostas, bem como procurar interações entre as variáveis. Para prosseguir com essa análise, primeiro transformamos a variável explicativa fibra em farelo e goma usando a função ifelse():
diet$bran <- factor ( ifelse ( test = diet$fiber == "bran" |
diet$fiber == "both", yes = "yes", no = "no"))
diet$gum <- factor ( ifelse ( test = diet$fiber == "gum" |
diet$fiber == "both", yes = "yes", no = "no"))
head (diet , n = 4)## fiber bloat count count2 bran gum
## 1 bran high 0 0.5 yes no
## 2 gum high 5 5.0 no yes
## 3 both high 2 2.0 yes yes
## 4 none high 0 0.5 no noA única linha vertical | dentro do argumento test significam “ou”, de modo que a primeira linha atribui a bran (farelo) um valor “sim” se fiber = “bran” ou fiber = “both”; caso contrário, bran é definido como “no”. Além disso, colocamos ifelse() dentro de factor() para que o resultado seja um tipo de classe de fator em vez de um caractere. Se isso não for feito, mensagens de aviso serão produzidas pelo multinom() enquanto ele realiza essa conversão internamente.
Em seguida, ajustamos o modelo multinomial que inclui efeitos para bran (farelo), gum (goma) e sua interação:
mod.fit.nom.inter <- multinom ( formula = bloat ~ bran + gum + bran :gum ,
weights = count , data = diet )## # weights: 20 (12 variable)
## initial value 66.542129
## iter 10 value 54.406806
## iter 20 value 54.196639
## final value 54.195746
## convergedsummary (mod.fit.nom.inter )## Call:
## multinom(formula = bloat ~ bran + gum + bran:gum, data = diet,
## weights = count)
##
## Coefficients:
## (Intercept) branyes gumyes branyes:gumyes
## low -0.4063036 -0.1532825 0.4076874 1.0674733
## medium -1.0973831 -0.8503013 1.5033271 0.8503872
## high -12.4349471 -2.0291752 13.3518350 1.1121073
##
## Std. Errors:
## (Intercept) branyes gumyes branyes:gumyes
## low 0.6456363 0.899801 1.190275 1.58416
## medium 0.8160908 1.345541 1.224613 1.86471
## high 204.7028906 561.368770 204.704601 561.37028
##
## Residual Deviance: 108.3915
## AIC: 132.3915Anova (mod.fit.nom.inter )## Analysis of Deviance Table (Type II tests)
##
## Response: bloat
## LR Chisq Df Pr(>Chisq)
## bran 2.5857 3 0.460010
## gum 16.2897 3 0.000989 ***
## bran:gum 0.4880 3 0.921514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logLik (mod.fit.nom.inter )## 'log Lik.' -54.19575 (df=12)O teste da razão de verossimilhanças (LRT) para a interação fornece um \(p\)-valor de 0.9215, indicando que não há evidência suficiente de uma interação. Além disso, vemos que o bran tem um \(p\)-valor grande, mas gum tem um \(p\)-valor pequeno. Assim, há evidências suficientes para indicar que a goma tem efeito sobre a gravidade do inchaço.
Como o modelo saturado está sendo ajustado aqui, os problemas com contagens zero na tabela que foram observados no exemplo anterior continuam sendo um problema neste modelo. No final, essas dificuldades não são uma preocupação para o LRT, porque a estatística \(-2\log(\Lambda )\) mudará muito pouco para um número maior de iterações. Para as duas ocorrências de contagem 0, a estimativa de \(\pi_j\) do modelo de regressão multinomial já está muito próxima de 0. A estimativa só pode ficar um pouco mais próxima de 0 se o critério de convergência for mais estrito do que o padrão usado aqui.
Como alternativa à estrutura multinomial deste capítulo, a análise de tabelas de contingência é frequentemente realizada assumindo que as contagens têm uma distribuição de Poisson ou variantes dela. Discutiremos esses modelos de regressão de Poisson no Capítulo 4.
3.4 Modelos de regressão de resposta ordinal
Muitas variáveis de resposta categórica têm uma ordenação natural para seus níveis. Por exemplo, uma variável de resposta pode ser medida usando uma escala Likert com categorias “discordo totalmente”, “discordo”, “neutro”, “concordo” ou “concordo totalmente”. Se os níveis de resposta puderem ser organizados de modo que a categoria 1 < categoria 2 < \(\cdots\) < categoria \(J\) em alguma escala conceitual de medida (por exemplo, quantidade de concordância), então os modelos de regressão podem incorporar essa ordenação por meio de uma variedade de transformações logit da resposta probabilidades. Nesta seção, focamos na modelagem de probabilidades acumuladas com base na ordenação das categorias.
A probabilidade acumulada para a categoria \(j\) de \(Y\) é \(P(Y\leq j) = \pi_1+\cdots +\pi_j\) para \(j = 1,\cdots,J\). Observe que \(P(Y\leq J) = 1\). Modelos de regressão para respostas multinomiais ordinais podem examinar os efeitos das variáveis explicativas \(x_1,\cdots, x_p\) nas probabilidades ou log-odds acumulados, também chamadas de logits cumulativos, \[ logit\big( P(Y\leq j)\big)=\log\left( \dfrac{P(Y\leq j)}{1-P(Y\leq j)}\right)=\log\left( \dfrac{\pi_1+\cdots+\pi_j}{\pi_{j+1}+\cdots+\pi_J}\right)\cdot \]
Em particular, o modelo de probabilidades proporcionais é um modelo especial que assume que o logit dessas probabilidades acumuladas muda linearmente à medida que as variáveis explicativas mudam e também que a inclinação dessa relação é a mesma independentemente da categoria \(j\).
Formalmente, o modelo é definido como \[ logit\big( P(Y\leq j)\big) =\beta_{j0}+\beta_1 x_1+\cdots+\beta_p x_p, \] para \(j=1,\cdots,J\).
Observe que não há subscritos \(j\) nos parâmetros \(\beta_1,\cdots,\beta_p\). O modelo assume que os efeitos das variáveis explicativas são os mesmos, independentemente de quais probabilidades acumuladas são usadas para formar os logaritmos. Assim, o nome “proporcional odds” deriva de cada odd ser um múltiplo de \(\exp(\beta_{j0})\).
Para um \(j\) fixo, aumentar \(x_r\) por \(c\) unidades muda cada log-odds na equação acima por \(c\, \beta_r\) ao manter as outras variáveis explicativas constantes. Por outro lado, a diferença no log-odds entre as categorias de resposta \(j\) e \(j\,'\) é constante, \(\beta_{j0}-\beta_{j\, '0}\) e não depende dos valores de \(x_1,\cdots,x_p\) quando eles são mantidos fixos.
Esses resultados estão diretamente relacionados às razões de chances, conforme detalhado na Seção 3.4.1. Além disso, observe que as chances devem aumentar à medida que \(j\) aumenta, porque colocamos progressivamente mais probabilidade no numerador, \(P(Y\leq j)\). Isso implica que \(\beta_{10} <\cdots < \beta_{J-1,0}\).
As probabilidades de observar uma determinada categoria de resposta \(j\) são encontradas observando que \[ \pi_j = P(Y=j) = P(Y\leq j)-P(Y\leq j-1) \]
onde \(P(Y\leq 0)=0\), \(P(Y\leq J)=1\) e \[ P(Y\leq j) =\dfrac{\exp\big( \beta_{j0}+\beta_1 x_1+\cdots+\beta_p x_p\big)}{1+\exp\big(\beta_{j0}+\beta_1 x_1+\cdots+\beta_p x_p\big)}\cdot \]
Por exemplo, a probabilidade para a categoria 1 é \[ \pi_1=P(Y\leq 1)-P(Y\leq 0) = \dfrac{\exp\big( \beta_{10}+\beta_1 x_1+\cdots+\beta_p x_p\big)}{1+\exp\big(\beta_{10}+\beta_1 x_1+\cdots+\beta_p x_p\big)} \]
e a probabilidade para a categoria \(J\) é \[ \pi_J =P(Y\leq J)-P(Y\leq J-1) = \dfrac{\exp\big( \beta_{J-1,0}+\beta_1 x_1+\cdots+\beta_p x_p\big)}{1+\exp\big(\beta_{J-1,0}+\beta_1 x_1+\cdots+\beta_p x_p\big)}\cdot \]
Exemplo 3.8: Gráficos de modelo de probabilidades proporcionais.
O objetivo deste exemplo é examinar as características do modelo de probabilidades proporcionais. Considere o modelo \[ logit\big(P(Y\leq j)\big) = \beta_{j0}+ \beta_1 x_1, \]
onde \(J = 4\), \(\beta_{10} = 0\), \(\beta_{20} = 2\), \(\beta_{30} = 4\) e \(\beta_1 = 2\). A figura mostra as probabilidades cumulativas e as probabilidades de cada categoria individual para este modelo. Observe que as curvas de probabilidade cumulativa no gráfico à esquerda têm exatamente a mesma forma, porque \(\beta_1\) é o mesmo para cada categoria de resposta \(j\).
O deslocamento horizontal é devido aos diferentes valores para \(\beta_{j0}\). Por exemplo, \(P(Y\leq j) = 0.5\) quando \(\beta_{j0} + \beta_1 x_1 = 0\), ou seja, quando \(x_1 = -\beta_{j0}/\beta_1\).
Probabilidades de resposta para categorias individuais \(\pi_j\), \(j = 1,2,3,4\), são encontrados em cada \(x_1\) das distâncias verticais entre curvas consecutivas no gráfico à esquerda, juntamente com \(P(Y\leq 0) = 0\) e \(P(Y\leq 4) = 1\). Estes são mostrados no gráfico à direita na figura, que também ilustra as ordenações entre as categorias para \(Y\) em relação a \(x_1\).
Por exemplo, os menores valores de \(x_1\) levam a \(Y = 4\) com a maior probabilidade. À medida que \(x_1\) cresce em tamanho, a probabilidade de observar \(Y = 3\) torna-se maior até que, eventualmente, tenha a maior probabilidade para um determinado intervalo de \(x_1\). Resultados semelhantes ocorrem para \(Y = 2\) e \(Y = 1\) nesta ordem.
# Plot of proportional odds model
# Remember that beta10 < beta20 < beta30
beta<-c(0, 2, 4, 2) #beta10, beta20, beta30, beta1
x.range<-c(-5, 3)
par(mfrow = c(1, 2),mar=c(4,4,2,2))
curve(expr = plogis(q = beta[1] + beta[4]*x), xlim = x.range, ylab = expression(P(Y < j)),
xlab = expression(x[1]), main = "Probabilidades acumuladas de Y", lwd = 2)
curve(expr = plogis(q = beta[2] + beta[4]*x), add = TRUE, lty = "dashed", col = "red", lwd = 2)
curve(expr = plogis(q = beta[3] + beta[4]*x), add = TRUE, lty = "dotted", col = "blue", lwd = 2)
legend(x = -5.5, y = 0.9, legend = c(expression(P(Y < 2)), expression(P(Y<3)), expression(P(Y<4))),
lty = c("solid", "dashed", "dotted", "dotdash"), col = c("black", "red", "blue"),
bty = "n", lwd = 2)
grid()
curve(expr = plogis(q = beta[1] + beta[4]*x), xlim = x.range, ylab = expression(pi[j]),
xlab = expression(x[1]), main = "Probabilidades de Y", lwd = 2)
curve(expr = plogis(q = beta[2] + beta[4]*x) - plogis(q = beta[1] + beta[4]*x), add = TRUE,
lty = "dashed", col = "red", lwd = 2)
curve(expr = plogis(q = beta[3] + beta[4]*x) - plogis(q = beta[2] + beta[4]*x), add = TRUE,
lty = "dotted", col = "blue", lwd = 2)
curve(expr = 1 - plogis(q = beta[3] + beta[4]*x), add = TRUE,
lty = "dotdash", col = "green", lwd = 2)
legend(x = -5.5, y = 0.9, legend = c(expression(pi[1]), expression(pi[2]), expression(pi[3]),
expression(pi[4])), lty = c("solid", "dashed", "dotted", "dotdash"),
col = c("black", "red", "blue", "green"), bty = "n", lwd = 2)
grid()
Os parâmetros do modelo de probabilidades proporcionais são estimados usando a máxima verossimilhança. Semelhante à Seção 3.3, a função de verossimilhança para uma amostra de tamanho \(m\) é simplesmente o produto de \(m\) distribuições multinomiais com parâmetros de probabilidade expressos como funções das variáveis explicativas. Procedimentos numéricos iterativos são então usados para ajustar o modelo.
A função polr() do pacote MASS executa os cálculos necessários. É importante garantir que os níveis da resposta categórica sejam ordenados da maneira desejada ao usar polr(); caso contrário, a ordenação dos níveis de \(Y\) não será considerada corretamente. Demonstraremos como verificar e, se necessário, alterar a ordenação no próximo exemplo. A matriz de covariância para as estimativas de parâmetros de regressão é encontrada usando procedimentos de verossimilhança padrão.
Os procedimentos de inferência baseados em Wald e LR também são realizados da maneira usual. As hipóteses associadas aos testes de parâmetros de regressão individuais são \[ H_0 \, : \, \beta_r = 0 \qquad \mbox{vs.} \qquad H_1 \, : \, \beta_r\neq 0\cdot \] Se a hipótese nula for verdadeira, isso diz que \(J-1\) log-odds comparando \(P(Y\leq j)\) a \(P (Y > j)\) não dependem de \(x_r\), mantendo todas as outras variáveis explicativas constantes. Se a hipótese alternativa for verdadeira, então as log-odds para cada probabilidade cumulativa aumentam ou diminuem com \(x_r\) dependendo do sinal de \(\beta_r\). Por sua vez, isso impõe uma ordenação nas probabilidades de categorias individuais, como a vista no gráfico à direita da figura no exemplo anterior.
Compare o teste anterior envolvendo \(\beta_r\) com um teste correspondente envolvendo o modelo de regressão multinomial com uma resposta nominal da Seção 3.3, \[ H_0 \, : \, \beta_{2r} = \cdots= \beta_{Jr} = 0 \qquad \mbox{vs.} \qquad H_1 \, : \, \mbox{Pelo menos um } \beta_{jr}\neq 0\cdot \] O modelo da hipótese alternativa coloca menos restrições sobre como \(x_r\) se relaciona com probabilidades de categorias individuais por meio do uso de mais parâmetros do que o modelo de probabilidades proporcionais. Assim, quando as hipóteses de probabilidades proporcionais são aplicáveis, o modelo de hipóteses alternativas para regressão multinomial descreve a relação entre a resposta e as variáveis explicativas de forma menos eficiente.
Exemplo 3.9: Grãos de trigo.
Grãos saudáveis são o tipo mais desejável de grãos de trigo. Pode-se também argumentar que a presença de doenças torna os grãos de sarna menos desejáveis do que os grãos germinados. Assim, usamos a ordenação de sarna (Y = 1) < sprout (broto) (Y = 2) < healthy (saudável) (Y = 3) e ajustamos um modelo de probabilidades proporcionais a esses dados.
Começamos criando a ordenação necessária dos níveis para a variável type. Existem várias maneiras de fazer isso e mostramos uma abordagem abaixo que é bastante direta:
levels ( wheat$type )## [1] "Healthy" "Scab" "Sprout"wheat$type.order <- factor ( wheat$type , levels = c("Scab", "Sprout", "Healthy"))
levels ( wheat$type.order )## [1] "Scab" "Sprout" "Healthy"A nova variável de trigo wheat chamada type.order contém o fator com a ordenação adequada de seus níveis. O modelo de probabilidades proporcionais é então estimado usando polr():
library ( package = MASS )
mod.fit.ord <- polr ( formula = type.order ~ class + density + hardness +
size + weight + moisture , data = wheat , method = "logistic")
summary (mod.fit.ord)## Call:
## polr(formula = type.order ~ class + density + hardness + size +
## weight + moisture, data = wheat, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## classsrw 0.17370 0.391764 0.4434
## density 13.50534 1.713009 7.8840
## hardness 0.01039 0.005932 1.7522
## size -0.29253 0.413095 -0.7081
## weight 0.12721 0.029996 4.2411
## moisture -0.03902 0.088396 -0.4414
##
## Intercepts:
## Value Std. Error t value
## Scab|Sprout 17.5724 2.2460 7.8237
## Sprout|Healthy 20.0444 2.3395 8.5677
##
## Residual Deviance: 422.4178
## AIC: 438.4178library ( package = car)
Anova (mod.fit.ord)## Analysis of Deviance Table (Type II tests)
##
## Response: type.order
## LR Chisq Df Pr(>Chisq)
## class 0.197 1 0.65749
## density 98.437 1 < 2.2e-16 ***
## hardness 3.084 1 0.07908 .
## size 0.499 1 0.47982
## weight 18.965 1 1.332e-05 ***
## moisture 0.195 1 0.65872
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O valor do argumento method = “logistic” instrui R a usar a transformação logit nas probabilidades cumulativas. A função polr() estima o modelo \[ logit\big(P(Y\leq j)\big) = \beta_{j0} - \eta_1 x_1 - \cdots -\eta_p x_p, \]
onde \(-\eta_r\) é \(\beta_r\) na equação \[ logit\big( P(Y\leq j)\big) =\beta_{j0}+\beta_1 x_1+\cdots+\beta_p x_p, \] para \(j=1,\cdots,J\).
Devido a essas diferenças de notação, precisamos, portanto, inverter o sinal de todos os valores de \(\eta_r\) estimados na saída para estabelecer o modelo estimado. Isso resulta no modelo \[ logit\big(\widehat{P}(Y\leq j)\big) = \widehat{\beta}_{j0}-0.17 \, \mbox{srw} - 13.51 \, \mbox{density} - 0.010 \, \mbox{hardness} + 0.29 \, \mbox{size} \\ - 0.13 \, \mbox{weight} + 0.039 \, \mbox{moisture}, \]
onde \(\widehat{\beta}_{10} = 17.57\) e \(\widehat{\beta}_{20} = 20.04\).
A coluna de t value na tabela Coefficients fornece as estatísticas Wald para testar \(H_0 \, : \, \beta_r = 0\) vs. \(H_1 \, : \, \beta_r\neq 0\) para \(r = 1,\cdots, 6\). A função Anova() fornece os LRTs correspondentes. Devido aos grandes valores dos testes estatísticos para density e weight, há evidências suficientes de que essas são variáveis explicativas importantes. Há evidências marginais de que a hardness também é importante. Assim como em outros testes como este nos capítulos anteriores, cada um desses testes está condicionado às outras variáveis que estão no modelo.
A função predict() encontra as probabilidades ou classes estimadas para cada categoria de resposta:
pi.hat.ord <- predict ( object = mod.fit.ord , type = "probs")
head (pi.hat.ord )## Scab Sprout Healthy
## 1 0.03661601 0.2738502 0.6895338
## 2 0.03351672 0.2576769 0.7088064
## 3 0.08379891 0.4362428 0.4799583
## 4 0.01694278 0.1526100 0.8304472
## 5 0.11408176 0.4899557 0.3959626
## 6 0.02874814 0.2308637 0.7403882head ( predict ( object = mod.fit.ord , type = "class"))## [1] Healthy Healthy Healthy Healthy Sprout Healthy
## Levels: Scab Sprout HealthyPor exemplo, a probabilidade estimada de estar saudável (healthy) para a primeira observação é \[ \widehat{\pi}_\mbox{Healthy} = 1-\dfrac{\exp\big(20.04 -0.17+\cdots + 0.04\times 12.02 \big)}{1+\exp\big(20.04 -0.17+\cdots + 0.04\times 12.02 \big)} = 0.6895\cdot \]
Essa probabilidade e outras são encontradas usando o argumento type = “probs” em predict(). O argumento type = “class” (o padrão) fornece a classificação prevista com base em sua maior probabilidade estimada. Para o primeiro kernel, Healthy tem a maior probabilidade estimada, então o kernel está previsto para esta condição, que é uma previsão correta.
Assim como os objetos de ajuste de modelo de multinom(), a função predict() não fornece as variações estimadas necessárias para formar intervalos de confiança Wald para \(\pi_\mbox{Scab}\), \(\pi_\mbox{Sprout}\) ou \(\pi_\mbox{Healthy}\). Além disso, a função deltaMethod.polr() não permite nenhuma função matemática dos termos de intercepto e outros problemas surgem porque polr() estima um modelo conforme definido na equação \[ logit\big(P(Y\leq j)\big) = \beta_{j0} - \eta_1 x_1 - \cdots -\eta_p x_p, \] em vez da equação \[ logit\big( P(Y\leq j)\big) =\beta_{j0}+\beta_1 x_1+\cdots+\beta_p x_p\cdot \] Por essas razões, construímos uma nova função deltaMethod.polr2() que calcula \(\widehat{\pi}_j\) e \(\widehat{\mbox{Var}}(\widehat{\pi}_j)\):
deltaMethod.polr2 <- function (object , g) {
# All beta -hat ’s where the slope estimates are adjusted
beta.hat <- c(- object$coefficients , object$zeta )
# Count the number of slopes and intercepts
numb.slope <- length ( object$coefficients )
numb.int <- length ( object$zeta )
# Name the corresponding parameters
names ( beta.hat) <- c( paste ("b", 1: numb.slope , sep ="") ,
paste ("b", 1: numb.int , "0" , sep ="") )
# Fix covariance matrix - All covariances between slopes and intercepts need to be multiplied by -1
cov.mat <- vcov ( object )
# Above diagonal
cov.mat [1: numb.slope , ( numb.slope + 1) :( numb.slope + numb.int)] <-
-cov.mat [1: numb.slope , ( numb.slope + 1) :( numb.slope + numb.int)]
# Below diagonal
cov.mat [( numb.slope + 1) :( numb.slope + numb.int), 1: numb.slope ] <-
-cov.mat [( numb.slope + 1) :( numb.slope + numb.int), 1: numb.slope ]
# deltaMethod . default () method function completes calculations
deltaMethod ( object = beta.hat , g = g, vcov. = cov.mat)
}A função é usada de maneira muito semelhante a deltaMethod() na Seção 3.3. O objeto de ajuste do modelo de polr() é fornecido no argumento object e a cadeia de caracteres para \(\widehat{\pi}_j\) é fornecida no argumento g.
No entanto, não há argumento parameterNames porque os nomes dos parâmetros são criados automaticamente dentro da função. Somente nomes da forma bj0 e br podem ser usados para \(\beta_{j0}\) e \(\beta_r\), respectivamente. Calculamos um intervalo para \(\pi_\mbox{Healthy}\) da seguinte forma:
x1 <- 0; x2 <- wheat [1 ,2]; x3 <- wheat [1 ,3]
x4 <- wheat [1 ,4]; x5 <- wheat [1 ,5]; x6 <- wheat [1 ,6]
g.healthy <- "1 - exp(b20 + b1*x1 + b2*x2 + b3*x3 + b4*x4 + b5*x5 + b6*x6) /
(1 + exp (b20 + b1*x1 + b2*x2 + b3*x3 + b4*x4 + b5*x5 + b6*x6))"
calc.healthy <- deltaMethod.polr2 ( object = mod.fit.ord , g = g.healthy )
calc.healthy$Estimate # pi - hat_Healthy## [1] 0.6895338calc.healthy$SE # sqrt (Var -hat(pi - hat_Healthy ))## [1] 0.08199814alpha <- 0.05
calc.healthy$Estimate + qnorm (p = c( alpha /2, 1- alpha /2)) * calc.healthy$SE## [1] 0.5288204 0.8502472O intervalo de confiança de 95% para \(\pi_\mbox{Healthy}\) é (0.5288, 0.8502). Os intervalos de confiança de 95% para \(\pi_\mbox{Scab}\) e \(\pi_\mbox{Sprout}\) são (0.0061, 0.0671) e (0.1380, 0.4097), respectivamente (consulte o programa para obter o código).
Um gráfico das probabilidades estimadas para cada categoria versus uma variável explicativa individual, mantendo as outras variáveis no modelo constantes, pode ser usado para interpretar o modelo. O Exercício 19 examina gráficos desta forma. Simplificamos o problema aqui usando novamente a densidade como a única variável explicativa no ajuste de um novo modelo de probabilidades proporcionais: \[ logit\big(\widehat{P}(Y\leq j)\big) = \widehat{\beta}_{j0} - 15.64\, \mbox{density}, \] com \(\widehat{\beta}_{10} = 17.41\) e \(\widehat{\beta}_{20} = 19.63\).
A Figura 3.5 apresenta um gráfico desse modelo junto com o modelo de regressão multinomial estimado anteriormente que trata a resposta como nominal. Esses modelos produzem curvas de probabilidade um pouco semelhantes com a mesma crosta < broto < ordenação saudável para densidade com ambos os modelos. Isso fornece alguma garantia de que a ordenação usada para o modelo de probabilidades proporcionais é apropriada e que a suposição de probabilidades proporcionais entre probabilidades cumulativas é razoável.
# Parallel coordinate plot
library(package = MASS) # Location of parcoord() function
# Reorder variables because class is binary (may distort plot)
# Create indicator variable for class
wheat2<-data.frame(kernel = 1:nrow(wheat), wheat[,2:6],
class.new = ifelse(test = wheat$class == "hrw", yes = 0, no = 1))
head(wheat2)## kernel density hardness size weight moisture class.new
## 1 1 1.349253 60.32952 2.30274 24.6480 12.01538 0
## 2 2 1.287440 56.08972 2.72573 33.2985 12.17396 0
## 3 3 1.233985 43.98743 2.51246 31.7580 11.87949 0
## 4 4 1.336534 53.81704 2.27164 32.7060 12.11407 0
## 5 5 1.259040 44.39327 2.35478 26.0700 12.06487 0
## 6 6 1.300258 48.12066 2.49132 33.2985 12.18577 0# Colors by condition:
wheat.colors<-ifelse(test = wheat$type=="Healthy", yes = "black",
no = ifelse(test = wheat$type=="Sprout", yes = "red", no = "green"))
# Line type by condition:
wheat.lty<-ifelse(test = wheat$type=="Healthy", yes = "solid",
no = ifelse(test = wheat$type=="Sprout", yes = "longdash", no = "dotdash"))
parcoord(x = wheat2, col = wheat.colors, lty = wheat.lty) # Plot
legend(x = 6.15, y = 0.75, legend = c("Healthy", "Sprout", "Scab"),
lty = c("solid", "longdash", "dotdash"), col=c("black", "red", "green"), cex=0.8, bty="n")
# Principal cordinate analysis in order to plot in 3D
save<-princomp(formula = ~ density + class.new + hardness + size + weight + moisture,
data = wheat2, cor = TRUE, scores = TRUE)
summary(save, loadings = TRUE, cutoff = 0.0)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.4718974 1.3126735 0.9593345 0.8454985 0.5333147
## Proportion of Variance 0.3610804 0.2871853 0.1533871 0.1191446 0.0474041
## Cumulative Proportion 0.3610804 0.6482656 0.8016527 0.9207973 0.9682014
## Comp.6
## Standard deviation 0.43679671
## Proportion of Variance 0.03179856
## Cumulative Proportion 1.00000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## density 0.286 0.306 0.622 0.650 0.044 0.114
## class.new 0.391 -0.536 0.155 0.004 -0.720 -0.132
## hardness -0.361 0.237 0.660 -0.525 -0.186 -0.260
## size 0.440 0.461 -0.087 -0.417 -0.235 0.598
## weight 0.558 0.327 -0.157 -0.135 0.158 -0.717
## moisture 0.360 -0.493 0.349 -0.331 0.604 0.174wheat.colors<-ifelse(test = wheat$type=="Healthy", yes = "black",
no = ifelse(test = wheat$type=="Sprout", yes = "red", no = "green"))
symbols(x = save$scores[,1], y = save$scores[,2], circles = save$scores[,3]-min(save$scores[,3]),
inches=0.25, xlab = "Principal component 1", ylab = "Principal component 2",
fg = wheat.colors, main = "Bubble plot for first three principal components \n Wheat data")
# Note that circles can not take on a negative value!
abline(h = 0, lty = 1, lwd = 2)
abline(v = 0, lty = 1, lwd = 2)
text(x = save$scores[,1], y = save$scores[,2], col = 2, cex = 0.5) # Put kernel number on plot
legend(2,-2, legend=c("Healthy", "Sprout", "Scab"), pch = c(1,1,1),
col=c("black", "red", "green"), cex=1, bty="n")
grid()
Transformações e termos de interação são incluídos em polr() usando a mesma sintaxe que em glm() e multinom().
Existem várias outras funções que podem estimar um modelo de probabilidades proporcionais. Mostramos como usar a função lrm() do pacote rms no programa correspondente a este exemplo. Também mostramos como usar a função vglm() do pacote VGAM em um exemplo posterior nesta seção.
O pacote mcprofile não pode ser usado para métodos de inferência baseados em LR com objetos resultantes de polr().
# Plot density only model
# Estimate model with density only
mod.fit.nom.density<-multinom(formula = type ~ density, data = wheat)## # weights: 9 (4 variable)
## initial value 302.118379
## iter 10 value 229.769334
## iter 20 value 229.712304
## final value 229.712290
## convergedsummary(mod.fit.nom.density)## Call:
## multinom(formula = type ~ density, data = wheat)
##
## Coefficients:
## (Intercept) density
## Scab 29.37827 -24.56215
## Sprout 19.12165 -15.47633
##
## Std. Errors:
## (Intercept) density
## Scab 3.676892 3.017842
## Sprout 3.337092 2.691429
##
## Residual Deviance: 459.4246
## AIC: 467.4246beta.hat<-coefficients(mod.fit.nom.density)
curve(expr = 1/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) + exp(beta.hat[2,1] + beta.hat[2,2]*x)),
ylab = expression(hat(pi)), xlab = "Density", xlim = c(min(wheat$density), max(wheat$density)),
col = "black", lty = "solid", lwd = 2, n = 1000, type = "n",
panel.first = grid(col = "gray", lty = "dotted"))
# Plot each pi_j
curve(expr = 1/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) + exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "black", lty = "solid", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Healthy"]),
max(wheat$density[wheat$type == "Healthy"]))) # Healthy
curve(expr = exp(beta.hat[1,1] + beta.hat[1,2]*x)/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) +
exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "green", lty = "dotdash", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Scab"]),
max(wheat$density[wheat$type == "Scab"]))) # Scab
curve(expr = exp(beta.hat[2,1] + beta.hat[2,2]*x)/(1 + exp(beta.hat[1,1] + beta.hat[1,2]*x) +
exp(beta.hat[2,1] + beta.hat[2,2]*x)),
col = "red", lty = "longdash", lwd = 2, n = 1000, add = TRUE,
xlim = c(min(wheat$density[wheat$type == "Sprout"]),
max(wheat$density[wheat$type == "Sprout"]))) # Sprout
legend(x = 1.4, y = 0.8, legend=c("Healthy", "Sprout", "Scab"), lty=c("solid","longdash","dotdash"),
col=c("black","red","green"), bty="n", lwd = c(2,2,2), seg.len = 4)
# Verify plot is correct
density.values<-seq(from = 0.8, to = 1.6, by = 0.1)
data.frame(density.values, round(predict(object = mod.fit.nom.density,
newdata = data.frame(density = density.values), type = "probs"), 2))## density.values Healthy Scab Sprout
## 1 0.8 0.00 0.95 0.05
## 2 0.9 0.00 0.89 0.11
## 3 1.0 0.01 0.76 0.24
## 4 1.1 0.05 0.54 0.41
## 5 1.2 0.27 0.25 0.48
## 6 1.3 0.69 0.05 0.25
## 7 1.4 0.92 0.01 0.07
## 8 1.5 0.98 0.00 0.02
## 9 1.6 1.00 0.00 0.00predict.data<-data.frame(class = "hrw", density = c(mean(wheat$density), mean(wheat$density)),
hardness = mean(wheat$hardness), size = mean(wheat$size),
weight = mean(wheat$weight), moisture = mean(wheat$moisture))
pi.hat<-predict(object = mod.fit, newdata = predict.data, type = "probs")
head(pi.hat)## Healthy Scab Sprout
## 1 0.2300126 0.2514559 0.5185315
## 2 0.2300126 0.2514559 0.5185315pi.hat[,1]## 1 2
## 0.2300126 0.2300126# Create plotting area first to make sure get the whole region with respect to x-axis
curve(expr = predict(object = mod.fit, newdata = data.frame(class = "hrw", density = x,
hardness = mean(wheat$hardness), size = mean(wheat$size), weight = mean(wheat$weight),
moisture = mean(wheat$moisture)), type = "probs")[,1], ylab = expression(hat(pi)),
xlab = "Density", xlim = c(min(wheat$density), max(wheat$density)), col = "black",
lty = "solid", lwd = 2, n = 1000, panel.first = grid(col = "gray", lty = "dotted"))
curve(expr = predict(object = mod.fit, newdata = data.frame(class = "hrw", density = x,
hardness = mean(wheat$hardness), size = mean(wheat$size), weight = mean(wheat$weight),
moisture = mean(wheat$moisture)), type = "probs")[,2], ylab = expression(hat(pi)),
xlab = "Density", xlim = c(min(wheat$density), max(wheat$density)), col = "green",
lty = "dotdash", lwd = 2, n = 1000, add = TRUE, panel.first = grid(col = "gray", lty = "dotted"))
curve(expr = predict(object = mod.fit, newdata = data.frame(class = "hrw", density = x,
hardness = mean(wheat$hardness), size = mean(wheat$size), weight = mean(wheat$weight),
moisture = mean(wheat$moisture)), type = "probs")[,3], ylab = expression(hat(pi)),
xlab = "Density", xlim = c(min(wheat$density), max(wheat$density)), col = "red", lty = "longdash",
lwd = 2, n = 1000, add = TRUE, panel.first = grid(col = "gray", lty = "dotted"))
legend(1.4,0.6, legend=c("Healthy", "Sprout", "Scab"),
lty=c("solid","longdash","dotdash"), col=c("black","red","green"), bty="n", lwd = c(2,2,2))
3.4.1 Odds ratios
As razões de probabilidades ou odds ratios baseadas nas probabilidades acumuladas são facilmente formadas porque o modelo de probabilidades proporcionais iguala os logaritmos das probabilidades ao preditor linear. Por exemplo, a razão de chances envolvendo uma variável explicativa, digamos \(x_1\), é \[ \dfrac{Odds_{x_1+c,\cdots,x_p} (Y\leq j)}{Odds_{x_1,\cdots,x_p} (Y\leq j)} = \dfrac{e^{\beta_{j0}+\beta_1(x_1+c)+\cdots +\beta_p x_p}}{e^{\beta_{j0}+\beta_1 x_1+\cdots +\beta_p x_p}} = e^{c\, \beta_1}, \] onde deixamos \(Odds_{x_1,\cdots,x_p} (Y\leq j)\) denotar as chances de observar a categoria \(j\) ou menor para \(Y\) . Observe que esta formulação da razão de chances assume que \(x_2,\cdots,x_p\) permanecem inalterados. A interpretação formal da razão de chances é:
As chances de \(Y\leq j\) vs. \(Y > j\) mudam em \(e^{c\, \beta_1}\) vezes para um aumento de \(c\) unidades em \(x_1\) enquanto as outras variáveis explicativas no modelo são mantidas constantes.
Curiosamente, a razão de chances permanece a mesma, não importa qual categoria de resposta seja usada para \(j\). Esta característica chave do modelo de probabilidades proporcionais ocorre devido à ausência de um subscrito \(j\) nos parâmetros de inclinação \(\beta_1,\cdots,\beta_p\).
Exemplo 3.10: Grãos de trigo.
As razões de chances estimadas para cada variável explicativa são calculadas como \[ \widehat{OR} = \exp\big(c \, \widehat{\beta}_r\big) \] para \(r = 1,\cdots,6\). Similarmente à Seção 3.3, definimos \(c\) igual a um desvio padrão para cada variável explicativa contínua e \(c = 1\) para a variável indicadora srw.
Lembre-se de que também precisamos multiplicar os parâmetros de regressão estimados produzidos por polr() por -1. Abaixo seguem os cálculos:
round (c.value , 2) # class = 1 is first value## density hardness size weight moisture
## 1.00 0.13 27.36 0.49 7.92 2.03round (exp(c.value * (-mod.fit.ord$coefficients )), 2)## density hardness size weight moisture
## 0.84 0.17 0.75 1.15 0.37 1.08round (1/ exp(c.value * (-mod.fit.ord$coefficients )), 2)## density hardness size weight moisture
## 1.19 5.89 1.33 0.87 2.74 0.92Na última linha do código, expressamos novamente as razões de chances em relação às reduções em suas respectivas variáveis explicativas. As interpretações das razões de chances estimadas incluem:
- As probabilidades estimadas do scab (\(Y\leq 1\)) vs. sprout ou resposta healthy (Y > 1) são 0.84 vezes maiores para o trigo de inverno vermelho macio em vez do duro, mantendo as outras variáveis constantes.
- As probabilidades estimadas da scab vs. sprout ou healthy (resposta saudável) mudam em 5.89 vezes para uma diminuição de 0.13 na densidade, mantendo as outras variáveis constantes.
- As probabilidades estimadas da scab vs. sprout ou healthy (resposta saudável) mudam em 2.74 vezes para uma diminuição de 7.92 no peso, mantendo as outras variáveis constantes./li>
Devido à suposição de chances proporcionais, cada uma dessas afirmações também se aplica às chances da scab ou sprout vs. healthy (resposta saudável). Por esta razão, é comum interpretar odds ratio, como para densidade, dizendo:
As probabilidades estimadas da qualidade do kernel estar abaixo de um determinado nível mudam 5.89 vezes para uma diminuição de 0.13 na densidade, mantendo as outras variáveis constantes.
No geral, vemos que quanto maior a densidade e o peso, maior a probabilidade de um kernel ser de uma categoria mais saudável. A função polr() cria objetos que têm um tipo de classe “polr”. A função method() com este tipo de classe mostra quais funções de método estão disponíveis:
class (mod.fit.ord)## [1] "polr"methods ( class = polr )## [1] anova Anova brief confint
## [5] Confint deltaMethod extractAIC linearHypothesis
## [9] logLik model.frame nobs poTest
## [13] predict print profile S
## [17] simulate summary vcov vif
## see '?methods' for accessing help and source codeJá usamos as funções summary.polr() e Anova.polr(). A função de método correspondente para confint() pode ser usada para construir intervalos de confiança LR perfilada para os parâmetros \(\beta_1,\cdots,\beta_p\). Os resultados desta função são usados para encontrar intervalos LR perfilada para as razões de chances:
conf.beta <- confint ( object = mod.fit.ord , level = 0.95)
ci <- exp (c.value*(- conf.beta ))
round ( data.frame (low = ci [,2] , up = ci [ ,1]) , 2)## low up
## classsrw 0.39 1.81
## density 0.11 0.26
## hardness 0.55 1.03
## size 0.77 1.72
## weight 0.23 0.58
## moisture 0.76 1.54round ( data.frame (low = 1/ ci [,1], up = 1/ ci [ ,2]) , 2)## low up
## classsrw 0.55 2.57
## density 3.87 9.36
## hardness 0.97 1.83
## size 0.58 1.29
## weight 1.73 4.40
## moisture 0.65 1.31A interpretação da razão de chances de densidade é:
- Com 95% de confiança, as chances da qualidade do kernel estar abaixo de um determinado nível mudam de 3.87 a 9.36 vezes quando a densidade diminui em 0.13, mantendo as outras variáveis constantes.
Os intervalos de confiança para density (densidade) e weight (peso) são os únicos intervalos que não contêm 1. Portanto, há evidências suficientes de que ambas as variáveis explicativas são importantes.
Embora geralmente prefiramos intervalos LR, os intervalos de confiança de Wald também podem ser calculados. No entanto, a função confint.default() não pode ser aplicada.
3.4.2 Tabelas de contingência
A Seção 3.2 discute como testar a independência em uma tabela de contingência usando um modelo de regressão multinomial para uma resposta nominal. Um teste semelhante para independência também pode ser realizado usando o modelo de chances proporcionais, mas a hipótese alternativa especifica o tipo de dependência que pode estar presente.
Isso ocorre porque o modelo de chances proporcionais assume uma estrutura específica para a associação entre uma variável explicativa categórica \(X\) e uma variável de resposta \(Y\) e, portanto, usa menos parâmetros do que o modelo nominal para resumir essa associação.
Em particular, para uma tabela de contingência \(I\times J\) com categorias ordinais de \(Y\), o modelo de chances proporcionais é \[ logit\big( P(Y\leq j) \big) = \beta_{j0}+\beta_2 x_2 + \cdots + \beta_I x_I, \] \(j=1,\cdots,J-1\) onde \(x_2,\cdots, x_I\) são variáveis indicadoras para as linhas \(2,\cdots,I\), respectivamente.
Qualquer \(\beta_i\neq 0\) para \(i = 2,\cdots,I\) significa que as probabilidades envolvendo probabilidades acumuladas para \(Y\) não são as mesmas nas linhas 1 e \(i\). Categorias inferiores de \(Y\) são mais prováveis de serem observadas na linha \(i\) do que na linha 1 se \(\beta_i > 0\) e menos prováveis se \(\beta_i < 0\).
Assim, a independência é testada como \(H_0 \, : \, \beta_2 = \beta_3 = \cdots = \beta_I = 0\) vs. \(H_1 \, : \, \mbox{qualquer } \beta_i\neq 0\); \(i = 2,\cdots,I\). Se a hipótese nula de independência for rejeitada, a associação é resumida com o auxílio dos sinais e valores de \(\widehat{\beta}_2,\cdots,\widehat{\beta}_I\), juntamente com intervalos de confiança para os parâmetros correspondentes.
Observe que a equação acima contém apenas \((I-1)+(J-1)\) parâmetros, em comparação com \(I(J-1)\) no modelo de regressão multinomial correspondente. A redução nos parâmetros vem da suposição de chances proporcionais, que especifica que a diferença nos logits de probabilidades acumuladas para quaisquer duas linhas é controlada por apenas um parâmetro \(\beta_i\), independentemente de \(j\).
Assim, o teste de independência usando o modelo de chances proporcionais é mais poderoso do que um teste de independência usando o modelo de regressão multinomial quando a associação realmente segue essa estrutura. Caso contrário, o teste usando o modelo de chances proporcionais pode falhar em detectar qualquer associação, mesmo quando for muito forte.
Exemplo 3.11: Biscoitos enriquecidos com fibras.
Ao examinar uma tabela \(4\times 4\) de contingência para os dados do biscoito enriquecido com fibra, descobrimos anteriormente que havia evidência moderada de uma associação entre a gravidade do inchaço e a fonte de fibra. Como a gravidade do inchaço é medida de maneira ordinal (nenhum < baixo < médio < alto), um modelo de chances proporcionais nos permite realizar uma análise potencialmente mais poderosa.
Usando a fibra como variável explicativa com variáveis indicadoras apropriadamente construídas (consulte o exemplo da Seção 3.3), nosso modelo é \[ logit\big( P(Y\leq j)\big) = \beta_{j0}+\beta_2 \, \mbox{bran} +\beta_3 \, \mbox{gum} +\beta_4 \, \mbox{both}, \] onde j corresponde aos níveis 1 (nenhum), 2 (baixo) e 3 (médio) de gravidade do inchaço.
Observe que esse modelo envolve um número menor de parâmetros do que o modelo de regressão multinomial da Seção 3.3. Estimamos o modelo usando a função polr():
library ( package = MASS )
levels ( diet$bloat )## [1] "none" "low" "medium" "high"mod.fit.ord <- polr ( formula = bloat ~ fiber , weights = count , data = diet , method = "logistic")
summary (mod.fit.ord)## Call:
## polr(formula = bloat ~ fiber, data = diet, weights = count, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## fiberbran -0.3859 0.7813 -0.494
## fibergum 2.4426 0.8433 2.896
## fiberboth 1.4235 0.7687 1.852
##
## Intercepts:
## Value Std. Error t value
## none|low 0.0218 0.5522 0.0395
## low|medium 1.6573 0.6138 2.7002
## medium|high 3.0113 0.7249 4.1539
##
## Residual Deviance: 112.2242
## AIC: 124.2242O modelo estimado é \[ logit\big( \widehat{P}(Y\leq j)\big) = \widehat{\beta}_{j0}+0.3859 \, \mbox{bran} -2.4426 \, \mbox{gum} -1.4235 \, \mbox{both}, \]
onde \(\widehat{\beta}_{10}=0.0218\), \(\widehat{\beta}_{20}=1.6573\) e \(\widehat{\beta}_{30}=3.0113\).
A estatística da razão de verossimilhanças (LRT) para o teste de independência \[ H_0 \, : \, \beta_2=\beta_3=\beta_4=0 \qquad \mbox{vs.} \qquad H_1 \ : \, \mbox{qualquer } \beta_i\neq 0, \] \(i=2,3,4\) produz:
library ( package = car)
Anova (mod.fit.ord)## Analysis of Deviance Table (Type II tests)
##
## Response: bloat
## LR Chisq Df Pr(>Chisq)
## fiber 15.048 3 0.001776 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como \(-2\log(\Lambda ) = 15.048\) é grande em relação a uma distribuição \(\chi^2_3\), \(p\)-valor = 0.0018, há fortes evidências de que existe associação na forma de uma tendência entre os log-odds para as probabilidades cumulativas.
Lembre-se de que com o modelo de regressão multinomial da Seção 3.3 e os testes de independência fornecidos na Seção 3.2, houve apenas evidência moderada de associação. Ao restringir o tipo de associação considerado por nosso teste, agora temos mais poder para rejeitar a hipótese nula.
As razões de chances estimadas para comparar cada fonte de fibra com o uso de nenhuma fibra são \(\widehat{OR} = \exp\big(\widehat{\beta}_r\big)\) para \(r = 2,3\) e 4. Abaixo estão essas razões de chances e os intervalos LR de perfil correspondentes:
round (exp(- coefficients (mod.fit.ord )), 2)## fiberbran fibergum fiberboth
## 1.47 0.09 0.24conf.beta <- confint ( object = mod.fit.ord , level = 0.95)
ci <- exp (- conf.beta )
round ( data.frame (low = ci [,2] , up = ci [ ,1]) , 2)## low up
## fiberbran 0.32 7.06
## fibergum 0.02 0.43
## fiberboth 0.05 1.05round ( data.frame (low = 1/ ci [,1], up = 1/ ci [ ,2]) , 2)## low up
## fiberbran 0.14 3.15
## fibergum 2.32 65.01
## fiberboth 0.95 19.82Observe que inverter os intervalos de confiança no último conjunto de código expressa o efeito de cada fonte de fibra em termos de chances de inchaço mais grave. Como o intervalo LR perfilada para <stronggum não contém 1, há evidências suficientes para indicar que o uso de gum (chiclete) como fonte de fibra aumenta a gravidade do inchaço; no entanto, não há evidências suficientes de que o bran (farelo) aumente a gravidade do inchaço porque seu intervalo de confiança contém 1. Ao combinar ambas as fontes de fibra, há evidência marginal de um efeito, porque o limite inferior do intervalo é próximo a 1. Formalmente, podemos escrever uma interpretação para os intervalos de confiança da fonte de chiclete vs. sem fibra como:
Com 95% de confiança, as chances de a gravidade do inchaço estar acima de um determinado nível são entre 2.32 e 65.01 vezes maiores ao usar chiclete como fonte de fibra do que ao não usar fibra.
Para comparar o chiclete diretamente com o farelo, a maneira mais simples é usar a função relevel() para alterar o nível de base da variável de origem da fibra. Um código semelhante ao apresentado acima pode ser usado para reestimar o modelo e calcular os intervalos LR perfilada.
Descobrimos que, com 95% de confiança, as chances de a gravidade do inchaço estar abaixo de um determinado nível são entre 0.01 e 0.30 vezes maior ao usar chiclete como fonte de fibra do que ao usar farelo.
Assim, isso indica que o chilcte causa mais inchaço do que o farelo. No geral, se for dada uma escolha entre as duas fontes de fibra e sua combinação, o farelo seria preferido se o único critério fosse a gravidade do inchaço.
3.4.3 Modelo de probabilidades não proporcionais
Um modelo de probabilidades proporcionais é uma das formas preferidas de explicar uma resposta multinomial ordenada, isso porque os parâmetros de regressão de inclinação são constantes nas categorias de resposta.
Embora isso possa simplificar bastante o modelo, impõe a suposição de que a associação afeta os logaritmos das probabilidades acumuladas da mesma forma para todo \(j = 1,\cdots, J - 1\). Isso pode não ser verdade em todas as situações. Um modelo alternativo que relaxa essa suposição é \[ logit\big(P(Y\leq j)\big) = \beta_{j0}+\beta_{j1} x_1 + \cdots +\beta_{jp} x_p, \] onde \(j = 1,\cdots, J- 1\). Observe que todos os parâmetros de regressão agora podem variar entre os níveis de \(Y\). A equação acima é chamada de modelo de probabilidade não proporcional.
Como o modelo de chances proporcionais é um caso especial da equação acima, podemos testar a suposição de chances proporcionais por meio das hipóteses \[ H_0 \, : \, \beta_{1r} = \cdots = \beta_{J-1,r} \quad \mbox{para} \; r = 1,\cdots,p \quad \mbox{vs.} \quad H_1 \, : \, \mbox{Nem todos são iguais}\cdot \]
O teste é conduzido como um LRT onde os graus de liberdade para a distribuição \(\chi^2\) de referência são encontrados a partir da diferença no número de parâmetros nos dois modelos, \[ (p + 1)(J - 1) − (p + J - 1) = p(J − 2)\cdot \]
Rejeitar a suposição de probabilidades proporcionais sugere que o modelo de probabilidades não proporcionais pode ser preferido. As probabilidades estimadas e as razões de chances têm uma forma diferente devido aos parâmetros extras. No entanto, o modelo de probabilidades proporcionais ainda pode ser preferido devido ao seu menor número de parâmetros.
Cada parâmetro estimado adiciona variabilidade à estimativa subsequente de probabilidades e razões de chances, de modo que se pode obter estimativas de razões de chances mais próximas da verdade usando um modelo menor com um defeito menor do que usando um modelo muito maior sem o defeito.
Além disso, um tamanho de amostra muito grande pode resultar na rejeição da hipótese nula, mesmo que os dados se desviem apenas ligeiramente da suposição de chances proporcionais. Deixar de rejeitar a hipótese das probabilidades proporcionais não é prova de que ela seja verdadeira. No entanto, oferece alguma garantia de que um modelo de probabilidades proporcionais fornece uma aproximação razoável para relações verdadeiras entre \(Y\) e as variáveis explanatórias. Examinaremos maneiras adicionais de avaliar o ajuste de um modelo no Capítulo 5.
Exemplo 3.12: Biscoitos enriquecidos com fibras.
Infelizmente, a função polr() não fornece uma maneira de testar a suposição de chances proporcionais. Em vez disso, usamos o vglm() do pacote VGAM para ajudar a realizar o teste. Primeiro, estimamos o modelo de probabilidades proporcionais e examinamos um resumo dos resultados:
library ( package = VGAM )
mod.fit.po <- vglm ( formula = bloat ~ fiber , family = cumulative ( parallel = TRUE ) ,
weights = count , data = diet [ diet$count != 0 ,])
summary ( mod.fit.po )##
## Call:
## vglm(formula = bloat ~ fiber, family = cumulative(parallel = TRUE),
## data = diet[diet$count != 0, ], weights = count)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.02182 0.55654 0.039 0.96872
## (Intercept):2 1.65734 0.62533 2.650 0.00804 **
## (Intercept):3 3.01128 0.72766 4.138 3.5e-05 ***
## fiberbran 0.38593 0.79412 0.486 0.62697
## fibergum -2.44258 0.82591 -2.957 0.00310 **
## fiberboth -1.42348 0.78090 -1.823 0.06832 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]),
## logitlink(P[Y<=3])
##
## Residual deviance: 112.2242 on 36 degrees of freedom
##
## Log-likelihood: -56.1121 on 36 degrees of freedom
##
## Number of Fisher scoring iterations: 6
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## fiberbran fibergum fiberboth
## 1.47098886 0.08693587 0.24087418O modelo estimado é praticamente o mesmo encontrado anteriormente com polr(). Dentro da chamada para vglm(), o argumento family = cumulative especifica um modelo logit cumulativo, onde o valor do argumento parallel = TRUE denota a versão das probabilidades proporcionais. Além disso, o argumento dos pesos é usado para denotar o número de indivíduos que cada linha do quadro de dados representa. Se o conjunto de dados tivesse sido formulado alternativamente de modo que cada linha representasse uma observação individual, ou seja, uma pessoa testando um biscoito, o argumento dos pesos não seria necessário. Observe que uma peculiaridade de vglm() requer a remoção de todas as contagens 0 do conjunto ou quadro de dados antes de usar a função.
O pacote VGAM usa uma estrutura de classe diferente da que vimos até agora. Essa estrutura de classe, geralmente referida como S4, produz objetos com componentes em slots. R e seus pacotes correspondentes são escritos de forma muito semelhante à linguagem de programação S. Essa linguagem foi desenvolvida pela primeira vez na década de 1970 nos Laboratórios Bell. As versões 3 e 4 do S são emuladas pelo R. A versão 3 (S3) é predominante e é ela que é usada principalmente neste texto.
Por exemplo, para acessar uma lista de todos os slots dentro do objeto mod.fit.po, usamos a função slotNames():
slotNames ( mod.fit.po )## [1] "extra" "family" "iter" "predictors"
## [5] "assign" "callXm2" "contrasts" "df.residual"
## [9] "df.total" "dispersion" "effects" "offset"
## [13] "qr" "R" "rank" "ResSS"
## [17] "smart.prediction" "terms" "Xm2" "Ym2"
## [21] "xlevels" "call" "coefficients" "constraints"
## [25] "control" "criterion" "fitted.values" "misc"
## [29] "model" "na.action" "post" "preplot"
## [33] "prior.weights" "residuals" "weights" "x"
## [37] "y"mod.fit.po@coefficients## (Intercept):1 (Intercept):2 (Intercept):3 fiberbran fibergum
## 0.02182272 1.65733908 3.01127560 0.38593487 -2.44258450
## fiberboth
## -1.42348056Para acessar slots específicos, o símbolo @ é usado na forma
mod.fit.npo <- vglm ( formula = bloat ~ fiber , family = cumulative ( parallel = FALSE ) ,
weights = count , data = diet [ diet$count != 0 ,])
summary ( mod.fit.npo )##
## Call:
## vglm(formula = bloat ~ fiber, family = cumulative(parallel = FALSE),
## data = diet[diet$count != 0, ], weights = count)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.680e-15 5.774e-01 0.000 1.0000
## (Intercept):2 1.609e+00 7.746e-01 2.078 0.0377 *
## (Intercept):3 1.802e+01 1.433e+03 NA NA
## fiberbran:1 3.365e-01 8.223e-01 0.409 0.6824
## fiberbran:2 7.885e-01 1.300e+00 0.606 0.5443
## fiberbran:3 4.083e-09 2.027e+03 0.000 1.0000
## fibergum:1 -1.609e+00 9.661e-01 -1.666 0.0957 .
## fibergum:2 -2.303e+00 9.874e-01 -2.332 0.0197 *
## fibergum:3 -1.768e+01 1.433e+03 -0.012 0.9902
## fiberboth:1 -1.609e+00 9.661e-01 -1.666 0.0957 .
## fiberboth:2 -1.273e+00 9.710e-01 -1.311 0.1899
## fiberboth:3 -1.641e+01 1.433e+03 -0.011 0.9909
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]),
## logitlink(P[Y<=3])
##
## Residual deviance: 108.3914 on 30 degrees of freedom
##
## Log-likelihood: -54.1957 on 30 degrees of freedom
##
## Number of Fisher scoring iterations: 16
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'
##
##
## Exponentiated coefficients:
## fiberbran:1 fiberbran:2 fiberbran:3 fibergum:1 fibergum:2 fibergum:3
## 1.400000e+00 2.200000e+00 1.000000e+00 2.000000e-01 1.000000e-01 2.089257e-08
## fiberboth:1 fiberboth:2 fiberboth:3
## 2.000000e-01 2.800000e-01 7.461632e-08round ( mod.fit.npo@coefficients , 2)## (Intercept):1 (Intercept):2 (Intercept):3 fiberbran:1 fiberbran:2
## 0.00 1.61 18.02 0.34 0.79
## fiberbran:3 fibergum:1 fibergum:2 fibergum:3 fiberboth:1
## 0.00 -1.61 -2.30 -17.68 -1.61
## fiberboth:2 fiberboth:3
## -1.27 -16.41Observe que agora há \(p(J − 2) = 6\) parâmetros adicionais em comparação com o modelo de chances proporcionais. As hipóteses para o teste da suposição de odds proporcionais são \(H_0 \, : \, \beta_{1i} = \beta_{2i} = \beta_{3i}\) para \(i\) = 2, 3, 4 vs. \(H_1 \, : \, \mbox{Nem todos iguais}\).
Nem Anova() nem anova() foram implementadas no VGAM, então calculamos \(-2\log(\Lambda)\) encontrando a diferença entre os dois desvios residuais do modelo, consulte a Seção 2.2.2 para um exemplo envolvendo regressão logística:
tran.LR <- deviance ( mod.fit.po ) - deviance ( mod.fit.npo )
df <- mod.fit.po@df.residual - mod.fit.npo@df.residual
p.value <- 1 - pchisq ( q = tran.LR , df = df )
data.frame ( tran.LR , df , p.value )## tran.LR df p.value
## 1 3.832754 6 0.6992972Usamos a função deviance() para extrair os deviances residuais, porque não há slot @deviance disponível. A estatística do teste da razão de verossimilhanças (LRT) resulta em uma estatística de \(-2\log(\Lambda) = 3.83\) e um \(p\)-valor de 0.6993. Portanto, não há evidências suficientes para indicar que a suposição de chances proporcionais foi violada.
Um problema com o uso de modelos de chances não proporcionais para uso geral é que o modelo não restringe adequadamente os parâmetros para evitar \(P(Y\leq j) < P(Y\leq j' )\) para \(j > j'\). Assim, as probabilidades cumulativas podem diminuir em algum ponto fazendo com que a probabilidade de uma categoria individual seja menor que 0.
Essa violação das regras de probabilidade ocorre porque o efeito que uma variável explicativa tem sobre \(logit\big(P (Y\leq j)\big)\) pode mudar para cada \(j\); isto é, os parâmetros \(\beta_{j1},\cdots,\beta_{jp}\) podem variar livremente sobre os níveis de \(Y\). Por esse motivo, é necessário ter cuidado com esses modelos para garantir que probabilidades sem sentido não ocorram. O próximo exemplo ilustra um caso em que um modelo de probabilidades não proporcional é inadequado.
Exemplo 3.13:
Gráficos de modelo de probabilidades não
proporcionais.
Considere o modelo \[ logit\big(P (Y\leq j)\big) = \beta_{j0} + \beta_{j1} x_1, \] onde \(\beta_{10} = 0\), \(\beta_{20} = 2\), \(\beta_{30} = 4\), \(\beta_{11} = 2\), \(\beta_{21} = 3\), \(\beta_{31} = 6\) e \(J = 4\).
A figura abaixo mostra os gráficos das probabilidades cumulativas e individuais das categorias.
# Plot of non-proportional odds model
beta<-c(0, 2, 4, 2, 3, 6) #beta10, beta20, beta30, beta11, beta21, beta31
#beta<-c(0, 2, 4, 2, 3, 4) #beta10, beta20, beta30, beta11, beta21, beta31 #Lines do not appear to cross
x.range<-c(-5, 3)
par(mfrow = c(1, 2))
curve(expr = plogis(q = beta[1] + beta[4]*x), xlim = x.range, ylab = expression(P(Y<=j)),
xlab = expression(x[1]), main = "Cumulative probabilities for Y", lwd = 2)
curve(expr = plogis(q = beta[2] + beta[5]*x), add = TRUE, lty = "dashed", col = "red", lwd = 2)
curve(expr = plogis(q = beta[3] + beta[6]*x), add = TRUE, lty = "dotted", col = "blue", lwd = 2)
legend(x = -5, y = 0.9, legend = c(expression(P(Y<=1)), expression(P(Y<=2)), expression(P(Y<=3))),
lty = c("solid", "dashed", "dotted", "dotdash"), col = c("black", "red", "blue", "green"),
bty = "n", lwd = 2)
curve(expr = plogis(q = beta[1] + beta[4]*x), xlim = x.range, ylab = expression(pi[j]),
xlab = expression(x[1]), main = "Probabilities for Y", lwd = 2)
curve(expr = plogis(q = beta[2] + beta[5]*x) - plogis(q = beta[1] + beta[4]*x), add = TRUE,
lty = "dashed", col = "red", lwd = 2)
curve(expr = plogis(q = beta[3] + beta[6]*x) - plogis(q = beta[2] + beta[5]*x), add = TRUE,
lty = "dotted", col = "blue", lwd = 2)
curve(expr = 1 - plogis(q = beta[3] + beta[6]*x), add = TRUE,
lty = "dotdash", col = "green", lwd = 2)
legend(x = -5, y = 0.9, legend = c(expression(pi[1]), expression(pi[2]), expression(pi[3]),
expression(pi[4])),
lty = c("solid", "dashed", "dotted", "dotdash"), col = c("black", "red", "blue", "green"),
bty = "n", lwd = 2)
Observe que a curva de probabilidade cumulativa para \(P (Y\leq 3)\) cruza as outras duas curvas. Assim, o modelo permite \(P (Y\leq 3) < P (Y\leq 2)\) e \(P (Y\leq 3) < P (Y\leq 1)\) para alguns valores de \(x_1\) . Como resultado, \(\pi_3\) é menor que 0 para alguns valores de \(x_1\).
3.5 Modelos de regressão adicionais
As seções 3.3 e 3.4 apresentam as duas formas mais comuns de modelar respostas multinomiais. Existem muitas outras abordagens. Por exemplo, não estamos limitados a usar apenas uma transformação logit de uma probabilidade (odds). As transformações probit e log-log complementares mostradas na Seção 2.3 podem ser usadas em vez disso. A função polr() estima esses modelos usando o valor de argumento method apropriado na chamada de função.
Outro modelo para respostas multinomiais é o modelo de categorias adjacentes: \[ \log(\pi_{j+1}/\pi_j) = \beta_{j0} +\beta_{j1} x_1 + \cdots + \beta_{jp} x_p, \] onde \(j = 1,\cdots, J-1\). Este modelo é semelhante ao modelo de regressão multinomial, mas também permite tirar vantagem de uma resposta ordinal devido às sucessivas comparações das categorias \(j + 1\) e \(j\). Para aproveitar ainda mais as propriedades ordinais em uma resposta, podemos usar uma versão simplificada do modelo: \[ \log(\pi_{j+1}/\pi_j) = \beta_{j0} + \beta_1 x_1 +\cdots +\beta_p x_p, \] onde os parâmetros de inclinação da regressão são forçados a ser os mesmos para todos os \(j\).
Este modelo assume que uma variável explicativa tem o mesmo efeito sobre as probabilidades para cada comparação de categoria sucessiva da mesma forma que o modelo de regressão de probabilidades proporcionais faz com probabilidades baseadas em probabilidades acumuladas. Ambos os modelos de categorias adjacentes podem ser estimados em R usando a função vglm(). Exploramos mais esses e outros modelos na Seção 3.6 Exercícios.
3.6 Exercícios
Discuta por que a distribuição marginal de \(n_{1+}\) é binomial com o parâmetro \(\pi_{1+}\).
A distribuição conjunta condicional de \((n_{11}, n_{12}, n_{21})\) dado \(n_{1+}\) pode ser mostrado como \[ f(n_{11}, n_{12},n_{21} \, | \, n_{1+}) = f(n_{11}, n_{12}, n_{21})/f(n_{1+}), \] onde \(f(n_{11}, n_{12}, n_{21})\) é a distribuição conjunta de \((n_{11}, n_{12}, n_{21})\) e \(f(n_{1+})\) é a distribuição marginal de \(n_{1+}\). Observe que não precisamos declarar \(n_{22}\) dentro de \(f(n_{11}, n_{12}, n_{21})\) porque \(n = n_{11}+n_{12}+n_{21}+n_{22}\) e \(n\) é conhecido. Da mesma forma, não precisamos declarar \(n_{2+}\) dentro de \(f(n_{1+})\).
Em geral para estruturas de tabelas de contingência \(I\times J\), estenda o resultado em (b).
Para um aumento de \(c\) unidades em \(x\) em um modelo de regressão multinomial, deduza a razão de chances que compara uma resposta de categoria \(j\) (\(j\neq 1\)) a uma resposta de categoria 1. Mostre a forma da variância que seria usada em um intervalo de confiança de Wald.
Para um aumento de \(c\) unidades em \(x\) em um modelo de regressão de probabilidades proporcionais, deduza a razão de chances que compara as probabilidades cumulativas. Mostre a forma da variância que seria usada em um intervalo de confiança de Wald.
Para um aumento de \(c\) unidades em \(x_1\), deduza a razão de chances correspondente que compara uma resposta de categoria \(j\) a uma resposta de categoria 1. Mostre a forma da variância que seria usada em um intervalo de confiança de Wald.
Repita este problema para um modelo de regressão de probabilidades proporcionais.
Mostre que \[ Odds_{x_1+c} (Y\leq j)/ Odds_{x_1} (Y\leq j) = \exp\big(c \, \beta_{j1}\big)\cdot \] Forneça uma interpretação formal da razão de chances.
Encontre \(\pi_j\) para \(j = 1,\cdots, J\). Por que um \(\pi_j\) pode ser menor que 0 com este modelo?
O que significaria independência e dependência entre as linhas e colunas no contexto de DAVFs? Qual conclusão você acha que os pesquisadores prefeririam para justificar seu sistema de classificação? Explique.
Cognard et ai. (1995) afirmam que “a tabela de contingência mostra uma diferença estatisticamente significativa (P = 0,0001)” sem mencionar o teste estatístico específico realizado. Realize um teste qui-quadrado de Pearson para independência para produzir um resultado de significância semelhante.
É improvável que pacientes com DAVF apresentem apenas um sintoma. Embora não declarado neste artigo, os autores podem ter resumido o sintoma “primário” ou “mais extremo” do paciente, de modo que cada paciente contribuiu para a contagem de apenas uma célula na tabela. No entanto, se este não for o caso, um paciente pode contribuir para nenhuma, uma ou mais de uma célula. Por exemplo, um paciente pode ter hipertensão intracraniana e convulsões enquanto possui uma classificação de nível 3. De fato, as contagens totais da tabela podem ser maiores ou menores que o número de pacientes. Esse cenário alternativo invalidaria o uso do teste qui-quadrado de Pearson nessa situação? Explique
Estime o modelo de regressão de chances proporcionais envolvendo farelo, goma e sua interação. Há evidências suficientes para indicar que uma interação é importante?
Avalie a suposição de probabilidades proporcionais estimando o modelo de probabilidades não proporcional correspondente e, em seguida, realizando um LRT.
Com o modelo mais apropriado encontrado através das investigações em (a) e (b), use as razões de chances para interpretar o efeito do farelo e da goma no inchaço.
Interprete os odds ratio estimados e os correspondentes intervalos de confiança para as outras variáveis explicativas não especificamente abordadas nos exemplos das Seções 3.3.1 e 3.4.1.
Calcule as estimativas de odds ratio comparando sarna vs. broto correspondente às variáveis referenciadas em (a). Calcule intervalos de confiança e interprete esses intervalos no contexto dos dados.
1- Para o caso de \(I = J = 2\), mostre que \[\sum_{i=1}^I \sum_{j=1}^J \dfrac{\big(n_{ij}-n_{i+}n_{+j}/n \big)^2}{n_{i+}n_{j+}/n} \] é equivalente à estatística \(X^2\) como indicada no Captulo 1.
2- Mostre que \(-2\log(\Lambda )\) simplifica para \[-2\sum_{i=1}Î \sum_{j=1}^J n_{ij}\log\left(\dfrac{n_{ij}}{n_{i+}n_{+j}/n}\right)\] ao testar a independência entre duas variáveis.
3- Considere uma estrutura de tabela de contingência \(2\times 2\).
4- Suponha que temos uma variável explicativa contínua \(x\) que deve ser representada em um modelo de regressão por um termo linear e um quadrático.
5- Suponha que um modelo de regressão multinomial tenha duas variáveis explicativas contínuas \(x_1\) e \(x_2\) e elas são representadas no modelo por seus termos lineares e de interação.
6- Considere o modelo de regressão de probabilidades não proporcionais \[ logit\big(P(Y\leq j)\big) = \beta_{j0} + \beta_{j1} x_1, \] onde \(j = 1,\cdots, J - 1\).
7- O Exercício 19 do Capítulo 2 examinou a prevalência da hepatite C entre os profissionais de saúde na área de Glasgow, na Escócia. Convirta os dados deste exercício para uma estrutura de tabela de contingência onde os grupos ocupacionais estão localizados nas linhas e a presença e ausência de hepatite estão localizadas nas colunas. Realize os testes qui-quadrado e LR de Pearson para independência usando esses dados. Seus resultados aqui são semelhantes aos encontrados para a análise de dados no Capítulo 2? Explique.
8- Os sistemas de classificação são frequentemente usados por médicos para ajudar a entender os sintomas dos pacientes e orientá-los na aplicação de tratamentos apropriados. Para isso, Cognard et al. (1995) propõem um novo sistema de classificação em relação às fístulas arteriovenosas durais (DAVFs). DAVFs são conexões anormais entre artérias e veias na cabeça de uma pessoa. Se não forem tratadas, podem levar à morte devido a derrames ou outras complicações. A tabela abaixo fornece os dados de 205 casos incluídos em seu estudo. As colunas da tabela fornecem a classificação de cada indivíduo, onde as classificações geralmente aumentam na gravidade da DAVFs da esquerda para a direita. Cada indivíduo recebe apenas uma classificação. As linhas representam diferentes sintomas experimentados pelos indivíduos. Suponha que cada indivíduo também tenha apenas um sintoma.
\[ \mbox{Classificação das DAVFs por sintomas. A fonte de dados é Cognard et al. (1995).} \\ \begin{array}{ccccc} \hline \mbox{} & & & & & \mbox{Classificação} & & & \\ \mbox{} & & 1 & 2a & 2b & 2a \; e \;2b & 3 & 4 & 5 \\\hline \mbox{} & \mbox{Hemorragia} & 0 & 0 & 2 & 1 & 10 & 19 & 5 \\ \mbox{} & \mbox{Hipertensão intracraniana} & 1 & 8 & 1 & 2 & 0 & 4 & 1 \\ \mbox{} & \mbox{Déficit neurológico focal} & 0 & 0 & 0 & 6 & 8 & 2 & 0 \\ \mbox{} \mbox{Sintomas} & \mbox{Convulsões} & 0 & 1 & 0 & 2 & 1 & 3 & 0 \\ \mbox{} & \mbox{Deficiência cardíaca} & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ \mbox{} & \mbox{Mielopatia} & 0 & 0 & 0 & 0 & 0 & 0 & 6 \\ \mbox{} & \mbox{Sintomas não agressivos} & 83 & 17 & 7 & 6 & 6 & 1 & 0 \\\hline \end{array} \]
9- A Seção 3.3.2 discute como o exemplo de crackers enriquecidos com fibras pode ser analisado separando a variável fonte de fibra em duas variáveis binárias separadas para farelo e goma. Esse tipo de análise pode ser preferido para avaliar separadamente os efeitos dessas duas fontes sobre o inchaço. Usando este formulário de dados, preencha o seguinte:
10- Para os dados do grão de trigo, preencha o seguinte:
11- Para maximizar as vendas, os itens dentro dos supermercados são colocados estrategicamente para chamar a atenção do cliente. Este exercício examina um tipo de item - cereal matinal. Normalmente, em grandes mercearias, as caixas de cereais são colocadas em conjuntos de prateleiras localizadas de um lado do corredor. Ao colocar caixas específicas de cereais em prateleiras específicas, os supermercados podem atrair mais clientes para eles. Para investigar isso ainda mais, uma amostra aleatória de tamanho 10 foi retirada de cada uma das quatro prateleiras de uma mercearia Dillons em Manhattan, KS. Esses dados são fornecidos no arquivo:
cereal = read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/cereal_dillons.csv")
head(cereal)## ID Shelf Cereal size_g sugar_g fat_g sodium_mg
## 1 1 1 Kellog's Razzle Dazzle Rice Crispies 28 10 0 170
## 2 2 1 Post Toasties Corn Flakes 28 2 0 270
## 3 3 1 Kellogg's Corn Flakes 28 2 0 300
## 4 4 1 Food Club Toasted Oats 32 2 2 280
## 5 5 1 Frosted Cheerios 30 13 1 210
## 6 6 1 Food Club Frosted Flakes 31 11 0 180A variável resposta é o número da prateleira, que é numerada de baixo (1) para cima (4), e as variáveis explicativas são o teor de açúcar, gordura e sódio dos cereais. Usando esses dados, conclua o seguinte:
- As variáveis explicativas precisam ser reformatadas antes de continuar. Primeiro, divida cada variável explicativa por seu tamanho de porção para contabilizar os diferentes tamanhos de porção entre os cereais. Em segundo lugar, redimensione cada variável para estar entre 0 e 1. O autor de multinom() sugere que isso pode ajudar na convergência de estimativas de parâmetros ao usar a função. Recomendamos o redimensionamento de 0-1 porque a convergência lenta ocorre aqui ao estimar o modelo na parte (d) sem redimensionar. Observe que a função vglm() também pode ser usada para estimar o modelo. Quando usamos esta função, a convergência é obtida no mesmo número de iterações com e sem o reescalonamento 0-1, então vglm() pode ser uma alternativa de ajuste de modelo mais conveniente. Abaixo está o código que usamos para reformatar os dados após o arquivo de dados é lido em um objeto chamado cereal:
stand01 <- function (x) { (x - min(x))/( max(x) - min(x)) }
cereal2 <- data.frame (Shelf = cereal$Shelf, sugar = stand01 (x = cereal$sugar_g / cereal$size_g ),
fat = stand01 (x = cereal$fat_g / cereal$size_g ),
sodium = stand01 (x = cereal$sodium_mg / cereal$size_g ))- Construa gráficos de box-plots (caixas) lado a lado com gráficos de pontos sobrepostos para cada uma das variáveis explicativas. Abaixo está o código que pode ser usado para parcelas envolvendo açúcar (sugar):
boxplot ( formula = sugar ~ Shelf , data = cereal2 , ylab = "Sugar", xlab = "Prateleira",
pars = list ( outpch = NA))
stripchart (x = cereal2$sugar ~ cereal2$Shelf , lwd = 2, col = "red",
method = "jitter", vertical = TRUE , pch = 1, add = TRUE)
Além disso, construa um gráfico de coordenadas paralelas para as variáveis explicativas e o número da prateleira. Discuta se existem possíveis diferenças de conteúdo entre as prateleiras.
A resposta tem valores de 1, 2, 3 e 4. Em que cenário seria desejável levar em conta a ordinalidade. Você acha que isso ocorre aqui?
Estime um modelo de regressão multinomial com formas lineares das variáveis açúcar, gordura e sódio. Realize testes da razão de verossimilhanças (LRTs) para examinar a importância de cada variável explicativa.
Mostre que não há interações significativas entre as variáveis explicativas (incluindo uma interação entre as três variáveis).
Kellogg’s Apple Jacks (http://www.applejacks.com) é um cereal comercializado para crianças. Para uma porção de 28 gramas, seu teor de açúcar é de 12 gramas, o teor de gordura é de 0.5 gramas e o teor de sódio é de 130 miligramas. Estime as probabilidades de prateleira para Apple Jacks.
Construa um gráfico onde a probabilidade estimada de uma prateleira está no eixo \(y\) e o teor de açúcar está no eixo \(x\). Use o conteúdo médio geral de gordura e sódio como os valores variáveis correspondentes no modelo. Interprete o gráfico em relação ao teor de açúcar.
Estime as razões de chances e calcule os intervalos de confiança correspondentes para cada variável explicativa. Relacione suas interpretações com os gráficos construídos para este exercício.
12- Semelhante ao Exercício 11, o arquivo
cereal_supersaver = read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/cereal_supersaver.csv")contém o teor de açúcar, gordura e sódio dos cereais colocados nas prateleiras de uma mercearia Super Saver em Lincoln, NE. Analise esses dados para entender melhor a colocação dos cereais nas prateleiras dessa mercearia.
13- Um exemplo da Seção 4.2.5 examina os dados da Pesquisa Social Geral dos EUA de 1991, que classifica as pessoas de acordo com: ideol, Ideologia política: Muito liberal (VL), Ligeiramente liberal (SL), Moderado (M), Ligeiramente conservador (SC) e Muito conservador (VC). party, Partido político: Democrata (D) ou Republicano (R) e gender, o sexo: Feminino (F) ou Masculino (M).
Considere a ideologia política como uma variável de resposta ordinal e o partido político e o gênero como variáveis explicativas. Os dados estão disponíveis no arquivo:
set1 = read.csv(file = "http://leg.ufpr.br/~lucambio/CE073/20222S/pol_ideol_data.csv")Utiliza a função factor() com a variável de ideologia para garantir que R coloque os níveis da variável de ideologia na ordem correta.
Como as duas variáveis explicativas são categóricas, podemos visualizar todo o conjunto de dados em uma estrutura de tabela de contingência tridimensional. As funções xtabs() e ftable() são úteis para este propósito. Abaixo está como essas funções podem ser usadas após o arquivo de dados ter sido lido no R, o arquivo de dados set1 contém os dados originais:
c.table <- xtabs ( formula = count ~ party + ideol + gender , data = set1 )
ftable (x = c.table , row.vars = c("gender", "party") , col.vars = "ideol")## ideol M SC SL VC VL
## gender party
## F D 118 23 47 32 44
## R 86 39 28 48 18
## M D 53 18 34 23 36
## R 62 45 18 51 12- Usando modelos de regressão de probabilidades multinomiais e proporcionais que incluem partido, gênero e sua interação, conclua o seguinte:
Estimar os modelos e realizar LRTs para testar a importância de cada variável explicativa.
Calcule as probabilidades estimadas para cada nível de ideologia, dadas todas as combinações possíveis dos níveis de partido e gênero.
Construa uma tabela de contingência com contagens estimadas do modelo. Essas contagens estimadas são encontradas tomando a probabilidade estimada para cada nível de ideologia multiplicada por seu número correspondente de observações para uma combinação de partido e gênero. Por exemplo, existem 264 observações para gender = “F” e party = “D”. Como o modelo de regressão multinomial resulta em \(\widehat{\pi}_{VL} = 0.1667\), a contagem estimada desse modelo é \(0.1667\times 264 = 44\).
As contagens estimadas são iguais às observadas? Explique.
Use as razões de chances calculadas a partir dos modelos estimados para ajudar a entender as relações entre as variáveis explicativas e a resposta.
- Compare os resultados para os dois modelos. Discuta qual modelo pode ser mais apropriado para esse cenário.
14- Continuando o Exercício 13, considere novamente as contagens encontradas usando o modelo de chances proporcionais estimadas. Encontre a combinação correta de contagens desta tabela que resulte em probabilidades estimadas iguais às probabilidades estimadas encontradas diretamente do modelo estimado \(\exp(\widehat{\beta}_{j0} + \widehat{\beta}_{1} x_1 + \widehat{\beta}_2 x_2 + \widehat{\beta}_3 x_1 x_2)\), onde \(x_1\) e \(x_2\) são variáveis indicadoras apropriadas representando party e gender, respectivamente.
15- Pesquisadores da Penn State University realizaram um estudo para determinar o teor de gordura ideal para sorvete. Detalhes do estudo e análise de dados correspondente estão disponíveis em https://onlinecourses.science.psu.edu/stat504/node/187.
Em resumo, 493 indivíduos foram convidados a provar e, em seguida, avaliar um determinado tipo de sorvete em uma escala de 9 pontos: 1 a 9, sendo 1 igual a não gostar e 9 igual a gostar muito. O sorvete oferecido aos indivíduos tinha um nível de proporção de gordura de 0, 0.04, … ou 0.28. Os dados estão disponíveis em
ice_cream <- read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/ice_cream.csv")
head(ice_cream)## fat rating count
## 1 0 1 4
## 2 0 2 17
## 3 0 3 8
## 4 0 4 16
## 5 0 5 5
## 6 0 6 6onde a coluna de contagem representa o número de indivíduos em um determinado grupo de gordura dando uma determinada classificação. Usando esses dados, conclua o seguinte:
Exiba os dados em uma estrutura de tabela de contingência 8\(\times\) 9.
Encontre as proporções de classificação observadas condicionadas ao teor de gordura. Construa um gráfico dessas proporções versus conteúdo de gordura. Recomendamos conectar as proporções observadas por linhas correspondentes à pontuação da avaliação.
Discuta por que o qui-quadrado de Pearson e os testes LR de independência da Seção 3.2 não são escolhas de análise ideais para esses dados.
Estime um modelo de chances proporcionais usando o conteúdo de gordura para prever a classificação. Investigue se um termo quadrático é ou não útil. Observe que a variável de classificação precisa ser convertida em um tipo de classe factor para usar polr().
Realize um teste da razão de verossimilhanças (LRT) para avaliar a suposição de chances proporcionais.
Use razões de chances e gráficos das probabilidades estimadas para interpretar o modelo resultante. Recomende um teor de gordura.
16- Para os dados de crackers enriquecidos com fibras, um modelo de categoria adjacente pode ser escrito como \[ \log(\pi_{j+1}/\pi_{j}) = \beta_{j0} + \beta_{j1} \, \mbox{bran} + \beta_{j2} \, \mbox{gum}+ \beta_{j3} \, \mbox{both}, \] para \(j=1,2\) e 3. Esse modelo pode ser estimado pela função vglm() do pacote VGAM de maneira semelhante ao modelo de probabilidades não proporcionais discutido na Seção 3.4.3.
- Estime o modelo de categoria adjacente usando o seguinte como uma continuação do código usado no exemplo:
library(VGAM)
mod.fit.full <- vglm (formula = bloat ~ fiber, weights = count, acat (parallel = FALSE),
data = diet [diet$count != 0,])
summary (mod.fit.full )##
## Call:
## vglm(formula = bloat ~ fiber, family = acat(parallel = FALSE),
## data = diet[diet$count != 0, ], weights = count)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -0.4055 0.6455 -0.628 0.530
## (Intercept):2 -0.6931 0.8660 -0.800 0.423
## (Intercept):3 -15.9091 1221.6440 NA NA
## fiberbran:1 -0.1542 0.8997 -0.171 0.864
## fiberbran:2 -0.6931 1.4142 -0.490 0.624
## fiberbran:3 0.5277 1803.6972 0.000 1.000
## fibergum:1 0.4055 1.1902 0.341 0.733
## fibergum:2 1.0986 1.2583 0.873 0.383
## fibergum:3 16.4199 1221.6442 0.013 0.989
## fiberboth:1 1.3218 1.0567 1.251 0.211
## fiberboth:2 0.1823 1.1328 0.161 0.872
## fiberboth:3 15.5036 1221.6443 0.013 0.990
##
## Names of linear predictors: loglink(P[Y=2]/P[Y=1]), loglink(P[Y=3]/P[Y=2]),
## loglink(P[Y=4]/P[Y=3])
##
## Residual deviance: 108.3914 on 30 degrees of freedom
##
## Log-likelihood: -54.1957 on 30 degrees of freedom
##
## Number of Fisher scoring iterations: 16
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):3'A função acat() usada no argumento famíly especifica um modelo de categoria adjacente.
Para aproveitar melhor a resposta ordinal da gravidade do inchaço, podemos estimar o modelo como \[ \log(\pi_j/\pi_{j+1}) = \beta_{j0} + \beta_1 \, \mbox{bran} + \beta_2 \, \mbox{gum} + \beta_3 \, \mbox{both}, \] onde os parâmetros são compartilhados entre as categorias de resposta para as fontes de fibra. Estime esse modelo usando acat(parallel = TRUE) no argumento vglm().
Examine as razões de chances usando os ajustes de ambos os modelos em (a) e (b). Compare seus valores.
Expresse por palavras o que é testado por \[ H_0 \, : \, \beta_{11} = \beta_{21} = \beta_{31}, \; \beta_{12} = \beta_{22} = \beta_{32}, \; \beta_{13} = \beta_{23} = \beta_{33} \quad \mbox{vs}.\\ H_1 \, : \, \mbox{Pelo menos uma desigualdade}\cdot \] Realize o teste usando um teste da razão de verossimilhanças (LRT).
17- Para os dados do grão de trigo, considere um modelo para estimar a condição do grão usando a variável explicativa densidade como termo linear.
Escreva uma função R que calcule a função log-verossimilhança para o modelo de regressão multinomial. Avalie a função nas estimativas dos parâmetros produzidas por multinom() e verifique se o valor calculado é o mesmo produzido por logLik(), use o objeto salvo de multinom() nesta função.
Maximize a função de log-verossimilhança usando optim() para obter os MLEs e a matriz de covariância estimada. Compare suas respostas com o que é obtido por multinom(). Observe que, para obter valores iniciais para optim(), uma abordagem é estimar modelos de regressão logística separados para \(\log(\pi_2/\pi_1)\) e \(\log(\pi_3/\pi_1)\). Esses modelos são estimados apenas para as observações que têm as respostas correspondentes, por exemplo, a \(Y = 1\) ou \(Y = 2\) para \(\log\big(\pi_2/\pi_1)\big)\).
Repita (a) e (b), mas agora para um modelo de regressão de chances proporcionais. Compare suas respostas com o que é obtido usando polr(). Observe que, para obter valores iniciais para optim(), uma abordagem é estimar modelos de regressão logística separados para \(logit\big(P(Y\leq 1)\big)\) e \(logit\big(P(Y\leq 2)\big)\). Os interceptos são usadas como estimativas iniciais para \(\beta_{10}\) e \(\beta_{20}\) e as duas inclinações são calculadas para obter uma estimativa inicial de \(\beta_1\).
18- Para os dados do grão de trigo, considere novamente os modelos de chances multinomiais e proporcionais estimados das Seções 3.3 e 3.4 que incluíram todas as variáveis explanatórias como termos lineares.
Construa gráficos dos modelos estimados com a probabilidade estimada no eixo \(y\) e density no eixo \(x\), onde as demais variáveis explicativas são definidas com seus valores médios. Interprete os gráficos e compare os modelos.
Como os gráficos de (a) mudam quando as variáveis explicativas restantes são definidas com valores diferentes de seus valores médios?
Construa gráficos adicionais como em (a), mas agora para variáveis explicativas diferentes de density no eixo \(x\).
19- As mudanças no habitat afetam diferentes espécies de vida selvagem de maneira diferente. Aves que caçam insetos no solo, por exemplo, podem ser menos afetadas pelo desmatamento da floresta tropical do que aquelas que se alimentam de néctares de plantas, como os beija-flores ou que extraem insetos das árvores, como os pica-paus. Root (1967) define uma guilda como “um grupo de espécies que exploram a mesma classe de recursos ambientais de maneira semelhante”. As aves capturadas em redes de neblina no exemplo da Seção 4.2.3 foram classificadas em guildas de acordo com seus hábitos alimentares. As guildas são rotuladas por fonte de alimento como Air-Arth, Leaf-Arth, Wood-Arth, Ground-Arth, Carnivore, Fruit, Seeds e Nectar, onde “Arth” significa artrópode (inseto e aranha). Queremos saber se a abundância relativa de aves de diferentes guildas está relacionada com a quantidade de “degradação” do habitat nas seis localidades. A degradação é definida aqui como um fator com três níveis que representam o quanto o habitat mudou em relação ao seu estado original de floresta. Os dois locais de floresta são “Low” (baixos), a borda da floresta e o fragmento de floresta são “Moderate” (Moderados) e os dois locais de pastagem são “High” (Altos). Os dados são apresentados na tabela abaixo. Realize uma análise para determinar quais guildas são mais afetadas ou menos afetadas pela degradação do habitat. \[ \mbox{Tabela cruzada de aves Equatorianas capturadas em redes de neblina} \\ \mbox{por degradação de habitat e guilda de fonte de alimento.} \\ \begin{array}{ccccc} \hline \mbox{} & & & \mbox{Degradação} & \\ \mbox{} & \mbox{Guilda} & \mbox{Alta} & \mbox{Moderada} & \mbox{Baixa} \\\hline \mbox{} & \mbox{Air-Arth} & 8 & 39 & 48 \\ \mbox{} & \mbox{Leaf-Arth} & 69 & 139 & 131 \\ \mbox{} & \mbox{Wood-Arth} & 2 & 39 & 50 \\ \mbox{} & \mbox{Ground-Arth} & 25 & 57 & 115 \\ \mbox{} & \mbox{Fruit} & 20 & 48 & 16 \\ \mbox{} & \mbox{Seeds} & 7 & 20 & 19 \\ \mbox{} & \mbox{Nectar} & 100 & 177 & 190 \\\hline \end{array} \]
20- A pneumoconiose (doença do pulmão negro) “é uma doença pulmonar que resulta da inalação de pó de carvão, grafite ou carbono artificial durante um longo período de tempo”. O quadro ou arquivo de dados pneumo no pacote VGAM contém contagens de mineiros de carvão agregados em oito grupos com base no tempo de trabalho na mina (em anos) e na gravidade da pneumoconiose (três níveis: normal, leve e grave). Analisar esses dados com os objetivos específicos de (1) entender a relação entre o tempo nas minas e a gravidade resultante da doença e (2) prever o status da doença dos mineiros após 5, 10, 15, 20 ou 25 anos de trabalho em uma mina.