Capítulo VIII. Modelos de regressão


O objetivo dos modelos de regressão é modelar a variação de uma variável de resposta quantitativa \(Y\) em termos da variação de uma ou várias variáveis explicativas \((X_1,\cdots,X_p)^\top\). Introduzimos esses modelos nos Capítulos III e VII onde os modelos lineares foram escritos como \[ Y = X\beta +\epsilon, \] onde \(Y_{n\times 1}\) é o vetor de observação para a variável de resposta, \(X_{n\times p}\) é a matriz de dados das \(p\) variáveis explicativas e \(\epsilon\) são os erros. Os modelos lineares não são restritos a lidar apenas com as relações lineares entre \(Y\) e \(X\). A curvatura é permitida, incluindo os termos de ordem superior apropriados na matriz de planejamento \(X\).


Exemplo VIII.1

Se \(Y\) representa a resposta e \(X_1\) e \(X_2\) dois fatores que explicam a variação de \(Y\) via um modelo quadrático: \[ y_i=\beta_0+\beta_1 x_{i1}+\beta_2 x_{i2}+\beta_3 x_{i1}^2+\beta_4 x_{i2}^2+\beta_5 x_{i1}x_{i2}+\epsilon_i, \qquad i=1,\cdots,n\cdot \]

O modelo acima pertence à classe dos modelos lineares porque é linear em \(\beta\). A matriz \(X\) é \[ X=\begin{pmatrix} 1 & x_{11} & x_{12} & x_{11}^2 & x_{12}^2 & x_{11}x_{12} \\ 1 & x_{21} & x_{22} & x_{21}^2 & x_{22}^2 & x_{21}x_{22} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_{n1} & x_{n2} & x_{n1}^2 & x_{n2}^2 & x_{n1}x_{n2}\end{pmatrix}\cdot \]


Observe também que os modelos não lineares puros às vezes podem ser reescritos como um modelo linear, escolhendo uma transformação apropriada das coordenadas das variáveis. Por exemplo, a função de produção de Cobb-Douglas \[ y_i=k x_{i1}^{\beta_1} x_{i2}^{\beta_2} x_{i3}^{\beta_3}, \]

onde \(y\) é o nível de produção de uma planta e \((x_1,x_2,X_3)^\top\) são três fatores de produção, por exemplo, trabalho, capital e energia; podem ser transformados em um modelo linear na escala logarítmica. Nós temos de fato \[ \log(y_i)=\beta_0+\beta_1\log(x_{i1})+\beta_2\log(x_{i2})+\beta_3\log(x_{i3}), \]

sendo \(\beta_0=\log(k)\) e os \(\beta_j\), \(j=1,2,3\) são as elasticidades \(\beta_j=\partial \log(y)/\partial \log(x_j)\).

Os modelos lineares são flexíveis e cobrem uma ampla classe de modelos. Se \(X\) tiver posto completo, eles poderão ser facilmente estimados por mínimos quadrados \(\widehat{\beta}=(X^\top X)^{-1}X^\top y\) e as restrições lineares nos \(\beta_j\) podem ser testadas usando as ferramentas desenvolvidas no Capítulo VII.

No Capítulo III, vimos que mesmo variáveis explicativas qualitativas podem ser usadas definindo a codificação apropriada dos valores nominais de \(X\). Neste capítulo, estenderemos nossa caixa de ferramentas mostrando como codificar esses fatores qualitativos de uma maneira que permita a introdução de vários fatores qualitativos, incluindo a possibilidade de interações. Isso abrange modelos ANOVA mais gerais do que os introduzidos no Capítulo III. Isso inclui os modelos ANCOVA, onde variáveis qualitativas e quantitativas estão presentes nas variáveis explicativas.

Quando a variável de resposta é qualitativa ou categórica, por exemplo, um indivíduo pode ser empregado ou desempregado, uma empresa pode ser falida ou não, a opinião de uma pessoa em relação a uma questão específica pode ser a favor, contra ou Indiferente, etc.; os modelos lineares precisam ser adaptados a essa situação em particular.

Os modelos mais úteis para esses casos serão apresentados na segunda parte do capítulo; Isso abrange os modelos log-lineares para tabelas de contingência, onde analisamos as relações entre várias variáveis categóricas e o modelo de logístico para respostas quantais ou binomiais, em que analisamos a probabilidade de estar em um estado em função das variáveis explicativas.


VIII.1. Modelos ANOVA gerais e modelos ANCOVA



VIII.1.1. Modelos ANOVA


Modelos de um fator


Na Seção III.5, introduzimos o exemplo de análise do efeito de um fator, três estratégias de marketing possíveis nas vendas de um produto, um pulôver. A maneira padrão de apresentar modelos ANOVA de um fator com \(p\) níveis é a seguinte \[ y_{kl}=\mu=\alpha_l +\epsilon_{kl}, \qquad k=1,\cdots,n_l, \quad l=1,\cdots,p, \]

sendo todos os \(\epsilon_{kl}\) independente.

Aqui \(l\) é o rótulo que indica o nível do fator e \(\alpha_l\) é o efeito do \(l\)-ésimo nível: ele mede o desvio da média global de \(y\), devido a esse nível do fator. Nesta notação, precisamos impor a restrição \(\sum_{l=1}^p \alpha_l=0\) para identificar \(\mu\) como a média de \(y\).

Esta apresentação é equivalente, mas um pouco diferente, à apresentada no Capítulo III, mas permite uma extensão mais fácil do caso de vários fatores. Observe também que aqui permitimos diferentes tamanhos de amostra para cada nível do fator, um design desequilibrado, mais geral do que o design equilibrado apresentado no Capítulo III.

Para simplificar a apresentação, suponha como no exemplo do pulôver, no qual \(p=3\). Nesse caso, pode-se tentar escrever o modelo acima sob a forma geral de um modelo linear usando três variáveis indicadoras \[ y_i = \mu+\alpha_1 x_{i1}+\alpha_2 x_{i2}+\alpha_3 x_{i3}+\epsilon_i \] onde \(x_{il}\) é igual a 1 ou 0 de acordo com a \(i\)-ésima observação, se pertencer ou não ao nível \(l\) do fator. Por simplicidade \(n_1=n_2=n_3\) e temos com \(\beta=(\mu,\alpha_1,\alpha_2,\alpha_3)^\top\) que \[ y=X\beta + \epsilon, \] onde a matriz de planejamento \(X\) é dada por \[ X=\begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{pmatrix}\cdot \]

Infelizmente, esse tipo de codificação não é útil porque a matriz \(X\) não é de posto completo, a soma de cada linha é igual à mesma constante 2 e, portanto, a matriz \(X^\top X\) não é invertível. Uma maneira de superar esse problema é alterar a codificação introduzindo a restrição adicional de que os efeitos somam zero. Há muitas maneiras de conseguir isso.

Observando que \(\alpha_3=-\alpha_1-\alpha_2\), não precisamos introduzir \(\alpha_3\) explicitamente no modelo. O modelo linear poderia realmente ser escrito como \[ y_i=\mu+\alpha_1 x_{i1}+\alpha_2 x_{i2}+\epsilon_i, \] com matriz de planejamento \[ X=\begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & -1 & -1 \\ 1 & -1 & -1 \end{pmatrix}, \] com automaticamente implicando que \(\alpha_3=-(\alpha_1+\alpha_2)\).

O modelo linear agora está correto com \(\beta=(\mu,\alpha_1,\alpha_2)^\top\). O estimador de mínimos quadrados \(\widehat{\beta}=(X^\top X)^{-1}x^\top y\) pode ser calculado, fornecendo o estimador dos parâmetros ANOVA \(mu\) e \(\alpha_l\)˛ \(l=1,2,3\). Qualquer restrição linear em \(\beta\) pode ser testada usando técnicas apropriadas. Por exemplo, a hipótese nula de nenhum fator de efeito \[ H_0 \, : \, \alpha_1=\alpha_2=\alpha_3=0, \] pode ser escrita matricialmente como \[ H_0 \, : \, A\beta=a, \] onde \[ A=\begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \] e \(a=(0,0)^\top\). 


pullover = read.table("http://leg.ufpr.br/~lucambio/MSM/pullover.txt", header = TRUE, sep = " ")
reg.pullover <- lm(Sales ~ Price+Advert+Asst.Hours, data= pullover)
summary(reg.pullover)
## 
## Call:
## lm(formula = Sales ~ Price + Advert + Asst.Hours, data = pullover)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.369  -9.406   1.599   5.151  19.729 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 65.66956   57.12507   1.150  0.29407   
## Price       -0.21578    0.32194  -0.670  0.52764   
## Advert       0.48519    0.08678   5.591  0.00139 **
## Asst.Hours   0.84373    0.37400   2.256  0.06491 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.7 on 6 degrees of freedom
## Multiple R-squared:  0.9067, Adjusted R-squared:  0.8601 
## F-statistic: 19.44 on 3 and 6 DF,  p-value: 0.001713


library(car)
linearHypothesis(reg.pullover, matrix(c(0,0,1,0,0,0,0,1), nrow=2, byrow=TRUE), c(0,0))
## Linear hypothesis test
## 
## Hypothesis:
## Advert = 0
## Asst.Hours = 0
## 
## Model 1: restricted model
## Model 2: Sales ~ Price + Advert + Asst.Hours
## 
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1      8 10080.8                                  
## 2      6   967.6  2    9113.2 28.255 0.0008843 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


restrito = lm(Sales ~ Price, data= pullover)
summary(restrito)
## 
## Call:
## lm(formula = Sales ~ Price, data = pullover)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.095  -8.898   2.036   9.805  64.725 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 210.7736    79.9837   2.635   0.0299 *
## Price        -0.3640     0.7571  -0.481   0.6435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.5 on 8 degrees of freedom
## Multiple R-squared:  0.02808,    Adjusted R-squared:  -0.09341 
## F-statistic: 0.2311 on 1 and 8 DF,  p-value: 0.6435
anova(restrito)
## Analysis of Variance Table
## 
## Response: Sales
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Price      1   291.3  291.27  0.2311 0.6435
## Residuals  8 10080.8 1260.10


Modelos de múltiplos fatores


A codificação acima pode ser estendida para situações mais gerais com muitas variáveis qualitativas, ou seja, fatores e com a possibilidade de interações entre os fatores. Suponha que, em um exemplo de marketing, as vendas de um produto possam ser explicadas por dois fatores: a estratégia de marketing com três níveis, como no exemplo do pulôver, mas também a localização da loja que pode ser em um grande shopping center ou em um local menos comercial, dois níveis para este fator. Podemos pensar também que há uma interação entre os dois fatores: a estratégia de marketing pode ter um efeito diferente em um shopping center do que em uma pequena área tranquila.

Para fixar a ideia, os dados são coletados conforme baixo:

dados.marketing = data.frame(cbind( Vendas=c(18,15,15,20,25,30,5,10,8,12,8,10,20,14,25), 
                        A = c("A1","A1","A1","A1","A1","A1","A2","A2","A2","A2","A2","A3","A3","A3","A3"),
                        B = c("B1","B2","B1","B2","B2","B2","B1","B2","B1","B2","B1","B1","B2","B1","B2")))
dados.marketing
##    Vendas  A  B
## 1      18 A1 B1
## 2      15 A1 B2
## 3      15 A1 B1
## 4      20 A1 B2
## 5      25 A1 B2
## 6      30 A1 B2
## 7       5 A2 B1
## 8      10 A2 B2
## 9       8 A2 B1
## 10     12 A2 B2
## 11      8 A2 B1
## 12     10 A3 B1
## 13     20 A3 B2
## 14     14 A3 B1
## 15     25 A3 B2

O modelo geral de dois fatores com interações pode ser escrito como \[ y_{ijk}=\mu+\alpha_i+\gamma_j+(\alpha\gamma)_{ij}+\epsilon_{ijk}, \]

\(i=1,\cdots,r\), \(j=1,\cdots,s\) e \(k=1,\cdots,n_{ij}\) onde as restrições de identificação sãO: \[ \sum_{i=1}^r \alpha_i =0, \qquad \sum_{j=1}^s \gamma_j=0, \qquad \sum_{i=1}^r (\alpha\gamma)_{ij} \] para \(j=1,\cdots,s\) e, para \(i=1,\cdots,r\) também exigimos \[ \sum_{i=1}^r (\alpha\gamma)_{ij}\cdot \]

Em nosso exemplo \(r=3\) e \(s=2\). Os \(\alpha\)’s medem o efeito da estratégia de marketing, em três níveis e os \(\gamma\)’s o efeito da localização, em dois níveis. Um valor positivo ou negativo de um desses parâmetros indicaria um efeito favorável ou desfavorável nas vendas esperadas; sendo a média global de vendas representada pelo parâmetro \(mu\). As interações são medidas pelos parâmetros \((\alpha\gamma)_{ij}\), \(i=1,\cdots,r\)˛ \(j=1,\cdots,s\), novamente restrições de identificação implica as \(r+s\) restrições nos termos das interações.

Por exemplo, um valor positivo de \((\alpha\gamma)_{11}\) indicaria que o efeito da estratégia de venda \(A_1\), anúncio em jornal local, se houver, é mais favorável nas vendas no local \(B_1\), em um grande centro comercial, do que em o local \(B_2\), não um centro comercial, com a relação \((\alpha\gamma)_{11}=-(\alpha\gamma)_{12}\). Como outro exemplo, um valor negativo de \((\alpha\gamma)_{31}\) indicaria que a estratégia de marketing \(A_3\), apresentação de luxo em vitrines, tem menos efeito, se houver, no local tipo \(B_1\) do que em \(B_2\): novamente \((\alpha\gamma)_{31}=-(\alpha\gamma)_{32}\), etc.

O bom é que é fácil estender a regra de codificação para o modelo de um fator para esta situação geral, a fim de apresentar o modelo um modelo linear padrão com a matriz de projeto apropriada \(X\). Para construir as colunas de \(X\) para o efeito de cada fator, precisaremos, como acima, de \(r-1\) ou \(s-1\) variáveis para codificar uma variável qualitativa com \(r\) ou \(s\), respectivamente níveis com a convenção definida acima no caso de um fator.

Para as interações precisaremos de \((r-1)\times (s-1)\) colunas adicionais que serão obtidas realizando o produto, elemento por elemento, das colunas de efeito principal correspondentes. Então, no final, para um modelo completo com todas as interações, temos os \(\big(1+r-1+s-1+(r-1)(s-1) \big)=rs\) parâmetros onde a primeira coluna de 1’s é para o intercepto, a constante.Ilustramos isso para nosso exemplo de marketing onde \(r=3\) e \(s=2\). Primeiro descrevemos um modelo sem interações. 

model.matrix(Vendas ~ A + B, data = dados.marketing)
##    (Intercept) AA2 AA3 BB2
## 1            1   0   0   0
## 2            1   0   0   1
## 3            1   0   0   0
## 4            1   0   0   1
## 5            1   0   0   1
## 6            1   0   0   1
## 7            1   1   0   0
## 8            1   1   0   1
## 9            1   1   0   0
## 10           1   1   0   1
## 11           1   1   0   0
## 12           1   0   1   0
## 13           1   0   1   1
## 14           1   0   1   0
## 15           1   0   1   1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$A
## [1] "contr.treatment"
## 
## attr(,"contrasts")$B
## [1] "contr.treatment"
ajuste.Vendas = lm(Vendas ~ A*B, data = dados.marketing)
summary(ajuste.Vendas)
## 
## Call:
## lm(formula = Vendas ~ A * B, data = dados.marketing)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -7.50  -2.00   1.00   1.75   7.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   16.500      2.963   5.569 0.000348 ***
## AA2           -9.500      3.825  -2.484 0.034777 *  
## AA3           -4.500      4.190  -1.074 0.310772    
## BB2            6.000      3.629   1.654 0.132611    
## AA2:BB2       -2.000      5.272  -0.379 0.713226    
## AA3:BB2        4.500      5.543   0.812 0.437812    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.19 on 9 degrees of freedom
## Multiple R-squared:  0.7851, Adjusted R-squared:  0.6658 
## F-statistic: 6.577 on 5 and 9 DF,  p-value: 0.007654
anova(ajuste.Vendas)
## Analysis of Variance Table
## 
## Response: Vendas
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## A          2 399.88 199.942 11.3891 0.003424 **
## B          1 153.65 153.648  8.7521 0.016001 * 
## A:B        2  23.80  11.901  0.6779 0.531819   
## Residuals  9 158.00  17.556                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Os contrastes podem ser usados ​​para fazer comparações específicas de tratamentos dentro de um modelo linear.

Um uso comum é quando um planejamento fatorial é usado, mas tratamentos de controle ou verificação são usados ​​além do planejamento fatorial. No primeiro exemplo abaixo, há dois tratamentos (D e C), cada um em dois níveis (1 e 2), e depois há um tratamento de Controle. A abordagem usada aqui é analisar o experimento como uma análise de variância unidirecional e, em seguida, usar contrastes para testar várias hipóteses.

Outro uso comum é quando existem vários tratamentos que podem ser pensados ​​como membros de um grupo. No segundo exemplo abaixo, há medidas para seis vinhos, alguns dos quais são tintos (Merlot, Cabernet, Syrah) e alguns são brancos (Chardonnay, Riesling, Gewürztraminer). Podemos comparar os tratamentos dentro do grupo do vinho tinto estabelecendo contrastes e realizando um teste F. Isso é análogo a testar o efeito principal do Vinho Tinto.

Os pacotes lsmeans e multcomp permitem testes ilimitados de contrastes de um grau, com uma correção de valor p para testes múltiplos. Eles também permitem um teste F para contrastes de várias linhas, por exemplo, ao testar dentro de grupos. A função aov no pacote de estatísticas nativo tem uma funcionalidade mais limitada.

Consulte os capítulos sobre Anova de uma via e Anova de duas vias para considerações gerais sobre a realização de análise de variância


VIII.1.2. Modelos ANCOVA


ANCOVA (ANálise de COVAriâncias) são modelos mistos onde algumas variáveis são qualitativas e outras são quantitativas. A mesma codificação da ANOVA será usada para a variável qualitativa. A matriz de planejamento \(X\) é concluída pelas colunas para as variáveis explicativas quantitativas \(x\). As interações entre uma variável qualitativa, um fator com \(r\) níveis e uma quantitativa \(x\) também são possíveis, isso corresponde a situações em que o efeito de \(x\) na resposta \(y\) é diferente de acordo com o nível do fator. Isso é conseguido adicionando à matriz de planejamento \(X\), uma nova coluna obtida pelo produto, elemento a elemento, da variável quantitativa com as variáveis codificadas para o fator, \(r-\) variáveis de interação se a variável categórica tiver \(r\) níveis.

Por exemplo, considere um modelo simples em que uma resposta \(y\) é explicada por uma variável explicativa \(x\) e um fator com dois níveis, por exemplo, o nível 1 de gênero para homens e nível 2 para mulheres, teríamos no caso \(n_1=n_2=3\) \[ X=\begin{pmatrix} 1 & x_1 & 1 & x_1 \\ 1 & x_2 & 1 & x_2 \\ 1 & x_3 & 1 & x_3 \\ 1 & x_4 & -1 & -x_4 \\ 1 & x_5 & -1 & -x_5 \\ 1 & x_6 & -1 & -x_6 \end{pmatrix}, \]

com \(\beta=(\beta_1,\beta_2,\beta_3,\beta_4)^\top\). O intercepto e a inclinação são \(\beta_1+\beta_3\) e \(\beta_1+\beta_4\) para os homens e \(\beta_1-\beta_3\) e \(\beta_1-\beta_4\) para as mulheres.


VIII.2. Respostas categóricas


VIII.2.1. Amostragem multinomial e tabelas de contingência


Em muitas aplicações, a variável de resposta de interesse é qualitativa ou categórica, no sentido de que a resposta pode assumir seu valor nominal em uma de \(K\) classes ou categorias. Freqüentemente, observamos as contagens \(y_k\), o número de observações na categoria \(k=1,\cdots,K\). Se o número total de observações \(n=\sum_{k=1}^K y_k\) for fixo e podemos assumir a independência das observações, obtemos uma amostra deum processo multinomial.

Denotemos por \(p_k\) a probalidade das observações na \(k\)-ésima categoria com \(\sum_{k=1}^K p_k=1\). Temos \(\mbox{E}(Y_k)=m_k=np_k\). A verossimilhança da amostra pode ser escrita como \[ L = \dfrac{n!}{\displaystyle \prod_{k=1}^K y_k!} \prod_{k=1}^K \left( \dfrac{m_k}{n}\right)^{y_k}\cdot \]

Nas tabelas de contingência, as categorias são definidas por várias variáveis qualitativas. Por exemplo, em uma tabela \(J\times K\) bidirecional, as observações contagens \(y_{jk}\), \(j=1,\cdots,J\) e \(k=1,\cdots,K\) são relatados para a linha \(j\) e a coluna \(k\). Aqui \(n=\sum_{j=1}^J\sum_{k=1}^K y_{jk}\).

Os modelos de log-linear introduzem uma estrutura linear nos logaritmos das frequências esperadas \(m_{jk}=\mbox{E}(Y_{jk})\), com \(\sum_{j=1}^J\sum_{K=1}^K p_{jk}=1\). Estruturas log-linear no \(m_{jk}\) imporão a mesma estrutura para o \(p_{jk}\), a estimativa do modelo será então obtida por máxima verossimilhança restrita. Tabelas de três vias \(J\times K\times L\) podem ser analisados da mesma maneira.

Às vezes, informações adicionais estão disponíveis sobre variáveis explicativas \(x\). Nesse caso, o modelo logit será apropriado quando a resposta categórica for binária \(K=2\). Introduziremos esses modelos quando a principal resposta de interesse for binária, por exemplo, tabelas \(2\times K\) ou \(2\times K\times L\). Além disso, mostraremos como eles podem ser adaptados ao caso de tabelas de contingência.


VIII.2.2. Modelo log-linear para tabelas de contingência


Tabelas com duas entradas (Two-way)


Considere um atable de dupla entrada \(J\times K\), onde \(y_{jk}\) é o número de observações com valor nominal \(j\) para o primeiro caractere qualitativo e valor nominal \(k\) para o segundo caractere. Como o número total de observações é fixo \[ n=\sum_{j=1}^J \sum_{k=1}^K y_{jk}, \]

existem \(JK-1\) células livres na tabela. A verossimilhança multinomial pode ser escrita como \[ L = \dfrac{n!}{\displaystyle \prod_{j=1}^J \prod_{k=1}^K y_{jk}!} \prod_{j=1}^J \prod_{k=1}^K \left( \dfrac{m_{jk}}{n}\right)^{y_{jk}}, \]

onde agora introduzimos uma estrutura log-linear para analisar o papel das linhas e das colunas para determinar os parâmetros \(m_{jk} = \mbox{E}(y_{jk})\) ou \(p_{jk}\).

  1. Modelo sem interações
    Suponha que não haja interação entre as linhas e as colunas: isso corresponde à hipótese de independência entre os dois caracteres qualitativos. Em outras palavras \(p_{jk}=p_j\, p_k\), para todo \(j,k\). Isso implica o modelo log-linear \[ \log\big( m_{jk}\big)=\mu + \alpha_j + \gamma_k , \qquad \mbox{para } \; j=1,\cdots,J,\; k=1,\cdots,K, \] onde, como nos modelos ANOVA para fins de identificação \[ \sum_{j=1}^J \alpha_j =\sum_{k=1}^K \gamma_k=0\cdot \]

    Usando os mesmos dispositivos de codificação acima, o modelo pode ser escrito como \[ \log(m)=X\beta\cdot \]

    Para o caso de uma tabela \(2\times 3\) temos: \[ \log(m)=\begin{pmatrix} \log(m_{11}) \\ \log(m_{12}) \\ \log(m_{13}) \\ \log(m_{21}) \\ \log(m_{22}) \\ \log(m_{23}) \end{pmatrix}, \qquad X=\begin{pmatrix} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & 0 \\ 1 & -1 & 0 & 1 \\ 1 & -1 & -1 & -1 \end{pmatrix}, \qquad \beta=\begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}, \]

    onde a primeira coluna de \(X\) é para o termo constante, a segunda coluna é a coluna codificada para o efeito da linha de 2 níveis e as duas últimas colunas são as colunas codificadas para o efeito da coluna de 3 níveis. A estimativa é obtida maximizando a log-verossimilhança, que é equivalente a maximizar a função \(L(\beta)\) em temos de \(\beta\): \[ L(\beta)=\sum_{j=1}^J \sum_{k=1}^K y_{jk}\log(m_{jk})\cdot \]

    A maximização está sob a restrição \(\sum_{jk} m_{jk}=n\). Em resumo, temos \(1+(J-1)+(K-1)-1\) parâmetros livres para \(JK-1\) células livres. O número de graus de liberdade no modelo é o número de células livres menos o número de parâmetros livres. É dado por \[ r=JK-1-(J-1)-(K-1)=(J-1)(K-1)\cdot \] No exemplo acima, temos, portanto, \((3-1)\times (2-1) = 2\) graus de liberdade. Os parâmetros originais do modelo podem ser estimados como: \[ \alpha_1=\beta_1, \quad \alpha_2=-\beta_1, \quad \gamma_1=\beta_2, \quad \gamma_2=\beta_3, \quad \gamma_3=-(\beta_1+\beta_3)\cdot \]


  2. Modelo com interações
    Nas tabelas bidirecionais, as interações entre as duas variáveis são de interesse. Este corresponde ao modelo geral (completo) \[ \log\big( m_{jk}\big)=\mu + \alpha_j + \gamma_k + (\alpha\gamma)_{jk} , \qquad \mbox{para } \; j=1,\cdots,J,\; k=1,\cdots,K, \]

    onde, além disso, temos as \(J+K\) restrições \[ \sum_{k=1}^K (\alpha\gamma)_{jk}=0, \qquad \mbox{para } \; j=1,\cdots,J \] e \[ \sum_{j=1}^J (\alpha\gamma)_{jk}=0, \qquad \mbox{para } \; k=1,\cdots,K\cdot \]

    Como no modelo ANOVA, as interações podem ser codificadas adicionando \((J-1)\times (K-1)\) colunas a \(X\), obtido pelo produto das variáveis codificadas correspondentes. Em nosso exemplo para tabela \(2\times 3\), a matriz de planejamento \(X\) é concluída com mais duas colunas: \[ X=\begin{pmatrix} 1 & 1 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & -1 & -1 & -1 & -1 \\ 1 & -1 & 1 & 0 & -1 & 0 \\ 1 & -1 & 0 & 1 & 0 & -1 \\ 1 & -1 & -1 & -1 & 1 & 1 \end{pmatrix}, \quad \beta=\begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \\ \beta_5 \end{pmatrix}\cdot \] Agora as interações são determinadas como \[ (\alpha\gamma)_{11}=\beta_4, \quad (\alpha\gamma)_{12}=\beta_5, \quad (\alpha\beta)_{13}=-\big((\alpha\gamma)_{11}+(\alpha\gamma)_{12} \big)=-(\beta_4+\beta_5), \]

    \[ (\alpha\beta)_{21}=-(\alpha\gamma)_{11}=-\beta_4, \]

    \[ (\alpha\beta)_{22}=-(\alpha\gamma)_{12}=-\beta_5, \quad (\alpha\beta)_{23}=-(\alpha\gamma)_{13}=\beta_4+\beta_5\cdot \]

    Temos novamente um modelo log-linear e a estimação passa pela maximização em \(\beta\) de \(L(\beta)\) sob a mesma restrição. O modelo com todos os termos de interação é chamado de modelo saturado. Neste modelo não há graus de liberdade, o número de parâmetros livres a serem estimados é igual ao número de células livres. Os parâmetros de interesse são as interações. Em particular, estamos interessados em testar seu significado. Essas questões serão abordadas a seguir.


Tabelas com três entradas (Three-way)


Os modelos apresentados acima para tabelas de dupla entrada podem ser estendidos a tabelas de ordem superior, mas a um custo de complexidade notacional. Mostramos como nos adaptar às tabelas de três entradas. Isso merece atenção especial devido à presença de interações de ordem superior no modelo saturado.

Uma tabela de três entradas \(J\times K\times L\) pode ser construída sob amostragem multinomial da seguinte forma: cada uma das \(n\) observações cai em uma, e apenas uma, categoria de cada uma das três variáveis categóricas com modalidades \(J\), \(K\) e \(L\), respectivamente. Acabamos com uma tabela tridimensional com \(JKL\) células contendo as contagens \(y_{jkl}\), onde \[ n=\sum_{j=1}^J \sum_{k=1}^K \sum_{l=1}^L y_{jkl}\cdot \] As contagens esperadas dependem das probabilidades desconhecidas \(p_{jkl}\) da maneira usual: \[ m_{jkl}=n p_{jkl}, \qquad j=1,\cdots,J, \; k=1,\cdots,K, \; l=1,\cdots,L\cdot \]

  1. Modelo saturado

    Um modelo log-linear saturado completo é lido da seguinte forma \[ \log\big(m_{jkl}\big)=\mu+\alpha_j +\beta_k+\gamma_l + (\alpha\beta)_{jk}+(\alpha\gamma)_{jl}+(\beta\gamma)_{kl}+(\alpha\beta\gamma)_{jkl}, \] para \(j=1,\cdots,J\), \(k=1,\cdots,K\) e \(l=1,\cdots,L\).

    As restrições são as seguintes, usando a notação \(\cdot\) para soma dos índices correspondentes: \[ \alpha_{(\bullet)}=\beta_{(\bullet)}=\gamma_{(\bullet)}=0, \]

    \[ (\alpha\beta)_{j\bullet}= (\alpha\gamma)_{j\bullet} = (\beta\gamma)_{k\bullet} =0, \]

    \[ (\alpha\beta)_{\bullet k} = (\alpha\gamma)_{\bullet l}=(\beta\gamma)_{\bullet l} =0 \] e \[ (\alpha\beta\gamma)_{jk \bullet} = (\alpha\beta\gamma)_{j \bullet l} = (\alpha\beta\gamma)_{\bullet kl} =0\cdot \]

    Os parâmetros \((\alpha\beta)_{jk}\), \((\alpha\gamma)_{jl}\) e \((\beta\gamma)_{kl}\) são chamados interações de primeira ordem. As interações de segunda ordem são os parâmetros \((\alpha\beta\gamma)_{jkl}\), eles permitem levar em consideração as heterogeneidades nas interações entre duas das três variáveis. Por exemplo, vamos \(l\) indicar as duas categorias de gênero \(L=2\), se supormos que \((\alpha\beta\gamma)_{jk1}= -(\alpha\beta\gamma)_{jk2}\neq 0\); queremos dizer que as interações entre as variáveis \(J\) e \(K\) não são as mesmas para ambas as categorias de gênero.

    A estimativa dos parâmetros do modelo saturao é obtida através da maximização da log-verossimilhança. No esquema de amostragem multinomial, corresponde à maximização da função: \[ L= \sum_{j,k,l} y_{jkl} \log\big( m_{jkl}\big), \] sob a restrição \(\sum_{j,k,l} m_{jkl}=n\).

    O número de graus de liberdade no modelo saturado é novamente zero. De fato, o número de parâmetros livres no modelo é \[ 1+(J-1)+(K-1)+(L-1)+(J-1)(K-1)+(J-1)(L-1)+(K-1)(L-1) + \\ \qquad \qquad \qquad \qquad (J-1)(K-1)(L-1)-1=JKL-1\cdot \]

    Isso é realmente igual ao número de células livres na tabela e, portanto, não há graus de liberdade.

  2. Modelos hierárquicos não saturados

    Como ilustrado acima, um modelo saturado não tem graus de liberdade. Modelos não saturados correspondem a modelos reduzidos em que alguns parâmetros são fixados para serem iguais a zero. Eles são, portanto, casos particulares do modelo saturado. Os modelos hierárquicos não saturados que consideraremos aqui são modelos em que uma vez que um conjunto de parâmetros é definido igual a zero, todos os parâmetros de ordem superior contendo os mesmos índices também são definidos iguais a zero.

    Por exemplo, se supormos \(\alpha_1=0\), consideramos apenas modelos não saturados onde também \((\alpha\beta)_{1k}=(\alpha\gamma)_{1l}=(\alpha\beta\gamma)_{1kl}=0\) para todos os valores de \(k\) e \(l\). Se apenas supormos que \((\alpha\beta)_{12}=0\), também assumimos isso para \((\alpha\beta\gamma)_{12l}\) para todo \(l\).

    Os modelos hierárquicos têm a vantagem de serem mais facilmente interpretáveis. De fato, sem essa hierarquia, os modelos seriam difíceis de interpretar. Qual seria, por exemplo, o significado do parâmetro \((\alpha\beta\gamma)_{12l}\), se soubermos que \((\alpha\beta)_{12}=0\)? A estimativa dos modelos não saturados será alcançada pela maneira usual, isto é, maximizando a função de log-verossimilhança como acima, mas as novas restrições do modelo reduzido.


VIII.2.3. Problemas de teste com dados de contagem


Um dos principais interesses práticos nos modelos de regressão para tabelas de contingência é testar restrições nos parâmetros de um modelo mais completo. Essas idéias de teste são criadas no mesmo espírito da Seção III.5 onde testamos restrições nos modelos ANOVA.

Nos modelos lineares, as estatísticas de teste são baseadas na comparação da bondade do ajuste para o modelo completo e para o modelo reduzido. A bondade do ajuste é medida pela soma residual de quadrados (RSS). A idéia aqui será a mesma, mas com uma medida mais apropriada para a qualidade do ajuste. Depois que um modelo for estimado, podemos calcular o valor previsto sob esse modelo para cada célula da tabela. Denotamos, como acima, o valor observado em uma célula por \(y_k\) e \(\widehat{m}_k\) indicará o valor esperado previsto pelo modelo. A bondade do ajuste pode ser apreciada medindo, de alguma forma, a distância entre a série de valores observados e previstos.

Duas estatísticas foram propostas: \(\chi^2\) de Pearson e o Deviance \(G^2\). Definidas como \[ \chi^2= \sum_{k=1}^K \dfrac{(y_k-\widehat{m}_k)^2}{\widehat{m}_k}, \] \[ G^2 = 2 \sum_{k=1}^K y_k \log\left( \dfrac{y_k}{\widehat{m}_k}\right), \] onde \(K\) é o número total de células da tabela. O desvio está diretamente relacionado à estatística da razão de log-verossimilhanças e geralmente é preferido porque pode ser usado para comparar modelos aninhados, como geralmente fazemos nesse contexto.

Sob a hipótese de que o modelo usado para calcular o valor previsto é verdadeiro, ambas as estatísticas para amostras grandes, são distribuídas aproximadamente como uma variável \(\chi^2\) com graus de liberdade d.f. dependendo do modelo. O d.f. pode ser calculado da seguinte forma: \[ d.f. = \mbox{No. de células livres} - \mbox{No. de parâmetros livres estimados}. \] Para modelos saturados, o ajuste é perfeito: \(\chi^2=G^2=0\) com \(d.f.=0\).


Exemplo VIII.4

Everitt e Dunn (1998) fornecem uma tabela de constagens tridimensional \(2\times 2\times 5\) de \(n=5.833\) pessoas entrevistadas. A contagem estava com medicamentos psicotrópicos prescritos na quinzena anterior à entrevista em função da idade e do sexo. Os dados estão resumidos na tabela abaixo, onde as categorias para os três fatores são “M” para masculino, “F” para feminino, “SIM” para sim, tendo tomado drogas, “NÃO” por não, não tomar drogas e as cinco categorias etárias: “A1” ( 16–29), “A2” (30–44), “A3” (45–64), “A4” (65–74), “A5” para mais de 74. A tabela fornece as frequências observadas \(y_{jkl}\) em cada uma das células da tabela de três vias: onde \(j\) significa gênero, \(k\) para drogas e \(l\) para categorias de idade.

O modelo saturado fornece as estimativas exibidas a seguir: 


contagens = c(21,32,70,43,19,683,596,705,295,99,46,89,169,98,51,738,700,847,336,196)
idades = rep(c("A1","A2","A3","A4","A5"),4)
sexo = c(rep("M",10),rep("F",10))
droga = rep(c(rep("SIM",5),rep("NÃO",5)),2)
tabela = glm(contagens ~ sexo*droga*idades, family = poisson)
summary(tabela)
## 
## Call:
## glm(formula = contagens ~ sexo * droga * idades, family = poisson)
## 
## Deviance Residuals: 
##  [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              6.603944   0.036811 179.404  < 2e-16 ***
## sexoM                   -0.077449   0.053096  -1.459 0.144656    
## drogaSIM                -2.775302   0.151968 -18.262  < 2e-16 ***
## idadesA2                -0.052863   0.052760  -1.002 0.316359    
## idadesA3                 0.137757   0.050355   2.736 0.006225 ** 
## idadesA4                -0.786833   0.065812 -11.956  < 2e-16 ***
## idadesA5                -1.325829   0.080356 -16.499  < 2e-16 ***
## sexoM:drogaSIM          -0.706670   0.268658  -2.630 0.008529 ** 
## sexoM:idadesA2          -0.083391   0.076978  -1.083 0.278671    
## sexoM:idadesA3          -0.106054   0.073609  -1.441 0.149646    
## sexoM:idadesA4          -0.052687   0.095839  -0.550 0.582497    
## sexoM:idadesA5          -0.605546   0.134247  -4.511 6.46e-06 ***
## drogaSIM:idadesA2        0.712858   0.189100   3.770 0.000163 ***
## drogaSIM:idadesA3        1.163500   0.173758   6.696 2.14e-11 ***
## drogaSIM:idadesA4        1.543159   0.190458   8.102 5.39e-16 ***
## drogaSIM:idadesA5        1.429013   0.218641   6.536 6.32e-11 ***
## sexoM:drogaSIM:idadesA2 -0.155391   0.343176  -0.453 0.650691    
## sexoM:drogaSIM:idadesA3  0.008769   0.308187   0.028 0.977299    
## sexoM:drogaSIM:idadesA4  0.013038   0.334669   0.039 0.968923    
## sexoM:drogaSIM:idadesA5  0.402278   0.399524   1.007 0.313986    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  5.7157e+03  on 19  degrees of freedom
## Residual deviance: -2.6201e-14  on  0  degrees of freedom
## AIC: 176.98
## 
## Number of Fisher Scoring iterations: 3


Vemos, por exemplo, que \(\widehat{\beta}_1<0\) (sexoM), então há menos homens do que mulheres no estudo embora não significativa a diferença, uma vez que \(\widehat{\beta}_{7}\) (sexoM:drogaSIM) também é negativo, parece que a tendência de os homens tomarem a droga é menos importante do que para as mulheres. Além disso, observe que de \(\widehat{\beta}_{12}\) a \(\widehat{\beta}_{15}\) formam uma sequência quase-crescente, de modo que o fator de idade parece aumentar a tendência de tomar o medicamento, até os 74 anos de idade. Note que, neste modelo saturado, não há graus de liberdade e o ajuste é perfeito, \(\widehat{m}_{jkl}= y_{jkl}\) para todas as células da tabela.

As interações de segunda ordem têm uma ordem mais baixa de magnitude, por isso queremos testar se forem significativamente diferentes de zero. Consideramos um modelo restrito em que \((\alpha\beta\gamma)_{jkl}\) estão todos definidos como zero. Isso pode ser alcançado testando \[ H_0 \, : \, \beta_{16}=\beta_{17}=\beta_{18}=\beta_{19}=0\cdot \]


tabela1 = glm(contagens ~ sexo*droga+sexo*idades+droga*idades, family = poisson)
summary(tabela1)
## 
## Call:
## glm(formula = contagens ~ sexo * droga + sexo * idades + droga * 
##     idades, family = poisson)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.03899  -0.64931  -0.02167   0.00051   1.03787   0.00686   0.15660   0.00684  
##        9        10        11        12        13        14        15        16  
## -0.00019  -0.41358   0.02646   0.41031   0.01397  -0.00034  -0.56902  -0.00660  
##       17        18        19        20  
## -0.14393  -0.00624   0.00018   0.30016  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        6.60419    0.03650 180.940  < 2e-16 ***
## sexoM             -0.07795    0.05217  -1.494  0.13512    
## drogaSIM          -2.77945    0.12883 -21.575  < 2e-16 ***
## idadesA2          -0.04767    0.05201  -0.917  0.35932    
## idadesA3           0.13773    0.04958   2.778  0.00547 ** 
## idadesA4          -0.78709    0.06411 -12.278  < 2e-16 ***
## idadesA5          -1.34759    0.07925 -17.004  < 2e-16 ***
## sexoM:drogaSIM    -0.69376    0.09274  -7.481 7.38e-14 ***
## sexoM:idadesA2    -0.09474    0.07484  -1.266  0.20556    
## sexoM:idadesA3    -0.10602    0.07113  -1.491  0.13608    
## sexoM:idadesA4    -0.05216    0.09046  -0.577  0.56418    
## sexoM:idadesA5    -0.54224    0.12452  -4.355 1.33e-05 ***
## drogaSIM:idadesA2  0.66776    0.15753   4.239 2.25e-05 ***
## drogaSIM:idadesA3  1.16636    0.14350   8.128 4.36e-16 ***
## drogaSIM:idadesA4  1.54735    0.15659   9.881  < 2e-16 ***
## drogaSIM:idadesA5  1.53331    0.18355   8.354  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 5715.6790  on 19  degrees of freedom
## Residual deviance:    2.3004  on  4  degrees of freedom
## AIC: 171.28
## 
## Number of Fisher Scoring iterations: 4


anova(tabela1,tabela, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: contagens ~ sexo * droga + sexo * idades + droga * idades
## Model 2: contagens ~ sexo * droga * idades
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         4     2.3004                     
## 2         0     0.0000  4   2.3004   0.6807


O modelo completo possui 0 graus de liberdade e o modelo reduzido possui 4 graus de liberdade. O deviance \(G^2\) é dado por 2.3004; possui 4 graus de liberdade. O \(p\)-valor do modelo restrito é de 0.6807, portanto, não rejeitamos a hipótese nula, o modelo restrito sem interação da 2ª ordem. Em outras palavras, a idade não interfere nas interações entre gênero e drogas, ou equivalentemente, o gênero não interfere nas interações entre idade e drogas. O leitor pode verificar se as interações de primeira ordem são significativas, tomando, por exemplo, o modelo sem interações da segunda ordem como o novo modelo completo e testando um modelo reduzido, onde todas as interações de primeira ordem são definidas como zero.


VIII.2.4. Modelo logístico


Os modelos logísticos são úteis para analisar como as variáveis explicativas influenciam uma resposta binária \(y\). A resposta \(y\) pode assumir os dois valores 1 e 0 para denotar a presença ou ausência de uma certa característica qualitativa, por exemplo, uma pessoa pode ser empregada ou desempregada, uma empresa pode ser falida ou não, um paciente pode ser afetado por uma certa doença ou não, etc.. Os modelos logísticos são projetados para estimar a probabilidade de \(Y=1\) como uma função logística de combinações lineares de \(x\).

Os modelos logísticos podem ser adaptados à análise de tabelas de contingência, onde uma das variáveis qualitativas é binária. Obtém-se a probabilidade de estar em um dos dois estados dessa variável binária em função das outras variáveis. Nós nos concentramos aqui em tabelas \(2\times K\) e \(2\times K\times L\).


VIII.2.4.1. Modelos logísticos para resposta binária


Considere o vetor \(n\times 1\) \(y\) de observações em uma variável de resposta binária: um valor de 1 indicando a presença de uma característica qualitativa particular e um valor de 0, sua ausência. O modelo logístico assume que a probabilidade de observar \(y_i=1\) dado um valor particular de \(x_i=(x_{i1},\cdots,x_{ip})^\top\) é dado pela função logística de uma pontuação, uma combinação linear de \(x\): \[ p(x_i)=P(y_i=1|x_i) = \dfrac{\exp\big( \beta_0+\sum_{j=1}^p \beta_j x_{ij}\big)}{1+\exp\big( \beta_0+\sum_{j=1}^p \beta_j x_{ij}\big)}\cdot \]

Isso implica a probabilidade da ausência do traço é: \[ 1-p(x_i)=P(y_i=0|x_i) = 1-\dfrac{\exp\big( \beta_0+\sum_{j=1}^p \beta_j x_{ij}\big)}{1+\exp\big( \beta_0+\sum_{j=1}^p \beta_j x_{ij}\big)}=\dfrac{1}{1+\exp\big( \beta_0+\sum_{j=1}^p \beta_j x_{ij}\big)}, \]

que implica \[ \log\left(\dfrac{p(x_i)}{1-p(x_i)} \right)=\beta_0+\sum_{j=1}^p \beta_j x_{ij}\cdot \]

Isso indica que o modelo logístico é equivalente a um modelo log-linear para a razão de chances \(p(x_i)/\big(1-p(x_i)\big)\). Um valor positivo de \(\beta_j\) indica uma variável explicativa \(x_j\) que favorecerá a presença do traço, pois melhora as chances. Um valor zero de \(\beta_j\) corresponde à ausência de efeito desta variável no aparecimento do traço qualitativo.

Para observações i.i.d., a função de verossimilhança é: \[ L(\beta_0,\beta)=\prod_{i=1}^n p(x_i)^{y_i} \big( 1-p(x_i)\big)^{1-y_i}\cdot \] Os estimadores de máxima verossimilhança dos \(\beta\)’s são obtidos como a solução do problema de maximização não linear \[ (\widehat{\beta}_0,\widehat{\beta}) = \arg \max_{\beta_0,\beta} \log \big(L(\beta_0,\beta) \big), \] onde \[ \log \big(L(\beta_0,\beta) \big) = \sum_{i=1}^n \Big(y_i \log\big(p(x_i)\big)+(1-y_i)\log\big(1-p(x_i) \big) \Big)\cdot \]

A teoria assintótica dos estimadores de máxima verossimilhança se aplica e, portanto, a inferência assintótica em \(\beta\) está disponível, como testes de hipóteses ou intervalos de confiança.


Exemplo VIII.5

Os arquivo de dados contém os indicadores de lucratividade (Profitability), alavancagem (Leverage) e falência (Bankruptcy) de 84 empresas.

library(aod); library(ggplot2)
Bankruptcy = read.csv(file = "http://leg.ufpr.br/~lucambio/MSM/bankruptcy.txt", header = TRUE, sep = "")
head(Bankruptcy)
##   Profitability Leverage Bankruptcy
## 1      0.022806   0.7816          1
## 2     -0.063584   0.7325          1
## 3     -0.227860   1.2361          1
## 4      0.021364   0.7350          1
## 5      0.042058   0.6339          1
## 6      0.021662   0.8614          1

O conjunto de dados contém informações sobre 42 das maiores empresas que entraram com pedido de proteção contra credores sob o Capítulo 11 do Código de Falências dos EUA em 2001–2002 após o crash da bolsa de 2000. As empresas falidas foram comparadas com 42 empresas sobreviventes com as capitalizações mais próximas e os mesmos códigos de classificação da indústria dos EUA disponíveis através da Divisão de Finanças Corporativas da Securities and Exchange Commission (SEC 2004).

As informações de cada empresa foram coletadas dos relatórios anuais de 1998-1999 (SEC 2004), ou seja, três anos anteriores à inadimplência das empresas falidas. O conjunto de dados a seguir contém índices de lucratividade e alavancagem calculados, respectivamente, como o índice de lucro líquido (NI) e ativo total (TA) e o índice de passivo total (TL) e ativo total (TA).

Com estas informações calculamos Profitability=NI/TA e Leverage = TL/TA. O modelo logístico pode ser usado para avaliar a probabilidade de falência em função desses índices financeiros. 

Resposta = factor(Bankruptcy$Bankruptcy, levels = c(-1,1), labels = c(0,1))
par(mfrow=c(1,2), pch=19)
plot(Bankruptcy$Profitability[Resposta==1], Bankruptcy$Leverage[Resposta==1], 
     xlim = c(-0.6,0.2), ylim = c(0.2,1.4), xlab = "Lucratividade", ylab = "Alavancagem",
     main = "Empresas falidas")
grid()
plot(Bankruptcy$Profitability[Resposta==0], Bankruptcy$Leverage[Resposta==0],
     xlim = c(-0.6,0.2), ylim = c(0.2,1.4), xlab = "Lucratividade", ylab = "Alavancagem",
     main = "Empresas não falidas")
grid()

sapply(Bankruptcy, mean)
## Profitability      Leverage    Bankruptcy 
##  -0.003009786   0.706240476   0.000000000
sapply(Bankruptcy, sd)
## Profitability      Leverage    Bankruptcy 
##     0.1100849     0.2138259     1.0060061


resultado = glm(Resposta ~ Bankruptcy$Profitability+Bankruptcy$Leverage, family = binomial)
summary(resultado)
## 
## Call:
## glm(formula = Resposta ~ Bankruptcy$Profitability + Bankruptcy$Leverage, 
##     family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5077  -0.9351  -0.1920   1.0327   1.8484  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                -1.565      1.088  -1.438   0.1504  
## Bankruptcy$Profitability  -13.211      5.305  -2.490   0.0128 *
## Bankruptcy$Leverage         2.458      1.495   1.644   0.1001  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.45  on 83  degrees of freedom
## Residual deviance:  94.92  on 81  degrees of freedom
## AIC: 100.92
## 
## Number of Fisher Scoring iterations: 5


library(car)
Anova(resultado, type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: Resposta
##                          LR Chisq Df Pr(>Chisq)   
## Bankruptcy$Profitability   9.7366  1   0.001806 **
## Bankruptcy$Leverage        2.8499  1   0.091380 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


## Intervalos confidenciais utilizando log-verossimilhança perfilada
confint(resultado)
##                               2.5 %     97.5 %
## (Intercept)               -3.799458  0.5225578
## Bankruptcy$Profitability -24.937439 -4.1420381
## Bankruptcy$Leverage       -0.388756  5.5457417
## Intervalos conficenciais utilizando erros padrão
confint.default(resultado)
##                                2.5 %     97.5 %
## (Intercept)               -3.6979570  0.5679058
## Bankruptcy$Profitability -23.6095400 -2.8133988
## Bankruptcy$Leverage       -0.4718278  5.3874836

Podemos testar um efeito geral de classificação usando a função wald.test da biblioteca aod. A ordem em que os coeficientes são dados na tabela de coeficientes é a mesma que a ordem dos termos no modelo. Isso é importante porque a função wald.test se refere aos coeficientes por sua ordem no modelo.

Usamos a opção b na função wald.test para fornecer os coeficientes, enquanto Sigma fornece a matriz de variâncias e covariâncias dos termos de erro, finalmente Terms informa ao R quais termos no modelo devem ser testados, neste caso, os termos 1 e 2, são os termos para as covariáveis. 

wald.test(b = coef(resultado), Sigma = vcov(resultado), Terms = 1:2)
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 11.9, df = 2, P(> X2) = 0.0026

A estatística do teste qui-quadrado de 11.9, com dois graus de liberdade está associada a um \(p\)-valor de 0.0026, indicando que o efeito geral da alavancagem (Profitability) e lucratividade (Leverage) é estatisticamente significativo.

Também podemos testar hipóteses adicionais sobre as diferenças nos coeficientes. Abaixo testamos que o coeficiente para posto=2 é igual ao coeficiente para o posto=3. A primeira linha de código abaixo cria um vetor l que define o teste que queremos realizar. Nesse caso, queremos testar a diferença (subtração) dos termos para posto=2 e posto=3, ou seja, o 2º e o 3º termos no modelo são iguais.

Para contrastar esses dois termos, multiplicamos um deles por 1 e o outro por -1. Os outros termos no modelo não estão envolvidos no teste, então eles são multiplicados por 0. A segunda linha de código abaixo usa L=l para dizer ao R que desejamos basear o teste no vetor l. 

l <- cbind(0, 1, -1)
wald.test(b = coef(resultado), Sigma = vcov(resultado), L = l)
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 9.3, df = 1, P(> X2) = 0.0023


A estatística do teste qui-quadrado é 9.3 com 1 grau de liberdade está associada a um \(p\)-valor de 0.0023, indicando que a diferença entre o coeficiente para posto=2 e o coeficiente para posto=3 é estatisticamente significativa.

Você também pode exponenciar os coeficientes e interpretá-los como odds-ratios ou razão de chances. R fará esse cálculo para você. Para obter os coeficientes exponenciados, você diz a R que deseja exponenciar (exp) e que o objeto que deseja exponenciar é chamado de coeficientes e faz parte de resultado, coef(resultado)). Podemos usar a mesma lógica para obter as razões de chances e seus intervalos de confiança, expondo os intervalos de confiança anteriores. Para colocar tudo em uma tabela, usamos cbind para vincular os coeficientes e intervalos de confiança em coluna. 

## odds ratios 
exp(coef(resultado))
##              (Intercept) Bankruptcy$Profitability      Bankruptcy$Leverage 
##             2.090827e-01             1.829497e-06             1.167942e+01


## odds ratios and 95% CI
exp(cbind(OR = coef(resultado), confint(resultado)))
##                                    OR        2.5 %       97.5 %
## (Intercept)              2.090827e-01 2.238290e-02   1.68633551
## Bankruptcy$Profitability 1.829497e-06 1.478454e-11   0.01589043
## Bankruptcy$Leverage      1.167942e+01 6.778997e-01 256.14448958



VIII.2.4.2. Modelos logísticos para tabelas de contingência


O modelo logístico pode conter variáveis explicativas quantitativas e qualitativas. Neste último caso, a variável pode ser codificada de acordo com as regras descritas nas seções ANOVA/ANCOVA acima. Isso permite uma revisita às tabelas de contingência onde uma das variáveis é binária e é a variável de interesse.

Como a probabilidade de tomar um dos dois valores nominais pode ser avaliada em função das outras variáveis? Mantemos as notações da Seção VIII.1 e suponha, sem perda de generalidade, que a primeira variável com \(J=2\) é a variável binária de interesse.


Tabelas \(2\times K\) com amostragem binomial


Seja \(p_k\) a probabilidade de cair na primeira linha da \(k\)-ésima coluna, \(k=1,\cdots,K\). Como estamos interessados principalmente nas probabilidades \(p_k\) em função de \(k\), supomos aqui que \(y_{\bullet k}\) são fixos para \(k=1,\cdots,K\) ou trabalhamos condicionalmente no valor observado desses totais das colunas, onde \(y_{\bullet k}=\sum_{j=1}^J y_{jk}\). Portanto, temos \(K\) processos binomiais independentes com parâmetros \((y_{\bullet k},p_k)\). Como a variável da coluna é nominal, podemos usar um modelo ANOVA para analisar o efeito da variável da coluna em \(p_k\) através dos logaritmos das probabilidades \[ \log\left( \dfrac{p_k}{1-p_k}\right) = \eta_0+\eta_k, \qquad k=1,\cdot,K, \]

onde \(\sum_{k=1}^K \eta_k=0\).

Assim como nos modelos ANOVA, um dos interesses será testar \(H_0 \, : \, \eta_1=\cdots=\eta_K=0\). O modelo log-linear para as probabilidades tem seu equivalente em uma forma logística para \(p_k\) \[ p_k=\dfrac{\exp\big(\eta_0+\eta_k\big)}{1+\exp\big(\eta_0+\eta_k\big)}, \qquad k=1,\cdots,K\cdot \]


VIII.3. Exercícios


  1. Para o modelo ANOVA de um fator, mostre que se o modelo for equilibrado, ou seja, se \(n_1=n_2=n_3\), temos \(\widehat{\mu}=\overline{y}\). Se o modelo não estiver equilibrado, mostre que \[ \overline{y}=\widehat{\mu}+n_1\widehat{\alpha}_1+n_2\widehat{\alpha}_2+n_3\widehat{\alpha}_3\cdot \]

  2. Refazer os cálculos do Exemplo VIII.2 e teste se os principais efeitos da estratégia de marketing e da localização são significativos.

  3. Refazer os cálculos do Exemplo VIII.3 com o conjunto de dados CarData.

  4. Calcule o intervalo de previsão para as vendas de pulôver, Exemplo III.2, correspondentes a \(\mbox{Price}=120\).

  5. Refazer os cálculos do exemplo de habitação de Boston (Boston Housing) na Seção VIII.1.3.

  6. Queremos analisar as variações no consumo de pacotes de cigarros por mês \(y\) em função da marca: A ou B (Brand)), do preço (Price) por pacote e em função do gênero do fumante: M ou F (Gender).
    Os dados estão em:

    cigarettes = read.csv(file = "http://leg.ufpr.br/~lucambio/MSM/cigarettes.txt", 
                      header = TRUE, sep = "")  
head(cigarettes)
    ##    y Price Gender Brand
## 1 30  3.50      M     A
## 2  4  4.00      F     B
## 3 20  4.10      F     B
## 4 15  3.75      M     A
## 5 24  3.25      F     A
## 6 11  5.00      F     B

    (a) Além dos efeitos da marca, preço e gênero, teste se houver uma interação entre a marca e o preço.
    (b) Como a matriz de planejamento de um modelo completo com todas as interações entre as variáveis apareceu? Qual seria o número de graus de liberdade desse modelo?
    (c) Gostaríamos de introduzir um termo de curvatura para a variável de preço. Como podemos proceder? Teste se esse coeficiente é significativo.

  7. No Exemplo VIII.4 de medicamentos, teste se as interações de primeira ordem são significativas.