No Capítulo 1, apresentamos a distribuição binomial para modelar dados de contagem resultantes da observação de uma resposta binária em \(n\) tentativas de Bernoulli independentes. As contagens observadas dessa maneira são restritas a situar-se entre 0 e \(n\), com variabilidade que diminui conforme a probabilidade de sucesso se aproxima de 0 ou 1.

Os dados de contagem podem surgir de outros mecanismos que nada têm a ver com as tentativas de Bernoulli. O número de carros que atravessam uma determinada ponte em um determinado dia, o número de ervas daninhas em um terreno de cultivo, o número de pessoas na fila de um balcão de atendimento e o número de pintas no corpo de uma pessoa são exemplos de contagens que não são limitados por nenhum número específico de tentativas.

Neste capítulo, usamos a distribuição de Poisson para modelar contagens desse tipo. Começamos explorando o modelo para uma única contagem e, em seguida, estendemos essa estrutura para modelos de regressão de Poisson que podem lidar com qualquer número de variáveis explicativas de qualquer tipo. Finalmente, consideramos alguns casos especiais de modelos de Poisson que são usados para descrever situações que não se encaixam perfeitamente na estrutura de regressão de Poisson padrão. Mostramos em todos os casos como a inferência é conduzida e como o R pode ser usado para auxiliar na análise.


4.1 Modelo de Poisson para dados de contagem



4.1.1 Distribuição de Poisson


Considere uma situação em que algum tipo de evento ocorre durante um período fixo de tempo, como carros passando por uma ponte. É fácil imaginar que a velocidade ou intensidade com que os carros cruzam a ponte pode variar de acordo com a hora do dia; por exemplo, mais carros na hora do rush, menos às 3 da manhã ou dia da semana; mais em um dia de trabalho, menos em fins de semana e feriados.

Suponha que contemos carros por exatamente uma hora no mesmo horário todas as semanas, digamos quarta-feira, das 7h às 8h. Mesmo que não haja diferenças práticas na intensidade populacional de semana para semana, não esperaríamos que exatamente o mesmo número de carros atravessasse a ponte a cada semana. Algumas pessoas que normalmente cruzam a ponte durante esse horário podem pegar carona ou pegar o trânsito; podem estar doentes, de férias, em reunião ou em teletrabalho; ou eles podem simplesmente chegar mais cedo ou mais tarde do que o normal e perder os cortes para a hora da contagem. Assim, esperamos que a contagem real de carros durante esse período varie aleatoriamente, mesmo que a taxa populacional subjacente ao processo não esteja mudando.

Essa é a natureza da distribuição de Poisson para contagens. Se pudermos assumir que:

  1. todas as contagens realizadas em algum processo têm a mesma intensidade subjacente e

  2. o período de observação é constante para cada contagem,

então as contagens seguem uma distribuição de Poisson.

Se \(Y\) é uma variável aleatória de Poisson, então a função de probabilidade (PMF) de \(Y\) é \[ P(Y=y)=\dfrac{e^{-\mu} \mu^y}{y!}, \qquad y=0,1,2,\cdots, \] onde \(\mu> 0\) é um parâmetro. Podemos abreviar essa distribuição escrevendo \(Y\sim Po(\mu )\).

Adaptações deste modelo serão dadas posteriormente para casos onde a intensidade não é constante para todas as contagens, Seção 4.2 e para quando o período de observação não é constante para todas as contagens, Seção 4.3.


Propriedades da distribuição de Poisson


A distribuição de Poisson tem várias propriedades atraentes que a tornam um modelo conveniente para se trabalhar. Em primeiro lugar, a forma da distribuição é bastante simples. Existe apenas um parâmetro \(\mu\), que representa a média e a variância da distribuição: \[ \mbox{E}(Y) = \mbox{Var}(Y) = \mu\cdot \]

Observe que, ao contrário da distribuição normal, a variância da distribuição de Poisson muda conforme a média muda. Isso faz sentido porque as contagens são limitadas abaixo por 0. Por exemplo, suponha que as médias para dois grupos diferentes de contagens sejam 5 e 50. Então as contagens no primeiro grupo podem estar principalmente na faixa de 0 a 10, enquanto as do segundo grupo pode variar naturalmente de, digamos, 20-80. Há mais “espaço” para variabilidade com uma média maior, então esperamos que a variância aumente com a média.

O modelo de Poisson é mais específico do que isso: requer que a variância seja igual à média. Veremos na Seção 5.3 que esse requisito às vezes é muito restritivo.

Outra propriedade útil é que somas de variáveis aleatórias de Poisson também são variáveis aleatórias de Poisson: Se \(Y_1, Y_2,\cdots, Y_m\) são independentes com \(Y_k\sim Po(\mu_k)\) para \(k = 1,\cdots,m\), então \[ \sum_{k=1}^m Y_k \sim Po\Big(\sum_{k=1}^m \mu_k\Big)\cdot \]

Assim, os totais de várias contagens podem ser modelados com distribuições de Poisson, desde que as contagens constituintes possam.

A forma da distribuição de Poisson se aproxima da distribuição normal à medida que \(\mu\) cresce. De fato, a distribuição normal às vezes é usada como modelo para contagens, embora não seja uma distribuição discreta. O Exercício 22 demonstra quando a distribuição normal é uma aproximação razoável da Poisson.


Exemplo 4.1: carros em um cruzamento.


A interseção das ruas 33rd e Holdrege em Lincoln, Nebraska, é uma típica interseção norte-sul/leste-oeste de 4 vias, exceto que há um quartel de bombeiros localizado a aproximadamente 45 metros ao norte da interseção no lado oeste da rua.

Um backup de veículos esperando para ir para o sul no cruzamento pode bloquear a entrada do corpo de bombeiros, o que impediria que os veículos de emergência saíssem da estação.

Para determinar a probabilidade de que isso pudesse acontecer, 40 ciclos consecutivos de semáforos foram observados e o número de veículos parados no lado norte da interseção foi contado. Observe que não havia veículos restantes na interseção por mais de um ciclo de semáforo.

Abaixo, lemos os dados em R e calculamos a média e a variância para as 40 contagens:

stoplight = read.csv( file = "http://leg.ufpr.br/~lucambio/ADC/Stoplight.csv")
head(stoplight)
##   Observation vehicles
## 1           1        4
## 2           2        6
## 3           3        1
## 4           4        2
## 5           5        3
## 6           6        3
# Summary statistics
mean ( stoplight$vehicles )
## [1] 3.875
var( stoplight$vehicles )
## [1] 4.317308


A média amostral é \(\overline{y} = 3.875\approx 3.9\). A variância da amostra ou variância amostral \(s^2 = 4.3\), é bastante próxima da média, como esperado se a variável aleatória de contagem de carros seguir uma distribuição de Poisson. Uma distribuição de frequência e um histograma dos dados, mostrados na figura abaixo, são produzidos pelo código a seguir.

# Frequencies
table ( stoplight$vehicles ) # Note that y = 0, 1, ... , 8 all have positive counts
## 
## 0 1 2 3 4 5 6 7 8 
## 1 5 7 3 8 7 5 2 2
rel.freq <- table ( stoplight$vehicles ) / length ( stoplight$vehicles )
rel.freq2 <- c( rel.freq , rep (0, times = 7))
# Poisson calculations
y <- 0:15
prob <- round ( dpois (x = y, lambda = mean ( stoplight$vehicles )), 4)
# Observed and Poisson
data.frame (y, prob , rel.freq = rel.freq2 )
##     y   prob rel.freq
## 1   0 0.0208    0.025
## 2   1 0.0804    0.125
## 3   2 0.1558    0.175
## 4   3 0.2013    0.075
## 5   4 0.1950    0.200
## 6   5 0.1511    0.175
## 7   6 0.0976    0.125
## 8   7 0.0540    0.050
## 9   8 0.0262    0.050
## 10  9 0.0113    0.000
## 11 10 0.0044    0.000
## 12 11 0.0015    0.000
## 13 12 0.0005    0.000
## 14 13 0.0001    0.000
## 15 14 0.0000    0.000
## 16 15 0.0000    0.000
plot (x = y -0.1 , y = prob , type = "h", ylab = "Probabilidade", xlab = "Número de veículos", 
      lwd = 2, xaxt = "n")
axis ( side = 1, at = 0:15)
lines (x = y+0.1 , y = rel.freq2 , type = "h", lwd = 2, lty = "solid", col = "red")
abline (h = 0)
legend (x = 9, y = 0.15 , legend = c("Poisson", "Observado") , lty = c("solid", "solid") , 
        lwd = c(2 ,2) , col = c("black", "red") , bty = "n")
grid()


O gráfico mostra concordância geralmente boa das observações com a função de probabilidade Poisson assumindo \(\mu=\overline{y} = 3.875\). O único grande desvio aparece em \(y = 3\). Não acreditamos que haja qualquer razão física que impeça grupos de três carros de se alinharem no cruzamento, então assumimos que esse desvio é apenas resultado do acaso.

No geral, parece que a distribuição de Poisson fornece uma aproximação satisfatória para a distribuição de carros parados no semáforo. Deixamos o uso desta aproximação para estimar a probabilidade de bloqueio do corpo de bombeiros como Exercício 4.



4.1.2 Verossimilhança e inferência de Poisson


O parâmetro \(\mu\) no modelo de Poisson é estimado por máxima verossimilhança. Dada uma amostra aleatória \(y_1,\cdots,y_n\); de uma distribuição de Poisson com média \(\mu\), a função de verossimilhança é \[ L(\mu; y_1,\cdots, y_n) = \prod_{i=1}^n \dfrac{e^{-\mu} \mu^{y_i}}{y_i!}\cdot \]

O estimador de máxima verossimilhança (MLE) é a média amostral, \(\widehat{\mu} =\overline{y}\). Usando a segunda derivada da log-verossimilhança, a variância de \(\widehat{\mu}\) é facilmente mostrada como estimada por \(\widehat{\mbox{Var}}(\widehat{\mu}) = \widehat{\mu}/n\). Observe que isso difere de nossa estimativa usual para a variância de \(\overline{y}\), \(s^2/n\), porque a média e a variância são as mesmas na distribuição de Poisson. Em grandes amostras de uma distribuição de Poisson, a média e a variância da amostra serão quase idênticas.

Assim como ocorre com o parâmetro da probabilidade de sucesso em uma distribuição binomial, há muitas maneiras de formar intervalos de confiança para o parâmetro de média da Poisson. Veja Swift (2009) para discussão do nível de confiança e propriedades de comprimento para alguns desses métodos.

Discutimos aqui Wald, razão de verossimilhança, escore e métodos exatos para desenvolver um intervalo de confiança para \(\mu\). Em particular, o intervalo de confiança Wald de \(100(1-\alpha)\%\) é \[ \widehat{\mu}\pm Z_{1-\alpha/2} \sqrt{\widehat{\mu}/n}\cdot \]

O intervalo escore é derivado dos testes de hipótese de \(H_0 \, : \, \mu = \mu_0\), que podem ser conduzidos usando a estatística escore \[ Z_0 = \dfrac{\widehat{\mu}-\mu_0}{\sqrt{\widehat{\mu}/n}}, \]

comparando \(Z_0\) com quantis apropriados da distribuição normal padrão. Invertendo os resultados deste teste encontramos o intervalo de confiança, \[ \left( \widehat{\mu}+\dfrac{Z^2_{1-\alpha/2}}{2n}\right) \pm Z_{1-\alpha/2}\sqrt{\dfrac{\widehat{\mu}+Z^2_{1-\alpha/2}/4n}{n}}\cdot \]

O intervalo da razão de verossimilhanças (LR) é o conjunto de valores de \(\mu_0\) para o qual \[ -2\log\left( \dfrac{L(\mu_0 \, | \, y_1,\cdots, y_n)}{L(\widehat{\mu} \, | \, y_1,\cdots, y_n)}\right) \leq \chi^2_{1,1-\alpha}\cdot \]

Isso não tem uma solução fechada e deve ser encontrado usando procedimentos numéricos iterativos. Portanto, geralmente não é usado para esse problema simples. O intervalo exato é desenvolvido usando uma lógica semelhante à que foi usada para o intervalo de Clopper-Pearson para a probabilidade binomial da Seção 1.1.2. Sua forma é \[ \chi^2_{2n\widehat{\mu},\alpha/2}/2n < \mu < \chi^2_{2(n\widehat{\mu}+1),1-\alpha/2}/2n\cdot \]

Alguns desses métodos são muito semelhantes aos testes e intervalos de confiança usados para inferência sobre as médias da distribuição normal. A principal vantagem desses métodos baseados na Poisson é o uso da função de verossimilhança de Poisson para obter a relação \(\mbox{Var}(\overline{Y}) = \mu/n\), o que permite intervalos de confiança mais curtos e testes mais poderosos, especialmente em amostras menores.

A utilização da função de verossimilhança de Poisson também pode ser sua desvantagem, pois é comum que os processos que geram contagens tenham intensidade que não permaneça rigidamente constante. Isso faz com que as contagens se comportem como se tivessem mais variabilidade do que o modelo de Poisson espera. Isso é chamado de superdispersão e é abordado na Seção 5.3.

Portanto, \(\widehat{\mbox{Var}}(\widehat{\mu}) = \widehat{\mu}/n\) geralmente subestima a verdadeira variabilidade da contagem média. Embora \(s^2\) seja uma estimativa menos precisa de \(\mbox{Var}(Y)\) quando a distribuição de Poisson é exata, ela estima a variabilidade presente devido a todas as fontes, incluindo quaisquer variações na intensidade do processo. Portanto, \(t\)-testes comuns e intervalos de confiança para uma média populacional são mais robustos a desvios da suposição de Poisson e às vezes são usados com dados de contagem extraídos de uma única população. No entanto, observe que os intervalos de confiança desenvolvidos especificamente para a distribuição de Poisson podem ser usados mesmo quando \(n = 1\), enquanto o intervalo baseado na disribuição \(t\) não pode.


Exemplo 4.2: carros em um cruzamento.


Em seguida, encontramos os intervalos de confiança para o número médio de carros parados no cruzamento. A função poisson.test() no pacote stats pode ser usada para encontrar um intervalo exato. O pacote epitools também contém funções que podem calcular intervalos de confiança para o parâmetro médio de Poisson.

Em particular, pois.exact() produz o intervalo exato e pois.approx() produz o intervalo de Wald. Exemplos de seu uso estão incluídos ainda neste exemplo.

alpha <- 0.05
n <- length ( stoplight$vehicles )
mu.hat <- mean ( stoplight$vehicles )
# Wald
mu.hat + qnorm (p = c( alpha /2, 1- alpha /2))* sqrt (mu.hat/n)
## [1] 3.264966 4.485034
# Score
(mu.hat + qnorm (p = c( alpha /2, 1- alpha /2)) /(2* n)) + qnorm (p =
c( alpha /2, 1- alpha /2)) * sqrt (( mu.hat + qnorm (p =
1- alpha /2) /(4* n))/n)
## [1] 3.239503 4.510497
# Exact
qchisq (p = c( alpha /2, 1- alpha /2) , df = c(2* n*mu.hat , 2*( n*mu.hat +1))) /(2* n)
## [1] 3.288979 4.535323
# Usual t- distribution based interval
t.test (x = stoplight$vehicles , conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  stoplight$vehicles
## t = 11.795, df = 39, p-value = 1.955e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3.210483 4.539517
## sample estimates:
## mean of x 
##     3.875


O intervalo de confiança escore de 95% é \(3.2 < \mu < 4.5\); há 95% de chance de que os intervalos produzidos dessa maneira contenham o verdadeiro número médio de carros no semáforo. Outros intervalos de confiança são bastante semelhantes e têm uma interpretação semelhante.

alpha <- 0.05
n <- length(stoplight$vehicles)
mu.hat <- mean(stoplight$vehicles)
# Wald
mu.hat + qnorm(p = c(alpha/2, 1 - alpha/2))*sqrt(mu.hat/n)
## [1] 3.264966 4.485034
# Score - original used with published book
(mu.hat + qnorm(p = c(alpha/2, 1 - alpha/2))/(2*n)) + 
  qnorm(p = c(alpha/2, 1 - alpha/2)) * sqrt((mu.hat + qnorm(p = 1 - alpha/2)/(4*n))/n)
## [1] 3.239503 4.510497
# Score - corrected after book published (missing ^2)
(mu.hat + qnorm(p = 1 - alpha/2)^2/(2*n)) + qnorm(p = c(alpha/2, 1 - alpha/2)) *
  sqrt((mu.hat + qnorm(p = 1 - alpha/2)^2/(4*n))/n)
## [1] 3.311097 4.534939
# Exact
qchisq(p = c(alpha/2, 1 - alpha/2), df = c(2*n*mu.hat, 2*(n*mu.hat + 1)))/(2*n)
## [1] 3.288979 4.535323
# Other code for exact
# qgamma(p = c(alpha/2, 1-alpha/2), shape = c(n*mu.hat, n*mu.hat+1), scale = 1)/n
# c(qchisq(p = alpha/2, df = 2*n*mu.hat),qchisq(p = 1-alpha/2, df = 2*(n*mu.hat+1)))/(2*n)
# Exact using poisson.test in stats package
poisson.test(x = mu.hat*n)$conf.int / n
## [1] 3.288979 4.535323
## attr(,"conf.level")
## [1] 0.95
# Usual t-distribution based interval
t.test(x = stoplight$vehicles, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  stoplight$vehicles
## t = 11.795, df = 39, p-value = 1.955e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3.210483 4.539517
## sample estimates:
## mean of x 
##     3.875
#######################################################################
# Repeating CI calculations using functions from epitools package.
#  For each function, x=sum of data (using n*mu.hat here) and 
#  pt = sample size in our applications.  See help(pois.conf.int).  
library(epitools)
pois.exact(x = mu.hat*n, pt = n, conf.level = 0.95)
##     x pt  rate    lower    upper conf.level
## 1 155 40 3.875 3.288979 4.535323       0.95
pois.approx(x = mu.hat*n, pt = n, conf.level = 0.95)
##     x pt  rate    lower    upper conf.level
## 1 155 40 3.875 3.264966 4.485034       0.95


Nossas recomendações sobre quais intervalos usar para estimar o parâmetro médio da distribuição de Poisson refletem aquelas para o parâmetro de probabilidade de sucesso do binomial. O intervalo de Wald não é bom, porque seu verdadeiro nível de confiança pode ser muito baixo e raramente atinge seu nível declarado. Por outro lado, o intervalo exato pode ser excessivamente conservador e, portanto, mais amplo do que o necessário. O intervalo escore é um bom compromisso, pois seu verdadeiro nível de confiança é geralmente melhor do que o Wald e geralmente é mais curto do que o intervalo exato.

Swift (2009) compara vários intervalos (excluindo Wald, infelizmente) e encontra outros intervalos que ele prefere ao intervalo escore. Por exemplo, há melhorias no intervalo exato, semelhantes àquelas descritas para intervalos de confiança para o parâmetro de probabilidade binomial na Seção 1.1.2, que podem torná-lo mais curto, mantendo o mesmo nível de confiança verdadeiro. As diferenças entre os vários intervalos geralmente são pequenas, no entanto, como mostra o próximo exemplo. Para mais detalhes sobre diferentes intervalos de confiança possíveis e suas propriedades, consulte Swift (2009).


Exemplo 4.3:
Comparação dos intervalos de confiança para a média de Poisson.


Neste exemplo, comparamos os verdadeiros níveis de confiança e os comprimentos médios dos intervalos de confiança Wald, escore, exato e \(t\) para a média da população quando os dados são amostrados de uma distribuição de Poisson.

O algoritmo para encontrar os verdadeiros níveis de confiança para Wald e intervalos escore segue aquele usado na Seção 1.1.2 para uma probabilidade binomial de sucesso. Encontramos a probabilidade de que uma amostra de tamanho \(n\) tenha uma soma observada igual a \(0, 1,\cdots, 40n\) quando a amostra é retirada de \(Po(\mu)\) para \(\mu= 0.01, 0.02,\cdots, 20\). Para fazer isso, usamos o fato de que a soma tem uma distribuição \(Po(n\mu )\).

Tecnicamente, existem infinitos valores possíveis da soma para qualquer \(\mu\), porque a distribuição \(Po(\mu )\) coloca probabilidade diferente de zero nas respostas \(y = 0, 1,2,\cdots\). Todos, exceto alguns, acontecem com probabilidades tão pequenas que são ignoráveis. O programa calcula probabilidades para o número relativamente pequeno de valores de \(y\) que possuem probabilidades não desprezíveis.

Como o intervalo baseado em \(t\) depende tanto da média quanto da variância estimadas é muito mais difícil encontrar exatamente seu verdadeiro nível de confiança. Precisaríamos identificar todas as combinações possíveis de média e variância amostrais para um determinado tamanho de amostra e encontrar a probabilidade de Poisson para cada um dos intervalos de confiança resultantes.

Em vez disso, estimamos o verdadeiro nível de confiança usando simulação de Monte Carlo como na Seção 1.1.3. Simulamos conjuntos de dados de respostas de contagem de \(Po(\mu )\), usamos esses conjuntos de dados para encontrar intervalos de confiança e estimamos o verdadeiro nível de confiança com a proporção observada desses intervalos que contém o \(\mu\) que foi usado para gerar os dados. Os comprimentos médios estimados são apenas os comprimentos médios de todos os \(t\) intervalos calculados a partir de cada \(\mu\).

# Initial settings and calculations
alpha <- 0.05
n <- 5
sum.y <- c(0:(40*n))
y.bar <- sum.y/n
# Wald
lower.wald <- y.bar - qnorm(p = 1 - alpha/2)*sqrt(y.bar/n)
upper.wald <- y.bar + qnorm(p = 1 - alpha/2)*sqrt(y.bar/n)
length.wald <- upper.wald - lower.wald
# Score
lower.score <- (y.bar + qnorm(p = 1 - alpha/2)^2/(2*n)) - qnorm(p = 1 - alpha/2) * 
  sqrt((y.bar + qnorm(p = 1 - alpha/2)^2/(4*n))/n)
upper.score <- (y.bar + qnorm(p = 1 - alpha/2)^2/(2*n)) + qnorm(p = 1 - alpha/2) * 
  sqrt((y.bar + qnorm(p = 1 - alpha/2)^2/(4*n))/n)
length.score <- upper.score - lower.score
# Exact
lower.exact <- qchisq(p = alpha/2, df = 2*sum.y)/(2*n)
upper.exact <- qchisq(p = 1 - alpha/2, df = 2*(sum.y + 1))/(2*n)
length.exact <- upper.exact - lower.exact
# All mu's
mu.seq <- seq(from = .01, to = 20, by = .01)
# Save true confidence levels in a matrix
save.true.conf <- matrix(data = NA, nrow = length(mu.seq), ncol = 9)
# Create counter for the loop
counter <- 1
# Function to simulate Poisson data and calculate mean and variance
mv.func <- function(mu, R, n) {
  y <- matrix(data = rpois(R*n, lambda = mu), nrow = n, ncol = R) 
  mean.y <- apply(X = y, MARGIN = 2, FUN = mean)
  var.y <- apply(X = y, MARGIN = 2, FUN = var)
  cbind(mean.y, var.y)
}
# Loop over each pi that the true confidence level is calculated on each
for(mu in mu.seq) {
  pmf <- dpois(x = sum.y, lambda = n*mu)
  # Wald
  save.wald <- ifelse(test = mu > lower.wald, 
                      yes = ifelse(test = mu < upper.wald, yes = 1, no = 0), no = 0)
  wald <- sum(save.wald * pmf)
  wald.l <- sum(length.wald * pmf)
  # Score
  save.score <- ifelse(test = mu > lower.score, 
                       yes = ifelse(test = mu < upper.score, yes = 1, no = 0), no = 0)
  score <- sum(save.score * pmf)
  score.l <- sum(length.score * pmf)
  # exact
  save.exact <- ifelse(test = mu > lower.exact, 
                       yes = ifelse(test = mu < upper.exact, yes = 1, no = 0), no = 0)
  exact <- sum(save.exact * pmf)
  exact.l <- sum(length.exact * pmf)
  # t: Generate data, compute mean and variance with mv.func
  R = 1825
  mean.var <- mv.func(mu = mu, R = R, n = n)
  # t intervals
  lowert <- mean.var[,1] - qt(p = 1 - alpha/2, df = n - 1) * sqrt(mean.var[,2]/n)
  uppert <- mean.var[,1] + qt(p = 1 - alpha/2, df = n - 1) * sqrt(mean.var[,2]/n)
  length.t <- uppert - lowert
  save.t <- ifelse(test = mu > lowert, yes = ifelse(test = mu < uppert, yes = 1, no = 0), no = 0)
  tt <- mean(save.t)
  t.l <- mean(length.t)
# Save all coverage estimates and lengths for this mu
  save.true.conf[counter,] <- c(mu, wald, score, exact, tt, wald.l, score.l, exact.l, t.l)
  counter <- counter + 1
}
###################
# CHRIS EDITED HERE
# Plots: Confidence level B/W
par(mfrow = c(1,4)) #1x4 plotting grid
par(mar =  c(5, 4, 4, 0) + 0.1) #Extends right margin to edge of plot
plot(x = save.true.conf[,1], y = save.true.conf[,2], main = paste("Wald,n=", n, sep=""),
  xlab = expression(mu), ylab = "True confidence level", type = "l", ylim = c(0.8,1))
abline(h = 1 - alpha, lty = "dotted") #Can add , col="red"
axis(side = 4, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
plot(x = save.true.conf[,1], y = save.true.conf[,3], main = paste("Score,n=", n, sep=""),
  xlab = expression(mu), ylab = "", type = "l", ylim = c(0.8,1),
  yaxt = "n")
axis(side = 2, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
axis(side = 4, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
abline(h = 1 - alpha, lty = "dotted") #Can add , col="red"
plot(x = save.true.conf[,1], y = save.true.conf[,4], main = paste("Exact,n=", n, sep=""),
  xlab = expression(mu), ylab = "", type = "l", ylim = c(0.8,1), yaxt = "n")
axis(side = 2, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
axis(side = 4, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
abline(h = 1 - alpha, lty = "dotted") #Can add , col="red"
plot(x = save.true.conf[,1], y = save.true.conf[,5], main = paste("t,n=", n, sep=""),
 xlab = expression(mu), ylab = NA, type = "l", ylim = c(0.8,1), yaxt = "n")
axis(side = 2, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
axis(side = 4, at = seq(from = 0.8, to = 1, by = 0.05), labels = FALSE)
abline(h = 1 - alpha + 1.96*sqrt(alpha*(1 - alpha)/R), lty = "dotted")
abline(h = 1 - alpha - 1.96*sqrt(alpha*(1 - alpha)/R), lty = "dotted")
abline(h = 1 - alpha, lty = "dotted") #Can add , col="red"


A partir da figura acima e das investigações para outros tamanhos de amostra, é óbvio que o verdadeiro nível de confiança para o intervalo exato é sempre de pelo menos 95%; para \(\mu\) pequenos é consideravelmente maior. Se for necessária a manutenção estrita do nível de confiança, como em algumas situações regulatórias, esse intervalo ou uma de suas variantes deve ser usado.

O nível de confiança para o intervalo escore é bem mantido perto de 95% em qualquer tamanho de amostra. O intervalo escore só se torna errático para \(\mu\) muito pequeno, por exemplo, \(\mu< 1\), mas ainda é o único intervalo entre os três intervalos aproximados cujo verdadeiro nível de confiança não cai para 0 à medida que \(\mu\) diminui em direção a 0.

Portanto, esse intervalo pode ser recomendado para uso geral, desde que se espere que a distribuição de Poisson se mantenha.

par(mar =  c(5, 4, 4, 2) + 0.1) #default
# Could also  reduce the second element of mar, but the "True confidence level"
# y-axis label on the first plot would be left off. 
# It may still be clear from the figure caption what is plotted.
###################
# Plots: average length: B/W
plot(x = save.true.conf[,1], y = save.true.conf[,6], 
     main = paste("Mean length of Intervals, n=", n, sep = ""),
     xlab = expression(mu), ylim = c(0, floor(30/sqrt(n))), 
     ylab = "Mean Length", type = "l", lty = 3, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,7],  type = "l", lty = 1, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,8],  type = "l", lty = 4, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,9],  type = "l", lty = 2, lwd = 2)
legend(x = 0, y = floor(30/sqrt(n)), legend = c("Wald", "Score", "Exact", "t"), 
       lty = c("solid", "dotted", "dotdash", "dashed"), lwd = 2, bty = "n")
grid()

# Plots: Average length at small mu: B/W
plot(x = save.true.conf[,1], y = save.true.conf[,6], 
     main = paste("Mean length of Intervals, n=", n, sep=""),
  xlab = expression(mu), ylab = "Mean Length", type = "l", lty = 3, lwd = 2, xlim = c(0,1),
  ylim= c(0,ceiling(4/sqrt(n))))
lines(x = save.true.conf[,1], y = save.true.conf[,7],  type = "l", lty = 1, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,8],  type = "l", lty = 4, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,9],  type = "l", lty = 2, lwd = 2)
legend(x = 0, y = floor(30/sqrt(n)), legend = c("Wald", "Score", "Exact", "t"), 
       lty = c("solid", "dotted", "dotdash", "dashed"), lwd = 2, bty = "n")
grid()

# Plots: average length: Color
plot(x = save.true.conf[,1], y = save.true.conf[,6], 
     main = paste("Mean length of Intervals, n=", n, sep = ""),
     xlab = expression(mu), ylim = c(0,floor(30/sqrt(n))), 
     ylab = "Mean Length", type = "l", lty = 3, lwd = 2)
lines(x = save.true.conf[,1], y = save.true.conf[,7],  type = "l", lty = 1, lwd = 2, col="blue")
lines(x = save.true.conf[,1], y = save.true.conf[,9],  type = "l", lty = 2, lwd = 2, col="orange")
lines(x = save.true.conf[,1], y = save.true.conf[,8],  type = "l", lty = 4, lwd = 2, col="green")
legend(x = 0, y = floor(30/sqrt(n)), legend=c("Wald", "Score", "Exact", "t"), 
       lty = c(3,1,4,2), col = c("black", "blue", "green", "orange"), lwd = 2, bty = "n")
grid()

# Plots: Average length at small mu: Color
plot(x = save.true.conf[,1], y = save.true.conf[,6], 
     main = paste("Mean length of Intervals, n=", n, sep = ""),
     xlab = expression(mu), ylab = "Mean Length", type = "l", lty = 3, 
     lwd = 2, xlim = c(0,1), ylim = c(0,ceiling(4/sqrt(n))))
lines(x = save.true.conf[,1], y = save.true.conf[,7],  type = "l", lty = 1, lwd = 2, col="blue")
lines(x = save.true.conf[,1], y = save.true.conf[,8],  type = "l", lty = 4, lwd = 2, col="green")
lines(x = save.true.conf[,1], y = save.true.conf[,9],  type = "l", lty = 2, lwd = 2, col="orange")
legend(x = 0, y = ceiling(4/sqrt(n)), legend = c("Wald", "Score", "Exact", "t"), 
       lty = c(3,1,4,2), col = c("black", "blue", "green", "orange"), lwd = 2, bty = "n")
grid()


O nível de confiança verdadeiro estimado do intervalo \(t\) parece bom - principalmente dentro de seus limites - para \(n = 5\), desde que \(\mu\) seja grande o suficiente para que a distribuição de \(Y\) seja razoavelmente aproximada por um normal, \(\mu > 7\) ou mais. É uma alternativa viável ao intervalo escore nesses casos, embora crie intervalos mais longos, como mostra a figura acima. No entanto, em \(n = 2\) (não mostrado), o intervalo \(t\) é interrompido e não atinge mais 95% para qualquer \(\mu< 20\). O intervalo escore mantém seu nível declarado mesmo para \(n = 1\).



4.2 Modelos de regressão de Poisson para respostas de contagem


Como vimos no Capítulo 2, os modelos de regressão binomial (binária) fornecem uma estrutura poderosa dentro da qual uma ampla variedade de problemas pode ser abordada. Modelos de regressão para dados de contagem também foram desenvolvidos com base na distribuição de Poisson. Estes são genericamente chamados de modelos de regressão de Poisson e estão intimamente relacionados aos modelos de regressão logística pelo fato de que ambas as famílias são classes de modelos lineares generalizados descritos na Seção 2.3.

Um tipo especial de modelo de regressão de Poisson, chamado de modelo loglinear, não apenas replica a análise clássica de tabelas de contingência descritas nas Seções 1.2.3, 2.2.6 e 3.2, mas também pode estender essas análises para qualquer número de variáveis e pode acomodar tanto variáveis contínuas e ordinais. Os modelos de regressão de Poisson podem até ser usados para replicar análises de regressão logística, embora seu uso dessa maneira seja um tanto complicado.

As etapas de modelagem com uma regressão de Poisson são praticamente as mesmas realizadas em qualquer outro processo de modelagem de regressão. Essas etapas incluem

  1. especificar o modelo, ou seja, a distribuição, os parâmetros a serem modelados; por exemplo, probabilidades ou médias e sua relação com as variáveis explicativas;

  2. selecionar variáveis explicativas;

  3. estimar os parâmetros do modelo;

  4. avaliar o ajuste do modelo e

  5. realizar inferências nos parâmetros do modelo e outras quantidades de interesse.

Assumimos por enquanto que as variáveis explicativas a serem usadas no modelo já são conhecidas ou escolhidas, deixando para o Capítulo 5 qualquer discussão sobre o processo pelo qual elas podem ter sido selecionadas. Também deixamos a avaliação do ajuste do modelo para o Capítulo 5.


4.2.1 Modelo para média: Ligação logarítmica


Temos \(n\) observações de uma variável aleatória de resposta \(Y\) e \(p\geq 1\) variáveis explanatórias fixas, \(x_1,\cdots,x_p\). Assumimos que para cada observação \(i = 1,\cdots,n\), \[ Y_i \sim Po(\mu_i), \] onde \[ \mu_i = exp\big(\beta_0 +\beta_1 x_{i1} + \cdots +\beta_p x_{ip}\big)\cdot \]

Assim, nosso modelo linear generalizado tem um componente aleatório de Poisson, um componente sistemático linear \(\beta_0 +\beta_1 x_{1} + \cdots +\beta_p x_{p}\) e uma ligação logarítmica, \[ \log(\mu ) = \beta_0 +\beta_1 x_{1} + \cdots +\beta_p x_{p}\cdot \]

Cada uma dessas três especificações é uma suposição que pode ser questionada. Por exemplo, a distribuição de Poisson pode fornecer um ajuste ruim aos dados. Veremos mais adiante neste capítulo e no Capítulo 5 que existem diferentes distribuições que podem servir como componentes aleatórios para dados de contagem neste caso.

Além disso, o preditor linear poderia ser substituído por algo que relacionasse as variáveis explicativas aos parâmetros de maneira não linear, ao custo da complexidade interpretativa e computacional. Por fim, a função de ligação pode assumir uma forma alternativa. Em particular, a ligação identidade \(\mu = \beta_0 +\beta_1 x_{1} + \cdots +\beta_p x_{p}\), pode parecer uma forma mais simples de usar. No entanto, a distribuição de Poisson exige que \(\mu> 0\) e a ligação identidade pode levar a um valor não positivo de \(\mu\) para determinados valores das variáveis explicativas. A função de ligação logaritmo garante \(\mu> 0\), por isso é quase universalmente usada.

Uma consequência do uso da função de ligação logaritmica com um preditor linear é que as variáveis explicativas afetam a média da resposta multiplicativamente. Considere um modelo de regressão de Poisson com uma variável explicativa: \(\mu(x) = \exp\big(\beta_0 + \beta_1 x\big)\), onde nossa notação enfatiza que \(\mu\) muda em função de \(x\). Quando aumentamos a variável explicativa em \(c\) unidades, o resultado é \[ \mu(x + c) = \exp\big(\beta_0 + \beta_1 (x + c)\big) = \mu(x)\exp\big(c \beta_1\big)\cdot \]

Assim, a razão das médias em \(x + c\) e em \(x\) é \[ \dfrac{\mu(x+c)}{\mu(x)} = \dfrac{\exp\big(\beta_0 + \beta_1 (x + c)\big)}{\exp\big(\beta_0 + \beta_1 x \big)}=\exp\big(c \beta_1\big)\cdot \]

Isso leva a uma maneira conveniente de interpretar o efeito de \(x\):
A variação percentual na resposta média que resulta de uma alteração de \(c\) unidades em \(x\) é \(PC = 100(e^{c\beta_1}-1)\%\).

Se houvesse variáveis adicionais no modelo, a mesma interpretação se aplicaria, desde que mantivéssemos constantes as variáveis explicativas adicionais. Quando há termos de interação ou transformações envolvendo a variável explicativa de interesse, a razão de médias é mais complicada, mas pode ser derivada de maneira semelhante à mostrada para razões de chances na Seção 2.2.5. Essas derivações são o assunto do Exercício 6.


Exemplo 4.4: Médias de regressão de Poisson.


A figura abaixo mostra como é a função média \[ \mu(x) = \exp\big(\beta_0 +\beta_1 x\big) \] para um modelo de regressão de Poisson quando \(\beta_1 > 0\) ou \(\beta_1 < 0\). Ela também mostra as regiões que provavelmente contêm a maioria dos dados. Observe em particular que, como a variância e a média são iguais, essas regiões se expandem à medida que a média cresce.

# set parameters for mean functions 
beta0a <- .1
beta1a <- 1
beta0b <- 5 
beta1b <- -1
# Parameter set "a"
curve(expr = exp(beta0a + beta1a*x), from = 0, to = 5, xlab = expression(x), 
      ylab = expression(mu = exp(beta[0] + beta[1]*x)), 
      lty = "solid", ylim = c(0,200), lwd = 2) 
curve(expr = exp(beta0a + beta1a*x) + 2*sqrt(exp(beta0a + beta1a*x)), 
      add = TRUE, lty = "dotted", lwd = 2)
curve(expr = exp(beta0a + beta1a*x) - 2*sqrt(exp(beta0a + beta1a*x)), 
      add = TRUE, lty = "dotted", lwd = 2)
# Parameter set "b"
curve(expr = exp(beta0b + beta1b*x), add = TRUE, lty = "dashed", lwd = 2)
curve(expr = exp(beta0b + beta1b*x) + 2*sqrt(exp(beta0b + beta1b*x)), 
      add = TRUE, lty = "dotted", lwd = 2)
curve(expr = exp(beta0b + beta1b*x) - 2*sqrt(exp(beta0b + beta1b*x)), 
      add = TRUE, lty = "dotted", lwd = 2)
legend(x = 1.5, y = 200, legend = c(expression(beta[1] > 0), expression(beta[1] < 0), 
      "95% probability range"), cex = 0.9, lty = c("solid", "dashed", "dotted"), lwd = 2, bty = "n")
grid()


figura acima mostrando funções médias usando ligação logaritmo para \(\beta_0 = 0.1\); \(\beta_1 = 1\) na linha sólida e \(\beta_0 = 5\); \(\beta_1 = -1\) na linha tracejada. As linhas pontilhadas representam regiões aproximadas nas quais 95% das observações devem cair de acordo com o modelo de regressão de Poisson, \(\pm\) 2 desvios padrão.



4.2.2 Estimação e inferência de parâmetros


O modelo de regressão Poisson assume que as observações \(y_1, y_2,\cdots, y_n\) são independentes. Por exemplo, eles não são serialmente correlacionados nem formam clusters e, portanto, a verossimilhança é formada pelo produto das funções de probabilidade individuais.

Seguindo as mesmas etapas apresentadas anteriormente, isso leva à log-verossimilhança \[ \begin{array}{rcl} \log\big(L(\beta_0,\cdots,\beta_p \, | \, y_1,\cdots,y_n) \big) & = & \displaystyle \sum_{i=1}^n -\exp\big(\beta_0+\beta_1 x_{i1} + \cdots + \beta_p x_{ip}\big) \\ & & \qquad \qquad + y_i \big(\beta_0+\beta_1 x_{i1} + \cdots + \beta_p x_{ip}\big) - \log(y_i!) \cdot \end{array} \]

Como de costume, diferenciamos a log-verossimilhança em relação a cada parâmetro \(\beta_j\), definimos cada uma das \(p+1\) equações resultantes iguais a 0 e resolvemos o sistema de equações para encontrar os estimadores de máxima verossimilhança (MLEs) \(\widehat{\beta}_0;\widehat{\beta}_1,\cdots,\widehat{\beta}_p\). Como foi o caso da regressão logística na Seção 2.2.1, essas equações normalmente não fornecem soluções de forma fechada; portanto, as soluções devem ser encontradas usando procedimentos numéricos iterativos.

Uma vez que temos as estimativas MLEs para os coeficientes de regressão, podemos calcular MLEs para qualquer função dos coeficientes tomando a mesma função dos MLEs. Por exemplo, o valor ajustado ou previsão é \[ \widehat{\mu}_i = \exp\big(\widehat{\beta}_0 +\widehat{\beta}_1 x_{i1} +\cdots+\widehat{\beta}_p x_{ip}\big), \] e a alteração estimada na média para uma alteração de \(c\) unidades em \(x_j\) é \[ \widehat{PC} = e^{c\, \widehat{\beta}_j}, \] mantendo as outras variáveis explicativas constantes.


Inferência


As abordagens padrão para inferência na estimativa de máxima verossimilhança são usadas aqui. Podemos derivar intervalos de confiança do tipo Wald e testes para parâmetros de regressão e várias funções deles apelando para a aproximação normal em amostras grandes para a distribuição dos MLEs \(\widehat{\beta}_0;\widehat{\beta}_1,\cdots,\widehat{\beta}_p\).

Embora os métodos de Wald sejam relativamente fáceis de executar, eles nem sempre funcionam bem. Métodos de razão de verossimilhança são escolhas melhores quando as rotinas computacionais estão disponíveis.


Exemplo 4.5: Consumo de álcool.


DeHart et ai. (2008) descrevem um estudo no qual “bebedores moderados a pesados”, estes definidos como àqueles que consomem pelo menos 12 bebidas alcoólicas/semana para mulheres, 15 para homens foram recrutados para manter um registro diário de cada bebida consumida durante um período de estudo de 30 dias.

Os participantes também completaram uma variedade de escalas de classificação que cobrem eventos diários em suas vidas e itens relacionados à auto-estima. Entre as hipóteses dos pesquisadores estava a de que eventos negativos, particularmente aqueles envolvendo relacionamentos românticos, podem estar relacionados à quantidade de álcool consumida, especialmente entre aqueles com baixa autoestima.

O aspecto de medidas repetidas deste estudo, os 30 dias de medições em cada participante, está além do nosso escopo atual. Em vez disso, focamos no primeiro sábado de cada participante no estudo, porque o sábado é um dia em que o consumo de álcool normalmente pode ser alto. Vamos modelar o número de bebidas consumidas - nossa resposta de contagem - como uma função das variáveis que medem o total de eventos positivos e negativos. Cada uma dessas variáveis explicativas é uma combinação de vários itens pontuados em uma escala de 0 a 3. Valores altos de eventos positivos ou negativos representam números maiores de eventos e/ou eventos mais extremos.

A planilha de dados DeHartSimplified.csv contém 13 variáveis, medidas nos primeiros 6 dias, sábado é codificado como “6” na variável dayweek. Dados gentilmente fornecidos pelo Dr. Steve Armeli, Escola de Psicologia, Fairleigh Dickinson University.

Nem todas as variáveis são necessárias em nossa análise atual, portanto, nossa primeira etapa é criar um quadro de dados que contenha o subconjunto necessário desses dados. As variáveis são:

id: identificador do participante,

numall número de bebidas alcoólicas ou “bebidas”, consumidas em um dia,

negevent um índice para combinar o número total e a intensidade de eventos negativos experimentados durante o dia e

posevent o mesmo que negevent, exceto com eventos positivos.

Observe que apenas 89/100 participantes tiveram dados de sábado nos primeiros 6 dias.

dehart <- read.table ( file = "http://leg.ufpr.br/~lucambio/ADC/DeHartSimplified.csv", 
                         header = TRUE , sep = ",", na.strings = " ")
head ( dehart )
##   id studyday dayweek numall     nrel      prel  negevent  posevent gender rosn
## 1  1        1       6      9 1.000000 0.0000000 0.4000000 0.5250000      2  3.3
## 2  1        2       7      1 0.000000 0.0000000 0.2500000 0.7000000      2  3.3
## 3  1        3       1      1 1.000000 0.0000000 0.2666667 1.0000000      2  3.3
## 4  1        4       2      2 0.000000 1.0000000 0.5333333 0.6083333      2  3.3
## 5  1        5       3      2 1.333333 0.3333333 0.6633333 0.6933333      2  3.3
## 6  1        6       4      1 1.000000 0.0000000 0.5900000 0.6800000      2  3.3
##        age  desired    state
## 1 39.48528 5.666667 4.000000
## 2 39.48528 2.000000 2.777778
## 3 39.48528 3.000000 4.222222
## 4 39.48528 3.666667 4.111111
## 5 39.48528 3.000000 4.222222
## 6 39.48528 4.000000 4.333333
# Reduce data to what is needed for examples
saturday <- dehart [ dehart$dayweek == 6, c(1 ,4 ,7 ,8)]
head ( round (x = saturday , digits = 3))
##    id numall negevent posevent
## 1   1      9    0.400    0.525
## 11  2      4    2.377    0.924
## 18  4      1    0.233    1.346
## 24  5      0    0.200    1.500
## 35  7      2    0.000    1.633
## 39  9      7    0.550    0.625
dim( saturday )
## [1] 89  4


Os modelos de regressão de Poisson são ajustados usando glm(). Detalhes sobre esta função são fornecidos na Seção 2.2.1. A principal adaptação para a regressão de Poisson é que agora devemos especificar family = poisson(link = “log”).

As mesmas funções de acompanhamento da regressão logística: summary(), confint(), anova(), Anova(), predict() e assim por diante, são usadas para fornecer estimativas de parâmetros, intervalos de confiança e testes.

Começamos com um modelo simplista usando o número total de eventos negativos para estimar o número médio de bebidas:

mod.neg <- glm( formula = numall ~ negevent , family = poisson ( link = "log") , data = saturday )
summary (mod.neg )
## 
## Call:
## glm(formula = numall ~ negevent, family = poisson(link = "log"), 
##     data = saturday)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9855  -1.3563  -0.2745   0.4736   5.8854  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.52049    0.07524  20.208   <2e-16 ***
## negevent    -0.26118    0.13597  -1.921   0.0547 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 250.34  on 88  degrees of freedom
## Residual deviance: 246.39  on 87  degrees of freedom
## AIC: 505.79
## 
## Number of Fisher Scoring iterations: 5
100*( exp(mod.neg$coefficients [2]) - 1)
##  negevent 
## -22.98564
beta1.int <- confint (mod.neg , parm = "negevent", level = 0.95)
100*( exp( beta1.int) - 1)
##       2.5 %      97.5 % 
## -41.5294841  -0.3479521
library (car)
Anova (mod.neg )
## Analysis of Deviance Table (Type II tests)
## 
## Response: numall
##          LR Chisq Df Pr(>Chisq)  
## negevent   3.9495  1    0.04688 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Matching confidence level with p- value of LRT to demonstrate equivalence
confint (mod.neg , parm = "negevent", level = 1 -0.04688)
##         2.3 %        97.7 % 
## -5.406268e-01 -5.932462e-06


O primeiro modelo que ajustamos é \(Y_i\sim Po(\mu_i)\) com \(\log(\mu_i) = \beta_0 + \beta_1 \, \mbox{negevent}_i\), onde \(Y_i\) é o número de bebidas consumidas, ou seja, numall para a pessoa \(i = 1,\cdots, 89\). As estimativas dos parâmetros são \(\widehat{\beta}_0 = 1.52\) e \(\widehat{\beta}_1 = -0.26\). O parâmetro “slope” negativo indica que o número de bebidas está diminuindo à medida que os eventos negativos aumentam. Assim, um aumento de 1 unidade em eventos negativos leva a uma diminuição estimada de \(\widehat{PC} = 23.0\%\) no número de bebidas alcoólicas, com intervalo de confiança LR perfilada de 95% \(-41.5\% < PC < -0.3\%\).

Esta alteração é marginalmente significativa de acordo com um teste de Wald, \(p\)-valor = 0.0547 e LRT \(p\)-valor = 0.047. Para fins de ilustração, reexecutamos a configuração dos cálculos do intervalo de confiança para o \(p\)-valor do teste, para mostrar a equivalência entre o intervalo de confiança LR e o LRT. Como esperado, isso gera um intervalo de confiança para 1 com um ponto final em zero, confirmando que as duas funções estão realizando os mesmos cálculos.

plot(x = saturday$negevent, y = saturday$numall, xlab = "Negative Event Index", 
     ylab = "Alcoholic Drinks Consumed", pch = 19)
grid()
curve(expr = exp(mod.neg$coefficients[1] + x*mod.neg$coefficients[2]), add = TRUE, lwd = 2)
# Function to find confidence interval
ci.mu <- function(newdata, mod.fit.obj, alpha) {
  lin.pred.hat <- predict(object = mod.fit.obj, newdata = newdata, type = "link", se = TRUE)
  lower <- exp(lin.pred.hat$fit - qnorm(1-alpha/2) * lin.pred.hat$se)
  upper <- exp(lin.pred.hat$fit + qnorm(1-alpha/2) * lin.pred.hat$se)
  list(lower = lower, upper = upper)
}
# Test
ci.mu(newdata = data.frame(negevent = 2), mod.fit.obj = mod.neg, alpha = 0.05)
## $lower
##        1 
## 1.748805 
## 
## $upper
##        1 
## 4.209492
# Add confidence interval bands
curve(expr = ci.mu(newdata = data.frame(negevent = x), mod.fit.obj = mod.neg, alpha = 0.05)$lower,
    col = "blue", add = TRUE, lty = "dotdash")
curve(expr = ci.mu(newdata = data.frame(negevent = x), mod.fit.obj = mod.neg, alpha = 0.05)$upper,
    col = "blue", add = TRUE, lty = "dotdash")


A figura acima mostra um gráfico dos dados com a curva média estimada \[ \widehat{\mu} = \exp(1.52 - 0.26 \, \mbox{negevent}) = 4.57(0.77)^\mbox{negevent}, \] usando o código fornecido no programa correspondente a este exemplo.

É evidente no gráfico que a média diminui à medida que os eventos negativos aumentam, mas a mudança é um pouco pequena, diminuindo de cerca de 4.5 drinques para cerca de 2.5 drinques em toda a gama de eventos negativos. O gráfico também mostra algumas outras características interessantes. Muito poucos desses bebedores “moderados a pesados” não consumiram álcool em seu primeiro sábado no estudo. Uma pessoa relatou consumir 21 bebidas!

Além disso, a relação inversa observada entre os eventos negativos e o consumo de álcool é algo surpreendente. Considerando que estamos usando um modelo muito simples e analisando apenas uma fração dos dias de dados de cada participante do estudo, veríamos esses resultados como preliminares e que precisam de verificação adicional. Finalmente, uma questão legítima que não abordamos aqui é se essas pessoas se comportaram de maneira diferente porque estavam participando de um estudo, ou seja, registrar seu consumo fez com que mudassem seu comportamento de beber?. Se assim for, isso pode explicar por que nenhuma das pessoas que experimentaram os eventos negativos mais extremos estava entre os consumidores mais pesados de álcool neste dia.

Uma investigação mais completa dos dados, por exemplo, comparando o primeiro sábado com os sábados posteriores, pode ajudar a resolver essa questão. Em seguida, consideramos um modelo estendido que também leva em consideração os efeitos de eventos positivos no comportamento de beber.

Agora modelamos a média como \[ \log(\mu) = \beta_0 + \beta_1 \, \mbox{negevent}+ \beta_2 \, \mbox{posevent} + \beta_3 \, \mbox{negevent}\times \mbox{posevent}\cdot \] As estimativas dos parâmetros e outras análises são dadas abaixo:

mod.negpos <- glm ( formula = numall ~ negevent * posevent , 
                    family = poisson ( link = "log") , data = saturday )
summary ( mod.negpos )
## 
## Call:
## glm(formula = numall ~ negevent * posevent, family = poisson(link = "log"), 
##     data = saturday)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0208  -1.2673  -0.5198   0.5172   5.9333  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         1.7571     0.1565  11.228  < 2e-16 ***
## negevent           -1.1988     0.3231  -3.710 0.000207 ***
## posevent           -0.1996     0.1199  -1.665 0.095941 .  
## negevent:posevent   0.7525     0.2253   3.341 0.000836 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 250.34  on 88  degrees of freedom
## Residual deviance: 233.07  on 85  degrees of freedom
## AIC: 496.47
## 
## Number of Fisher Scoring iterations: 5
confint ( mod.negpos )
##                        2.5 %      97.5 %
## (Intercept)        1.4490794  2.06217629
## negevent          -1.8492097 -0.58240139
## posevent          -0.4388871  0.03049091
## negevent:posevent  0.3137587  1.19645560
Anova ( mod.negpos )
## Analysis of Deviance Table (Type II tests)
## 
## Response: numall
##                   LR Chisq Df Pr(>Chisq)    
## negevent            4.0153  1  0.0450894 *  
## posevent            1.9081  1  0.1671724    
## negevent:posevent  11.4129  1  0.0007293 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


O modelo estimado é \[ \log(\widehat{\mu}) = 1.76 - 1.20 \, \mbox{negevent} - 0.20 \, \mbox{posevent} + 0.75 \, \mbox{negevent} \times \mbox{posevent}\cdot \]

Uma característica interessante aqui é a força do termo de interação. O intervalo de confiança LR perfilada para \(\beta_3\) de confint() é 0.31 < \(\beta_3\) < 1.20, excluindo claramente 0. O LRT de Anova() coloca o \(p\)-valor para testar \(H_0 \, : \, \beta_3 = 0\) em 0.0007.

# Prep work to make plots: creating grid of points: x1=neg events, y1=pos events
#  Then creating matrix of estimated means.  Matrix form required by contour()
x1 <- seq(from = 0, to = 2.5, by = .01)
y1 <- seq(from = 0, to = 3.5, by = .01)
xy1 <- data.frame(expand.grid(negevent = x1, posevent = y1))
surface = matrix(predict(object = mod.negpos, newdata = xy1, type = "response"), nrow = length(x1))
colorpal <- cm.colors(max(saturday$numall))
drink.col2 <- colorpal[saturday$numall]
plot(x = saturday$negevent, y = saturday$posevent, xlim = c(0,2.5), ylim = c(0,3.5), 
     xlab = "Negative Event Index", ylab = "Positive Event Index", pch = 21, bg = drink.col2, 
     cex = 1.5, main = "Number of Drinks vs. Negative and Positive Events")
grid()
contour(x = x1, y = y1, z = surface, xlim = c(0,2.5), ylim = c(0,3.5), labcex = 1, 
        levels = c(1,2,3,4,5,6,7,8,9,10,15,20,30,40,50,60,80), add = TRUE, color = "gray10")


A figura acima mostra um gráfico dos dados sombreados de acordo com o número de bebidas, com tons mais escuros indicando mais bebidas. Os contornos da média estimada do modelo são sobrepostos. Vemos que para pequenos valores do índice de evento positivo, por exemplo, < 1, as sombras parecem ficar mais claras para valores maiores de eventos negativos.

No entanto, para valores maiores do índice de eventos positivos, por exemplo, > 2, as sombras tendem a ser mais claras para valores menores de eventos negativos. Isso descreve claramente a natureza da interação. Os contornos também mostram isso: eles estão aumentando de esquerda para a direita na metade superior do gráfico e diminuindo na metade inferior.

Sem interação, os contornos seriam linhas paralelas. Observe também que as médias estimadas no canto superior direito do gráfico, onde não temos dados, são obviamente absurdas. Isso serve como um lembrete dos perigos de extrapolar uma função de regressão além do intervalo dos dados.

# Plot of data and estimated mean at quartiles of POSEVENT
posev.colors <- gray(1 - saturday$posevent/max(saturday$posevent))
posev.quart <- summary(saturday$posevent)
plot(x = saturday$negevent, y = saturday$numall, xlab = "Negative Event Index", 
     ylab = "Alcoholic Drinks Consumed", pch = 21, bg = posev.colors, cex = 1.5, 
     main = "Number of Drinks vs. Negative Events")
curve(expr = exp(mod.negpos$coefficients[1] + x*mod.negpos$coefficients[2] + 
      posev.quart[2]*mod.negpos$coefficients[3] + x*posev.quart[2]*mod.negpos$coefficients[4]), 
      add = TRUE, lwd = 2, col = "red")
curve(expr = exp(mod.negpos$coefficients[1] + x*mod.negpos$coefficients[2] + 
      posev.quart[3]*mod.negpos$coefficients[3] + x*posev.quart[3]*mod.negpos$coefficients[4]), 
      add = TRUE, lty = "dashed", lwd = 2, col = "red")
curve(expr = exp(mod.negpos$coefficients[1] + x*mod.negpos$coefficients[2] + 
      posev.quart[5]*mod.negpos$coefficients[3] + x*posev.quart[5]*mod.negpos$coefficients[4]), 
      add = TRUE, lty = "dotted", lwd = 2, col = "red")
legend(x = 1.0, y = 20, legend = c("1st", "2nd", "3rd"), lty = c("solid", "dashed", "dotted"), 
       col = c("red","red","red"), lwd = 2, title = "Quartile of Positive Event Index", bty = "n")
grid()


O gráfico acima mostra as curvas médias estimadas para bebidas contra eventos negativos nos três quartis de eventos positivos, 0.68, 1.03 e 1.43; uma técnica sugerida por Milliken and Johnson (2001) para examinar interações envolvendo variáveis contínuas.

As três curvas mostram que, quando há eventos positivos baixos, há uma diminuição mais aguda nas bebidas à medida que os eventos negativos aumentam do que quando há eventos positivos mais altos. Também podemos calcular a razão das médias correspondentes a um aumento de 1 unidade em eventos em cada um dos três quartis de eventos positivos. Para um valor fixo de eventos positivos, digamos \(a\), a razão das médias é \[ \dfrac{\mu(\mbox{negevent}+1, a)}{\mu(\mbox{negevent}, a)} = \exp\big(\beta_1+a \, \beta_3\big)\cdot \]

Essas quantidades são calculadas abaixo, juntamente com as variações percentuais correspondentes, \(100\big(\exp(\beta_1 + a\beta_3) -1\big)\) e seus intervalos de confiança. O processo computacional usado por mcprofile() pode às vezes produzir intervalos de confiança ligeiramente diferentes em diferentes execuções do mesmo programa. As diferenças são geralmente insignificantes, mas podem ser eliminadas usando set.seed() antes de mcprofile().

# Estimate percentage change per unit negevent at three quartiles of posevent
posev.quart
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.6833  1.0333  1.1583  1.4333  3.4000
mean.ratio <- exp(mod.negpos$coefficients[2] + posev.quart[c(2,3,5)]*mod.negpos$coefficients[4])
mean.ratio
##   1st Qu.    Median   3rd Qu. 
## 0.5042619 0.6561968 0.8866492
100*(mean.ratio - 1)
##   1st Qu.    Median   3rd Qu. 
## -49.57381 -34.38032 -11.33508
# Profile LR interval using mcprofile
library(package = mcprofile)
K <- matrix(data = c(0, 1, 0, 1*posev.quart[2],
                     0, 1, 0, 1*posev.quart[3],
                     0, 1, 0, 1*posev.quart[5]), nrow = 3, ncol = 4, byrow = TRUE)
linear.combo <- mcprofile(object = mod.negpos, CM = K)  #Calculate -2log(Lambda)
ci.beta <- confint(object = linear.combo, level = 0.95)
# ci.beta$confint
100*(exp(ci.beta$estimate) - 1) #Verifies got same answer as above
##     Estimate
## C1 -49.57381
## C2 -34.38032
## C3 -11.33508
100*(exp(ci.beta$confint) - 1)
##       lower      upper
## 1 -68.16128 -23.137341
## 2 -54.17361  -8.861619
## 3 -36.66359  21.986574


Vemos que o consumo de álcool diminui no primeiro quartil de eventos positivos em quase 50% por unidade de aumento em eventos negativos e o intervalo de confiança correspondente exclui claramente 0, -68% < \(PC\) < -23%; enquanto no terceiro quartil de eventos positivos o aumento é de apenas 11% por unidade e o intervalo de confiança não exclui 0, -37% < \(PC\) < 22%.

O próximo passo nesta análise seria uma avaliação do ajuste do modelo. Isso é abordado na Seção 5.2.



4.2.3 Variáveis explicativas categóricas


As variáveis explicativas categóricas podem ser incluídas em um modelo de regressão de Poisson de maneira semelhante aos modelos de regressão para respostas binárias ou multicategorias. Uma variável explanatória categórica \(X\) com \(I\) níveis é convertida em \(I-1\) variáveis indicadoras, cada uma indicando um dos níveis de \(X\). Para conveniência de notação posterior e para emular a maneira como R parametriza suas variáveis indicadoras, denotamos essas variáveis como \(x_2, x_3,\cdots,x_I\).

Suponha, por enquanto, que \(X\) tenha quatro níveis e que seja a única variável explicativa do modelo. Então podemos escrever o modelo para a média como \[ \log(\mu)=\beta_0+\beta_2 x_2 +\beta_3 x_3 +\beta_4 x_4\cdot \] Isso especifica uma estrutura análoga a um modelo ANOVA, conforme descrito na tabela abaixo.

\[ \mbox{Tabela 4.1: Variáveis indicadoras para uma variável explicativa categórica de 4 níveis.}\\ \begin{array}{c|ccc|c} \hline \mbox{Níveis} & x_2 & x_3 & x_4 & \mbox{Logaritmo da média} \\\hline \mbox{1} & 0 & 0 & 0 & \log(\mu_1) = \beta_0 \\ \mbox{2} & 1 & 0 & 0 & \log(\mu_2) = \beta_0+\beta_1 \\ \mbox{3} & 0 & 1 & 1 & \log(\mu_3) = \beta_0+\beta_3 \\ \mbox{4} & 0 & 0 & 1 & \log(\mu_4) = \beta_0+\beta_4 \\\hline \end{array} \]

A tabela acima também demonstra que podemos escrever o modelo como \(\log(\mu_i) =\beta_0 +\beta_i x_i\) para o nível \(i = 2,\cdots,I\) de \(X\). Isso leva imediatamente à interpretação dos parâmetros e torna clara nossa escolha de rótulos de índice: \(\beta_0\) é a média logarítmica para o nível 1 de \(X\); \(\beta_2\) é a diferença em médias logarítmicas entre os níveis 2 e 1 e outros parâmetros são interpretados de forma semelhante a \(\beta_2\).

Os parâmetros de regressão \(\beta_2,\cdots,\beta_q\) são freqüentemente chamados de parâmetros de efeito para a variável categórica \(X\). Observe que \(\beta_2 = \log(\mu_2) - \log(\mu_1)\), o que implica que \(e^{\beta_2} =\mu_2/\mu_1\); com uma interpretação semelhante para \(\beta_3\) e \(\beta_4\). Assim, os parâmetros de efeito são muito úteis para comparar médias em diferentes níveis contra a média no primeiro nível. Além disso, a diferença entre quaisquer dois parâmetros de efeito também mede uma razão entre duas médias. Por exemplo, \(\log(\mu_2)-\log(\mu_3) =\beta_2 - \beta_3\), de modo que \(\mu_2/\mu_3 = e^{\beta_2-\beta_3}\). Se escolhermos escrever \(\beta_1 = 0\) por conveniência, então todas as razões entre pares de médias são da forma \(\mu_i/\mu_j= e^{\beta_i-\beta_j}\).

As inferências normalmente se concentram em estimar as médias para cada nível de \(X\) e comparar essas médias de alguma maneira escolhida. As médias são estimadas como \(\widehat{\mu}_i = e^{\widehat{\beta}_0+\widehat{\beta}_i}\); \(i = 1,\cdots,q\). Intervalos de confiança para médias individuais são encontrados exponenciando intervalos de confiança para as somas correspondentes de parâmetros de regressão, \(\beta_0 +\beta_i\).

É comum testar \(H_0 \, : \, \mu_1 = \mu_2 = \cdots = \mu_q\), que é equivalente a \(H_0 \, : \, \beta_2 = \cdots = \beta_q = 0\). Esta hipótese é facilmente testada usando um LRT. Outras comparações entre as médias são realizadas por meio de testes ou intervalos de confiança em combinações de parâmetros.

Por exemplo, conforme observado acima, \(\mu_i/\mu_j = e^{\beta_i-\beta_j}\), portanto, um intervalo de confiança para a diferença entre dois parâmetros de efeito pode ser exponenciado em um intervalo de confiança para a razão de suas médias. Os métodos de LR perfilada são preferidos para intervalos de confiança. A função do método confint.glm() pode encontrar intervalos de confiança para parâmetros individuais, permitindo assim comparações diretas entre \(\mu_1\) e qualquer outra média. Para comparar outras médias, pode-se reordenar os níveis da variável explicativa categórica, reestimar o modelo e usar confint() novamente (consulte a Seção 2.2.6).

Mais geralmente, os intervalos de confiança LR para combinações lineares de parâmetros estão disponíveis no pacote mcprofile. Os testes para médias ou razões individuais são realizados usando a mesma abordagem dos intervalos de confiança, criando combinações lineares dos parâmetros do modelo correspondente.

Vários pacotes em R têm a capacidade de produzir esses testes, por exemplo, a função deltaMethod() em car, a função glht() em multcomp ou a função summary() em mcprofile. Os testes também podem ser realizados informalmente “invertendo” o intervalo de confiança, por exemplo, determinando se um intervalo de confiança para \(\beta_2-\beta_3\) contém 0, embora isso não produza um \(p\)-valor.

Observe que trabalhamos continuamente com proporções de médias em vez de diferenças. Isso ocorre porque a ligação logaritmo leva muito naturalmente ao cálculo de diferenças na escala logaritmo, o que leva a razões das médias. Os estimadores de parâmetros de regressão são aproximadamente normalmente distribuídos e têm variâncias que são relativamente fáceis de estimar. As combinações lineares de variáveis aleatórias normais também são variáveis aleatórias normais, de modo que as combinações lineares de estimadores de parâmetros também são aproximadamente distribuídas normalmente.

Portanto, os intervalos de confiança para certas funções de combinações lineares de parâmetros, como \(e^{\beta_i-\beta_j}\), são fáceis de desenvolver e calcular. Por outro lado, funções de médias que não podem ser reduzidas a combinações lineares de parâmetros, como \[ \dfrac{\mu_2+\mu_3}{2} = \dfrac{e^{\beta_0+\beta_2}+e^{\beta_0+\beta_3}}{2} \] e \[ \mu_2-\mu_3 = e^{\beta_0+\beta_2}-e^{\beta_0+\beta_3}, \] são menos convenientes de se trabalhar. A função deltaMethod() em car permite estimar funções não lineares dos parâmetros, aplicando o método delta para encontrar um erro padrão aproximado para a função. Os resultados podem então ser usados para realizar inferências Wald manualmente.


Exemplo 4.6: Contagem de pássaros equatorianos.


O desmatamento de florestas para terras agrícolas geralmente resulta em uma mudança na composição da vida vegetal e animal na área. Em um estudo de espécies de aves que habitam regiões de florestas nubladas do Equador, os pesquisadores usaram redes de neblina para capturar e identificar aves em seis locais diferentes. Dois locais estavam dentro de floresta não perturbada, um estava em um fragmento de floresta cercado por terra desmatada, um estava na borda entre a floresta e a terra desmatada e dois estavam em pastagens que já haviam sido florestadas.

Várias visitas de redes foram feitas a cada local ao longo de três anos. Para cada ave capturada, os pesquisadores registraram a espécie, o sexo e várias outras medições. Ver Becker et al. (2008) para detalhes sobre o estudo. As contagens totais de todas as aves de cada visita estão listadas na tabela abaixo.

\[ \mbox{Tabela 4.2: Dados sobre a contagem de aves em 6 locais no Equador.}\\ \begin{array}{l|l} \hline \mbox{Locais} & \mbox{Contagens} \\\hline \mbox{Floresta (ForA) A} & 155, 84\\ \mbox{Floresta (ForB) B} & 77, 57, 38, 40, 51, 67 \\ \mbox{Beira (Edge)} & 53, 49, 44, 47 \\ \mbox{Fragmento (Frag)} & 35, 50, 66, 51, 49, 75 \\ \mbox{Pasto A (PasA)} & 40, 39 \\ \mbox{Pasto B (PasB)} & 33, 31, 38, 50 \\\hline \end{array} \]


Estamos interessados em comparar a abundância geral de aves em diferentes tipos de habitat. Para fazer isso, ajustamos um modelo de regressão de Poisson para estimar a contagem média de capturas de aves em cada local como um índice da abundância geral e comparamos essas médias em diferentes tipos de habitats. Assim, usamos glm() para estimar as seis médias usando 5 variáveis indicadoras que indicam os últimos cinco locais em ordem alfabética.

Os dados são inseridos e essas variáveis indicadoras são exibidas abaixo. Observe que os resultados de contrasts() mostram uma forma idêntica à Tabela 4.1.

alldata <- read.table(file = "http://leg.ufpr.br/~lucambio/ADC/BirdCounts.csv", 
                      sep = ",", header = TRUE)
head(alldata)
##    Loc Birds
## 1 ForA   155
## 2 ForA    84
## 3 ForB    77
## 4 ForB    57
## 5 ForB    38
## 6 ForB    40
library(contrast)
contrasts(as.factor(alldata$Loc))
##      ForA ForB Frag PasA PasB
## Edge    0    0    0    0    0
## ForA    1    0    0    0    0
## ForB    0    1    0    0    0
## Frag    0    0    1    0    0
## PasA    0    0    0    1    0
## PasB    0    0    0    0    1


Em seguida, ajusta-se o modelo e obtém-se o resumo.

#####################################################################
# Model using regular likelihood
# Fit Poisson Regression and test whether all means are equal
M1 <- glm(formula = Birds ~ Loc, family = poisson(link = "log"), data = alldata)
summary(M1)
## 
## Call:
## glm(formula = Birds ~ Loc, family = poisson(link = "log"), data = alldata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4322  -0.7594  -0.1302   0.8874   3.1038  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.87640    0.07198  53.853   <2e-16 ***
## LocForA      0.90692    0.09678   9.371   <2e-16 ***
## LocForB      0.13094    0.09062   1.445   0.1485    
## LocFrag      0.11874    0.09082   1.307   0.1911    
## LocPasA     -0.20010    0.13356  -1.498   0.1341    
## LocPasB     -0.23881    0.10844  -2.202   0.0277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 216.944  on 23  degrees of freedom
## Residual deviance:  67.224  on 18  degrees of freedom
## AIC: 217.87
## 
## Number of Fisher Scoring iterations: 4


Como \(\widehat{\beta}_0 = 3.88\), a contagem média estimada para o habitat de beira (edge) é \(\widehat{\beta}_1 = e^{3.88} = 48.3\). Os \(\widehat{\beta}_i\)’s restantes, \(i = 2,\cdots,6\) medem a diferença estimada em médias logarítmicas entre os locais fornecidos e a borda (edge), \(\log(\widehat{\mu}_i)-\log(\widehat{\mu}_1)\); \(i = 2,\cdots,6\). Por exemplo, a diferença em médias logarítmicas entre o fragmento e a aresta ou beira é representada por \(\widehat{\beta}_4 = 0.12\). Assim, a razão dessas duas médias é \(\widehat{\mu}_4/\widehat{\mu}_1 = e^{0.11874} = 1.126\).

Fornecemos o código no programa correspondente para este exemplo para mostrar que todas as médias previstas deste modelo são as mesmas que suas respectivas médias amostrais. Isso acontece porque estamos permitindo que o modelo estime uma média separada para cada habitat e os MLEs neste caso podem ser mostrados como as respectivas médias amostrais. Assim, de certa forma, parece que estamos fazendo um grande esforço extra para fazer algo que deveria ser mais fácil de fazer. A vantagem da abordagem baseada em modelo reside nos numerosos procedimentos de inferência que

Para testar se todas as médias são iguais, \[ H_0 \,: \, \mu_1=\mu_2=\cdots=\mu_6 \\ H_1 \, : \, \mbox{Nem todas as médias são iguais}, \] usamos anova() para executar o LRT; Anova() também pode ser usado:

library(car)
Anova(M1, test = "F") #Also could use Anova() from car package
## Analysis of Deviance Table (Type II tests)
## 
## Response: Birds
## Error estimate based on Pearson residuals 
## 
##            Sum Sq Df F value    Pr(>F)    
## Loc       149.720  5  8.0238 0.0003964 ***
## Residuals  67.174 18                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Get predicted means and CIs
pred.data <- data.frame(Loc = c("ForA", "ForB", "Frag", "Edge", "PasA", "PasB"))
means <- predict(object = M1, newdata = pred.data, type = "link", se.fit = TRUE)
alpha <- 0.05
# Wald CI for log means
lower.logmean <- means$fit + qnorm(alpha/2)*means$se.fit
upper.logmean <- means$fit + qnorm(1 - alpha/2)*means$se.fit
# Combine means and confidence intervals in count scale
mean.wald.ci <- data.frame(pred.data, round(cbind(exp(means$fit), exp(lower.logmean), 
                                                  exp(upper.logmean)), digits = 2))
colnames(mean.wald.ci) <- c("Location", "Mean", "Lower", "Upper")
mean.wald.ci
##   Location   Mean  Lower  Upper
## 1     ForA 119.50 105.27 135.65
## 2     ForB  55.00  49.37  61.27
## 3     Frag  54.33  48.74  60.56
## 4     Edge  48.25  41.90  55.56
## 5     PasA  39.50  31.68  49.25
## 6     PasB  38.00  32.41  44.55

Esta hipótese nula é fortemente rejeitada, \(-2\log(\Lambda) = 67.2\) e o \(p\)-valor é inferior a \(2.2\times 10^{-16}\), então claramente nem todos os locais têm a mesma média. Para entender a causa desse resultado, um gráfico das médias e intervalos de confiança correspondentes pode ajudar.

Mostramos a seguir como calcular médias estimadas e intervalos de confiança LR correspondentes usando mcprofile(). Começamos especificando a matriz de coeficientes \(K\) para estimar cada \(\log(\widehat{\mu}_i) = \widehat{\beta}_0 + \widehat{\beta}_i\). As fileiras estão dispostas em ordem decrescente de habitat florestal: forest, fragment, edge, pastures.

######################################################################
# Get LR intervals for means.  Each log-mean is a combination of intercept and a parameter.
#   Want means ordered by level of forestation: Forests, Fragment, Edge, Pastures 
#  
library(mcprofile)
K <- matrix(data = c(1, 1, 0, 0, 0, 0,
                     1, 0, 1, 0, 0, 0,
                     1, 0, 0, 1, 0, 0,
                     1, 0, 0, 0, 0, 0,
                     1, 0, 0, 0, 1, 0,
                     1, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)
K
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    0    0    0    0
## [2,]    1    0    1    0    0    0
## [3,]    1    0    0    1    0    0
## [4,]    1    0    0    0    0    0
## [5,]    1    0    0    0    1    0
## [6,]    1    0    0    0    0    1
# Profile LR
linear.combo <- mcprofile(object = M1, CM = K)
ci.log.mu <- confint(object = linear.combo, level = 0.95, adjust = "none")
mean.LR.ci1 <- data.frame(Loc = pred.data, Estimate = exp(ci.log.mu$estimate), 
                          Lower = exp(ci.log.mu$confint[,1]), Upper = exp(ci.log.mu$confint[,2]))
mean.LR.ci1
##     Loc  Estimate     Lower     Upper
## C1 ForA 119.50000 104.98310 135.29716
## C2 ForB  55.00000  49.27736  61.14941
## C3 Frag  54.33333  48.64676  60.44668
## C4 Edge  48.25000  41.75903  55.38108
## C5 PasA  39.50000  31.41702  48.86316
## C6 PasB  38.00000  32.27465  44.36544
mean.LR.ci1$Loc2 <- factor(mean.LR.ci1$Loc, levels = levels(as.factor(mean.LR.ci1$Loc))[c(2,3,4,1,5,6)])
# StripChart in color
stripchart(Lower ~ Loc2, data = mean.LR.ci1, vertical = FALSE, xlim = c(20,150), 
           col = "red", pch = "(", main = "", xlab = "Bird Count", ylab = "Location")
stripchart(Upper ~ Loc2, data = mean.LR.ci1, vertical = FALSE, col = "red", pch = ")", add = TRUE)
stripchart(Estimate ~ Loc2, data = mean.LR.ci1, vertical = FALSE, col = "red", pch = "+", add = TRUE)
grid(nx = NA, ny = NULL, col = "gray", lty = "dotted")
abline(v = mean(alldata$Birds), col = "darkblue", lwd = 4)


Gráfico acima mostra as MLEs para as contagens médias de aves, denotado por “+” e intervalos de confiança LR, denotados por parênteses, para capturas de aves com redes de neblina em 6 locais no Equador. A linha vertical representa a média geral de todas as contagens.

Este gráfico mostra que há uma tendência geral de aumento nas contagens médias à medida que a quantidade de florestamento aumenta. Também é aparente que um local de floresta tem uma média que é drasticamente diferente dos outros locais. Esperava-se que as florestas tivessem contagens médias mais altas e não se esperava que essas médias fossem necessariamente idênticas em todas as florestas. Isso mostra que a variabilidade de um local para outro pode ser muito grande.

######################################################################
# Alternative parameterization for model to get LR intervals for means
# "-1" in formula removes intercept and causes model to estimate log-means directly 
M1.rp <- glm(formula = Birds~ Loc - 1, family = poisson(link = "log"), data = alldata)
LR.ci <- exp(confint(M1.rp))
mean.LR.ci2 <- data.frame(Estimate = exp(M1.rp$coefficients), Lower = LR.ci[,1], Upper = LR.ci[,2])
mean.LR.ci2
##          Estimate     Lower     Upper
## LocEdge  48.25000  41.75893  55.38074
## LocForA 119.50000 104.98282 135.29653
## LocForB  55.00000  49.27724  61.14921
## LocFrag  54.33333  48.64664  60.44649
## LocPasA  39.50000  31.41768  48.86164
## LocPasB  38.00000  32.27461  44.36506
#####################################################################
# Plots based on Wald Intervals
#Redefining location group factor so that plot does not alphabetize
mean.wald.ci$Loc2 <- factor(mean.wald.ci$Location, levels = c("ForA", "ForB", "Frag", "Edge", "PasA", "PasB"))
mean.wald.ci
##   Location   Mean  Lower  Upper Loc2
## 1     ForA 119.50 105.27 135.65 ForA
## 2     ForB  55.00  49.37  61.27 ForB
## 3     Frag  54.33  48.74  60.56 Frag
## 4     Edge  48.25  41.90  55.56 Edge
## 5     PasA  39.50  31.68  49.25 PasA
## 6     PasB  38.00  32.41  44.55 PasB
# StripChart in color
# NOTE: vertical = FALSE shows intervals horizontally
stripchart(Lower ~ Loc2, data = mean.wald.ci, vertical = FALSE, xlim = c(20,150), 
           col = "red", pch = "(", main = "", xlab = "Bird Count", ylab = "Location")
stripchart(Upper ~ Loc2, data = mean.wald.ci, vertical = FALSE, col = "red", pch = ")", add = TRUE)
stripchart(Mean ~ Loc2, data = mean.wald.ci, vertical = FALSE, col = "red", pch = "+", add = TRUE)
grid(nx = NA, ny = NULL, col = "gray", lty = "dotted")
abline(v = mean(alldata$Birds), col = "darkblue", lwd = 4)


Comparar as duas florestas ou as duas pastagens entre si não é de fato interessante. Em vez disso, queremos combinar as médias para cada par de locais do mesmo habitat antes de compará-los com outros habitats. Seria inapropriado combinar seus dados e estimar uma única média para os dois locais do mesmo habitat, porque os dados de pares de locais não vêm necessariamente da mesma distribuição de Poisson. Em vez disso, estimamos suas médias separadamente e depois combinamos as médias para fazer comparações. Para realizar inferências, precisamos desenvolver uma matriz de coeficientes para as combinações lineares de parâmetros que representam as comparações que estão sendo feitas.

Por exemplo, para comparar os dois habitats florestais com o fragmento, podemos comparar a razão da média geométrica das duas médias da floresta contra a média do fragmento, ou seja, \((\mu_2\mu_3)^{1/2}/\mu_4\). A média geométrica de \(a_1,\cdots,a_n\) é \[ \Big(\prod_{i=1}^n a_i \Big)^{1/n}, \] de modo que o logaritmo da média geométrica seja \[ \dfrac{1}{n}\sum_{i=1}^n \log(a_i)\cdot \] O logaritmo dessa quantidade em termos dos parâmetros de regressão é \((\beta_2 +\beta_3)/2-\beta_4\). Os coeficientes são então (0, 0.5, 0.5, -1, 0, 0). Conjuntos semelhantes de coeficientes comparando cada par de habitats são dados abaixo:

#####################################################################
# Confidence intervals for linear combinations of parameters
#Create coefficients for Lin Combos
contr.mat <- rbind(c(0,.5,.5,0,0,0), c(0,.5,.5,-1,0,0), c(0,.5,.5,0,-.5,-.5), 
                   c(0,0,0,1,0,0), c(0,0,0,1,-.5,-.5), c(0,0,0,0,-.5,-.5))
rownames(contr.mat) <- c("For-Edge", "For-Frag", "For-Past", "Frag-Edge", "Frag-Past", "Edge-Past")
contr.mat
##           [,1] [,2] [,3] [,4] [,5] [,6]
## For-Edge     0  0.5  0.5    0  0.0  0.0
## For-Frag     0  0.5  0.5   -1  0.0  0.0
## For-Past     0  0.5  0.5    0 -0.5 -0.5
## Frag-Edge    0  0.0  0.0    1  0.0  0.0
## Frag-Past    0  0.0  0.0    1 -0.5 -0.5
## Edge-Past    0  0.0  0.0    0 -0.5 -0.5


# Wald inferences using multcomp package
library(multcomp)
loc.test <- glht(model = M1, linfct = contr.mat)
# Defaults use multiplicity adjustment for simultaneous confidence level
summary(loc.test)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = Birds ~ Loc, family = poisson(link = "log"), data = alldata)
## 
## Linear Hypotheses:
##                Estimate Std. Error z value Pr(>|z|)    
## For-Edge == 0   0.51893    0.08358   6.209   <0.001 ***
## For-Frag == 0   0.40019    0.06979   5.734   <0.001 ***
## For-Past == 0   0.73838    0.08132   9.080   <0.001 ***
## Frag-Edge == 0  0.11874    0.09082   1.307    0.553    
## Frag-Past == 0  0.33819    0.08875   3.811   <0.001 ***
## Edge-Past == 0  0.21945    0.09995   2.196    0.122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(loc.test)$confint)
##           Estimate       lwr      upr
## For-Edge  1.680227 1.3564214 2.081332
## For-Frag  1.492103 1.2478380 1.784184
## For-Past  2.092546 1.6990713 2.577142
## Frag-Edge 1.126079 0.8923470 1.421033
## Frag-Past 1.402413 1.1172376 1.760380
## Edge-Past 1.245395 0.9640842 1.608789
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.561489
# Options to get unadjusted (univariate) tests and CIs
summary(loc.test, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = Birds ~ Loc, family = poisson(link = "log"), data = alldata)
## 
## Linear Hypotheses:
##                Estimate Std. Error z value Pr(>|z|)    
## For-Edge == 0   0.51893    0.08358   6.209 5.33e-10 ***
## For-Frag == 0   0.40019    0.06979   5.734 9.81e-09 ***
## For-Past == 0   0.73838    0.08132   9.080  < 2e-16 ***
## Frag-Edge == 0  0.11874    0.09082   1.307 0.191077    
## Frag-Past == 0  0.33819    0.08875   3.811 0.000139 ***
## Edge-Past == 0  0.21945    0.09995   2.196 0.028124 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
exp(confint(loc.test, calpha = qnorm(0.975))$confint)
##           Estimate       lwr      upr
## For-Edge  1.680227 1.4263561 1.979283
## For-Frag  1.492103 1.3013404 1.710831
## For-Past  2.092546 1.7842491 2.454112
## Frag-Edge 1.126079 0.9424543 1.345482
## Frag-Past 1.402413 1.1785034 1.668865
## Edge-Past 1.245395 1.0238272 1.514912
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964


Testes e intervalos de confiança para essas combinações lineares são facilmente obtidos usando mcprofile() para intervalos LR e o pacote multcomp ou mcprofile para Wald. Mostramos o uso de mcprofile() acima. Isso é seguido por summary() para testes e confint() para intervalos de confiança. O padrão é ajustar os resultados para considerar múltiplas inferências, o que garante que a probabilidade de nenhum erro tipo I ou falha de cobertura em toda a família de inferências seja mantida em pelo menos \(1-\alpha\).

Esses padrões podem ser contornados, se desejado, especificando argumentos adicionais conforme mostrado abaixo. No entanto, é considerado uma boa prática apresentar resultados de uma análise ajustados e não ajustados à multiplicidade. Os resultados ajustados são mais conservadores, mas quaisquer achados estatisticamente significativos têm maior probabilidade de resistir a investigações futuras do que aqueles encontrados apenas sem ajustes. Veja Westfall and Young (1993) para uma discussão sobre a importância de taxas de erro simultâneas e coberturas em inferência múltipla.

# Do the same thing with LR in mcprofile
linear.combo <- mcprofile(object = M1, CM = contr.mat)
summary(linear.combo)  
## 
##    mcprofile - Multiple Testing
## 
## Adjustment:   single-step 
## Margin:       0 
## Alternative:  two.sided 
## 
##                Estimate Statistic Pr(>|z|)    
## For-Edge == 0      0.52      6.45   <0.001 ***
## For-Frag == 0      0.40      5.81   <0.001 ***
## For-Past == 0      0.74      9.54   <0.001 ***
## Frag-Edge == 0     0.12      1.31    0.549    
## Frag-Past == 0     0.34      3.86    0.001 ***
## Edge-Past == 0     0.22      2.19    0.123    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(confint(linear.combo)$confint)
##       lower    upper
## 1 1.3608975 2.089856
## 2 1.2492084 1.787121
## 3 1.7045785 2.587743
## 4 0.8940097 1.424983
## 5 1.1191925 1.765186
## 6 0.9633350 1.609484
summary(linear.combo, adjust = "none")
## 
##    mcprofile - Multiple Testing
## 
## Adjustment:   none 
## Margin:       0 
## Alternative:  two.sided 
## 
##                Estimate Statistic Pr(>|z|)    
## For-Edge == 0      0.52      6.45   <2e-16 ***
## For-Frag == 0      0.40      5.81   <2e-16 ***
## For-Past == 0      0.74      9.54   <2e-16 ***
## Frag-Edge == 0     0.12      1.31    0.189    
## Frag-Past == 0     0.34      3.86   <2e-16 ***
## Edge-Past == 0     0.22      2.19    0.028 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(confint(linear.combo, adjust = "none")$confint)
##       lower    upper
## 1 1.4292924 1.983755
## 2 1.3023075 1.712299
## 3 1.7878645 2.459671
## 4 0.9436183 1.347467
## 5 1.1798823 1.671273
## 6 1.0235287 1.515020


A partir dos testes ou intervalos ajustados, os dois primeiros conjuntos de resultados, vemos que em quatro das seis comparações as proporções mostram diferenças que refletem as expectativas dos pesquisadores.

As florestas têm uma média maior do que os habitats de borda (edge), fragmento e pastagem; linhas 1, 2 e 3, respectivamente, com uma proporção especialmente grande em comparação com a pastagem. Essa razão estimada é \(\exp(0.74) = 2.09\), com um intervalo de confiança de 1.71 a 2.59, sugerindo aproximadamente o dobro de pássaros na floresta do que no pasto.

Continuando a examinar outras comparações, vemos que o fragmento tinha mais aves do que os pastos. As razões das médias para fragmento/borda e borda/pasto não são significativamente diferentes de 1; seus logaritmos não são significativamente diferentes de 0, de acordo com os resultados ajustados. Os resultados não ajustados sugerem que a relação borda/pasto pode realmente ser diferente de 1. No entanto, esta conclusão não pode ser feita tão fortemente quanto as outras porque, sem os ajustes de multiplicidade, ela vem de um conjunto de testes cuja chance de fazer qualquer tipo I erros é inflado em relação a \(\alpha\).

As inferências de Wald fornecem intervalos de confiança quase idênticos e a mesma interpretação. Mostramos os resultados de glht() de cálculos multcomp e cálculos manuais no programa correspondente a este exemplo. Para os cálculos manuais, precisamos estimar as combinações lineares dos parâmetros, estimar as variâncias correspondentes dessas combinações lineares e combiná-las da maneira usual para criar pontos finais dos intervalos.

Os cálculos são realizados facilmente usando métodos matriciais. O programa também mostra como usar deltaMethod() do pacote car para estimar a diferença, em vez da razão, entre as médias da floresta e do fragmento.

################################
# Using deltaMethod
library(package = car)
# My original comparison of (mu2 + mu3)/2 - mu4
lin.comb <- deltaMethod(object = M1, g = "exp(b0 + b2)/2 + exp(b0 + b3)/2 - exp(b0 + b4)",
  parameterNames = c("b0", "b2", "b3", "b4", "b5", "b6"))  #Remember there is no beta1
# "For-Frag" comparison as done in the example
lin.comb <- deltaMethod(object = M1, g = "(exp(b0 + b2)*exp(b0 + b3))^0.5 / exp(b0 + b4)",
    parameterNames = c("b0", "b2", "b3", "b4", "b5", "b6"))  #Remember there is no beta1
names(lin.comb)
## [1] "Estimate" "SE"       "2.5 %"    "97.5 %"
lin.comb
##                                                Estimate      SE   2.5 % 97.5 %
## (exp(b0 + b2) * exp(b0 + b3))^0.5/exp(b0 + b4)  1.49210 0.10414 1.28800 1.6962
z0 <- lin.comb$Estimate/lin.comb$SE
pvalue <- 2*(1-pnorm(q = abs(z0)))
data.frame(z0, pvalue)
##         z0 pvalue
## 1 14.32805      0
lin.comb$Estimate + qnorm(p=0.975)*lin.comb$SE
## [1] 1.696211
################################
# Wald CIs by manual computation
# Get out coefficients and variances
beta <- matrix(coef(M1), ncol = 1)
v.beta <- vcov(M1)
# Estimate Lin Combos and standard errors as matrix computations
log.contrasts <- contr.mat %*% beta
SElog.contrasts <- matrix(sqrt(diag(contr.mat %*% v.beta %*% t(contr.mat))), ncol = 1)
# Compute confidence intervals in linear predictor scale
alpha <- 0.05
lower.log <- log.contrasts + qnorm(alpha/2)*SElog.contrasts
upper.log <- log.contrasts + qnorm(1-alpha/2)*SElog.contrasts
# Combine Lin Combo coefficients, estimates of contrasts, and confidence intervals in mean scale
wald.ci <- round(cbind(exp(log.contrasts), exp(lower.log), exp(upper.log)), digits = 2)
colnames(wald.ci) <- c("Estimate", "Lower", "Upper")
wald.ci
##           Estimate Lower Upper
## For-Edge      1.68  1.43  1.98
## For-Frag      1.49  1.30  1.71
## For-Past      2.09  1.78  2.45
## Frag-Edge     1.13  0.94  1.35
## Frag-Past     1.40  1.18  1.67
## Edge-Past     1.25  1.02  1.51




4.2.4 Regressão de Poisson para tabelas de contingência: modelos loglineares


Um uso especial para variáveis explicativas categóricas em um modelo de regressão Poisson é modelar contagens de tabelas de contingência. Praticamente qualquer análise que se possa normalmente conduzir usando métodos tradicionais para tabelas de contingência - como as descritas nas Seções 1.2 e 3.2 - pode ser realizada usando a regressão Poisson. As ferramentas disponíveis na estrutura do modelo linear generalizado para realizar análises adicionais o tornam uma alternativa atraente aos métodos tradicionais baseados em tabelas. Como a regressão Poisson em tabelas de contingência é tão comum, ela tem um nome alternativo - o modelo loglinear - que deriva do uso do ligação logaritmica para a média. Observe que o modelo loglinear pode ser amplamente equivalente aos modelos de regressão multinomial discutidos na Seção 3.3.2. Compararemos as vantagens e desvantagens das duas abordagens no final da próxima seção.

Usar modelos de regressão Poisson para uma estrutura de tabela de contingência bidimensional requer nada mais do que tratar a variável de linha e a variável de coluna como variáveis explicativas categóricas. A regressão Poisson é então aplicada, usando os conjuntos resultantes de variáveis indicadoras juntas em um modelo. Esse processo é muito semelhante à modelagem ANOVA de respostas contínuas. No entanto, algumas explicações sobre os detalhes são necessárias para ver como o modelo reproduz as estatísticas tradicionais.

Seja \(X\) representar a variável de linha com \(I\) níveis e \(Z\) representar a variável de coluna com \(J\) níveis. Usamos as variáveis indicadoras \(x_2,\cdots,x_I\) para \(X\) e \(z_2,\cdots,z_J\) para \(Z\) dentro do nosso modelo de regressão Poisson: \[ \log(\mu)=\beta_0+\beta_2^X x_2+\beta_3^X x_3+\cdots+\beta_I^X x_I+\beta_2^Z z_2+\beta_3^Z z_3 +\cdots+\beta_J^Z z_J, \] onde \(\beta_i^X\), \(i=2,\cdots,I\) e \(\beta_j^Z\), \(j=2,\cdots,J\) são os parâmetros de efeito de regressão comuns com o sobrescrito adicionado por conveniência. Assim, o modelo para a média na célula \((i,j)\) da tabela de contingência pode ser escrito como \[ \log(\mu_{ij})=\beta_0+\beta_i^X +\beta_j^Z, \qquad i=1,\cdots,I, \quad j=1,\cdots,J, \] sendo implícito que \(\beta_1^X=\beta_1^Z=0\).

\[ \mbox{Tabela 4.3: Variáveis indicadoras e logaritmo das médias de um }\\ \mbox{modelo loglinear para uma tabela de contingência } 2\times 2.\\ \begin{array}{cccc|l} \hline X & Z & x_2 & z_2 & \mbox{Logaritmo das médias} \\\hline 1 & 1 & 0 & 0 & \log(\mu_{11})=\beta_0 \\ 1 & 2 & 0 & 1 & \log(\mu_{12})=\beta_0+\beta_2^Z \\ 2 & 1 & 1 & 0 & \log(\mu_{21})=\beta_0+\beta_2^X \\ 2 & 2 & 1 & 1 & \log(\mu_{22})=\beta_0+\beta_2^X+\beta_2^Z \\\hline \end{array} \]

Para ver o que esses parâmetros medem, avaliamos as equações acima para uma tabela de contingência \(2\times 2\) na Tabela 4.3. Vemos que \[ \log(\mu_{21})-\log(\mu_{11}) = \log(\mu_{22})-\log(\mu_{12})=\beta_2^X, \] então \(\beta_2^X\) mede a diferença em média logarítmica entre as duas contagens na mesma coluna das linhas 2 e 1.

Assim, \(\exp(\beta_2^X)\) mede a razão entre quaisquer duas dessas médias. De forma similar, \[ \log(\mu_{12}) - \log(\mu_{11}) = \log(\mu_{22}) - \log(\mu_{21}) = \beta_2^Z, \] então \(\beta_2^Z\) mede a diferença em média logarítmica entre as duas contagens na mesma linha das colunas 2 e 1 e \(\exp(\beta_2^Z)\) mede a razão dessas médias.

Em tabelas maiores, \(\beta^X_i\) mede a diferença em média logaritmica, ou seja, o logaritmo da razão de médias entre as contagens nas linhas \(i\) e 1 na mesma coluna. Uma interpretação análoga vale para \(\beta^Z_j\). Portanto, pode-se mostrar que esses parâmetros também medem diferenças em média logaritmica ou log-ratios de médias, nas margens da tabela: \[ \beta^X_i = \log\big(\mu_{i+}/\mu_{1+}\big)\cdot \]

Não há parâmetros nas equações anteriores que permitam que células individuais em uma tabela \(I\times J\) se desviem do que determinam as margens de linha e coluna. Isso corresponde à independência entre \(X\) e \(Z\), ver Seção 3.2. Outra maneira de ver isso é considerar a razão de chances entre quaisquer duas linhas e colunas.

Por exemplo, suponha que temos uma tabela \(2\times 2\) com probabilidades de célula \(\pi_{ij}\), \(i = 1,2\); \(j = 1,2\), como na Seção 3.2.1. Usando a equivalência entre os modelos Poisson e multinomial discutidos na Seção 4.2.5, podemos denotar a média da contagem do modelo de regressão Poisson na célula \((i,j)\) por \(\mu_{ij} = n\pi_{ij}\). É então razoavelmente fácil mostrar que \[ OR = \dfrac{\mu_{11}\mu_{22}}{\mu_{21}\mu_{12}}\cdot \]

As esquações anteriores implicam que \[ \begin{array}{rcl} \log(OR) & = & \log(\mu_{11})+\log(\mu_{22})-\log(\mu_{21})-\log(\mu_{12}) \\ & = & \beta_0+\big(\beta_0+\beta_2^X+\beta_2^Z\big)-\big(\beta_0+\beta_2^X\big)-\big(\beta_0+\beta_2^Z\big)=0\cdot \end{array} \]

Assim, \(OR = 1\). A mesma coisa ocorre em modelos para tabelas maiores independentemente de quais duas linhas ou colunas são escolhidas. Ou seja, o modelo loglinear representa a independência entre \(X\) e \(Z\). Os parâmetros desse modelo servem apenas para fazer com que a tabela de contagens previstas tenha as mesmas distribuições marginais da tabela de contagens observadas.

Na verdade, as contagens previstas são exatamente as mesmas que as contagens esperadas que calculamos na Seção 3.2 para testar a independência em uma tabela de contingência. Como a inferência nas distribuições marginais normalmente não é nossa principal preocupação em uma tabela de contingência, precisamos aumentar o modelo loglinear para permitir a associação entre \(X\) e \(Z\).

Isso pode ser feito adicionando ao modelo termos de produtos cruzados na forma \(\beta^{XZ}_{ij}\), \(i = 2,\cdots,I\), \(j = 2,\cdots,J\), para criar uma interação entre \(X\) e \(Z\). Como existem \(I-1\) indicadores \(x_i\) e \(J-1\) indicadores \(z_j\), existem \((I-1)(J-1)\) novos termos a serem adicionados ao modelo. Isso cria um modelo saturado, porque o número total de parâmetros é igual ao número de médias de células que estão sendo estimadas.

Para ver isso, observe primeiro que existem \(IJ\) contagens amostrais independentes na tabela. O intercepto \(\beta_0\) representa 1 grau de liberdade (df). Os parâmetros de efeito para \(X\), \(\beta^X_2,\cdots,\beta^X_I\), respondem por \(I-1\) df, enquanto os de \(Z\) representam um adicional de \(J-1\). Com os termos de interação, o novo modelo consome \[ 1 + (I-1) + (J-1) + (I-1)(J-1) = IJ, \] graus de liberdade (df).

O novo modelo pode ser escrito como \[ \log(\mu_{ij})=\beta_0+\beta_i^X+\beta_j^Z+\beta_{j}^{XZ}\cdot \] Os parâmetros de efeito \(\beta^X_i\) e \(\beta^Z_j\) operam novamente nas contagens marginais, mas agora os \(XZ\) parâmetros de interação \(\beta^{XZ}_{ij}\), \(i = 1,\cdots,I\), \(j=1,\cdots,J\) permitem que as proporções de médias entre as células em duas linhas mudem nas colunas e vice-versa. Para ver isso, considere novamente a tabela \(2\times 2\) e lembre-se de que definimos \(\beta^{XZ}_{ij}= 0\) quando \(i = 1\) ou \(j = 1\). Então \[ \begin{array}{rcl} \log(OR) & = & \log(\mu_{11})+\log(\mu_{22})-\log(\mu_{21})-\log(\mu_{12}) \\ & = & \big(\beta_0+\beta_1^X+\beta_1^Z+\beta_{11}^{XZ}\big)+\big(\beta_0+\beta_2^X+\beta_2^Z+\beta_{22}^{XZ}\big) \\ & & \qquad \qquad \qquad - \big(\beta_0+\beta_2^X+\beta_1^Z+\beta_{21}^{XZ}\big) - \big(\beta_0+\beta_1^X+\beta_2^Z+\beta_{12}^{XZ}\big)\\ & = & \beta_{11}^{XZ} + \beta_{22}^{XZ} - \beta_{21}^{XZ} - \beta_{12}^{XZ} \, = \, \beta_{22}^{XZ} \end{array} \]

Assim, a razão de chances em uma tabela \(2\times 2\) é apenas \(\exp\big(\beta_{XZ}^{22}\big)\). De forma mais geral, a razão de chances entre as linhas \(i\) e \(i'\) e as colunas \(j\) e \(j'\) é \[ OR_{i \, i', j \, j'} = \exp\big(\beta_{i \, j'}^{XZ} + \beta_{i' \, j'}^{XZ}-\beta_{i'\, j}^{XZ}-\beta^{XZ}_{i\, j'}\big), \] onde os subscritos extras nos lembram que há mais de uma razão de chances possível.

Como o modelo loglinear com interações é um modelo saturado, suas médias estimadas \(\widehat{\mu}_{ij}\) correspondem exatamente às contagens observadas \(y_{ij}\). Portanto, as razões de chance estimadas acima são exatamente as mesmas baseadas nas contagens correspondentes na tabela observada.

Os testes de independência em tabelas de contingência descritos na Seção 3.2.3 podem ser replicados exatamente a partir do modelo de regressão Poisson por meio de testes de comparação de modelos. A hipótese nula de independência implica no modelo loglinear sem interações, enquanto a alternativa implica o modelo irrestrito com interações.

Ajustar ambos os modelos às contagens tabuladas \(y_{ij}\), \(i = 1,\cdots,I\), \(j = 1,\cdots,J\) e comparar desvios residuais fornece a estatística LRT. Nesse caso, a estatística também é apenas o deviance residual do modelo nulo, pois o modelo alternativo está saturado, então seu deviance residual é zero. Essa estatística pode ser exatamente igual à estatística LRT apresentada na Seção 3.2 para testar a independência.

Os graus de liberdade para a distribuição \(\chi^2\) amostral \((I-1)(J-1)\), vêm do número de parâmetros adicionais no modelo completo em comparação com o modelo nulo, sem contar aqueles definidos como zero. As estimativas da razão de chances e os intervalos de confiança são encontrados exponenciando estimativas e intervalos de confiança para combinações lineares dos parâmetros das interação \(XZ\) conforme determinado pela acima.

Deve-se notar que a estatística de Pearson apresentada na Seção 3.2.3 é geralmente preferida à estatística LR para testar a independência em uma tabela \(I\times j\) de contingência, veja também a discussão em 1.2.3. No entanto, a estatística LR é mais fácil de adaptar para testes em outros modelos, como discutimos nas Seções 4.2.5 e 4.2.6.


Exemplo 4.7: Arremesso de lance livre de Larry Bird.


Aqui, usamos um modelo loglinear para reanalisar os dados do arremesso de lance livre de Larry Bird que foi usado anteriormente no Capítulo 1. Mostramos que a estimativa da razão de chances e o intervalo de confiança de Wald são exatamente os mesmos fornecidos na Seção 1.2.5.

Além disso, fornecemos um intervalo de confiança LR para a razão de chances. Existem duas variáveis, First e Second, cada uma com níveis “made” e “missed”. Para este exemplo, usamos o objeto c.table criado no Capítulo 1 e convertemos o objeto em um arquivo de dados:

c.table <- array ( data = c(251 , 48, 34, 5) , dim = c(2 ,2) , 
                    dimnames = list ( First = c("made", "missed") , Second = c("made", "missed")))
all.data <- as.data.frame (as.table (c.table ))
all.data
##    First Second Freq
## 1   made   made  251
## 2 missed   made   48
## 3   made missed   34
## 4 missed missed    5


As contagens estão contidas na variável Freq. Estimamos o modelo saturado \[ \log(\mu_{ij}) = \beta_0 + \beta^F_i + \beta^S_j + \beta^{FS}_{ij}, \] onde \(\mu_{ij}\) é a contagem média na célula \((i, j)\), \(i = 1,2\), \(j =1,2\) e o índice sobrescrito \(F\) significa “First” e \(S\) para “Second”. O modelo é estimado usando o valor do argumento formula como Freq ~ First*Second em glm().

É importante observar no contexto da Tabela 4.3 que o nível 1 para ambas as variáveis é representado por “made” e o nível 2 por “missed”. Assim, \(\beta_0\) é a contagem média logarítmica de sequências “made”, \(\beta^F_2\) é a razão logarítmica das contagens médias para a primeira perdida para a primeira feita quando a segunda é feita e \(\beta^S_2\) é a razão logarítmica das contagens médias para o segundo perdido para segundo feito quando o primeiro é feito. O parâmetro de interação \(\beta^{FS}_{22}\) é o logaritmo da razão das chances de perder o segundo dado que o primeiro foi perdido, para as chances de perder o segundo dado que o primeiro foi feito. Isso é equivalente ao logaritmo da razão das chances de fazer o segundo dado que o primeiro foi feito, para as chances de fazer o segundo dado que o primeiro foi perdido. O resumo do ajuste do modelo está abaixo:

M1 <- glm ( formula = Freq ~ First *Second , family = poisson ( link = "log") , data = all.data )
summary (M1)
## 
## Call:
## glm(formula = Freq ~ First * Second, family = poisson(link = "log"), 
##     data = all.data)
## 
## Deviance Residuals: 
## [1]  0  0  0  0
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               5.52545    0.06312  87.540   <2e-16 ***
## Firstmissed              -1.65425    0.15754 -10.501   <2e-16 ***
## Secondmissed             -1.99909    0.18275 -10.939   <2e-16 ***
## Firstmissed:Secondmissed -0.26267    0.50421  -0.521    0.602    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4.0206e+02  on 3  degrees of freedom
## Residual deviance: 3.8636e-14  on 0  degrees of freedom
## AIC: 29.926
## 
## Number of Fisher Scoring iterations: 3


Os parâmetros estimados são \(\widehat{\beta}_0 = 5.53\), \(\widehat{\beta}^F_2 = -1.65\), \(\widehat{\beta}^S_2 = -2.00\) e \(\widehat{\beta}^{FS}_{22} = -0.26\). A Tabela 4.3 mostra como estimar cada contagem média de acordo com os parâmetros do modelo. Como o modelo está saturado, eles correspondem exatamente às contagens observadas. Por exemplo, \(\widehat{\mu}_{12} = \exp(5.53 - 2.00) = 34\) é exatamente o mesmo que a contagem observada para o primeiro resultado, segundo omisso (Secondmissed).

Seguimos com as funções de análise usuais para estimar a razão de chances e encontrar intervalos de confiança:

round (exp( M1$coefficients [4]) , digits = 2)
## Firstmissed:Secondmissed 
##                     0.77
inter <- confint (M1 , parm = "Firstmissed:Secondmissed")
round (exp( inter ), digits = 2)
##  2.5 % 97.5 % 
##   0.25   1.91
round (exp( confint.default (M1 , parm = "Firstmissed:Secondmissed")), digits = 2)
##                          2.5 % 97.5 %
## Firstmissed:Secondmissed  0.29   2.07


A razão de chances estimada é \(\widehat{OR} = 0.77\), portanto, há uma probabilidade estimada menor de errar o segundo arremesso após um primeiro arremesso errado do que após um primeiro arremesso feito. Da mesma forma, a probabilidade de sucesso é maior no segundo tiro após um primeiro tiro errado.

O intervalo de confiança de LR perfilada obrido de exp(confint()) é \(0.25 < OR < 1.91\), portanto, com 95% de confiança, as chances de sucesso no segundo tiro estão entre 1/4 a quase duas vezes maiores após um primeiro tiro de sucesso em comparação com um erro primeiro. Para fins de comparação, o intervalo de confiança de Wald calculado por confint.default() é \(0.29 < OR < 2.07\), o que corresponde ao que foi fornecido na Seção 1.2.5.

Finalmente, o LRT para \(H_0 \, : \, \beta^{FS}_{22} = 0\), \(H_1 \, : \, \beta^{FS}_{22} \neq 0\), que é equivalente a \(H_0 \, : \, OR = 1\), \(H_1 \, : \, OR \neq 1\) é fornecido usando a função Anova():

library (car)
Anova (M1)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Freq
##              LR Chisq Df Pr(>Chisq)    
## First         174.958  1     <2e-16 ***
## Second        226.812  1     <2e-16 ***
## First:Second    0.286  1      0.593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


A estatística do teste é \(-2\log(\Lambda ) = 0.29\) com um \(p\)-valor de 0.59, indicando que não há evidências suficientes para concluir que as chances e, portanto, a probabilidade de sucesso para o segundo lance livre dependem do resultado do primeiro lance livre. Observe que esse \(p\)-valor é exatamente o mesmo do LRT realizado na Seção 1.2.3.



4.2.5 Grandes modelos loglineares


Os modelos loglineares podem ser estendidos para modelar contagens classificadas de forma cruzada por muitas variáveis categóricas arbitrárias. Uma tabela de contingência de contagens correspondente à classificação cruzada de \(p\) variáveis categóricas é chamada de tabela \(p\)-way. Por exemplo, considere a tabela de contingência de três vias que consiste em três variáveis categóricas, \(X\), \(Z\) e \(W\) com níveis \(I\), \(J\) e \(K\), respectivamente. Podemos estender nosso modelo loglinear adicionando um subscrito para a variável \(W\) às médias e contagens observadas: \(y_{ijk}\) tem média \(\mu_{ijk}\), \(i = 1,\cdots,I\); \(j= 1,\cdots,J\) e \(k = 1,\cdots, K\).

O modelo de independência é expandido para incluir variáveis indicadoras para \(W\), \(w_2,w_3,\cdots, w_K\) definidas analogamente às de \(X\) e \(Z\). Seus correspondentes parâmetros de regressão, \(\beta_2^W,\cdots,\beta_K^W\), afetam as contagens marginais de uma tabela de três vias ajustada, por exemplo, \(K\) tabelas de duas vias diferentes resumindo as contagens para \(X\) e \(Z\) em cada valor de \(W\), de modo que as razões das médias dos níveis \(k\) e \(k'\) de \(W\) são constantes.

Assim, \(\log\big(\mu_{ijk}/\mu_{ijk'}\big) = \exp\big(\beta_k^W-\beta_{k'}^W\big)\), em todos os níveis de \(X\) e \(Z\). Esse modelo também representa independência entre todas as três variáveis, porque qualquer razão de chance entre dois níveis de uma variável e dois níveis de outra variável é 1. Por causa disso, às vezes é chamado de modelo de independência mútua.

Modelos que permitem associações entre as três variáveis são desenvolvidos dependendo de quais associações são de interesse ou esperadas. Alguns deles estão resumidos na tabela a seguir. Existem agora três interações emparelhadas a serem consideradas: \(XZ\), \(XW\) e \(ZW\). Cada uma é formada pela adição de produtos das variáveis indicadoras ao modelo de regressão e cada conjunto de parâmetros de interação permite que as razões de chances para o par correspondente de variáveis sejam algo diferente de 1. Se todas as três interações pareadas forem incluídas no modelo, temos \[ \log(\mu_{ijk}) = \beta_0+\beta_i^X + \beta_j^Z + \beta_k^W + \beta_{ij}^{XZ}+\beta_{ik}^{XW}+\beta_{jk}^{ZW}\cdot \]


\[ \mbox{Modelos comuns e suas razões de chances associadas para um modelo loglinear } \\ \mbox{de três vias com variáveis } X, Z, W.\\ \begin{array}{ccc} \hline \mbox{Modelo} & \mbox{Notação} & \mbox{Odds Ratio de } XZ \mbox{ no nível } k \mbox{ de } W \\\hline \mbox{Independência Mútua} & X, Z, W & 1 \mbox{ para todo } k \\\hline \mbox{Homogeneous Association} & XZ, XW, ZW & \exp\big(\beta_{ij}^{XZ}+\beta_{i'j'}^{XZ}-\beta_{i'j}^{XZ}-\beta_{ij'}^{XZ}\big) \mbox{ para todo } k \\\hline \mbox{Saturada} & XZW & \big(\beta_{ij}^{XZ}+\beta_{i'j'}^{XZ}-\beta_{i'j}^{XZ}-\beta_{ij'}^{XZ}\big) \\ \mbox{(Associação Heterogênea)} & & +\big(\beta_{ijk}^{XZW}+\beta_{i'j'k}^{XZW}-\beta_{i'jk}^{XZW}-\beta_{ij'k}^{XZW}\big) \\\hline \end{array} \]


Para ver as razões de chances (odds ratio) correspondentes para este modelo, considere formar \(OR\) para os níveis \(i\), \(i'\) de \(X\) e \(j\), \(j'\) de \(Z\) no nível \(k\) de \(W\): \[ OR_{ii',jj'(k)}=\dfrac{\mu_{ijk}\mu_{i'j'k}}{\mu_{i'jk}\mu_{ij'k}}\cdot \]

Aplicando a expressão de \(\log(\mu_{ijk})\) acima temos \[ \begin{array}{rcl} \log\big(OR_{ii',jj'(k)}\big) & = & \log(\mu_{ijk}) + \log(\mu_{i'j'k})-\log(\mu_{i'jk})-\log(\mu_{ij'k}) \\ & = & \beta_{ij}^{XZ}+\beta_{i'j'}^{XZ}-\beta_{i'j}^{XZ}-\beta_{ij'}^{XZ}\cdot \end{array} \]

Todos os outros parâmetros se cancelam, exceto aqueles envolvidos na interação a partir da qual a razão de chances de interesse é calculada. O mesmo acontece com qualquer par de variáveis categóricas em quaisquer pares de níveis. Portanto, esse modelo permite que todas as razões de chances correspondentes a qualquer par de variáveis variem livremente entre os níveis dessas duas variáveis. Observe, no entanto, que uma determinada razão de chances é forçada a ser constante entre os níveis do terceiro fator. Por exemplo, a razão de chances acima é \[ \exp\big(\beta_{ij}^{XZ}+\beta_{i'j'}^{XZ}-\beta_{i'j}^{XZ}-\beta_{ij'}^{XZ}\big), \] independentemente do nível \(k\) de \(W\).

Assim, esse modelo é chamado de modelo de associação homogênea. Podem ser consideradas versões deste modelo em que apenas uma ou duas das interações estão presentes, permitindo modelar associações entre alguns pares de variáveis e assumindo independência entre outras.

Em todos os casos, testes e estimativas são obtidos usando as mesmas ferramentas de antes. Finalmente, interações de ordem superior são possíveis. A interação de três variáveis permite que as associações entre qualquer par de variáveis variem entre os níveis da terceira. Termos os termos \(\beta_{ijk}^{XZW}\) adicionados e esses termos entram em qualquer razão de chances. Por exemplo, a razão de chances para os níveis \(i\) e \(i'\) de \(X\) e os níveis \(j\) e \(j'\) de \(Z\) no nível \(k\) de \(W\) torna-se \[ \begin{array}{rcl} \log\big(OR_{ii',jj'(k)}\big) & = & \log(\mu_{ijk}) + \log(\mu_{i'j'k})-\log(\mu_{i'jk})-\log(\mu_{ij'k}) \\ & = & \beta_{ij}^{XZ}+\beta_{i'j'}^{XZ}-\beta_{i'j}^{XZ}-\beta_{ij'}^{XZ}+\big(\beta_{ijk}^{XZW}+\beta_{i'j'k}^{XZW}-\beta_{i'jk}^{XZW}-\beta_{ij'k}^{XZW}\big), \end{array} \] o que mostra que as razões de chances \(XZ\) têm o potencial de variar dependendo do nível \(k\) de \(W\). Embora essa fórmula possa parecer complicada, é previsível e não é difícil de estimar em R.

Em modelos com mais variáveis, termos de ordem ainda superior são tecnicamente possíveis. Conjuntos de parâmetros de interação, indexados pelos nomes e níveis das variáveis, podem ser criados para corresponder a qualquer combinação de variáveis. No entanto, as interações de ordem superior são difíceis de interpretar. Uma interpretação geral pode ser desenvolvida recursivamente da seguinte forma: uma interação de 3 variáveis permite que interações de 2 variáveis, ou seja, razões de chances variem dependendo do nível da terceira variável; uma interação de 4 variáveis permite que as interações de 3 variáveis variem dependendo do nível da quarta variável; e assim por diante.

Compreender o que essas declarações implicam no contexto de um problema específico é uma questão diferente. A tabulação cuidadosa ou a representação gráfica das razões de chances entre os níveis de outras variáveis pode ajudar a extrair o significado subjacente de uma interação de ordem superior.

Freqüentemente, o objetivo de uma análise é determinar quais associações, se houver, são importantes e subseqüentemente estimar essas associações. Geralmente é aceito que os termos de efeito principal que controlam as contagens marginais devem permanecer no modelo quando termos de interação de 2 variáveis e de ordem superior são o foco do estudo. De fato, deixar os efeitos principais fora do modelo em glm() quando a interação correspondente permanece no modelo resulta em uma interpretação completamente diferente do termo de interação.

Para ver isso, tente executar \[ \mbox{M2 = glm(formula = Freq ~ First + First:Second, family = poisson("log"), data = all.data) } \] com o exemplo dados dados de Larry Bird e examine os resultados. Em particular, observe summary(M2), model.matrix(M2) e predict(M2).

Além disso, o princípio da marginalidade determina que todos os termos de ordem inferior que são subconjuntos de quaisquer termos de ordem superior também devem permanecer no modelo. Por exemplo, a interação \(XZW\) contém \(XZ\), \(XW\) e \(ZW\), bem como \(X\), \(Z\) e \(W\). Se um modelo contém \(XZW\), ele também contém esses outros seis termos. Seguindo essa convenção, podemos usar a notação abreviada que identifica qualquer modelo listando apenas os termos de ordem mais alta no modelo. Nosso exemplo atual é denotado como “Modelo \(XZW\)”. Isso também corresponde aos termos que precisam ser especificados no argumento formula de glm() para ajustar o modelo. Por exemplo, formula = count ~ X * Z * W caberia neste modelo, onde count denota uma contagem de células. Como outro exemplo, “Modelo \(X, ZW\)” consistiria em termos para \(X\), \(Z\), \(W\) e \(ZW\) , especificando a independência entre \(X\) e \(W\) e entre \(X\) e \(Z\), mas permitindo a associação entre \(Z\) e \(W\) que é homogênea entre os níveis de \(X\). Este modelo pode ser ajustado usando formula = count ~ X + Z * W.

Na maioria dos casos, a natureza exata da estrutura de associação não é conhecida com antecedência, portanto, os termos do modelo para associações precisam ser selecionados por meio de análise estatística. Como selecionar esses termos é objeto de algum debate. Procedimentos análogos à seleção direta ou métodos de eliminação regressiva que já foram populares na análise de regressão (Kutner et al., 2004) às vezes são recomendados (por exemplo, Agresti, 2007). Embora os testes de significância possam nos dizer se uma hipótese específica deve ou não ser rejeitada o uso múltiplo e sequencial de testes de significância para selecionar ou eliminar associações repetidamente de acordo com seu tamanho estimado viola os fundamentos de probabilidade por trás da definição de erro tipo I (Burnham and Anderson, 2002). Este procedimento deve, portanto, ser usado como último recurso quando não houver métodos melhores disponíveis e o nível de \(\alpha\) usado não precisa ser o tradicional 0.05. Discutimos métodos de seleção de modelos no Capítulo 5.


Exemplo 4.8: Ideologia política.


Agresti (2007) apresenta um resumo do US General Social Survey de 1991 classificando pessoas de acordo com:

  1. ideologia política: uma variável de 5 níveis com níveis 1: “Muito Liberal” (“VL”), 2: “Ligeiramente Liberal” (“SL”), 3: “Moderado” (“M”), 4: “Ligeiramente Conservador” (“SC”), e 5: “Muito Conservador” (“VC”)

  2. partido político: uma variável de 2 níveis com níveis “Republicano” (“R”) e “Democrata” (“D”)

  3. gênero: uma variável de 2 níveis com níveis “Feminino” (“F”) e “Masculino” (“M”).

A fonte desses dados é uma pesquisa em larga escala com um desenho complexo de amostragem baseado em probabilidade. Deixar de levar em conta o projeto de pesquisa na análise pode levar a estimativas de parâmetros enviesadas e deturpação das verdadeiras relações entre as variáveis. As análises que apresentamos aqui são usadas apenas como ilustração das técnicas de análise atuais e como comparação com trabalhos publicados anteriormente.

Ajustamos a esses dados vários modelos que diferem na interpretação de quais variáveis estão associadas e se as associações são homogêneas. Observe que a variável ideologia é ordinal. Na próxima seção, exploramos como essa estrutura pode ser explorada para fornecer uma descrição mais parcimoniosa da associação. Por enquanto, tratamos a ideologia como nominal. Abaixo segue parte dos dados:

alldata <- read.csv ( file = "http://leg.ufpr.br/~lucambio/ADC/PolIdeolData.csv" )
head ( alldata )
##   gender party ideol count
## 1      F     D    VL    44
## 2      F     D    SL    47
## 3      F     D     M   118
## 4      F     D    SC    23
## 5      F     D    VC    32
## 6      F     R    VL    18


Os modelos que consideramos são todos regressões Poisson com as variáveis gender \((G)\), party \((P)\) e ideol \((I)\) presentes em todos os modelos como efeitos principais. Os modelos diferem em como as interações são tratadas. O modelo saturado \((GPI)\) é \[ \log(\mu_{ijk})=\beta_0+\beta_i^G+\beta_j^P +\beta_k^I + \beta_{ij}^{GP}+ \beta_{ik}^{GI}+ \beta_{jk}^{PI}+\beta_{ijk}^{GPI}, \] onde \(i=F,M\), \(j=D,R\) e \(k=1,\cdots,5\). O modelo de associação homogênea \((GP, GI, PI)\) é dado por \[ \log(\mu_{ijk})=\beta_0+\beta_i^G+\beta_j^P +\beta_k^I + \beta_{ij}^{GP}+ \beta_{ik}^{GI}+ \beta_{jk}^{PI}\cdot \]

Consideramos também o modelo que assume independência mútua entre todas as variáveis \((G, P, I)\), \[ \log(\mu_{ijk})=\beta_0+\beta_i^G+\beta_j^P +\beta_k^I \cdot \]

Considera-se mais um modelo intermediário aos dois últimos. Temos quase certeza de que uma associação de \(IP\) deve existir porque ideologias diferentes fazem parte do que define as duas partes. A verdadeira questão é se as outras duas associações, \(GP\) e \(GI\), existem. Portanto, esperamos o modelo \(G, P, I\) para fornecer um ajuste ruim e, portanto, consideramos um aumento dele para permitir uma associação para \(PI\), modelo \(G, PI\).

Comparando este modelo com o modelo \(GP, GI, PI\) permitirá testar se as outras duas associações, cuja importância é incerta, podem ser retiradas simultaneamente do modelo.

Ajustamos todos os modelos usando glm() e comparamos modelos usando anova():

# Saturated Model : GPI
mod.sat <- glm( formula = count ~ gender * party *ideol , family = poisson ( link = "log"), 
                data = alldata )
# Homogeneous association model in all 3 associations : GP ,GI ,PI
mod.homo <- glm( formula = count ~ ( gender + party + ideol )^2, family = poisson ( link = "log"), 
                 data = alldata )
anova (mod.homo , mod.sat , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: count ~ (gender + party + ideol)^2
## Model 2: count ~ gender * party * ideol
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         4     3.2454                     
## 2         0     0.0000  4   3.2454   0.5176
# Model assuming only PI association : G,PI
mod.homo.PI <- glm( formula = count ~ gender + party *ideol , family = poisson ( link = "log"), 
                    data = alldata )
anova (mod.homo.PI , mod.homo , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: count ~ gender + party * ideol
## Model 2: count ~ (gender + party + ideol)^2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1         9    17.5186                       
## 2         4     3.2454  5   14.273  0.01396 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model assuming mutual independence : G,P,I
mod.indep <- glm( formula = count ~ gender + party + ideol , family = poisson ( link = "log"), 
                  data = alldata )
anova (mod.indep , mod.homo.PI , test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: count ~ gender + party + ideol
## Model 2: count ~ gender + party * ideol
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        13     79.851                          
## 2         9     17.519  4   62.333 9.376e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


As comparações de modelos estão resumidas na Tabela 4.5. Por exemplo, a comparação dos dois últimos modelos testa \(H_0 \, : \, \beta^{PI}_{jk} = 0\), \(j = D,R\), \(k = 1,\cdots,5\) contra \(H_1 \, : \, \mbox{Nem todos os } \beta^{PI}_{jk} = 0\). O modelo \(GP, GI, PI\) parece ajustar-se razoavelmente bem em relação ao modelo saturado, \(-2\log(\Lambda ) = 3.24\) e \(p\)-valor = 0.52, indicando que todas as associações de duas variáveis parecem ser razoavelmente consistentes entre os níveis da terceira variável.

As associações \(GP\) e \(GI\) parecem contribuir com algo para o ajuste do modelo, conforme indicado pelo ajuste um tanto ruim de \(G, PI\) em relação ao \(GP, GI, PI\), \(-2\log(\Lambda ) = 14.27\) e \(p\)-valor=0.01. Claramente, o modelo de independência \(G, P, I\) não se ajusta bem, como esperado \(-2\log(\Lambda) = 62.33\) e \(p\)-valor \(\approx\) 0, mostrando a importância da associação \(PI\). Diante desses resultados, consideramos \(GP, GI, PI\) como o modelo de trabalho a partir do qual examinar essas associações mais a fundo.

\[ \mbox{Tabela 4.5: Análise de deviance da associação partido-ideologia. } \\ \mbox{O teste em cada linha usa o modelo dessa linha como hipótese nula} \\ \mbox{e o modelo diretamente acima dela como alternativa. } \\ \begin{array}{cccccc} \hline & \mbox{Deviance} & \mbox{df} & \mbox{Estatística LRT} & \mbox{df} & \\ \mbox{Modelo} & \mbox{residual} & \mbox{residual} & -2\log(\Lambda) & \mbox{LRT} & p-\mbox{valor} \\\hline \mbox{GPI} & 0 & 0 & & & \\ \mbox{GP, GI, PI} & 3.25 & 4 & 3.25 & 4 & 0.52 \\ \mbox{(G, PI} & 17.52 & 9 & 14.27 & 5 & 0.01 \\ \mbox{G, P, I} & 79.85 & 13 & 62.33 & 4 & \approx 0 \\\hline \end{array} \]

Nosso próximo passo é obter um resumo das estimativas dos parâmetros do modelo para que possamos ver a ordem em que estão listados. Embora possamos prever essa ordenação com base no conhecimento dos níveis dos fatores e da ordem em que aparecem no modelo ajustado, examiná-los pode ajudar a formar as combinações lineares necessárias para a estimativa da razão de chances.

Também obtemos testes de cada termo no modelo usando Anova():

# Summary and LR tests of homogeneous association model 
round ( summary (mod.homo ) $coefficients , digits = 3)
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)        4.740      0.089  53.345    0.000
## genderM           -0.704      0.137  -5.149    0.000
## partyR            -0.244      0.125  -1.962    0.050
## ideolSC           -1.638      0.197  -8.302    0.000
## ideolSL           -0.832      0.157  -5.308    0.000
## ideolVC           -1.306      0.176  -7.438    0.000
## ideolVL           -0.899      0.162  -5.534    0.000
## genderM:partyR     0.276      0.147   1.878    0.060
## genderM:ideolSC    0.534      0.216   2.476    0.013
## genderM:ideolSL    0.236      0.216   1.093    0.274
## genderM:ideolVC    0.448      0.201   2.230    0.026
## genderM:ideolVL    0.372      0.227   1.636    0.102
## partyR:ideolSC     0.826      0.222   3.716    0.000
## partyR:ideolSL    -0.437      0.217  -2.015    0.044
## partyR:ideolVC     0.702      0.203   3.458    0.001
## partyR:ideolVL    -0.861      0.243  -3.548    0.000
library (car)
Anova (mod.homo )
## Analysis of Deviance Table (Type II tests)
## 
## Response: count
##              LR Chisq Df Pr(>Chisq)    
## gender         20.637  1  5.551e-06 ***
## party           0.528  1    0.46736    
## ideol         154.131  4  < 2.2e-16 ***
## gender:party    3.530  1    0.06026 .  
## gender:ideol    8.965  4    0.06198 .  
## party:ideol    60.555  4  2.218e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Additional sub-models of homogeneous association
mod.homo.PIGI <- glm(formula = count ~ gender*ideol + party*ideol, 
                     family = poisson(link = "log"), data = alldata)
Anova(mod.homo.PIGI)
## Analysis of Deviance Table (Type II tests)
## 
## Response: count
##              LR Chisq Df Pr(>Chisq)    
## gender         20.637  1  5.551e-06 ***
## ideol         154.131  4  < 2.2e-16 ***
## party           0.528  1    0.46736    
## gender:ideol   10.743  4    0.02961 *  
## ideol:party    62.333  4  9.376e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.homo.PIGP <- glm(formula = count ~ gender*party + party*ideol, 
                     family = poisson(link = "log"), data = alldata)
Anova(mod.homo.PIGP)
## Analysis of Deviance Table (Type II tests)
## 
## Response: count
##              LR Chisq Df Pr(>Chisq)    
## gender         20.637  1  5.551e-06 ***
## party           0.528  1    0.46736    
## ideol         154.131  4  < 2.2e-16 ***
## gender:party    5.308  1    0.02123 *  
## party:ideol    62.333  4  9.376e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Observe que os níveis da variável ideol foram ordenados alfabeticamente, de modo que M (Moderado) seja a linha de base e todos os parâmetros ideol representem comparações a esse nível. Isso pode ser alterado usando a função relevel() ou factor() conforme descrito na Seção 2.2.6, se desejado.

Os LRTs ou testes de razão de verossimilhanças para os termos do modelo mostram que a associação \(PI\) é altamente significativa, mas há apenas evidência moderada para a existência das outras duas, ambos os \(p\)-valores são 0.06. Acontece que ambos os termos são significativos em \(\alpha= 0.05\) se o outro termo for primeiro removido do modelo. Este é um exemplo da dificuldade em usar testes de hipóteses como uma ferramenta de seleção de modelos. Consulte o Capítulo 5 para métodos alternativos de seleção de modelos.

Com a ordem das estimativas de parâmetros agora clara, podemos usá-las para estimar as razões de chances para as três associações, \(GP\), \(GI\) e \(PI\). Como as associações são homogêneas, não precisamos considerar uma terceira variável ao estimar uma razão de chances entre as outras duas. Para \(GP\), construímos a razão das chances de ser republicano para homens versus mulheres: \[ \log(OR) \, = \, \log\left( \dfrac{\mu_{MR}/\mu_{MD}}{\mu_{FR}/\mu_{FD}}\right)=\beta_{MR}^{GP}-\beta_{MD}^{GP}-\beta_{FR}^{GP}+\beta_{FD}^{GP} \, = \, \beta_{MR}^{GP}\cdot \]

A última igualdade resulta do fato de que \(\beta^{GP}_{MD} = \beta^{GP}_{FR} = \beta^{GP}_{FD} = 0\) porque “F” representa nível 1 de <strong.gender e “D” representa nível 1 de party. Assim, a razão de chances estimada é \(\exp\big(\widehat{\beta}^{GP}_{MR}\big) = \exp(0.276) = 1.3\).

Para \(GI\), construiremos a proporção de chances de estar em uma ideologia mais conservadora em relação a uma ideologia menos conservadora para homens versus mulheres, ou seja, quantidades como “Probabilidades de ser Muito Conservador em comparação com Ligeiramente Conservador dado Homem, dividido por Probabilidades de sendo Muito Conservador em comparação com Ligeiramente Conservador dado Feminino”. Por exemplo, \[ \log(OR) \, = \, \beta_{M1}^{GI}-\beta_{M2}^{GI}-\beta_{F1}^{GI}+\beta_{F2}^{GI} \, = \, \beta_{M1}^{GI}-\beta_{M2}^{GI}, \] onde, como lembrete, \(1=VC\) e \(2=SC\). Usamos um tratamento semelhante de \(PI\), comparando as chances de estar em uma ideologia mais conservadora para republicanos versus democratas. Aqui, esperamos que os republicanos sejam mais conservadores, portanto, esperamos razões de chances positivas de magnitude crescente com o aumento da diferença entre as duas ideologias comparadas.

Há uma razão de chances para \(GP\), dez para \(GI\) e dez para \(PI\). Isso leva a 21 conjuntos de coeficientes de combinação linear nos 16 parâmetros listados na saída acima. Os primeiros sete parâmetros não contribuem para as razões de chances, então seus coeficientes são todos 0. Outros coeficientes são 0, 1 ou -1 conforme determinado pelos cálculos de \(\log(OR)\). Isso leva à matriz de contraste mostrada abaixo para as razões de chances \(GI\):

#####################################################################
# Confidence intervals: LR
# Already did math to determine which coefficients contribute to ORs.
# Develop contrast matrices for the different sets of ORs
# We split them into sets according to the interaction to better define a "family" for simultaneous inference
contr.mat.GP <- rbind(c(rep(0, 7), 1,0,0,0,0,0,0,0,0))
row.names(contr.mat.GP) <- c("GP Rep | M:F")
library(mcprofile)
LRCI <- mcprofile(mod.homo, CM = contr.mat.GP)
exp(confint(LRCI, adjust = "none"))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##              Estimate lower upper
## GP Rep | M:F     1.32 0.988  1.76
exp(confint(LRCI)) #Same as unadjusted here because only one OR in family
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   single-step 
## 
##              Estimate lower upper
## GP Rep | M:F     1.32 0.988  1.76
contr.mat.GI <- rbind (c(rep (0 ,7) ,0 , -1 ,0 ,1 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 , -1 ,1 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 ,0 ,1 , -1 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,1 , -1 ,0 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,1 ,0 ,0 , -1 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 , -1 ,0 ,0 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 ,0 ,0 , -1 ,0 ,0 ,0 ,0) ,
                       c(rep (0 ,7) ,0 ,0 ,1 ,0 , -1 ,0 ,0 ,0 ,0))
row.names(contr.mat.GI) <- c("GI VC:SC | M:F", "GI VC:M | M:F", "GI VC:SL | M:F", "GI VC:VL | M:F", 
                             "GI SC:M | M:F", "GI SC:SL | M:F", "GI SC:VL | M:F", "GI M:SL | M:F", 
                             "GI M:VL | M:F", "GI SL:VL | M:F")
LRCI <- mcprofile(mod.homo, CM = contr.mat.GI)
exp(confint(LRCI, adjust = "none"))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##                Estimate lower upper
## GI VC:SC | M:F    0.917 0.571  1.47
## GI VC:M | M:F     1.565 1.056  2.32
## GI VC:SL | M:F    1.236 0.764  2.01
## GI VC:VL | M:F    1.079 0.652  1.79
## GI SC:M | M:F     1.706 1.118  2.61
## GI SC:SL | M:F    1.348 0.813  2.24
## GI SC:VL | M:F    1.177 0.694  2.00
## GI M:SL | M:F     0.790 0.518  1.21
## GI M:VL | M:F     0.690 0.442  1.08
## GI SL:VL | M:F    0.873 0.519  1.47
exp(confint(LRCI)) 
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   single-step 
## 
##                Estimate lower upper
## GI VC:SC | M:F    0.917 0.475  1.77
## GI VC:M | M:F     1.565 0.906  2.71
## GI VC:SL | M:F    1.236 0.634  2.42
## GI VC:VL | M:F    1.079 0.536  2.18
## GI SC:M | M:F     1.706 0.948  3.08
## GI SC:SL | M:F    1.348 0.668  2.74
## GI SC:VL | M:F    1.177 0.565  2.46
## GI M:SL | M:F     0.790 0.440  1.43
## GI M:VL | M:F     0.690 0.372  1.29
## GI SL:VL | M:F    0.873 0.424  1.80
contr.mat.PI <- rbind(c(rep(0, 7), 0, 0, 0, 0, 0,-1, 0, 1, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0, 0, 1, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0,-1, 1, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0, 0, 1,-1),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 1, 0, 0, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 1,-1, 0, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 1, 0, 0,-1),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0,-1, 0, 0),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0, 0, 0,-1),
                      c(rep(0, 7), 0, 0, 0, 0, 0, 0, 1, 0,-1))
row.names(contr.mat.PI) <- c("PI VC:SC | R:D", "PI VC:M | R:D", "PI VC:SL | R:D", "PI VC:VL | R:D", 
                             "PI SC:M | R:D", "PI SC:SL | R:D", "PI SC:VL | R:D", "PI M:SL | R:D", 
                             "PI M:VL | R:D", "PI SL:VL | R:D")
LRCI <- mcprofile(mod.homo, CM = contr.mat.PI)
exp(confint(LRCI, adjust = "none"))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##                Estimate lower upper
## PI VC:SC | R:D    0.884 0.535  1.45
## PI VC:M | R:D     2.019 1.360  3.02
## PI VC:SL | R:D    3.124 1.922  5.13
## PI VC:VL | R:D    4.775 2.823  8.24
## PI SC:M | R:D     2.284 1.485  3.55
## PI SC:SL | R:D    3.535 2.112  6.00
## PI SC:VL | R:D    5.403 3.108  9.60
## PI M:SL | R:D     1.548 1.016  2.38
## PI M:VL | R:D     2.365 1.484  3.85
## PI SL:VL | R:D    1.528 0.880  2.68
exp(confint(LRCI)) 
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   single-step 
## 
##                Estimate lower upper
## PI VC:SC | R:D    0.884 0.439  1.76
## PI VC:M | R:D     2.019 1.169  3.54
## PI VC:SL | R:D    3.124 1.597  6.25
## PI VC:VL | R:D    4.775 2.314 10.25
## PI SC:M | R:D     2.284 1.260  4.24
## PI SC:SL | R:D    3.535 1.735  7.40
## PI SC:VL | R:D    5.403 2.521 12.08
## PI M:SL | R:D     1.548 0.864  2.82
## PI M:VL | R:D     2.365 1.244  4.68
## PI SL:VL | R:D    1.528 0.712  3.35


Por exemplo, a primeira linha multiplicada pelas estimativas de parâmetro, conforme ordenadas pelo R, fornece \(-\widehat{\beta}^{GI}_{M2}+\widehat{\beta}^{GI}_{M1}\), que é a segunda razão de chances descrita no parágrafo anterior.

Com esta matriz, podemos obter os intervalos de confiança de verossimilhança perfilada usando mcprofile(). Alternativamente, podemos formar intervalos de Wald usando glht() do pacote multcomp ou por cálculo direto. Os pacotes mcprofile e multcomp permitem cálculos com ou sem ajuste para níveis de confiança simultâneos. Observe que, para inferência simultânea, precisamos definir uma “família” sobre a qual queremos aplicar nossas declarações de confiança.

Embora possamos ter 95% de confiança de que todos os intervalos de confiança em uma análise ampla cobrirão seus parâmetros, a desvantagem de escolher um tamanho de família maior é que os intervalos se tornam mais amplos. Por exemplo, se escolhermos todos os 21 intervalos para serem uma família aqui, o nível de confiança para cada intervalo usando o ajuste de Bonferroni seria \(\alpha ' = 0.0024\), portanto, \(Z_{1-\alpha '} = 3.03\).

Em vez disso, parece sensato para este problema fazer uma declaração de confiança separada sobre as razões de chances para cada um dos três efeitos de interação. Assim, a família \(GP\) contém 1 intervalo, enquanto as famílias \(GI\) e \(PI\) contêm 10 cada uma. Um ajuste de Bonferroni nas últimas famílias usaria \(\alpha ' = 0.005\), portanto, \(Z_{1-\alpha '} = 2.81\). Os intervalos LR ajustados a 95%, usando o ajuste de passo único preferido em vez de Bonferroni são mostrados abaixo para a família \(PI\), os detalhes são fornecidos no programa abaixo.

#####################################################################
# Wald Confidence intervals using multcomp
library(multcomp)
wald.GP <- glht(mod.homo, linfct = contr.mat.GP)
wald.GI <- glht(mod.homo, linfct = contr.mat.GI)
wald.PI <- glht(mod.homo, linfct = contr.mat.PI)
# Defaults use multiplicity adjustment for simultaneous confidence level
summary(wald.GP)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                   Estimate Std. Error z value Pr(>|z|)  
## GP Rep | M:F == 0   0.2756     0.1467   1.878   0.0604 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(wald.GP)$confint)
##              Estimate       lwr      upr
## GP Rep | M:F 1.317256 0.9880084 1.756223
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
# Options to get unadjusted (univariate) tests and CIs
summary(wald.GP, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                   Estimate Std. Error z value Pr(>|z|)  
## GP Rep | M:F == 0   0.2756     0.1467   1.878   0.0604 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
exp(confint(wald.GP, calpha = qnorm(0.975))$confint)
##              Estimate       lwr      upr
## GP Rep | M:F 1.317256 0.9880084 1.756223
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
# Defaults use multiplicity adjustment for simultaneous confidence level
summary(wald.GI)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)  
## GI VC:SC | M:F == 0 -0.08632    0.24141  -0.358   0.9964  
## GI VC:M | M:F == 0   0.44809    0.20090   2.230   0.1653  
## GI VC:SL | M:F == 0  0.21198    0.24586   0.862   0.9089  
## GI VC:VL | M:F == 0  0.07634    0.25720   0.297   0.9983  
## GI SC:M | M:F == 0   0.53441    0.21580   2.476   0.0941 .
## GI SC:SL | M:F == 0  0.29830    0.25842   1.154   0.7737  
## GI SC:VL | M:F == 0  0.16266    0.26944   0.604   0.9740  
## GI M:SL | M:F == 0  -0.23610    0.21592  -1.093   0.8068  
## GI M:VL | M:F == 0  -0.37175    0.22729  -1.636   0.4694  
## GI SL:VL | M:F == 0 -0.13565    0.26457  -0.513   0.9858  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(wald.GI)$confint)
##                 Estimate       lwr      upr
## GI VC:SC | M:F 0.9173007 0.4756188 1.769149
## GI VC:M | M:F  1.5653164 0.9061790 2.703898
## GI VC:SL | M:F 1.2361283 0.6332155 2.413101
## GI VC:VL | M:F 1.0793255 0.5360899 2.173038
## GI SC:M | M:F  1.7064375 0.9486172 3.069657
## GI SC:SL | M:F 1.3475715 0.6671040 2.722138
## GI SC:VL | M:F 1.1766321 0.5652865 2.449135
## GI M:SL | M:F  0.7896987 0.4388552 1.421025
## GI M:VL | M:F  0.6895255 0.3715148 1.279748
## GI SL:VL | M:F 0.8731501 0.4250738 1.793550
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.720792
# Options to get unadjusted (univariate) tests and CIs
summary(wald.GI, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)  
## GI VC:SC | M:F == 0 -0.08632    0.24141  -0.358   0.7207  
## GI VC:M | M:F == 0   0.44809    0.20090   2.230   0.0257 *
## GI VC:SL | M:F == 0  0.21198    0.24586   0.862   0.3886  
## GI VC:VL | M:F == 0  0.07634    0.25720   0.297   0.7666  
## GI SC:M | M:F == 0   0.53441    0.21580   2.476   0.0133 *
## GI SC:SL | M:F == 0  0.29830    0.25842   1.154   0.2484  
## GI SC:VL | M:F == 0  0.16266    0.26944   0.604   0.5460  
## GI M:SL | M:F == 0  -0.23610    0.21592  -1.093   0.2742  
## GI M:VL | M:F == 0  -0.37175    0.22729  -1.636   0.1019  
## GI SL:VL | M:F == 0 -0.13565    0.26457  -0.513   0.6082  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
exp(confint(wald.GI, calpha = qnorm(0.975))$confint)
##                 Estimate       lwr      upr
## GI VC:SC | M:F 0.9173007 0.5715123 1.472305
## GI VC:M | M:F  1.5653164 1.0558347 2.320643
## GI VC:SL | M:F 1.2361283 0.7634643 2.001420
## GI VC:VL | M:F 1.0793255 0.6519627 1.786826
## GI SC:M | M:F  1.7064375 1.1178864 2.604852
## GI SC:SL | M:F 1.3475715 0.8120491 2.236255
## GI SC:VL | M:F 1.1766321 0.6938993 1.995193
## GI M:SL | M:F  0.7896987 0.5172104 1.205745
## GI M:VL | M:F  0.6895255 0.4416506 1.076519
## GI SL:VL | M:F 0.8731501 0.5198583 1.466536
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
# Defaults use multiplicity adjustment for simultaneous confidence level
summary(wald.PI)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)    
## PI VC:SC | R:D == 0  -0.1235     0.2547  -0.485  0.98840    
## PI VC:M | R:D == 0    0.7024     0.2031   3.458  0.00472 ** 
## PI VC:SL | R:D == 0   1.1392     0.2503   4.551  < 0.001 ***
## PI VC:VL | R:D == 0   1.5634     0.2728   5.731  < 0.001 ***
## PI SC:M | R:D == 0    0.8259     0.2223   3.716  0.00173 ** 
## PI SC:SL | R:D == 0   1.2627     0.2660   4.748  < 0.001 ***
## PI SC:VL | R:D == 0   1.6869     0.2872   5.875  < 0.001 ***
## PI M:SL | R:D == 0    0.4368     0.2167   2.015  0.25313    
## PI M:VL | R:D == 0    0.8610     0.2427   3.548  0.00349 ** 
## PI SL:VL | R:D == 0   0.4242     0.2833   1.497  0.55765    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
exp(confint(wald.PI)$confint)
##                 Estimate       lwr       upr
## PI VC:SC | R:D 0.8837949 0.4422589  1.766145
## PI VC:M | R:D  2.0186145 1.1620383  3.506601
## PI VC:SL | R:D 3.1242912 1.5820460  6.169982
## PI VC:VL | R:D 4.7750185 2.2746883 10.023704
## PI SC:M | R:D  2.2840304 1.2482516  4.179282
## PI SC:SL | R:D 3.5350861 1.7155196  7.284576
## PI SC:VL | R:D 5.4028579 2.4751935 11.793370
## PI M:SL | R:D  1.5477404 0.8586360  2.789890
## PI M:VL | R:D  2.3654930 1.2229526  4.575449
## PI SL:VL | R:D 1.5283526 0.7075086  3.301531
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.718453
# Options to get unadjusted (univariate) tests and CIs
summary(wald.PI, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
##     data = alldata)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)    
## PI VC:SC | R:D == 0  -0.1235     0.2547  -0.485 0.627644    
## PI VC:M | R:D == 0    0.7024     0.2031   3.458 0.000545 ***
## PI VC:SL | R:D == 0   1.1392     0.2503   4.551 5.34e-06 ***
## PI VC:VL | R:D == 0   1.5634     0.2728   5.731 9.97e-09 ***
## PI SC:M | R:D == 0    0.8259     0.2223   3.716 0.000202 ***
## PI SC:SL | R:D == 0   1.2627     0.2660   4.748 2.06e-06 ***
## PI SC:VL | R:D == 0   1.6869     0.2872   5.875 4.24e-09 ***
## PI M:SL | R:D == 0    0.4368     0.2167   2.015 0.043876 *  
## PI M:VL | R:D == 0    0.8610     0.2427   3.548 0.000388 ***
## PI SL:VL | R:D == 0   0.4242     0.2833   1.497 0.134340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
exp(confint(wald.PI, calpha = qnorm(0.975))$confint)
##                 Estimate       lwr      upr
## PI VC:SC | R:D 0.8837949 0.5364995 1.455907
## PI VC:M | R:D  2.0186145 1.3556182 3.005865
## PI VC:SL | R:D 3.1242912 1.9128328 5.103005
## PI VC:VL | R:D 4.7750185 2.7975605 8.150244
## PI SC:M | R:D  2.2840304 1.4774593 3.530923
## PI SC:SL | R:D 3.5350861 2.0989763 5.953775
## PI SC:VL | R:D 5.4028579 3.0775081 9.485231
## PI M:SL | R:D  1.5477404 1.0120592 2.366957
## PI M:VL | R:D  2.3654930 1.4701130 3.806209
## PI SL:VL | R:D 1.5283526 0.8771217 2.663098
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
#####################################################################
# Wald CIs by manual calculation
# Not doing multiplicity adjustments, so put all contrasts together into one matrix
contr.mat <- rbind(contr.mat.GP, contr.mat.GI, contr.mat.PI)
# Get out coefficients and variances
beta <- matrix(coef(mod.homo), ncol = 1)
v.beta <- vcov(mod.homo)
# Estimate Lin Combos and standard errors as matrix computations
log.contrasts <- contr.mat %*% beta
SElog.contrasts <- matrix(sqrt(diag(contr.mat %*% v.beta %*% t(contr.mat))), ncol = 1)
# Compute confidence intervals in linear predictor scale
alpha = 0.05
lower.log <- log.contrasts + qnorm(alpha/2)*SElog.contrasts
upper.log <- log.contrasts + qnorm(1 - alpha/2)*SElog.contrasts
# Combine Lin Combo coefficients, estimates of contrasts, and confidence intervals in mean scale
wald.ci <- round(data.frame(exp(log.contrasts), exp(lower.log), exp(upper.log)), digits = 2)
# Attach contrast names to and columns.                  
colnames(wald.ci) <- c("Estimate", "Lower CI", "Upper CI")
wald.ci
##                Estimate Lower CI Upper CI
## GP Rep | M:F       1.32     0.99     1.76
## GI VC:SC | M:F     0.92     0.57     1.47
## GI VC:M | M:F      1.57     1.06     2.32
## GI VC:SL | M:F     1.24     0.76     2.00
## GI VC:VL | M:F     1.08     0.65     1.79
## GI SC:M | M:F      1.71     1.12     2.60
## GI SC:SL | M:F     1.35     0.81     2.24
## GI SC:VL | M:F     1.18     0.69     2.00
## GI M:SL | M:F      0.79     0.52     1.21
## GI M:VL | M:F      0.69     0.44     1.08
## GI SL:VL | M:F     0.87     0.52     1.47
## PI VC:SC | R:D     0.88     0.54     1.46
## PI VC:M | R:D      2.02     1.36     3.01
## PI VC:SL | R:D     3.12     1.91     5.10
## PI VC:VL | R:D     4.78     2.80     8.15
## PI SC:M | R:D      2.28     1.48     3.53
## PI SC:SL | R:D     3.54     2.10     5.95
## PI SC:VL | R:D     5.40     3.08     9.49
## PI M:SL | R:D      1.55     1.01     2.37
## PI M:VL | R:D      2.37     1.47     3.81
## PI SL:VL | R:D     1.53     0.88     2.66


Os nomes das linhas foram escolhidos para descrever a razão de chances que está sendo estimada. Por exemplo, “PI VC:SC | R:D” refere-se à proporção de [as chances de ser muito conservador versus ligeiramente conservador para um republicano] e [as chances de ser muito conservador versus ligeiramente conservador para um democrata].

Todas as razões de chances comparam uma ideologia mais conservadora a uma menos conservadora para republicanos vs. democratas e, portanto, espera-se que todas sejam \(\geq 1\). Vemos que isso geralmente é verdade: sete das dez razões de chances \(PI\) têm intervalos de confiança inteiramente acima de 1 e as estimativas e intervalos correspondentes tendem a estar mais distantes de 1 conforme as ideologias comparadas estão mais distantes. Temos 95% de confiança de que esse procedimento produz um conjunto de intervalos que cobrem todos os seus parâmetros.



4.2.5.1 Modelo loglinear vs. modelo logístico e multinomial


Vimos que o modelo loglinear e o modelo de regressão multinomial podem ser aplicados a tabelas de contingência. Existe uma clara preferência por um ou outro?

Em muitos casos a resposta é não. As duas abordagens para modelar contagens em uma tabela de contingência fornecem respostas idênticas devido a uma relação matemática entre as distribuições de Poisson e multinomial. Por um lado, a distribuição multinomial assume que a contagem total em toda a tabela foi fixada por design, por exemplo, especificando um tamanho de amostra total antes da coleta de dados. Por outro lado, o modelo Poisson assume que o total é uma quantidade aleatória, pois é uma soma de variáveis aleatórias independentes. No entanto, verifica-se que a distribuição de probabilidade conjunta de qualquer conjunto de variáveis aleatórias Poisson independentes, condicionadas à sua soma, é uma distribuição multinomial. Ou seja, se tratarmos a contagem total de uma tabela de contingência como uma quantidade fixa, os dois modelos são idênticos.

Além disso, nunca há mal algum em tratar essa soma como fixa, porque não é uma quantidade que tentamos modelar. Portanto, os dois modelos podem ser usados de forma intercambiável na maioria das configurações. Uma característica única dos modelos loglineares é que eles podem ser usados não apenas nos casos em que há uma variável específica que é considerada uma variável de “resposta”, mas também nos casos em que há múltiplas variáveis de resposta ou nenhuma. A regressão logística e os modelos multinomiais requerem uma variável de resposta cujas probabilidades devem ser modeladas. Em uma análise de modelo loglinear com uma variável de resposta, o interesse principal está em modelar associações entre essa variável e as outras. É típico, nesse caso, incluir em qualquer modelo todos os efeitos principais e interações envolvendo as variáveis explicativas e não explorar mais esses efeitos. Por exemplo, se a variável de resposta for \(X\) e as variáveis explicativas forem \(A\), \(B\), \(C\), então o modelo mínimo a ser ajustado seria o Modelo \(ABC,X\).

Para ver por que isso é necessário, considere o caso em que a variável de resposta \(X\) tem dois níveis e há duas variáveis explicativas, \(A\) e \(B\), cada uma com dois níveis. As estruturas de tabela e os parâmetros usados pelos dois modelos são apresentados na Tabela 4.6.

\[ \mbox{Tabela 4.6: Tabela e estruturas de parâmetros para regressão logística (parte superior) } \\ \mbox{e modelo loglinear (parte inferior) quando } X \mbox{ é uma variável de resposta de 2 níveis e } \\ A \mbox{ e } B \mbox{ são variáveis explicativas de 2 níveis cada.} \\ \begin{array}{cc|cc|cc|cc|c} \hline & & B = 1 & & & & B = 2 & & \\ & & X = 0 & X = 1 & & & X = 0 & X = 1 & \\\hline A & 1 & \pi_{11} & 1-\pi_{11} & 1 & 1 & \pi_{12} & 1-\pi_{12} & 1 \\ & 2 & \pi_{21} & 1-\pi_{12} & 1 & 2 & \pi_{22} & 1-\pi_{22} & 1 \\ \hline & & & & & & & & \\ & & B = 1 & & & & B = 2 & & \\ & & X = 0 & X = 1 & & & X = 0 & X = 1 & \\\hline A & 1 & \mu_{111} & 1-\mu_{112} & \mu_{11+} & 1 & \mu_{121} & 1-\mu_{122} & \mu_{12+} \\ & 2 & \mu_{211} & 1-\mu_{212} & \mu_{21+} & 2 & \mu_{221} & 1-\mu_{222} & \mu_{22+} \\\hline \end{array} \]

O modelo de regressão logística saturado escrito em um formato ANOVA especifica \[ \log\left(\dfrac{\pi_{ij}}{1-\pi_{ij}}\right)=\lambda+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB}, \qquad i=1,2, \quad j=1,2, \] onde usamos \(\lambda\)’s para representar os parâmetros de regressão aqui para evitar confusão com o modelo loglinear saturado, que especifica \[ \log(\mu_{ijk})=\beta_0+\beta_i^A+\beta_j^B+\beta_{ij}^{AB}+\beta_k^X+\beta_{ik}^{AX}+\beta_{jk}^{BX}+\beta_{ijk}^{ABX}, \] \(i=1,2\), \(j=1,2\) e \(k=1,2\).

Os parâmetros desses modelos têm algumas interpretações comuns:

  1. Lembre-se de que os efeitos principais no modelo logito, \(\lambda^A_i\) e \(\lambda^B_j\), representam associações entre essa variável e a variável resposta \(X\). No modelo loglinear, as interações de 2 variáveis com \(X\), \(\beta^{AX}_{ik}\) e \(\beta^{BX}_{jk}\), representam as mesmas associações. As razões de chances no modelo logístico são encontradas a partir de diferenças exponenciais nos parâmetros de efeito principal. No modelo loglinear, as mesmas razões de chances são encontradas a partir de contrastes exponenciados (1, -1, -1, 1) entre os parâmetros de interação com \(X\). Pode-se mostrar que as razões de chances calculadas a partir desses dois modelos são idênticas.

  2. No modelo logito, os parâmetros \(\lambda^{AB}_{ij}\) controlam a heterogeneidade da associação. No modelo loglinear, os parâmetros \(\beta_{ijk}^{ABX}\) servem ao mesmo propósito. Em ambos os modelos, as razões de chances previstas entre \(A\) e \(X\) ou \(B\) e \(X\) calculadas separadamente para cada nível da outra variável explicativa são as mesmas que as razões de chances observadas correspondentes nas tabelas. Além disso, em ambos os modelos, a ausência da interação de ordem mais alta resulta em um modelo com associação homogênea.

O resultado final é que qualquer modelo de regressão logística para o modelo de variáveis explicativas categóricas pode ser escrito como um modelo loglinear equivalente com (a) todos os efeitos principais e interações entre as variáveis explicativas, (b) um efeito principal para a variável resposta e (c) uma interação com a resposta para cada termo no modelo logístico, por exemplo, um efeito \(A\) no modelo logístico produz um efeito \(AX\) no modelo loglinear. Se isso for feito, então os parâmetros de associação loglinear \(\beta^{AX}_{ik}\) e \(\beta^{BX}_{jk}\), \(i = 1,\cdots,I\), \(j = 1,\cdots,J\), \(k = 1,2\) são idênticos aos parâmetros de regressão logística \(\lambda^A_i\) e \(\lambda^B_j\), respectivamente e os parâmetros de heterogeneidade \(\beta^{ABX}_{ijk}\) são idênticos aos parâmetros de interação \(\lambda^{AB}_{ij}\). Por exemplo, o código \[ \mbox{logit.homo = glm(formula = party ~ gender + ideol, family = binomial(link = "logit"),} \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \mbox{weights = count, data = alldata)} \] estima o mesmo modelo que mod.homo no exemplo de ideologia política.

Para ajustar um modelo baseado na distribuição Poisson como alternativa a uma regressão multinomial, cada conjunto observado de \(J\) contagens multinomiais é convertido em \(J\) observações com valores idênticos para as variáveis explicativas. Uma nova variável explicativa categórica é criada consistindo nos níveis das categorias da resposta. A variável resposta para cada observação é a contagem observada para essas variáveis explicativas de observações e categoria da resposta. O modelo contém todos os efeitos principais das variáveis explicativas e a nova variável categórica para a categoria da resposta, bem como as interações entre a nova variável e as variáveis explicativas. Os principais interesses residem então neste conjunto de interações.

Como a interpretação das razões de chances e as probabilidades é um pouco mais simples no modelo de regressão logística, tendemos a preferi-la ao modelo loglinear quando há uma variável de resposta de 2 níveis. Os modelos de regressão multinomial não são tão fáceis de usar e interpretar quanto a regressão logística, portanto, a escolha entre os modelos loglinear e multinomial quando uma variável resposta tem mais de dois níveis não é tão clara. A vantagem da formulação loglinear Poisson é que ferramentas muito mais simples de criação, diagnóstico e correção de modelos podem ser usadas, conforme discutido no Capítulo 5. Além disso, quando há uma variável resposta ordinal, o modelo loglinear é mais flexível do que o modelo de regressão multinomial em termos dos tipos de associações que podem ser modeladas, conforme descrito a seguir.


4.2.6 Variáveis categóricas ordinais


Na Seção 3.4, vimos uma abordagem para modelar contagens quando uma variável resposta é ordinal e surge de uma distribuição multinomial. Aqui fornecemos uma abordagem alternativa usando um modelo de regressão Poisson. Este modelo também funciona quando uma ou mais variáveis ordinais são explicativas ou quando não há nenhuma variável resposta designada.

Existem duas abordagens gerais para analisar contagens em diferentes níveis de uma variável ordinal: ignorar a ordinalidade ou atribuir escores numéricos para representar os níveis ordenados. No primeiro caso, aplicam-se as análises das seções anteriores, mas perde-se a informação relativa às contagens médias e odds ratio (razões de chance) à ordenação dos níveis. Este não é um tipo de análise eficiente ou informativo se as contagens médias ou as razões de chance variarem de alguma maneira sistemática com a ordenação. De fato, no exemplo da seção anterior, vimos que analisar as associações entre a variável de ideologia ordinal e as variáveis nominais era um tanto trabalhoso. Portanto, tentamos incorporar a estrutura ordinal das variáveis em uma análise sempre que possível. Em um modelo de regressão Poisson, atribuir escores numéricos para representar os níveis de uma variável ordinal nos permite tratar a variável ordinal como uma variável numérica. Mais uma vez, os métodos das seções anteriores são usados para a análise.

Assim, a presença de variáveis ordinais não representa nenhuma dificuldade em uma regressão Poisson. Portanto, nos concentramos na seleção de escores e na interpretação dos resultados. Existe também um procedimento disponível para verificar se a escolha dos escores falhou em revelar algum aspecto importante da relação entre a variável e as contagens.


4.2.6.1 Escolhendo Escores


A atribuição de escores é subjetiva, principalmente porque normalmente não existe um único sistema de pontuação que seja obviamente melhor do que todos os outros. Por exemplo, dados categóricos ordenados geralmente surgem de questionários que pedem ao respondente que avalie algo em uma escala como “Excelente”, “Muito Bom”, “Bom”, “Regular” e “Ruim”. Se alguém atribui escores numéricas para representar esses níveis, quais devem ser os escores? Devemos escolher 1-2-3-4-5 ou 5-4-3-2-1 ou alguma outra sequência? Isso importa? A resposta varia de acordo com o problema. Escores diferentes às vezes podem influenciar fortemente os resultados da análise. Outras vezes, todas as escolhas razoáveis de escores levam essencialmente aos mesmos resultados.

Os escores devem sempre ser escolhidas em consulta com o especialista no assunto cujas perguntas estão sendo respondidas. A característica importante em qualquer conjunto de escores são os tamanhos relativos dos espaçamentos ou lacunas, entre os escores. Quanto mais semelhantes forem as duas categorias, menor deve ser a diferença entre dois escores. Para o exemplo acima, se pensarmos que a noção por trás de “Regular” está mais próxima de “Bom” do que de “Ruim”, então a diferença entre “Bom” e “Regular” deve ser menor do que a diferença entre “Regular” e “Pobre”. Quanto mais perto é uma questão de julgamento. Se esses três escores devem ser 4-3-1 lacunas de 1 e 2, 5-4-1 lacunas de 1 e 3, 6-4-1 lacunas de 2 e 3 ou algo totalmente diferente é muitas vezes uma questão de adivinhação.

Quando há alguma incerteza na definição dos escores, muitas vezes é uma boa prática duplicar uma análise em vários conjuntos de escores razoáveis. Se os resultados forem muito semelhantes, a escolha dos escores não importa muito e as conclusões podem ser tiradas de acordo. Se os resultados mudarem consideravelmente, mais exploração pode ser necessária para entender a causa da diferença. Em qualquer caso, a menos que um conjunto de escores tenha sido fortemente preferido antes de ver os resultados da análise, ambos os conjuntos de resultados devem ser relatados com consideração aproximadamente igual.

Ao escolher escores, os tamanhos relativos das diferentes lacunas são importantes, não os números reais atribuídos. Por exemplo, escores 1-2-3 fornecem os mesmos resultados de análise que 10-20-30 ou 3-2-1. Claro, as estimativas de proporções de médias dependerão dos valores numéricos que definem uma “mudança de 1 unidade” no escore. Por exemplo, com escores 1-2-3, uma mudança de 1 unidade equivale a mover-se para uma categoria adjacente; com 3-2-1, o movimento é para uma categoria adjacente na direção oposta; e com 10-20-30, uma alteração de 10 unidades é necessária para alterar as categorias, portanto, os parâmetros de inclinação relacionados a uma alteração de 1 unidade podem precisar ser multiplicados por 10 para serem significativos. No entanto, os resultados do teste serão idênticos em todos os casos.


Exemplo 4.9: Ideologia política.


Vimos que a variável ideologia neste exemplo é ordinal, então tentaremos usar esta estrutura para construir uma descrição mais significativa e parcimoniosa da associação. Como esse exemplo surge de “dados interessantes” e não de um problema trazido a nós por um cientista político, estamos agindo como especialistas no assunto.

Escolhemos abordar dois conjuntos diferentes de questões usando a ordenação das categorias para a variável ideologia de maneiras diferentes. O primeiro conjunto de questões faz uso do fato de que as cinco categorias são naturalmente ordenadas de mais liberais para menos liberais, equivalentemente, de menos conservadores para mais conservadores. Como padrão, escolhemos 1-2-3-4-5 para representar essa progressão, em que os escores mais altos representam ideologias mais conservadoras. No entanto, suspeitamos que os grupos “ligeiramente” liberais e conservadores possam estar mais próximos do meio do que dos extremos de suas respectivas ideologias. Portanto, consideramos escores com lacunas maiores entre “muito” e “levemente” do que entre “levemente” e “moderado”.

De forma um tanto arbitrária, optamos por dobrar os tamanhos das lacunas, criando escores de 0-2-3-4-6. Certamente esperamos que qualquer conjunto de escores seja associado aos partidos políticos, já que cada partido busca deliberadamente apelar para uma ideologia diferente. Seria interessante ver se existe uma associação entre gênero e ideologia; os homens são mais liberais ou conservadores que as mulheres? ou se a associação entre partido e ideologia é diferente para homens e mulheres.

O segundo conjunto de questões refere-se a se um partido ou gênero tem opiniões mais extremas do que o outro. Em outras palavras, um partido ou gênero tende a ter mais pessoas com pontos de vista “muito” fortes em oposição a pontos de vista “ligeiramente” fortes ou centrais? Para examinar isso, criamos escores para medir o extremo da ideologia: 2-1-0-1-2. Isso basicamente cria três categorias que medem a força da ideologia, baixa, média e alta, com escores de 0-1-2 nesses níveis. Os escores são criados e adicionados aos dados conforme mostrado abaixo:

lin.score1 <- c( rep(x = c(1 ,2 ,3 ,4 ,5) , times = 4))
lin.score2 <- c( rep(x = c(0 ,2 ,3 ,4 ,6) , times = 4))
extrm.score <- c(rep(x = c(2 ,1 ,0 ,1 ,2) , times = 4))
alldata <- data.frame ( alldata , lin.score1 , lin.score2 , extrm.score )
head ( alldata )
##   gender party ideol count lin.score1 lin.score2 extrm.score
## 1      F     D    VL    44          1          0           2
## 2      F     D    SL    47          2          2           1
## 3      F     D     M   118          3          3           0
## 4      F     D    SC    23          4          4           1
## 5      F     D    VC    32          5          6           2
## 6      F     R    VL    18          1          0           2


Observe que a função ifelse() poderia ter sido usada para criar os escores.



4.2.6.2 Modelos e interpretação


Para uma única variável categórica ordenada \(X\) com \(I\) níveis, seja \(s_i\), \(i = 1,\cdots,I\) o conjunto de escores numéricos atribuídas aos níveis. O modelo ordinal de regressão Poisson trata a contagem \(y_i\) como \(Po(\mu_i)\), com \(\log(\mu_i) = \beta_0 + \beta_1 s_i\). A interpretação dos parâmetros de regressão é a mesma de qualquer variável explicativa numérica em uma regressão Poisson conforme descrito na Seção 4.2.1.

No entanto, as “unidades” aplicadas aos escores geralmente não têm significado, exceto nos casos em que as categorias ordenadas representam intervalos de um continuum. Por exemplo, se substituirmos categorias como “Bom”, “Regular” e “Ruim” por valores 3-2-1, uma alteração de 1 unidade no escore é apenas uma alteração na categoria. No entanto, se substituirmos as faixas etárias “<25”, “25-45” e “>45” por escores, como 20-35-60, uma alteração de 1 unidade no escore é aproximadamente equivalente a uma alteração de 1 ano na idade.

Examinar proporções de médias quando \(s\) muda em 1 unidade pode não atender adequadamente a todos os interesses se algumas das lacunas entre os escores forem maiores que 1. Portanto, frequentemente desejaremos examinar mudanças para uma mudança de \(c\) unidades em \(s\), onde \(c\) representa a diferença nos escores entre quaisquer duas categorias de interesse. Conforme descrito na Seção 4.2.1, isso significa simplesmente analisar \(c \,\beta\).

Quando há duas ou mais variáveis explicativas categóricas, pelo menos uma das quais é ordinal, então a forma como os escores entram no modelo depende do objetivo da análise; ou seja, se o foco está na modelagem e previsão de médias ou na modelagem de associações entre variáveis. Se o problema for mais parecido com uma regressão típica e as médias forem o foco principal, quaisquer conjuntos de escores serão tratados como variáveis numéricas e as análises usuais serão realizadas conforme descrito na Seção 4.2.1. Adicionar termos de interação ao modelo é opcional e depende dos objetivos da análise e do ajuste dos modelos resultantes. O procedimento para verificar o ajuste dos escores descrito posteriormente nesta seção pode ser aplicado separadamente para cada variável ordinal na regressão.

Quando as associações são o foco principal, são necessários modelos especiais, especialmente quando há variáveis categóricas ou ordinais adicionais. As associações entre variáveis categóricas são medidas pelas razões de chances, portanto, precisamos considerar as suposições que diferentes modelos fazem sobre como as razões de chances são estruturadas. Para começar, considere um problema com uma variável nominal \(X\) e uma variável ordinal \(Z\), com escores \(s_j\), \(j=1,\cdots,J\), como uma tabela de contingência com \(X\) como variável de linha e \(Z\) como variável de coluna. Os parâmetros de associação no modelo nominal usual de duas variáveis, \[ \log(\mu_{ij}) = \beta_0+ \beta^X_i + \beta^Z_j + \beta^{XZ}_{ij}, \]

não assume nenhuma estrutura particular para as razões de chances entre \(X\) e \(Z\).

Em vez disso, usamos um modelo que assume que, para qualquer par de linhas, o \(\log(OR)\) muda linearmente com a diferença nos escores entre as duas colunas que estão sendo comparadas. Especificamente, o modelo é \[ \log(\mu_{ij}) = \beta_0+ \beta^X_i + \beta^Z_j + \beta^{XZ}_{ij} s_j\cdot \]

Para ver a estrutura que esse modelo impõe às associações, considere novamente a razão de chances envolvendo os níveis \(i\) e \(i'\) de \(X\) e \(j\) e \(j'\) de \(Z\): \[ \begin{array}{rcl} \log(OR_{i\, i', \, j\, j'}) & = & \log(\mu_{ij})+\log(\mu_{i' \, j'})-\log(\mu_{i' \, j})-\log(\mu_{i \,j'}) \\ & = & \beta^{XZ}_{i} s_j + \beta^{XZ}_{i'} s_{j'}- \beta^{XZ}_{i'} s_j-\beta^{XZ}_{i} s_{j'} \, = \, \big(\beta^{XZ}_{i}-\beta^{XZ}_{i'} \big)(s_j-s_{j'})\cdot \end{array} \]

Por exemplo, suponha que \(J = 3\) e os escores sejam 1-2-3 e queremos estimar as razões de chances entre dois níveis de \(Z\) nos níveis 1 e 2 de \(X\). Então, temos \[ \begin{array}{rcl} \log(OR_{12,12}) & = & -1 \, \big(\beta^{XZ}_{1}-\beta^{XZ}_{2} \big) \\ \log(OR_{12,13}) & = & -2 \, \big(\beta^{XZ}_{1}-\beta^{XZ}_{2} \big) \\ \log(OR_{12,23}) & = & -1 \, \big(\beta^{XZ}_{1}-\beta^{XZ}_{2} \big)\cdot \end{array} \] Esse modelo às vezes é chamado de modelo de associação linear, embora a linearidade esteja tecnicamente no logaritmo das razões de chances. Observe que a equação \(\log(\mu_{ij}) = \beta_0+ \beta^X_i + \beta^Z_j + \beta^{XZ}_{ij} s_j\) especifica apenas \((I-1)\) parâmetros para descrever a associação, em oposição aos \((I-1)(J-1)\) parâmetros usuais para associação em uma tabela nominal. Assim, usar a ordenação nos permite modelar um conjunto de dados com menos parâmetros, o que é sempre desejável, mas ao custo de mais suposições sobre a estrutura dos dados. Mais adiante nesta seção, veremos como podemos testar se as suposições adicionadas sobre a estrutura de associação se ajustam bem aos dados.

Observe também que os escores não são usadas na definição dos efeitos de \(Z\) nas contagens marginais: o efeito principal \(\beta^Z_j\) é retido em vez de \(\beta^Z \, s_j\). Isso ocorre porque estamos examinando se as associações são afetadas por \(X\) e \(Z\) e os termos de efeito principal não entram nas razões de chances. Portanto, mantemos um modelo nominal para cada uma das contagens marginais. Então, qualquer diferença no ajuste geral do modelo devido ao uso de associação linear, em vez de nominal, é causada diretamente pela alteração na estrutura da associação e não também por uma alteração no ajuste das médias marginais.

Um teste para o termo de interação, \(H_0\, : \, \beta^{XZ}_i = 0\), \(i = 1,\cdots,I\) para este modelo testa a hipótese nula de independência contra a alternativa de uma mudança linear em \(\log(OR)\). Deixar de rejeitar \(H_0\) aqui não significa que não haja associação entre \(X\) e \(Z\). Significa apenas que não há evidência suficiente de uma tendência linear para \(\log(OR)\). Pode existir um padrão de associação de modo que as razões de chances não fiquem continuamente mais extremas à medida que a diferença entre seus escores aumentam. Da mesma forma, rejeitar a hipótese nula não implica que a associação seja inteiramente linear para \(\log(OR)\). Pode haver uma tendência geral de aumento ou diminuição com o aumento da distância entre os escores, mas a tendência não precisa ser linear. Testar a qualidade do ajuste de um modelo conforme descrito após o exemplo pode ajudar a avaliar a possibilidade de padrões alternativos.


Exemplo 4.10: Ideologia política.


Nosso trabalho com modelos nominais sugeriu que a interação de três variáveis não é necessária, então começamos com um modelo de associação homogênea e consideramos adicionar escores. Para começar, consideramos os escores 1-2-3-4-5 contidos em lin.score.1 e os aplicamos ao termo \(PI\), pois acreditamos que as razões de chance dessa associação devem aumentar à medida que aumenta a distância entre as ideologias.

O modelo que ajustamos é \[ \log(\mu_{ijk})=\beta_0+\beta^G_i+\beta^P_j+\beta^I_k+\beta^{GP}_{ij}+\beta^{GI}_{ik}+\beta^{PI}_j \, s_k^I, \] que assume associação linear entre partido e ideologia que é homogênea entre os gêneros.

Para as ideologias \(k\) e \(k'\) temos \[ \log(OR^{PI}_{RD,k \, k'}) = \big(\beta^{PI}_{R}-\beta^{PI}_D \big)\big(s_k^I-s_{k'}^I \big), \] para ambos os sexos, onde o valor de \(s_k^I-s_{k'}^I\) é a diferença de escores entre as duas ideologias comparadas. Observe que isso significa que as proporções de chances envolvendo quaisquer duas ideologias consecutivas são as mesmas. Além disso, por exemplo, a razão de chances comparando as chances do republicano para ideologias muito conservadoras versus ligeiramente liberais é a mesma que compara as chances do republicano entre ligeiramente conservador e muito liberal; ambas as diferenças de escore são 3. Nós ajustamos e analisamos este modelo em detalhes abaixo.

Mais modelos são ajustados no programa que acompanha este exemplo. Os termos da primeira equação neste exemplo são especificados concisamente em glm() usando apenas os termos de interação conforme mostrado abaixo:

# Homogeneous, Linear in PI: Score Set 1-2-3-4-5
mod.homo.lin1.PI <- glm(formula = count ~ gender*party + gender*ideol + party*lin.score1, 
                        family = poisson(link = "log"), data = alldata)
summary(mod.homo.lin1.PI)
## 
## Call:
## glm(formula = count ~ gender * party + gender * ideol + party * 
##     lin.score1, family = poisson(link = "log"), data = alldata)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.24363  -0.54155   0.02178   0.52057   0.90547  
## 
## Coefficients: (1 not defined because of singularities)
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        4.72852    0.08169  57.881  < 2e-16 ***
## genderM           -0.71139    0.13656  -5.209 1.90e-07 ***
## partyR            -1.51082    0.20862  -7.242 4.42e-13 ***
## ideolSC           -1.40587    0.14897  -9.437  < 2e-16 ***
## ideolSL           -0.83121    0.13685  -6.074 1.25e-09 ***
## ideolVC           -1.41136    0.15424  -9.150  < 2e-16 ***
## ideolVL           -0.89370    0.14906  -5.996 2.03e-09 ***
## lin.score1              NA         NA      NA       NA    
## genderM:partyR     0.28723    0.14555   1.973  0.04845 *  
## genderM:ideolSC    0.55859    0.21413   2.609  0.00909 ** 
## genderM:ideolSL    0.23706    0.21543   1.100  0.27115    
## genderM:ideolVC    0.43621    0.20135   2.166  0.03028 *  
## genderM:ideolVL    0.37450    0.22682   1.651  0.09872 .  
## partyR:lin.score1  0.43058    0.06005   7.170 7.48e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 255.1481  on 19  degrees of freedom
## Residual deviance:   8.3986  on  7  degrees of freedom
## AIC: 142.83
## 
## Number of Fisher Scoring iterations: 4


Como o escore linear lin.score1 aparece apenas em uma interação e não como um efeito principal no argumento formula, um efeito principal é criado automaticamente para ele.

Da mesma forma, os efeitos principais são criados para gender, party e ideol. Isso significa que o principal efeito da ideologia é incluído no modelo tanto como uma série de variáveis indicadoras quanto como uma variável numérica. Essa duplicação leva à observação sobre “singularidades” na listagem dos coeficientes do modelo estimado e os resultados “NA” para o efeito principal lin.score.1.

A causa exata da “singularidade” está além do escopo deste texto. Os leitores que estão familiarizados com matrizes de modelo, por exemplo, a matriz \(X\) na especificação do modelo linear \(Y = X\beta + \epsilon\), podem apreciar o fato de que a variável lin.score1 pode ser escrita como uma combinação linear de (Intercept) e os quatro indicadores que representam a variável ideol. Como resultado, lin.score1 não contém nenhuma informação nova que já não tenha sido contabilizada pelo termo ideol, então as colunas da matriz que representa este modelo não são linearmente independentes.

Essa redundância e a resultante falta de estimativa não criam dificuldade de interpretação, porque não estamos modelando os principais efeitos do escore da ideologia. O parâmetro para a interação entre party e lin.score.1 é estimado corretamente, então podemos ignorar o problema, exceto por alguns detalhes computacionais observados no programa. Em particular, glht() e mcprofile() não funcionarão por causa do parâmetro extra. Podemos coagir glht() a funcionar removendo o NA do vetor de estimativas de parâmetros, os detalhes estão no programa. Infelizmente, mcprofile() ainda não funcionará após essa alteração.

Começamos com os testes de efeitos usando Anova().

library(car)
Anova(mod.homo.lin1.PI)
## Analysis of Deviance Table (Type II tests)
## 
## Response: count
##                  LR Chisq Df Pr(>Chisq)    
## gender             20.637  1  5.551e-06 ***
## party               0.528  1    0.46736    
## ideol             163.007  3  < 2.2e-16 ***
## lin.score1                 0               
## gender:party        3.901  1    0.04826 *  
## gender:ideol        9.336  4    0.05323 .  
## party:lin.score1   55.402  1  9.825e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Vemos que existe uma associação linear muito forte entre partido e ideologia. Os testes para as outras associações produzem resultados muito semelhantes aos do caso nominal, embora se usarmos um nível de significância rígido de 0.05, \(GP\) agora é significativo.

Na prática, os dois testes estão nos dizendo a mesma coisa: há evidências moderadas para sugerir que as interações \(GP\) e \(GI\) estão presentes. As razões de chances para \(GP\) e \(GI\) são estimadas como no caso nominal, construindo combinações lineares dos parâmetros do modelo estimado. Para \(PI\), \(\beta^P_R\) é estimado na saída pelo coeficiente partyR:lin.score1, enquanto \(\beta^P_D = 0\), então os coeficientes necessários para a combinação linear são apenas as diferenças possíveis nos escores, 1, 2, 3 ou 4.

Ou seja, uma diferença de \(c\) unidades é resultado em escores em \(\widehat{OR} = \exp\big (c \, \widehat{\beta}^P_R\big)\). Novamente, como os republicanos são mais conservadores do que os democratas, esperamos razões de chances positivas de magnitude crescente à medida que a diferença nos escores aumenta para as duas ideologias comparadas.

Os cálculos para os intervalos de confiança são mostrados no programa abaixo.

# Print out model matrix to understand the singularity on linear association model
model.matrix(mod.homo.lin1.PI)
##    (Intercept) genderM partyR ideolSC ideolSL ideolVC ideolVL lin.score1
## 1            1       0      0       0       0       0       1          1
## 2            1       0      0       0       1       0       0          2
## 3            1       0      0       0       0       0       0          3
## 4            1       0      0       1       0       0       0          4
## 5            1       0      0       0       0       1       0          5
## 6            1       0      1       0       0       0       1          1
## 7            1       0      1       0       1       0       0          2
## 8            1       0      1       0       0       0       0          3
## 9            1       0      1       1       0       0       0          4
## 10           1       0      1       0       0       1       0          5
## 11           1       1      0       0       0       0       1          1
## 12           1       1      0       0       1       0       0          2
## 13           1       1      0       0       0       0       0          3
## 14           1       1      0       1       0       0       0          4
## 15           1       1      0       0       0       1       0          5
## 16           1       1      1       0       0       0       1          1
## 17           1       1      1       0       1       0       0          2
## 18           1       1      1       0       0       0       0          3
## 19           1       1      1       1       0       0       0          4
## 20           1       1      1       0       0       1       0          5
##    genderM:partyR genderM:ideolSC genderM:ideolSL genderM:ideolVC
## 1               0               0               0               0
## 2               0               0               0               0
## 3               0               0               0               0
## 4               0               0               0               0
## 5               0               0               0               0
## 6               0               0               0               0
## 7               0               0               0               0
## 8               0               0               0               0
## 9               0               0               0               0
## 10              0               0               0               0
## 11              0               0               0               0
## 12              0               0               1               0
## 13              0               0               0               0
## 14              0               1               0               0
## 15              0               0               0               1
## 16              1               0               0               0
## 17              1               0               1               0
## 18              1               0               0               0
## 19              1               1               0               0
## 20              1               0               0               1
##    genderM:ideolVL partyR:lin.score1
## 1                0                 0
## 2                0                 0
## 3                0                 0
## 4                0                 0
## 5                0                 0
## 6                0                 1
## 7                0                 2
## 8                0                 3
## 9                0                 4
## 10               0                 5
## 11               1                 0
## 12               0                 0
## 13               0                 0
## 14               0                 0
## 15               0                 0
## 16               1                 1
## 17               0                 2
## 18               0                 3
## 19               0                 4
## 20               0                 5
## attr(,"assign")
##  [1] 0 1 2 3 3 3 3 4 5 6 6 6 6 7
## attr(,"contrasts")
## attr(,"contrasts")$gender
## [1] "contr.treatment"
## 
## attr(,"contrasts")$party
## [1] "contr.treatment"
## 
## attr(,"contrasts")$ideol
## [1] "contr.treatment"
# Now estimate associations for Linear Association model and compute Wald CIs
contr.mat <- rbind(c(rep(0, 7), 1, 0, 0, 0, 0, 0), 
                   c(rep(0, 7), 0,-1, 0, 1, 0, 0),
                   c(rep(0, 7), 0, 0, 0, 1, 0, 0),
                   c(rep(0, 7), 0, 0,-1, 1, 0, 0),
                   c(rep(0, 7), 0, 0, 0, 1,-1, 0),
                   c(rep(0, 7), 0, 1, 0, 0, 0, 0),
                   c(rep(0, 7), 0, 1,-1, 0, 0, 0),
                   c(rep(0, 7), 0, 1, 0, 0,-1, 0),
                   c(rep(0, 7), 0, 0,-1, 0, 0, 0),
                   c(rep(0, 7), 0, 0, 0, 0,-1, 0),
                   c(rep(0, 7), 0, 0, 1, 0,-1, 0),
                   c(rep(0, 7), 0, 0, 0, 0, 0, 1),
                   c(rep(0, 7), 0, 0, 0, 0, 0, 2),
                   c(rep(0, 7), 0, 0, 0, 0, 0, 3),
                   c(rep(0, 7), 0, 0, 0, 0, 0, 4))
# Usual method for Wald intervals creates error message:
#  "Error in modelparm.default(model, ...) : 
#    dimensions of coefficients and covariance matrix don't match"
#  because of extra parameter with missing value in coefficients that is not in variance.
# We can remove this parameter using the code below: 
# First, we use is.na(mod.homo.lin1.PI$coefficients) to identify which elements of the 
#  coefficient vector are NA. The result is TRUE if the element is NA and FALSE otherwise. 
# Second, we reverse the TRUE and FALSE responses by using !, so that the positions of 
#  the non-NA elements are now given by TRUE. 
# Finally, we have R return only the non-NA elements and put them back in the original 
#  component of the mod.homo.lin1.PI object. 
# Note that the estimated covariance matrix given by vcov(mod.homo.lin1.PI) does not contain 
#  lin.score1 so we do not need to perform the same steps for it. 
# With the model fit object amended in this way, glht() works However, mcprofile() still 
#  will not work. 
check.na <- is.na(mod.homo.lin1.PI$coefficients)
mod.homo.lin1.PI$coefficients <- mod.homo.lin1.PI$coefficients[!check.na]
library(multcomp)
wald <- glht(mod.homo.lin1.PI, linfct = contr.mat)
wald.ci <- round(exp(confint(wald, calpha = qnorm(0.975))$confint), 2)
row.names(wald.ci) <- c("GP Rep | M:F", "GI VC:SC | M:F", "GI VC:M | M:F", 
                        "GI VC:SL | M:F", "GI VC:VL | M:F", "GI SC:M | M:F", 
                        "GI SC:SL | M:F", "GI SC:VL | M:F", "GI M:SL | M:F", 
                        "GI M:VL | M:F", "GI SL:VL | M:F", "PI REP | 1 Cat Ideol", 
                        "PI REP | 2 Cat Ideol", "PI REP | 3 Cat Ideol", "PI REP | 4 Cat Ideol")
colnames(wald.ci) <- c("Estimate", "Lower CI", "Upper CI")
wald.ci
##                      Estimate Lower CI Upper CI
## GP Rep | M:F             1.33     1.00     1.77
## GI VC:SC | M:F           0.88     0.55     1.42
## GI VC:M | M:F            1.55     1.04     2.30
## GI VC:SL | M:F           1.22     0.75     1.98
## GI VC:VL | M:F           1.06     0.64     1.76
## GI SC:M | M:F            1.75     1.15     2.66
## GI SC:SL | M:F           1.38     0.83     2.28
## GI SC:VL | M:F           1.20     0.71     2.03
## GI M:SL | M:F            0.79     0.52     1.20
## GI M:VL | M:F            0.69     0.44     1.07
## GI SL:VL | M:F           0.87     0.52     1.46
## PI REP | 1 Cat Ideol     1.54     1.37     1.73
## PI REP | 2 Cat Ideol     2.37     1.87     2.99
## PI REP | 3 Cat Ideol     3.64     2.56     5.18
## PI REP | 4 Cat Ideol     5.60     3.50     8.96
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
# This fix does not help with making mcprofile work.  See error messages.
#library(mcprofile)
#LRCI <- mcprofile(mod.homo.lin1.PI, CM = contr.mat)
#exp(confint(LRCI, adjust = "none")) # Unadjusted for multiplicity
# Wald CIs by manual calculation
# Get out coefficients and variances
# NOTE THAT THE PARAMETER ESTIMATES HAVE ALREADY HAD THE MISSING VALUE REMOVED.
beta <- matrix(coef(mod.homo.lin1.PI), ncol = 1)
# Note that the variance-covariabce matrix does NOT include the singularity, so no deletions are needed
v.beta <- vcov(mod.homo.lin1.PI)
# Estimate Lin Combos and standard errors as matrix computations
log.contrasts <- contr.mat %*% beta
SElog.contrasts <- matrix(sqrt(diag(contr.mat %*% v.beta %*% t(contr.mat))), ncol = 1)
# Compute confidence intervals in linear predictor scale
alpha = 0.05
lower.log <- log.contrasts + qnorm(alpha/2)*SElog.contrasts
upper.log <- log.contrasts + qnorm(1 - alpha/2)*SElog.contrasts
# Combine Lin Combo coefficients, estimates of contrasts, and confidence intervals in mean scale
wald.ci.1 <- round(data.frame(exp(log.contrasts), exp(lower.log), exp(upper.log)), digits = 2)
# Attach contrast names to rows and columns.                  
rownames <- c("GP Rep | M:F", "GI VC:SC | M:F", "GI VC:M | M:F", "GI VC:SL | M:F", 
              "GI VC:VL | M:F", "GI SC:M | M:F", "GI SC:SL | M:F", "GI SC:VL | M:F", 
              "GI M:SL | M:F", "GI M:VL | M:F", "GI SL:VL | M:F", "PI REP | 1 Cat Ideol", 
              "PI REP | 2 Cat Ideol", "PI REP | 3 Cat Ideol", "PI REP | 4 Cat Ideol")
colnames(wald.ci.1) <- c("Estimate", "Lower CI", "Upper CI")
wald.ci.1
##    Estimate Lower CI Upper CI
## 1      1.33     1.00     1.77
## 2      0.88     0.55     1.42
## 3      1.55     1.04     2.30
## 4      1.22     0.75     1.98
## 5      1.06     0.64     1.76
## 6      1.75     1.15     2.66
## 7      1.38     0.83     2.28
## 8      1.20     0.71     2.03
## 9      0.79     0.52     1.20
## 10     0.69     0.44     1.07
## 11     0.87     0.52     1.46
## 12     1.54     1.37     1.73
## 13     2.37     1.87     2.99
## 14     3.64     2.56     5.18
## 15     5.60     3.50     8.96
library(xtable)
print(xtable(cbind(rownames, 1/wald.ci.1)), floating = FALSE, include.rownames = FALSE)
## % latex table generated in R 4.2.2 by xtable 1.8-4 package
## % Tue Nov  7 20:38:06 2023
## \begin{tabular}{lrrr}
##   \hline
## rownames & Estimate & Lower CI & Upper CI \\ 
##   \hline
## GP Rep $|$ M:F & 0.75 & 1.00 & 0.56 \\ 
##   GI VC:SC $|$ M:F & 1.14 & 1.82 & 0.70 \\ 
##   GI VC:M $|$ M:F & 0.65 & 0.96 & 0.43 \\ 
##   GI VC:SL $|$ M:F & 0.82 & 1.33 & 0.51 \\ 
##   GI VC:VL $|$ M:F & 0.94 & 1.56 & 0.57 \\ 
##   GI SC:M $|$ M:F & 0.57 & 0.87 & 0.38 \\ 
##   GI SC:SL $|$ M:F & 0.72 & 1.20 & 0.44 \\ 
##   GI SC:VL $|$ M:F & 0.83 & 1.41 & 0.49 \\ 
##   GI M:SL $|$ M:F & 1.27 & 1.92 & 0.83 \\ 
##   GI M:VL $|$ M:F & 1.45 & 2.27 & 0.93 \\ 
##   GI SL:VL $|$ M:F & 1.15 & 1.92 & 0.68 \\ 
##   PI REP $|$ 1 Cat Ideol & 0.65 & 0.73 & 0.58 \\ 
##   PI REP $|$ 2 Cat Ideol & 0.42 & 0.53 & 0.33 \\ 
##   PI REP $|$ 3 Cat Ideol & 0.27 & 0.39 & 0.19 \\ 
##   PI REP $|$ 4 Cat Ideol & 0.18 & 0.29 & 0.11 \\ 
##    \hline
## \end{tabular}
# Test fit of scores
# Fit same model with nominal association for PI
mod.homo <- glm(formula = count ~ (gender + party + ideol)^2, family = poisson(link = "log"), 
                data = alldata)


Comparando essas estimativas e intervalos de confiança com aqueles obtidos para a análise nominal, há concordância muito próxima para todos os odds ratios \(GP\) e \(GI\). Por exemplo, dentre as razões de chances fornecidas acima, quatro delas estimam como as chances de ser republicano mudam à medida que as ideologias diferem progressivamente em mais categorias. Para categorias adjacentes de ideologia, a chance de ser republicano é 1.54 vezes maior para a ideologia mais conservadora do que para a mais liberal, \(1.37 < OR < 1.73\).

Para uma diferença de 2 categorias, as chances de ser republicano aumentam para 2.37 vezes mais na ideologia mais conservadora, \(1.87 < OR < 2.99\); para 3 categorias para 3.64, \(2.56 < OR < 5.18\) e, finalmente, as chances de ser republicano são 5.6 vezes maiores, \(3.50 < OR < 8.96\) para uma ideologia muito conservadora em comparação com uma ideologia muito liberal.

Observe que tanto a estimativa quanto a interpretação da associação \(PI\) são muito mais fáceis quando a ordinalidade é usada. Isso destaca o apelo de usar a ordinalidade quando ela está disponível e quando a associação pode ser razoavelmente descrita pela tendência linear ao longo do logaritmos das razões de chance. Quando vistas nominalmente, as estimativas baseadas em ideologias consecutivas foram 0.88, 2.28, 1.55 e 1.53. Nossa estimativa de 1.54 está no meio desses números. No entanto, observe que as duas primeiras dessas estimativas, VC:SC e SC:M possuem intervalos de confiança que não se sobrepõem na análise nominal, 0.54 a 1.46 para VC:SC e 1.48 a 3.52 para SC:M. Isso sugere que as verdadeiras razões de chances que esses dois intervalos estão estimando podem não ter o mesmo valor. Assim, resta saber se o padrão linear para \(\log(OR)\) realmente se ajusta a esses dados.

Esses intervalos de confiança também demonstram a outra grande vantagem de usar a ordinalidade explicitamente: os intervalos de confiança para razões de chances comparáveis são geralmente muito mais curtos na análise ordinal do que na análise nominal. Por exemplo, o intervalo de confiança comparando as duas ideologias mais extremas foi \(2.3 < OR < 10.3\) na análise nominal, que é consideravelmente mais longo do que o intervalo de confiança correspondente da análise ordinal. A capacidade de representar a associação \(PI\) usando menos parâmetros permite que esses parâmetros sejam estimados com mais precisão. Essa precisão é transferida para inferências sobre as razões de chances.


Quando ambas as variáveis em um modelo de duas variáveis são ordinais, podemos definir escores \(s^X_i\), \(i = 1,\cdots,I\) para \(X\) e \(s^Z_j\), \(j = 1,\cdots,J\) para \(Z\). Os \(I-1\) parâmetros de associação são substituídos por um único parâmetro usando o modelo \[ \log(\mu_{ij})=\beta_0+\beta^X_i+\beta^Z_j+\beta^{XZ} s^X_{i} s_j^Z\cdot \]

Pode-se ver que para este modelo \[ \log(OR_{ii', jj'}) = \beta^{XZ} \big(s^X_i-s^X_{i'}\big)\big(s^X_j-s^X_{j'}\big)\cdot \] Essa estrutura às vezes é chamada de associação linear-por-linear. A redução dos parâmetros de associação vem à custa de suposições mais estruturadas sobre a natureza dos logaritmos das razões de chance.

Portanto, quando um teste para \(H_0 \, : \, \beta^{XZ} = 0\) é rejeitado, isso não implica que a associação entre \(X\) e \(Z\) seja totalmente linear; pode haver estrutura adicional não capturada por \(\beta^{XZ}\). Da mesma forma, deixar de rejeitar \(H_0 \, : \, \beta^{XZ} = 0\) significa apenas que não há evidência suficiente de uma associação linear-por-linear, não que não haja nenhuma associação.

Com mais de duas variáveis, as associações continuam a ser definidas por interações de duas variáveis e de ordem superior. Os termos de associação nominal, linear e linear-por-linear para diferentes pares de variáveis podem aparecer juntos em qualquer combinação — de fato, o exemplo da ideologia política continha dois termos de associação nominal e um termo de associação linear. A colocação de estruturas lineares em interações de três variáveis pressupõe que os \(\log(OR)\) entre duas variáveis mudam linearmente nos escores da terceira variável. A modelagem de interações envolvendo mais de três variáveis pode ser feita. A principal dificuldade então está interpretando os modelos resultantes.


4.2.6.3 Verificando a bondade de ajuste para um conjunto de escores


Para uma única variável categórica ordenada \(X\) com \(I\) níveis, \(s_i\) representa o conjunto de escores atribuídos aos níveis. Como na Seção 4.2.3, sejam \(x_2,\cdots,x_I\) são as variáveis indicadoras que seriam utilizadas para representar os níveis de \(X\) caso a ordenação não fosse considerada. Podemos considerar o modelo nominal \[ \log(\mu_i) =\beta_0 +\beta_2 x_2 + \cdots + \beta_I x_I \] e o modelo ordinal \[ \log(\mu_i) = \gamma_0 +\gamma_1 s_i \]

como duas explicações alternativas para as mudanças na contagem média.

O modelo nominal cria um parâmetro separado para cada média e, portanto, representa um modelo saturado, de modo que o desvio residual é zero quando há uma observação por nível de \(X\). O modelo ordinal assume que as médias da população mudam de forma log-linear com os escores. Quando isso é verdade, os dados devem ajustar-se razoavelmente ao modelo e, portanto, o desvio residual deve ser relativamente pequeno. Se o modelo ordinal estiver errado, as médias previstas por ele podem estar longe das médias observadas e o desvio residual será grande. A comparação dos desvios residuais dos modelos nominal e ordinal representa um teste para a suposição de que as médias seguem um padrão linear com os escores escolhidos.

Formalmente, este é um teste da razão de verossimilhanças (LRT). A estatística de teste \(-2\log(\Lambda)\) é a diferença entre os desvios do modelo nominal e ordinal. Os graus de liberdade para o teste são derivados da diferença no número de parâmetros nos dois modelos. No modelo nominal, a variável \(X\) usa \(I-1\) df, os parâmetros \(\beta_2,\cdots,\beta_I\); enquanto no modelo ordinal \(X\) utiliza apenas um parâmetro \(\gamma_1\). Portanto, a estatística de teste é comparada a uma distribuição \(\chi^2_{I-2}\).

Observe que o uso de diferentes conjuntos de escores pode afetar o desvio residual do modelo ordinal e, portanto, a estatística e a conclusão do teste. Portanto, este teste é freqüentemente usado para ajudar a determinar se um determinado conjunto de escores é apropriado. Esse tipo particular de teste é chamado de teste de qualidade de ajuste.

As hipóteses para o teste de qualidade do ajuste são \[ H_0 \, : \, \log(\mu_i) = \gamma_0+\gamma_1 s_i \qquad \mbox{vs.} \qquad H_1 \, : \, \log(\mu_i) = \beta_0 + \beta_2 x_2 +\cdots + \beta_I x_I\cdot \] Em palavras, estamos testando a hipótese nula de que os escores fornecem uma explicação razoável para o padrão de contagens médias, contra a alternativa de que o padrão de médias não se ajusta à forma implícita nos escores. Deixar de rejeitar \(H_0\) implica que os escores parecem fazer um trabalho razoavelmente bom em capturar a tendência nas médias entre os níveis de \(X\).

Observe que isso não é o mesmo que dizer que os escores estão corretas. É bem possível que modelos baseados em vários conjuntos diferentes de escores possam fornecer ajustes “adequados” com base na falta de significância no teste. Observe também que este teste é diferente do LRT para a significância da variável de escores \(s\) no modelo ordinal, para o qual as hipóteses são \(H_0 \, : \, \log(\mu_i) =\gamma_0\) vs. \(H_1 \, : \, \log(\mu_i) =\gamma_0 +\gamma_1 s_i\). Essas duas questões são diferentes uma da outra; suas hipóteses nulas podem ser verdadeiras ou falsas em qualquer combinação.

Este procedimento se estende a modelos de associação para duas ou mais variáveis. A hipótese nula para o teste de ajuste é sempre representada pelo modelo que usa escores, enquanto a alternativa é o modelo em que um ou mais conjuntos de escores foram substituídos por variáveis indicadoras. Rejeitar \(H_0\) significa que a estrutura linear imposta pelos escores não é uma descrição adequada da associação para a qual são utilizados.

Quando existem múltiplas variáveis ordinais, é possível criar uma série de modelos afrouxando sucessivamente as suposições em uma determinada associação. Por exemplo, se \(X\) e \(Z\) são ambos ordinais, então em cada variável temos a opção de usar ou não escores em uma estrutura linear. Suponha que usamos um subscrito \(L\) para denotar que uma variável aparece em uma interação em sua forma de escores lineares. No termo de interação, poderíamos usar escores para ambas as variáveis \(X_L Z_L\), escores apenas para \(X\), \(X_L Z\), escores apenas para \(Z\), \(XZ_L\) ou deixar ambas as variáveis nominais \(XZ\). Os dois modelos intermediários que usam escores para apenas uma variável são generalizações do modelo linear-por-linear e são versões restritas do modelo totalmente nominal. Assim, os testes de comparação de modelos podem ser usados de várias maneiras diferentes, comparando \(H_0 \, : \, X_L Z_L\) com uma das várias alternativas possíveis, \(H_1 \, : \, X_L Z\) ou \(XZ_L\) ou \(XZ\). Também podemos testar \(H_0 \, : \, X_L Z\) ou \(H_0 \, : \, XZ_L\) contra \(H_1 \, : \, XZ\). Somente \(X_LZ\) e \(XZ_L\) não são diretamente comparáveis usando LRTs.

Testar pares de modelos por vez pode terminar com resultados ambíguos. Pode acontecer, por exemplo, que \(X_LZ_L\) seja rejeitado como modelo reduzido para \(X_LZ\) e/ou \(XZ_L\), mas não como modelo reduzido para \(XZ\). Este é um problema toda vez quando séries de testes são usados para tentar selecionar um único modelo. Uma alternativa é usar critérios de informação ou alguma outra medida para avaliar o ajuste de cada modelo individualmente, em vez de em comparação com algum outro modelo. Essas medidas são discutidas com mais detalhes no Capítulo 5.


Exemplo 4.11: Ideologia política.


Testamos o ajuste do modelo que assume uma associação linear na interação \(PI\) contra aquele que assume uma associação nominal. A hipótese nula é \[ H_0 \, : \, \mbox{A associação } PI \mbox{ é adequadamente descrita pelos escores 1-2-3-4-5} \] e a alternativa é simplesmente \[ H_1 \, : \, \mbox{Modelo } GP, GI, PI\cdot \] O teste é mostrado abaixo.

anova(mod.homo.lin1.PI, mod.homo, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: count ~ gender * party + gender * ideol + party * lin.score1
## Model 2: count ~ (gender + party + ideol)^2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         7     8.3986                     
## 2         4     3.2454  3   5.1532   0.1609


A estatística de teste é \(-2\log(\Lambda) = 5.15\) e o \(p\)-valor é 0.16, portanto, usando qualquer nível comum de, não rejeitaríamos \(H_0\): Podemos concluir que não há evidências suficientes de que os parâmetros extras fornecidos pela versão nominal das variáveis do modelo são necessárias.

Menos formalmente, isso significa que a suposição de razões de chances loglineares usando o conjunto de escores 1-2-3-4-5 fornece uma explicação razoável para a associação entre partido e ideologia.

Outras simplificações deste modelo podem ser consideradas. Podemos tentar usar os mesmos escores para descrever a associação \(GI\); ou seja, podemos especificar que as probabilidades formadas a partir de duas ideologias devem ter uma razão para homens versus mulheres que muda de forma log-linear com a distância entre as ideologias, conforme medido pelos escores. Podemos tentar os escores alternativas 0-2-3-4-6 para uma ou ambas as associações ordinais. E, finalmente, e talvez o mais interessante, podemos usar os escores de “extremismo”, 2-1-0-1-2, nas associações com partido, gênero ou ambos. Na verdade, podemos até usar esses escores além daquelas que medem as ideologias linearmente. Deixamos essas buscas como exercícios.



4.3 Regressão da taxa Poisson


Considere a seguinte afirmação: “A maioria dos acidentes automobilísticos ocorre a menos de oito quilômetros de casa”. Se essa afirmação for verdadeira, parece implicar que as estradas mais perigosas estão perto de sua casa - e perto de nossas casas e perto da casa de todos! Uma explicação que se pode inventar para essa afirmação é que nos tornamos tão familiarizados com as estradas perto de nossas próprias casas que relaxamos demais quando dirigimos perto de casa e, assim, nos expomos a graves perigos.

Uma explicação muito mais simples para esse resultado é que talvez a maior parte do trânsito ocorra a 8 km de casa! Se for assim, então a taxa com que colidimos com outras pessoas pode não ser diferente ou até menor, perto de casa do que em qualquer outro lugar, mas a exposição a acidentes em potencial é muito maior. Nesse contexto, contar o número de acidentes em diferentes distâncias de casa pode ser uma medida de risco menos significativa do que estimar a taxa de acidentes por milha percorrida.

Isso destaca um ponto importante levantado no início deste capítulo: uma suposição por trás do uso do modelo de Poisson para contagens é que tanto a intensidade ou taxa de ocorrência do evento quanto a oportunidade ou exposição para contagem são constantes para todas as observações disponíveis. Na Seção 4.2, investigamos modelos que permitem que uma contagem média mude. Implicitamente, estávamos assumindo que a exposição era constante para todas as observações.

Por exemplo, ao contar as bebidas consumidas em um sábado, a duração da observação “sábado” era constante para todos os participantes, de modo que variava apenas a taxa ou a intensidade do consumo de bebida. Em seguida, permitimos que as contagens médias, ou seja, taxa de consumo variassem usando um modelo de regressão de Poisson para relacioná-las com variáveis explicativas.

Em problemas onde a exposição não é constante entre todas as observações, modelar as contagens diretamente deixa as interpretações confusas. Por exemplo, se medimos alguns participantes no estudo de álcool por uma hora e outros por um dia e outros por uma semana ou um mês, deveríamos esperar que o número de bebidas consumidas entre os participantes fosse diferente simplesmente por causa da maneira como conduzimos nosso estudo. Isso pode obscurecer quaisquer efeitos que as variáveis explicativas possam ter. Teríamos, portanto, de desconsiderar o efeito da exposição, para podermos chegar ao cerne da questão que relaciona as variáveis explicativas à taxa de consumo.

Geralmente, “taxa” é definida como contagem média por unidade de exposição, por exemplo, bebidas por dia, animais presos por visita, acidentes por milha percorrida, árvores por hectare e assim por diante. Ou seja, \(R =\mu/t\), onde \(R\) é a taxa, \(t\) é a exposição e \(\mu\) é a contagem média ao longo de uma duração de exposição de \(t\). Quando todas as nossas contagens observadas têm a mesma exposição, modelando a contagem média em função das variáveis explicativas \(x_1,\cdots,x_p\) é o mesmo que modelar a taxa.

Quando as exposições variam, ainda podemos usar um modelo de regressão de Poisson para as médias, mas precisamos levar em conta a exposição no modelo. Isso é feito usando um modelo de regressão da taxa Poisson escrevendo o modelo como: \[ \log(r_i)=\log(\mu_i/t_i)+\beta_0+\beta_1 x_{1i}+\cdots+\beta_p x_{pi} \] ou \[ \log(\mu_i)=\log(t_i)+\beta_0+\beta_1 x_{1i}+\cdots+\beta_p x_{pi}, \] \(i=1,\cdots,n\).

Observe que o termo \(\log(t_i)\) é conhecido e não é multiplicado por um parâmetro. Ele simplesmente ajusta a média para cada observação diferente para explicar sua exposição. O termo \(\log(t)\) é chamado de offset (deslocamento), portanto, esse modelo às vezes é chamado de regressão de Poisson com deslocamentos. Uma consequência direta da equação acima é que \[ \mu = t\times \exp\big(\beta_0 +\beta_1 x_1 +\cdots+ \beta_p x_p \big), \] que mostra a influência direta da exposição na contagem média.

Existe uma semelhança entre taxas e proporções binomiais que às vezes causa confusão. Ambos são expressos como \(\{\mbox{contagens de eventos}\}/\{\mbox{oportunidades de eventos}\}\). Pode-se usar incorretamente o número de tentativas em um modelo binomial como uma exposição e aplicar uma regressão de taxa de Poisson ou, inversamente, pensar incorretamente na exposição em um problema de taxa de Poisson como um número de tentativas e aplicar uma regressão logística. No entanto, as duas situações são diferentes e a distinção é facilmente feita fazendo a seguinte pergunta: cada unidade de exposição pode produzir mais de um evento ou cada unidade de exposição deve produzir exatamente 0 ou 1 evento? Se mais de um evento for possível em uma única unidade de exposição, isso viola as suposições do modelo binomial e o modelo de taxa de Poisson é mais apropriado.

Por exemplo, não há nada que impeça que vários acidentes aconteçam em uma única milha percorrida e, de fato, um motorista muito bêbado ou descuidado pode passar por esse infortúnio. Um modelo binomial seria inadequado neste caso. Por outro lado, cada chute ou lance livre deve ser um sucesso ou uma falha e não se pode contar mais de um evento bem-sucedido por tentativa. Para estes problemas deve ser utilizado o modelo de regressão logística.

Os parâmetros \(\beta_0,\cdots,\beta_p\) são estimados usando a estimativa de máxima verossimilhança como antes, mas seu interpreto é em termos de taxa e não de média. Por exemplo, \(\beta_0\) é a verdadeira taxa de evento quando todas as variáveis explicativas são definidas como 0 e \(\beta_1\) é a mudança na taxa de ocorrência por unidade de aumento em \(x_1\). Observe que as unidades de medida para a taxa são unidades de contagem por unidade de exposiçã, por exemplo, acidentes por milha percorrida.

Às vezes é útil converter a exposição para outras unidades para criar uma taxa mais interpretável. Por exemplo, se estivermos analisando acidentes e a exposição for dada como milhas percorridas, o número médio de acidentes em 1 milha deve ser extremamente pequeno. Se preferirmos, podemos optar por reescrever a taxa como uma quantidade mais interpretável, como “acidentes por 1.000 milhas”, dividindo as milhas de exposição por 1.000 antes de aplicar o modelo.

Dado um conjunto de estimativas \(\widehat{\beta}_0,\cdots,\widehat{\beta}_p\), taxas previstas para valores dados \(x_1,\cdots, x_p\) são encontrados a partir de \[ \widehat{R}=\exp\big( \widehat{\beta}_0+\widehat{\beta}_1 x_1+\cdots +\widehat{\beta}_p x_p \big)\cdot \] Uma contagem prevista para uma determinada exposição \(t\) é encontrada simplesmente multiplicando a taxa pela exposição, \(\widehat{\mu} = t\times \widehat{R}\).


Exemplo 4.12: Resposta de postura de besouros à aglomeração.


Tauber et ai. (1996) descrevem um experimento examinando os efeitos da aglomeração nas propriedades reprodutivas de uma certa espécie de besouro da folha, Galerucella nymphaeae. Gaiolas de tamanho fixo contendo 1 fêmea e 1 macho ou 5 fêmeas e 5 machos foram mantidas em câmaras de temperatura a 21²C ou 24²C. Entre as medições feitas em cada gaiola está o número de massas de ovos produzidas pelas fêmeas na gaiola ao longo da vida das fêmeas.

Usamos a regressão Poisson para examinar se a aglomeração e/ou a temperatura influenciam o número de massas de ovos produzidas por cada fêmea. Uma característica complicada é que em gaiolas com 5 fêmeas, não há uma maneira fácil de identificar quais fêmeas colocaram uma determinada massa de ovos, então devemos considerar “massas de ovos por fêmea” como a taxa de resposta de interesse.

Os dados são inseridos usando o código abaixo. A variável de temperatura Temp assume os valores 21 e 24; TRT, que é a variável de tratamento de aglomeração, assume valores “I” para gaiolas com fêmeas individuais e “G” para gaiolas com grupos de 5; Unit identifica a gaiola dentro de uma combinação de temperatura e nível de aglomeração, existem 35 gaiolas individuais e 7 gaiolas de grupo em cada temperatura; NumEggs conta as massas de ovos colocadas em cada caixa ao longo da vida das fêmeas. Criamos uma variável, “females” para representar o número de fêmeas na gaiola que servirá como nossa variável de exposição.

# Read in beetle egg mass data
eggdata <- read.table(file = "http://leg.ufpr.br/~lucambio/ADC/BeetleEggCrowding.txt", 
                      header = TRUE, stringsAsFactors = TRUE)
# NOTE 2021-10-19: stringsAsFactors = TRUE was added due to changes in R version 4.0.0
eggdata$females <- ifelse(test = eggdata$TRT == "I", yes = 1, no = 5)
head(eggdata)
##   Temp TRT Unit NumEggs females
## 1   21   I    1       8       1
## 2   21   I    2       0       1
## 3   21   I    3       3       1
## 4   21   I    4       5       1
## 5   21   I    5       1       1
## 6   21   I    6       3       1


Como resumo inicial, as taxas médias amostrais, massas de ovos por fêmea para cada grupo foram calculadas, produzindo os resultados 3.54 para TRT=I, Temp=21, 4.34 para TRT=I, Temp=24, 0.51 para TRT=G, Temp=21 e 1.46 para TRT=G, Temp=24.

O modelo que vamos ajustar é \[ \log(\mu ) = \log(t)+ \beta_0+\beta_1 Temp+\beta_2 TRTI+ \beta_3 Temp\times TRTI, \] onde TRTI é o indicador para TRT=I e \(t\) é o número de fêmeas na gaiola, 1 ou 5. Em glm(), o deslocamento é incorporado ao modelo usando o argumento deslocamento conforme mostrado abaixo:

# Getting mean rate per female for each combination of TRT and Temp
aggregate( NumEggs ~ TRT + Temp, data = eggdata, FUN = mean, subset = TRT == "I")
##   TRT Temp  NumEggs
## 1   I   21 3.542857
## 2   I   24 4.342857
aggregate( NumEggs ~ TRT + Temp, data = eggdata, FUN = mean, subset = TRT == "G")$NumEggs/5
## [1] 0.5142857 1.4571429
# Fit model using temperature, crowding level, and interaction
eggmod1 <- glm(NumEggs ~ Temp*TRT, family = poisson(link = "log"), offset = log(females), data = eggdata)
round(summary(eggmod1)$coefficients, digits = 3)
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -7.955      2.125  -3.743    0.000
## Temp           0.347      0.091   3.799    0.000
## TRTI           7.795      2.314   3.369    0.001
## Temp:TRTI     -0.279      0.100  -2.796    0.005


A partir dos coeficientes positivos nos termos de efeito principal, fica evidente que a taxa de postura é maior na temperatura mais alta e nas gaiolas individuais. O coeficiente de interação negativo indica que a razão de taxas entre gaiolas individuais e coletivas é maior na temperatura mais baixa.

As taxas podem ser calculadas a partir desses coeficientes. Por exemplo, a taxa média para TRT=I, Temp=21 é \[ \exp\big( -7.955 + 0.347(21) + 7.795(1) - 0.279(21)(1)\big) = 3.54, \] correspondendo à estimativa anterior.

Uma maneira fácil de obter taxas previstas é usar predict() no conjunto de dados original. Isso produz uma previsão para cada gaiola, que abreviamos para uma para cada combinação de TRT e Temp na saída abaixo. Observe que isso retorna a contagem média prevista, que devemos dividir pela exposição para formar uma taxa. As contagens médias para outros níveis de exposição podem ser estimadas multiplicando essas taxas pela nova exposição.

# Get the true egg-laying rate per female:
# Predicted means using observed t. 
newdata <- data.frame(TRT = c("G", "I", "G", "I"), Temp = c(21, 21, 24, 24), females = c(5, 1, 5, 1))
mu.hat <- round(predict(object = eggmod1, newdata = newdata, type = "response"), 2)
# Convert means to rates
data.frame(newdata, mu.hat, rate = mu.hat/newdata$females)
##   TRT Temp females mu.hat  rate
## 1   G   21       5   2.57 0.514
## 2   I   21       1   3.54 3.540
## 3   G   24       5   7.29 1.458
## 4   I   24       1   4.34 4.340


A razão pela qual as taxas estimadas de postura rate são idênticas às nossas estimativas amostrais é que o modelo que estamos usando está saturado. As quatro taxas observadas são totalmente explicadas pelos 4 df representado por \(\beta_0\), \(\beta_1\), \(\beta_2\) e \(\beta_3\). Essas taxas previstas geralmente não correspondem às taxas médias observadas em um modelo reduzido.

Abaixo está um resumo dos testes da razão de verossimilhanças (LRTs) para os termos do modelo:

# Type II ANOVA Tests
library(car)
Anova(eggmod1)
## Analysis of Deviance Table (Type II tests)
## 
## Response: NumEggs
##          LR Chisq Df Pr(>Chisq)    
## Temp       10.842  1   0.000992 ***
## TRT       132.994  1  < 2.2e-16 ***
## Temp:TRT    8.450  1   0.003650 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O termo de interação Temp:TRT é significativo, indicando que as razões das médias entre os dois tratamentos de aglomeração não são as mesmas em ambas as temperaturas \[ -2\log(\Lambda) = 8.4, \] df = 1, \(p\)_valor = 0.004.

Os intervalos de confiança de Wald para a verdadeira taxa de postura por fêmea são produzidos usando a codificação usual. Eles são representados junto com as taxas estimadas na figura abaixo. O efeito da aglomeração é aparente: menos massas de ovos são colocadas por fêmea nas gaiolas com a condição lotada. O efeito de interação é explorado com mais detalhes no Exercício 20.

# Matrix of linear combination coefficients for predicted rates at each combination
coef.mat <- as.matrix(rbind(c(1,21,1,21), c(1,24,1,24), c(1,21,0,0), c(1,24,0,0)))
# Wald inferences using multcomp package
library(multcomp)
wald <- glht(eggmod1, linfct = coef.mat)
# Defaults use multiplicity adjustment for simultaneous confidence level
summary(wald)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = NumEggs ~ Temp * TRT, family = poisson(link = "log"), 
##     data = eggdata, offset = log(females))
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)    
## 1 == 0  1.26493    0.08980  14.086   <2e-16 ***
## 2 == 0  1.46853    0.08111  18.105   <2e-16 ***
## 3 == 0 -0.66498    0.23570  -2.821   0.0190 *  
## 4 == 0  0.37648    0.14003   2.689   0.0284 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
round(exp(confint(wald)$confint), digits = 2)
##   Estimate  lwr  upr
## 1     3.54 2.83 4.43
## 2     4.34 3.55 5.32
## 3     0.51 0.29 0.93
## 4     1.46 1.03 2.07
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 2.490844
# Options to get unadjusted (univariate) tests and CIs
summary(wald, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = NumEggs ~ Temp * TRT, family = poisson(link = "log"), 
##     data = eggdata, offset = log(females))
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)    
## 1 == 0  1.26493    0.08980  14.086  < 2e-16 ***
## 2 == 0  1.46853    0.08111  18.105  < 2e-16 ***
## 3 == 0 -0.66498    0.23570  -2.821  0.00478 ** 
## 4 == 0  0.37648    0.14003   2.689  0.00718 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
wald.ci <- exp(confint(wald, calpha = qnorm(0.975))$confint)
round(wald.ci, digits = 2)
##   Estimate  lwr  upr
## 1     3.54 2.97 4.22
## 2     4.34 3.70 5.09
## 3     0.51 0.32 0.82
## 4     1.46 1.11 1.92
## attr(,"conf.level")
## [1] 0.95
## attr(,"calpha")
## [1] 1.959964
# Use mcprofile to get LR confidence intervals.
library(mcprofile)
LRCI <- mcprofile(eggmod1, CM = coef.mat)
exp(confint(LRCI, level = 0.95))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   single-step 
## 
##    Estimate lower upper
## C1    3.543 2.808 4.396
## C2    4.343 3.524 5.280
## C3    0.514 0.268 0.878
## C4    1.457 1.006 2.026
exp(confint(LRCI, level = 0.95, adjust = "none"))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##    Estimate lower upper
## C1    3.543 2.955  4.20
## C2    4.343 3.689  5.07
## C3    0.514 0.312  0.79
## C4    1.457 1.093  1.89
wald.calc <- wald(LRCI)
exp(confint(wald.calc, level = 0.95, adjust = "none"))
## 
##    mcprofile - Confidence Intervals 
## 
## level:        0.95 
## adjustment:   none 
## 
##    Estimate lower upper
## C1    3.543 2.971 4.225
## C2    4.343 3.705 5.091
## C3    0.514 0.324 0.816
## C4    1.457 1.107 1.917
# Plot mean rates and confidence intervals: Color
plot(x = c(21,24,21,24), y = wald.ci[,1], xlim = c(20,25.5), ylim = c(0,5), 
     xlab = "Temperature", ylab = "Mean Egg Masses per Female")
segments(x0 = 21, x1 = 24, y0 = wald.ci[1,1], y1 = wald.ci[2,1], col = "blue", lwd = 2)
segments(x0 = 21, x1 = 24, y0 = wald.ci[3,1], y1 = wald.ci[4,1], col = "red", lwd = 2)
segments(x0 = 21, x1 = 21, y0 = wald.ci[1,2], y1 = wald.ci[1,3], col = "blue", lwd = 2)
segments(x0 = 21, x1 = 21, y0 = wald.ci[3,2], y1 = wald.ci[3,3], col = "red", lwd = 2)
segments(x0 = 24, x1 = 24, y0 = wald.ci[2,2], y1 = wald.ci[2,3], col = "blue", lwd = 2)
segments(x0 = 24, x1 = 24, y0 = wald.ci[4,2], y1 = wald.ci[4,3], col = "red", lwd = 2)
legend(x = 24, y = 1, legend = c("Individual", "Crowded"), cex = 0.9, col = c("blue", "red"), 
       lwd = 2, title = "TRT Level", bty = "n")
grid()


A figura acima mostra o gráfico das taxas médias de produção de massa de ovos e intervalos de confiança para todas as combinações de temperatura de tratamento no exemplo de postura de besouros. As gaiolas individuais possuem uma fêmea e um macho; gaiolas lotadas têm cinco fêmeas e cinco machos. As linhas verticais representam intervalos de confiança Wald de 95% para a taxa.


Observe que às vezes acontece que uma taxa pode depender da quantidade de exposição. Por exemplo, as taxas de criminalidade como crimes por 10.000 habitantes, costumam ser mais altas em cidades com populações maiores. Isso também estava implícito no exemplo do ovo de besouro acima, onde o tratamento de aglomeração estava relacionado ao número de fêmeas em cada caixa, que também era a exposição para a taxa de ovos/fêmea.

Como fica claro no exemplo, certamente é possível usar a variável de exposição como outra variável explicativa na regressão, além de usá-la como compensação. Ou seja, podemos definir \(x_{p+1} = t\) ou \(\log(t)\), se preferir e adicionar \(\beta_{p+1} x_{p+1}\) à equação. O uso duplo do deslocamento não causa problemas para a estimativa do parâmetro do modelo como aqueles vistos no exemplo de Ideologia Política na Seção 4.2.6, porque nenhum parâmetro é estimado para o deslocamento, em essência, o parâmetro é fixado em 1.


4.4 Inflação de zeros


Muitas populações consistem em subgrupos de diferentes tipos de unidades. Uma situação particularmente comum é que uma população consiste em duas classes: indivíduos “suscetíveis” e “imunes”. Por exemplo, máquinas de ressonância magnética (MRI) são equipamentos médicos grandes e muito caros. Alguns hospitais, mas não todos, têm uma máquina de ressonância magnética. Se um hospital tem uma, certamente não fica parada por longos períodos de tempo. Suponha que os administradores do hospital sejam solicitados em uma pesquisa a relatar o número de ressonâncias magnéticas que eles realizam em um mês.

Hospitais sem máquina de ressonância magnética não estão realizando nenhuma ressonância magnética, enquanto o restante relatará alguma contagem diferente de 0. Como um exemplo diferente, suponha que registremos o número de peixes capturados em vários lagos em viagens de pesca de 4 horas em Minnesota. Alguns lagos em Minnesota são muito rasos para que os peixes sobrevivam ao inverno, então pescar nesses lagos não produzirá capturas. Por outro lado, mesmo em um lago onde os peixes são abundantes, podemos ou não pegar algum peixe devido às condições ou à nossa própria competência. Assim, o número de peixes capturados será zero se o lago não comportar peixes e será zero, um ou mais se comportar.

Essas situações são exemplos em que as respostas vêm de misturas de distribuições e, em particular, são casos em que uma das distribuições colocam toda a sua massa em zero. Ou seja, os “suscetíveis” na população retornam uma contagem de acordo com alguma distribuição, enquanto todos os “imunes” retornam um zero.

Os dados resultantes desse tipo de mistura têm a propriedade de haver mais contagens de zero do que seria previsto por um único modelo distributivo. Isso é conhecido como inflação zero. Em geral, seja \(\pi\) a proporção de imunes na população e as contagens de suscetíveis tenham uma distribuição \(Po(\mu )\). Essa mistura é chamada de modelo de Poisson inflado de zero. A figura abaixo mostra histogramas das distribuições de Poisson infladas de zero com valores variáveis de \(\pi\) e \(\mu\), observe que \(\pi = 0\) é a distribuição de Poisson ordinária.

# Zero Inflated Poisson with probability of immune=pi, mean for Poisson=mu
# Create grid of pi, mu values, add pmf calculation
all <- expand.grid(y = c(0:20), pi = c(0.5,0.2,0), mu = c(10,2))
all <- data.frame(all, p = (1 - all$pi) * dpois(x = all$y, lambda = all$mu) + ifelse(all$y == 0, 
                                                                                     yes = all$pi, no = 0))
head(all)
##   y  pi mu            p
## 1 0 0.5 10 0.5000227000
## 2 1 0.5 10 0.0002269996
## 3 2 0.5 10 0.0011349982
## 4 3 0.5 10 0.0037833275
## 5 4 0.5 10 0.0094583187
## 6 5 0.5 10 0.0189166374
# Make plots
par(mfrow = c(2,3))
plot(x = all$y[1:21], y = all$p[1:21], type = "h", xlim = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0.5, ", mu, "=10")))
grid()
plot(x = all$y[22:42], y = all$p[22:42], type = "h", xlim = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0.2, ", mu, "=10")))
grid()
plot(x = all$y[43:63], y = all$p[43:63], type = "h", xlim = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0, ", mu, "=10")))
grid()
plot(x = all$y[64:84], y = all$p[64:84], type = "h", xlim = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0.5, ", mu, "=2")))
grid()
plot(x = all$y[85:105], y = all$p[85:105], type = "h", xlim = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0.2, ", mu, "=2")))
grid()
plot(x = all$y[106:126], y = all$p[106:126], type = "h", xlim  = c(0,20), xlab = "y", ylab = "P(Y=y)", 
     lty = "solid", lwd = 3, main = expression(paste(pi, "=0, ", mu, "=2")))
grid()


A inflação de zeros é óbvia quando \(\mu\gg 0\) como é o caso na linha superior. Dificilmente alguém seria tentado a usar uma única distribuição de Poisson como modelo para dados com um histograma como os mostrados para \(\pi> 0\). No entanto, quando a contagem média entre os suscetíveis pode colocar uma massa apreciável em zero ou um, linha inferior, então o a presença de inflação zero pode passar despercebida se a natureza do problema não o sugerir antecipadamente. A presença de variáveis explanatórias também pode obscurecer a distinção entre as classes, porque a contagem média de susceptíveis e as probabilidades de imunidade podem variar muito entre os membros da população. Na prática, a inflação zero é muitas vezes descoberta apenas ao ajustar um modelo e encontrando o ajuste ruim usando as técnicas descritas no Capítulo 5.


Modelos corretivos


Quando se sabe ou acredita-se que uma população consiste em dois subgrupos, conforme descrito acima, ou quando um excesso de contagens zero aparece inesperadamente, os modelos de Poisson comuns provavelmente produzirão ajustes ruins aos dados. Vários modelos foram propostos para explicar os zeros em excesso. O modelo principal é o modelo de Poisson inflado de zero (ZIP) de Lambert (1992).

O modelo ZIP especifica que a imunidade é uma variável aleatória binária que pode depender de variáveis explanatórias denotadas coletivamente por \(z\), enquanto as contagens entre os suscetíveis seguem uma distribuição de Poisson cuja média pode depender de variáveis explanatórias coletivamente denotadas por \(x\).

Desta forma, \[ \begin{array}{ccccc} Y & = & 0 & \mbox{com probabilidade} & \pi(z) \\ Y &\sim & Po(\mu(x)) & \; \mbox{com probabilidade} & 1-\pi(x) \end{array} \]

Normalmente, um modelo logístico é assumido para \(\pi(z)\) \[ logit\big(\pi(z)\big) = \gamma_0+\gamma_1 \, z_1 + \cdots + \gamma_r \, z_r, \] enquanto o modelo com ligação logaritmo usual é usado para \(\mu(x)\) \[ \log(\big(\mu(x)\big) = \beta_0+\beta_1 \, x_1 + \cdots + \beta_p \, x_p, \] onde \(\gamma_0,\gamma_1,\cdots,\gamma_r\) e \(\beta_0,\beta_1,\cdots,\beta_p\) são parâmetros de regressão desconhecidos.

É admissível que \(x = z\) de forma que as mesmas variáveis explicativas sejam utilizadas tanto para distinguir as classes quanto para modelar a média dos suscetíveis, mas isso não é obrigatório. Geralmente, \(x\) e \(z\) podem ser conjuntos de variáveis arbitrariamente semelhantes ou diferentes. Pode-se mostrar que

  1. \(P(Y = 0 \, | \, x, z) = \pi(z) + \big(1-\pi(z)\big)e^{-\mu(x)}\),

  2. \(\mbox{E}(Y \, | \, x, z) = \mu(x)\big(1-\pi(z)\big)\), e

  3. \(\mbox{Var}(Y \, | \, x, z) = \mu(x)\big(1-\pi(z)\big)\big(1 + \mu(x)\pi(z)\big)\).

Assim, a média de \(Y\) quando \(\pi(z) > 0\) é sempre menor que \(\mu(x)\), mas a variância de \(Y\) é maior que sua média por um fator de \(1 + \mu(x)\pi (z)\).

Um parente do modelo ZIP, chamado de modelo hurdle (Mullahy, 1986), difere em que uma contagem é sempre maior que zero sempre que algum obstáculo real ou conceitual é ultrapassado. O exemplo de ressonância magnética que iniciou esta seção é deste tipo, porque se um hospital gastou a instalação de um aparelho de ressonância magnética, o “obstáculo”, certamente o usará. Neste caso, todos os zeros correspondem a imunes, àquelas unidades que não ultrapassaram o obstáculo de instalação de uma máquina, enquanto o modelo de contagem para susceptíveis assume valores começando em 1 em vez de 0.

Isso é feito impondo probabilidade 0 em \(Y = 0\) e, em seguida, aumentando todas as outras probabilidades proporcionalmente para que elas mais uma vez somem 1. A distribuição resultante é chamada de “Poisson truncado à esquerda”, porque a cauda esquerda da distribuição \((Y = 0)\) é cortado da distribuição normal de Poisson.

O modelo logístico para a probabilidade de imune e o logaritmo da média da distribuição de Poisson de contagens truncada à esquerda para susceptíveis são usados no modelo hurdle como no modelo ZIP. Ambos os modelos podem ser estimados por máxima verossimilhança usando funções do pacote pscl. A sintaxe para o ajuste de ambos os modelos é muito semelhante e é descrita em detalhes por Zeileis et al. (2008).

O ajuste do modelo pode ser mais fácil com o modelo hurdle, porque fica imediatamente claro que todos os zeros são imunes. Então a modelagem da função média é baseada no subconjunto de dados com \(Y > 0\) e é independente da probabilidade de imunes. No entanto, o modelo ZIP parece ser o preferido na maioria das áreas. Nossa preferência geralmente é determinada por qual modelo para os imunes concorda melhor com a fonte dos dados.

Se uma contagem zero define um imune, como no exemplo da ressonância magnética, um obstáculo faz mais sentido. Se os suscetíveis às vezes podem fornecer contagens zero, como no exemplo da pesca, isso aponta para um modelo ZIP. Se não estiver claro se os suscetíveis podem produzir contagens zero, então considerações empíricas, como medidas de ajuste, podem prevalecer. Cconsulte o Capítulo 5 para obter detalhes.


Exemplo 4.13: Resposta de postura de besouro à aglomeração.


Um dos objetivos do experimento de Tauber et al. (1996) foi observar uma resposta particular chamada diapausa nas fêmeas. A diapausa ocorre em insetos quando eles cessam temporariamente o desenvolvimento biológico de alguma forma como resposta às condições ambientais. Nesse caso, os pesquisadores suspeitaram que os ciclos reprodutivos das fêmeas poderiam parar em resposta à temperatura ou aglomeração percebida. Portanto, eles esperavam que as fêmeas entrassem em diapausa e não produzissem ovos com maior probabilidade para algumas condições do que para outras. O fato de algumas fêmeas nas gaiolas coletivas poderem entrar em diapausa enquanto outras não cria uma complicação porque aquelas gaiolas de fêmeas mistas ainda produzirão ovos, mas em uma taxa reduzida.

Assim, estimar a probabilidade de que uma determinada fêmea produza zero massas de ovos em uma gaiola durante toda a duração do estudo é um exemplo de um problema de teste de grupo; consulte Bilder, 2009 para uma introdução ao teste de grupo. Isso está além do nosso escopo atual, portanto, para este exemplo, nos limitaremos às gaiolas individuais.

Usaremos um modelo ZIP para determinar se a temperatura afeta a probabilidade de diapausa ou o número médio de ovos postos por fêmeas que não estão em diapausa.

Primeiro, examinamos a extensão da inflação de zeros. Usamos o seguinte código para reduzir os dados apenas para as gaiolas individuais e mostrar esses dados:

eggdata <- read.table( file = "http://leg.ufpr.br/~lucambio/ADC/BeetleEggCrowding.txt", 
                       header = TRUE, stringsAsFactors = TRUE)
# NOTE 2021-10-19: stringsAsFactors = TRUE was added due to changes in R version 4.0.0
eggdata2 <- eggdata[eggdata$TRT == "I",]
# Plot histograms of egg mass counts for each group  
stripchart(x = NumEggs ~ Temp, method = "stack", data = eggdata2, ylab = "Temperature", 
           xlab = "Number of egg masses", at = c(1.1,1.8))
grid()


Em nosso exemplo anterior, as médias estimadas foram 3.54 para TRT = I, Temp = 21 e 4.34 para TRT = I, Temp = 24. A figura acima mostra que há consideravelmente mais zeros do que seria esperado das distribuições de Poisson com essas médias; as probabilidades estimadas são \(P(Y = 0) = e^{-\widehat{\mu}} = 0.029\) e 0.013, respectivamente. Além disso, há alguns zeros a mais em Temp = 24 do que em Temp = 21, 16 vs. 11, e as médias das contagens diferentes de zero parecem diferir entre as temperaturas.

No geral, o modelo ZIP é preferido para esta configuração porque uma fêmea que não está em diapausa pode, no entanto, produzir uma contagem zero; por exemplo, devido a algum outro processo biológico, como infertilidade ou morte precoce. Nosso modelo ZIP usa uma regressão logística para descrever a probabilidade de “imunidade” à postura de ovos, ou seja, diapausa e uma regressão de Poisson nas contagens de todas as fêmeas que não são “imunes”. Em ambos os modelos, permitiremos que as probabilidades ou as contagens médias variem de acordo com a temperatura.

Ou seja, estamos ajustando \[ \log\big(\pi/(1+\pi)\big) = \gamma_0+\gamma_1 \, \mbox{Temp} \]

e \[ \log(\mu)=\beta_)+\beta_1 \, \mbox{Temp}\cdot \]

O modelo ZIP é ajustado usando a função zeroinfl() do pacote pscl. O principal argumento necessário para zeroinfl() é a fórmula do modelo. Em geral, a sintaxe é \[ \mbox{formula} = y \sim x_1 + x_2 + \cdots \, | \, z_1 + z_2 + \cdots \] ajusta o modelo de contagem às variáveis \(x\) e o modelo de probabilidade de imunidade às variáveis \(z\). Se as mesmas variáveis forem usadas em ambos os modelos, como é o caso dos besouros, o atalho \(y \sim x_1 + x_2 \cdots\) pode ser usado.

Primeiro ajustamos o modelo no qual tanto a probabilidade de imunidade quanto as médias podem variar de acordo com a temperatura:

library(pscl)
zip.mod.tt <- zeroinfl(NumEggs ~ Temp | Temp, dist = "poisson", data = eggdata2)
summary(zip.mod.tt)
## 
## Call:
## zeroinfl(formula = NumEggs ~ Temp | Temp, data = eggdata2, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1688 -0.9659 -0.5211  0.8134  3.2601 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4631     0.9249  -1.582 0.113671    
## Temp          0.1476     0.0407   3.626 0.000288 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -5.1843     3.7897  -1.368    0.171
## Temp          0.2088     0.1671   1.249    0.212
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 10 
## Log-likelihood: -187.2 on 4 Df


A saída da chamada zeroinfl() fornece as estimativas dos parâmetros para o modelo de média de Poisson primeiro, seguidas por aquelas para o modelo de probabilidade logística. A função média estimada é \[ \widehat{\mu} = \exp(-1.46+0.15 \, \mbox{Temp}), \] e a probabilidade estimada de imunidade \[ \widehat{\pi} = \exp( -5.20 + 0.21 \, \mbox{Temp})/(1 + \exp( -5.20 + 0.21\, \mbox{Temp}))\cdot \] Os testes de Wald para efeitos de temperatura sugerem significância no modelo médio, mas não no modelo de probabilidade. Exploramos isso mais tarde.

Não há método anova() para objetos produzidos por zeroinfl(), então devemos usar outras ferramentas para testes da razão de verossimilhanças (LRTs). O pacote lmtest fornece uma variedade de ferramentas de teste e diagnóstico para modelos lineares, algumas das quais foram estendidas para lidar com objetos da classe zeroinfl. Em particular, a função lrtest() do pacote lmtest pode produzir LRTs para comparar modelos aninhados.

Usamos aqui esses desenvolvimentos para comparar o modelo completo fornecido em zip.mod.tt com os ajustes de três modelos reduzidos: (1) sem efeito de temperatura em \(\pi\),

  1. sem efeito de temperatura em \(\mu\) e

  2. sem efeito de temperatura em qualquer lugar.

# Fit ZIP Models: Temp in mean only
zip.mod.t0 <- zeroinfl(NumEggs ~ Temp | 1, dist = "poisson", data = eggdata2)
# Fit ZIP Models: Temp in probability only
zip.mod.0t <- zeroinfl(NumEggs ~ 1 | Temp, dist = "poisson", data = eggdata2)
# Fit ZIP Models: no explanatories
zip.mod.00 <- zeroinfl(NumEggs ~ 1 | 1, dist = "poisson", data = eggdata2)
# LRTs of each reduced model against largest model.
library(lmtest)
lrtest(zip.mod.t0, zip.mod.tt)
## Likelihood ratio test
## 
## Model 1: NumEggs ~ Temp | 1
## Model 2: NumEggs ~ Temp | Temp
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   3 -187.97                     
## 2   4 -187.18  1 1.5866     0.2078
lrtest(zip.mod.0t, zip.mod.tt)
## Likelihood ratio test
## 
## Model 1: NumEggs ~ 1 | Temp
## Model 2: NumEggs ~ Temp | Temp
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   3 -193.82                         
## 2   4 -187.18  1 13.294  0.0002662 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(zip.mod.00, zip.mod.tt)
## Likelihood ratio test
## 
## Model 1: NumEggs ~ 1 | 1
## Model 2: NumEggs ~ Temp | Temp
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -194.58                         
## 2   4 -187.18  2 14.808  0.0006087 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


A primeira LRT mostra que reduzir do modelo completo zip.mod.tt, com efeitos de temperatura tanto na média quanto na probabilidade, para zip.mod.t0, com efeitos de temperatura apenas na média, não causa uma queda significativa no valor observado da função de verossimilhança, \(-2\log(\Lambda ) = 1.59\), \(df=1\), \(p\)-valor = 0.21. No entanto, remover o efeito da temperatura da média resulta em um ajuste significativamente pior, \(-2\log(\Lambda ) = 13.29\), \(df=1\), \(p\)-valor = 0.0002. Portanto, consideraríamos o zip.mod.t0 como um bom modelo de trabalho, dependendo da avaliação de diagnóstico, conforme discutido no Capítulo 5.

Resumimos este modelo de várias maneiras usando o código abaixo para este exemplos. Primeiro, estimamos o efeito da temperatura como a razão entre as contagens médias de massa de ovos para cada aumento de 1℃ na temperatura. Também encontramos um intervalo de confiança de 95% para essa quantidade. A razão estimada é \(\exp(0.1476) = 1.16\) e o intervalo de confiança vai de 1.07 a 1.25. Assim, com 95% de confiança estimamos que o número médio de posturas de fêmeas fora da diapausa aumenta entre 7% e 25% para cada 1℃ de aumento na temperatura.

# Summary of chosen model
summary(zip.mod.t0)
## 
## Call:
## zeroinfl(formula = NumEggs ~ Temp | 1, data = eggdata2, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0998 -1.0318 -0.6531  0.8632  3.2015 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.44897    0.92249  -1.571 0.116250    
## Temp         0.14700    0.04061   3.620 0.000295 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.4720     0.2466  -1.914   0.0556 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 5 
## Log-likelihood:  -188 on 3 Df
# Estimated ratio for means 1C apart (Temp-mean is 2nd parameter)
exp(coef(zip.mod.t0)[2])
## count_Temp 
##   1.158352
#Confidence interval for all parameters
confint(zip.mod.t0) 
##                         2.5 %     97.5 %
## count_(Intercept) -3.25701992 0.35908153
## count_Temp         0.06740321 0.22659318
## zero_(Intercept)  -0.95529354 0.01123516
# Confidence interval for ratio of means 1C apart
round(exp(confint(zip.mod.t0)[2,]), digits = 2)
##  2.5 % 97.5 % 
##   1.07   1.25


Em seguida, calculamos estimativas e intervalos individuais de confiança Wald de 95% para a contagem média condicional da massa de ovos em cada grupo de temperatura, condicionada em insetos suscetíveis e para a probabilidade geral de imunidade (diapausa).

O código usa no pacote multcomp a função glht. Os resultados são: \(\widehat{\mu} = 5,14\), 4.31 \(< \mu < 6.15\) a 21℃, \(\widehat{\mu} = 8.00\), 6.82 \(< \mu < 9.38\) a 24℃ e \(\widehat{\pi} = 0.38\), 0.28 \(<\pi <\) 0.50. Assim, mais ovos são postos pelas fêmeas fora da diapausa a 24℃ do que a 21℃ e aproximadamente 38% das fêmeas entram em diapausa independentemente da temperatura.

# Wald Confidence Intervals using glht() from multcomp.
# Matrix of linear combination coefficients for predicted rates at each combination
coef.mat <- rbind(c(1,21,0), c(1,24,0), c(0,0,1))
library(multcomp)
wald <- glht(zip.mod.t0, linfct = coef.mat)
# Options to get unadjusted (univariate) tests and CIs
# Tests
summary(wald, test = univariate())
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: zeroinfl(formula = NumEggs ~ Temp | 1, data = eggdata2, dist = "poisson")
## 
## Linear Hypotheses:
##        Estimate Std. Error z value Pr(>|z|)    
## 1 == 0  1.63799    0.09081  18.037   <2e-16 ***
## 2 == 0  2.07899    0.08123  25.594   <2e-16 ***
## 3 == 0 -0.47203    0.24657  -1.914   0.0556 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)
# Confidence Intervals for means
wald.ci <- exp(confint(wald, calpha = qnorm(0.975))$confint)
round(wald.ci[c(1,2),], digits = 2)
##   Estimate  lwr  upr
## 1     5.14 4.31 6.15
## 2     8.00 6.82 9.38
# Confidence Intervals for P(Immune)
round(wald.ci[3,]/(1 + wald.ci[3,]), digits = 2)
## Estimate      lwr      upr 
##     0.38     0.28     0.50


Para comparar o tratamento de zeros nesses modelos ZIP com o modelo de Poisson ajustado na Seção 4.3, calculamos o número de contagens de zero previstas por todos os quatro modelos ZIP e pelo modelo de Poisson. Isso é feito calculando a probabilidade estimada de uma resposta zero para cada observação com base em suas variáveis explicativas e, em seguida, somando essas probabilidades em todo o conjunto de dados.

O número de zeros observados foi 27. Cada modelo ZIP se aproxima muito desse total, com totais esperados variando de 27.04 a 27.11. O modelo de Poisson que não leva em conta a inflação zero prevê apenas 1.47 zeros. Claramente, o modelo de Poisson não é adequado para esses dados. No entanto, os diagnósticos feitos no Exercício 33 do Capítulo 5 revelam que o modelo ZIP final fornecido em zip.mod.t0 também apresenta problemas!

# Computing ratio of means and confidence interval
# Making this flexible to choosing number of degrees for comparison (1C here)
deg <- 1
# Summary of chosen model
summary(zip.mod.t0)
## 
## Call:
## zeroinfl(formula = NumEggs ~ Temp | 1, data = eggdata2, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0998 -1.0318 -0.6531  0.8632  3.2015 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.44897    0.92249  -1.571 0.116250    
## Temp         0.14700    0.04061   3.620 0.000295 ***
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -0.4720     0.2466  -1.914   0.0556 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Number of iterations in BFGS optimization: 5 
## Log-likelihood:  -188 on 3 Df
# Estimated ratio for means "deg" degrees apart (Temp-mean is 2nd parameter)
exp(deg*coef(zip.mod.t0)[2])
## count_Temp 
##   1.158352
# Confidence interval for all parameters
confint(zip.mod.t0) 
##                         2.5 %     97.5 %
## count_(Intercept) -3.25701992 0.35908153
## count_Temp         0.06740321 0.22659318
## zero_(Intercept)  -0.95529354 0.01123516
# Confidence interval for ratio of means "deg" degrees apart
round(deg*exp(confint(zip.mod.t0)[2,]), digits = 2)
##  2.5 % 97.5 % 
##   1.07   1.25
#############################################################
# Demonstrate zero handling of all 4 models and regular Poisson model
# Fit regular Poisson model
poi.mod <- glm(NumEggs ~ Temp, family = poisson(link = "log"), data = eggdata2)
# Compute predicted number of zeroes by summing the probability of 0 across all observations
mu.poi <- exp(predict(poi.mod)) 
zero.poi <- sum(exp(-mu.poi))
zero.zip.tt <- sum(predict(object = zip.mod.tt, type = "prob")[,1])
zero.zip.t0 <- sum(predict(object = zip.mod.t0, type = "prob")[,1])
zero.zip.0t <- sum(predict(object = zip.mod.0t, type = "prob")[,1])
zero.zip.00 <- sum(predict(object = zip.mod.00, type = "prob")[,1])
data.frame(observed = sum(eggdata2$NumEggs == 0), Poi = round(zero.poi, digits = 2), 
           ZIP.tt = round(zero.zip.tt, digits = 2), ZIP.t0 = round(zero.zip.t0, digits = 2), 
           ZIP.0t = round(zero.zip.0t, digits = 2), ZIP.00 = round(zero.zip.00, digits = 2))
##   observed  Poi ZIP.tt ZIP.t0 ZIP.0t ZIP.00
## 1       27 1.47     27  27.02     27     27



Exemplo 4.14: Simulação comparando os modelos ZIP e Poisson.


Uma dificuldade potencial com o modelo ZIP é o fato de que a contagem zero pode surgir de duas fontes distintas:

  1. membros imunes da população e

  2. membros suscetíveis que por acaso registram um zero.

Observações individuais de zero geralmente não fornecem nenhuma indicação direta sobre se eles são imunes ou suscetíveis. Portanto, é difícil separar as duas fontes com precisão usando apenas dados de contagens observadas, particularmente quando a média de susceptíveis é próxima de zero, de modo que uma fração apreciável da população suscetível produz contagens de zero.

O objetivo desta simulação de Monte Carlo é explorar as propriedades das estimativas de \(\pi\) e \(\mu\) sob várias condições quando seus valores verdadeiros são conhecidos. Geramos um grande número de conjuntos de dados a partir de um modelo ZIP com valores fixos para \(\pi\) e \(\mu\).

Em seguida, usamos o modelo ZIP para estimar esses parâmetros e encontrar os intervalos de confiança de Wald para eles. Calculamos o valor médio das estimativas de parâmetro para que possamos verificar se há viés nas estimativas. Estimamos o verdadeiro nível de confiança dos intervalos de Wald e calculamos os pontos finais medianos para ver como seria um intervalo “típico”.

# Data generation logic: 
#   Generate n Bernoullis, 1=susceptible, 0=immune [P(Susc.) = 1-P(Immune)}
#   Generate n Poisson counts 
#   Multiply count by Bernoulli! If Immune, count=0.  Otherwise, count=Poisson
# Analysis function
#   Function to fit ZIP models with no explanatories
#   Get out the estimated mean and P(0) with corresponding CIs
# Note on limitations:
#   Pilot runs highlighted a problem with the use of large values of pi:
#   some data sets were generated that contained 35 counts of 0. The symptom that we 
#    observed were NA values for some of our summary statistics, which led us to 
#    discover the source of the problem. When there are no non-zero counts, there 
#    is no way to distinguish between very small mu and a very large pi, so some 
#    estimates of infinity are produced. We can calculate that 
#    P(All Y's are 0)=P(Y=0)^{35}=[pi+(1-pi)exp(-mu)]^{35}
#    =.009 when pi=0.8 and 0.102 when pi=0.9.  
#   Thus, in generating 1000 data sets from each of these cases, we expect roughly 
#    9 and 102 sets of all-zero data to be generated, respectively. Indeed, we 
#    observed 9 and 101 such sets, respectively, in our pilot study, and none for 
#    other values of pi, for which the probability of all-zero data is <1/1000. 
#   We therefore added code to omit these cases from the calculation of the summary 
#    statistics.
library(pscl)
zipmod <- function(y) {
  mod.fit <- zeroinfl(y ~ 1 | 1, dist = "poisson")
  log.mu.stats <- c(coef(mod.fit)[1], confint(mod.fit)[c(1,3)])
  logit.pi.stats <- c(coef(mod.fit)[2], confint(mod.fit)[c(2,4)])
  mu.stats <- exp(log.mu.stats)
  pi.stats <- exp(logit.pi.stats)/(1 + exp(logit.pi.stats))
  c(mu.stats, pi.stats)
}
# Check the function
suscept <- rbinom(n = 10000, size = 1, prob = 0.2)
count <- rpois(n = 10000, lambda = 1)
y = suscept*count
zipmod(y)
## count_(Intercept)                                      zero_(Intercept) 
##         0.9528860         0.8886449         1.0217713         0.7940992 
##                                     
##         0.7801133         0.8074152
# Function to create and analyze data
sim.zip <- function(n, pi, mu, sets, seed) {
  # Set seed number to reproduce results and simulate responses
  set.seed(seed)
  suscept <- rbinom(n = n*sets, size = 1, prob = 1 - pi)
  count <- rpois(n = n*sets, lambda = mu)
  y = suscept*count
  # Restructure y for apply() function
  y.mat <- matrix(data = y, nrow = n, ncol = sets)
  # Get rid of all=zero data sets, which create "INF" estimates
  y.mat <- y.mat[,which(apply(X = y.mat, MARGIN = 2, FUN = max) > 0)]
  # y.mat[1:3,1:5] # Check
  # Analyze data
  save.results <- apply(X = y.mat, MARGIN = 2, FUN = zipmod)
  cat("Legitimate data sets (not all 0) = ", ncol(save.results), "\n")
  # Calculate the true confidence level and investigate problems
  cat("Average estimated mean = ", mean(save.results[1,]), "\n")
  mu.cover <- ifelse(test = mu > save.results[2,], 
                     yes = ifelse(test = mu < save.results[3,], yes = 1, no = 0), no = 0)
  cat("Est. true conf. level for mean:", sum(mu.cover)/sets, "\n")
  cat("Median confidence interval for mean = ", 
      median(save.results[2,]), ",", median(save.results[3,]), "\n\n")
  cat("Average estimated P(Immune) = ", mean(save.results[4,]), "\n")
  pi.cover <- ifelse(test = pi > apply(X = save.results[c(5,6),], MARGIN = 2, FUN = min), 
      yes = ifelse(test = pi < apply(X = save.results[c(5,6),], MARGIN = 2, FUN = max), 
                   yes = 1, no = 0), no = 0)
  cat("Est. true conf. level for P(Immune):", sum(pi.cover)/sets, "\n")  
  cat("Median confidence interval for P(Immune) = ", median(save.results[5,]), ",", 
      median(save.results[6,]), "\n")
  save.results
}


A geração de dados é realizada da seguinte forma. Dado o tamanho da amostra \(n\) para um único conjunto de dados e o número de conjuntos de dados que desejamos criar, digamos \(R\), geramos \(nR\) observações independentes de duas variáveis aleatórias independentes. A primeira é uma variável aleatória de Bernoulli, digamos \(S\), que cria os membros imunes e suscetíveis da amostra. Definimos \(S = 1\) para suscetível por motivos que serão aparentes momentaneamente.

Então \(P(S = 1) = 1-\pi\). A segunda variável aleatória é uma variável de contagem, digamos \(C\), de \(Po(\mu )\). Por fim, criamos observações ZIP tomando \(Y = S\times C\). Dessa forma, os imunes são gerados com probabilidade \(\pi\) e suas contagens são definidas como 0 porque \(S = 0\), enquanto os suscetíveis são gerados com probabilidade \(1-\pi\) e suas contagens são distribuídas como \(C \sim Po(\mu )\). Formamos uma matriz de dados dessas contagens e aplicamos a função de análise ZIP a elas, um conjunto de \(n\) por vez. Em seguida, calculamos os intervalos de confiança para ambos \(\mu\) e \(\pi\).

# Keeping mu=1 and varying pi
save.0 <- sim.zip(n = 35, pi = 0, mu = 1, sets = 1000, seed = 23661882)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.072022 
## Est. true conf. level for mean: 0.942 
## Median confidence interval for mean =  0.7166866 , 1.590983 
## 
## Average estimated P(Immune) =  0.06426987 
## Est. true conf. level for P(Immune): 0 
## Median confidence interval for P(Immune) =  1.769848e-98 , 1
save.1 <- sim.zip(n = 35, pi = 0.1, mu = 1, sets = 1000, seed = 985125)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.050734 
## Est. true conf. level for mean: 0.924 
## Median confidence interval for mean =  0.6465357 , 1.701634 
## 
## Average estimated P(Immune) =  0.1256381 
## Est. true conf. level for P(Immune): 0.885 
## Median confidence interval for P(Immune) =  0.001851984 , 0.9053821
save.2 <- sim.zip(n = 35, pi = 0.2, mu = 1, sets = 1000, seed = 329482929)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.006759 
## Est. true conf. level for mean: 0.917 
## Median confidence interval for mean =  0.5785375 , 1.755724 
## 
## Average estimated P(Immune) =  0.187234 
## Est. true conf. level for P(Immune): 0.93 
## Median confidence interval for P(Immune) =  0.01874614 , 0.7754983
save.5 <- sim.zip(n = 35, pi = 0.5, mu = 1, sets = 1000, seed = 601299438)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  0.9931811 
## Est. true conf. level for mean: 0.904 
## Median confidence interval for mean =  0.4689279 , 2.008926 
## 
## Average estimated P(Immune) =  0.4418719 
## Est. true conf. level for P(Immune): 0.965 
## Median confidence interval for P(Immune) =  0.1881219 , 0.8136951
save.8 <- sim.zip(n = 35, pi = 0.8, mu = 1, sets = 1000, seed = 9287467)
## Legitimate data sets (not all 0) =  991 
## Average estimated mean =  0.9601104 
## Est. true conf. level for mean: 0.801 
## Median confidence interval for mean =  0.3004691 , 3.084581 
## 
## Average estimated P(Immune) =  0.6448954 
## Est. true conf. level for P(Immune): 0.984 
## Median confidence interval for P(Immune) =  0.4524859 , 0.951394
save.9 <- sim.zip(n = 35, pi = 0.9, mu = 1, sets = 1000, seed = 012839299)
## Legitimate data sets (not all 0) =  899 
## Average estimated mean =  0.9064185 
## Est. true conf. level for mean: 0.57 
## Median confidence interval for mean =  0.1374885 , 3.959458 
## 
## Average estimated P(Immune) =  0.5830128 
## Est. true conf. level for P(Immune): 0.891 
## Median confidence interval for P(Immune) =  0.469282 , 0.9850391
# Plot estimates of mean and P(Immune) in Color
par(mfrow = c(3,2))
line.width <- 4
# mu 1
hist(x = save.1[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.1, ", mu, "=1")))
abline(v = 1, col = "red", lwd = line.width)
abline(v = mean(save.1[1,]), col = "blue", lwd = line.width)
# pi 0.1
hist(x = save.1[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)), 
     main = expression(paste(pi, "=0.1, ", mu, "=1")))
abline(v = 0.1, col = "red", lwd = line.width)
abline(v = mean(save.1[4,]), col = "blue", lwd = line.width)
legend(x = 0.4, y = 500, legend = c("True Value", "Mean Estimate"), 
       col = c("red", "blue"), lwd = line.width, bty = "n")
# mu 1
hist(x = save.5[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.5, ", mu, "=1")))
abline(v = 1, col = "red", lwd = line.width)
abline(v = mean(save.5[1,]), col = "blue", lwd = line.width)
# pi 0.5
hist(x = save.5[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)),
     main = expression(paste(pi, "=0.5, ", mu, "=1")))
abline(v = 0.5, col = "red", lwd = line.width)
abline(v = mean(save.5[4,]), col = "blue", lwd = line.width)
# mu 1
hist(x = save.9[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.9, ", mu, "=1")))
abline(v = 1, col = "red", lwd = line.width)
abline(v = mean(save.9[1,]), col = "blue", lwd = line.width)
# pi 0.9
hist(x = save.9[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)), 
     main = expression(paste(pi, "=0.9, ", mu, "=1")))
abline(v = 0.9, col = "red", lwd = line.width)
abline(v = mean(save.9[4,]), col = "blue", lwd = line.width)

# Plot estimates of mean and P(Immune) in B/W
par(mfrow = c(3,2))
# mu 1
hist(x = save.1[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.1, ", mu, "=1")))
abline(v = 1, lty = "dotted", lwd = line.width)
abline(v = mean(save.1[1,]), lty = "dashed", lwd = line.width)
# pi 0.1
hist(x = save.1[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)), 
     main = expression(paste(pi, "=0.1, ", mu, "=1")))
abline(v = 0.1, lty = "dotted", lwd = line.width)
abline(v = mean(save.1[4,]), lty = "dashed", lwd = line.width)
legend(x = 0.4, y = 500, legend = c("True Value", "Mean Estimate"), 
       lty = c("dotted", "dashed"), lwd = line.width, bty = "n")
# mu 1
hist(x = save.5[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.5, ", mu, "=1")))
abline(v = 1, lty = "dotted", lwd = line.width)
abline(v = mean(save.5[1,]), lty = "dashed", lwd = line.width)
# pi 0.5
hist(x = save.5[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)), 
     main = expression(paste(pi, "=0.5, ", mu, "=1")))
abline(v = 0.5, lty = "dotted", lwd = line.width)
abline(v = mean(save.5[4,]), lty = "dashed", lwd = line.width)
# mu 1
hist(x = save.9[1,], breaks = c(0:12)/2, xlab = expression(hat(mu)), 
     main = expression(paste(pi, "=0.9, ", mu, "=1")))
abline(v = 1, lty = "dotted", lwd = line.width)
abline(v = mean(save.9[1,]), lty = "dashed", lwd = line.width)
# pi 0.9
hist(x = save.9[4,], breaks = c(0:10)/10, xlab = expression(hat(pi)), 
     main = expression(paste(pi, "=0.9, ", mu, "=1")))
abline(v = 0.9, lty = "dotted", lwd = line.width)
abline(v = mean(save.9[4,]), lty = "dashed", lwd = line.width)

# Increasing sample size for difficult cases.  
save.9.100 <- sim.zip(n = 100, pi = 0.9, mu = 1, sets = 1000, seed = 927849027)
## Legitimate data sets (not all 0) =  998 
## Average estimated mean =  0.9853552 
## Est. true conf. level for mean: 0.912 
## Median confidence interval for mean =  0.3557912 , 2.668612 
## 
## Average estimated P(Immune) =  0.8287459 
## Est. true conf. level for P(Immune): NA 
## Median confidence interval for P(Immune) =  0.7356356 , NA
save.9.500 <- sim.zip(n = 500, pi = 0.9, mu = 1, sets = 1000, seed = 12674726)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  0.9839414 
## Est. true conf. level for mean: 0.962 
## Median confidence interval for mean =  0.6281473 , 1.510906 
## 
## Average estimated P(Immune) =  0.8961218 
## Est. true conf. level for P(Immune): 0.946 
## Median confidence interval for P(Immune) =  0.8462975 , 0.9350361
save.0.100 <- sim.zip(n = 100, pi = 0, mu = 1, sets = 1000, seed = 67728891)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.045762 
## Est. true conf. level for mean: 0.947 
## Median confidence interval for mean =  0.8194361 , 1.326069 
## 
## Average estimated P(Immune) =  0.03873342 
## Est. true conf. level for P(Immune): 0 
## Median confidence interval for P(Immune) =  1.378535e-77 , 1
save.0.500 <- sim.zip(n = 500, pi = 0, mu = 1, sets = 1000, seed = 2356432)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.017809 
## Est. true conf. level for mean: 0.948 
## Median confidence interval for mean =  0.9110115 , 1.128744 
## 
## Average estimated P(Immune) =  0.01915619 
## Est. true conf. level for P(Immune): 0 
## Median confidence interval for P(Immune) =  6.680828e-56 , 1
# Trying different combinations of pi and mu that have the same overall mean count.  
save.0.2 <- sim.zip(n = 35, pi = 0.2, mu = 5/4, sets = 1000, seed = 29950288)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  1.26443 
## Est. true conf. level for mean: 0.936 
## Median confidence interval for mean =  0.791857 , 1.995204 
## 
## Average estimated P(Immune) =  0.1901421 
## Est. true conf. level for P(Immune): 0.933 
## Median confidence interval for P(Immune) =  0.03428871 , 0.6669797
save.0.5 <- sim.zip(n = 35, pi = 0.5, mu = 2, sets = 1000, seed = 884377129)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  2.000111 
## Est. true conf. level for mean: 0.954 
## Median confidence interval for mean =  1.355286 , 2.99212 
## 
## Average estimated P(Immune) =  0.4918913 
## Est. true conf. level for P(Immune): 0.971 
## Median confidence interval for P(Immune) =  0.2994437 , 0.6939398
save.0.8 <- sim.zip(n = 35, pi = 0.8, mu = 5, sets = 1000, seed = 448423156)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  4.937007 
## Est. true conf. level for mean: 0.951 
## Median confidence interval for mean =  3.500033 , 6.949617 
## 
## Average estimated P(Immune) =  0.7989219 
## Est. true conf. level for P(Immune): 0.964 
## Median confidence interval for P(Immune) =  0.6336043 , 0.9009084
# Values from Beetle Egg 21C
save.0.38 <- sim.zip(n = 35, pi = 0.38, mu = 5.14, sets = 1000, seed = 1200201)
## Legitimate data sets (not all 0) =  1000 
## Average estimated mean =  5.138207 
## Est. true conf. level for mean: 0.947 
## Median confidence interval for mean =  4.229646 , 6.173225 
## 
## Average estimated P(Immune) =  0.3774121 
## Est. true conf. level for P(Immune): 0.964 
## Median confidence interval for P(Immune) =  0.2263839 , 0.5381585


Escolhemos \(n = 35\) porque esse era o tamanho da amostra no exemplo de postura de besouro acima. Definimos \(\mu= 1\) como um caso em que a distribuição de contagem produzirá naturalmente um número razoavelmente grande de zeros \[ P(C = 0) = \exp(-1)\approx 0.37 \] e exploramos valores de \(\pi= 0\), 0.1, 0.2, 0.5, 0.8 e 0.9. Finalmente, usamos \(R = 1000\) para que possamos obter uma estimativa razoavelmente boa das distribuições de \(\widehat{\pi}\) e \(\widehat{\mu}\) e para o nível verdadeiro dos intervalos de confiança correspondentes.

As execuções piloto destacaram um problema com o uso de grandes valores de \(\pi\): alguns conjuntos de dados foram gerados contendo 35 contagens de 0. Isso criou problemas descritos no programa para este exemplo. Portanto, adicionamos código para omitir esses casos do cálculo das estatísticas resumidas. Os resultados das simulações são apresentados na tabela abaixo e na figura acima.

\[ \begin{array}{c|cccc|cccc}\hline & & & \widehat{\pi} & & & & \widehat{\mu} & & \\\hline & & \mbox{Nível} & \mbox{I.C.} & \mbox{mediano} & & \mbox{Nível} & \mbox{I.C.} & \mbox{mediano} \\ \pi & \mbox{Média} & \mbox{Conf.} & \mbox{Mínimo} & \mbox{Máximo} & \mbox{Média} & \mbox{Conf.} & \mbox{Mínimo} & \mbox{Máximo} \\\hline 0 & 0.06 & 0 & 0 & 1 & 1.07 & 94.2 & 0.72 & 1.59 \\ 0.1 & 0.13 & 88.5 & 0.00 & 0.91 & 1.05 & 92.4 & 0.65 & 1.70 \\ 0.2 & 0.19 & 93.0 & 0.02 & 0.77 & 1.01 & 91.7 & 0.58 & 1.76 \\ 0.5 & 0.44 & 96.5 & 0.19 & 0.81 & 0.99 & 90.4 & 0.47 & 2.01 \\ 0.8 & 0.64 & 98.4 & 0.45 & 0.95 & 0.96 & 80.1 & 0.30 & 3.08 \\ 0.9 & 0.58 & 89.1 & 0.47 & 0.99 & 0.91 & 57.0 & 0.14 & 3.96 \\\hline \end{array} \]


É evidente a partir dos resultados que às vezes pode ser difícil estimar os parâmetros de um modelo ZIP. Quando a probabilidade de imune é muito pequena, há grande incerteza nas estimativas de \(\pi\), conforme indicado pelos intervalos de confiança excessivamente amplos. Os níveis de confiança verdadeiros estimados dos intervalos de confiança correspondentes são muito baixos, apesar de suas vastas larguras. Com \(\mu = 1\) os susceptíveis geram em média 37% de zeros e os imunes acrescentam muito pouco ao total. O modelo estima uma probabilidade próxima de zero de imune, mas com um imenso erro padrão em \(logit(\pi)\). Assim, os pontos finais do intervalo são transformados em quase (mas não totalmente) 0 e 1.

Por outro lado, com valores muito grandes de \(\pi\) há poucas contagens diferentes de zero para fornecer estimativas confiáveis de \(\mu\). As estimativas \(\widehat{\mu}\) e \(\widehat{\pi}\) são muito pequenas, indicando que o procedimento tende a atribuir muitos zeros a um \(\mu\) pequeno em vez de um \(\pi\) grande. A maior parte dessa dificuldade está associada ao pequeno tamanho da amostra. Não há informações suficientes para estimar com segurança uma pequena média e uma grande probabilidade de imunidade com \(n = 35\). Simulações adicionais usando \(n = 100\) e \(n = 500\) com \(\pi= 0.9\) fornecem estimativas médias de 0.99 e 0.98 para \(\mu\) e 0.83 e 0.90 para \(\pi\), espectivamente.

A partir deste exemplo, podemos concluir que é necessário algum cuidado ao usar os modelos de Poisson inflacionados por zero com tamanhos de amostra pequenos, especialmente quando há muito poucas observações diferentes de zero. No Exercício 21 as simulações podem ser repetidas para o modelo hurdle.



4.5 Exercícios