Modelo multinomial log-linear

---

Introdução

O valor estético e ornamental de pimenteiras do gênero Capsicum é em grande parte devido a coloração dos frutos, considerada a alta variabilidade dessa característica (Figura).

Considere um estudo experimental em que foi determinada a cor predominante dos frutos em estágio intermediário de maturação de quatro famílias (G1-G4) de pimenteiras. Para cada família, foram tomadas três repetições, cada uma constituída de dez plantas, sendo o experimento conduzido em casa-de-vegetação. Foram encontrados frutos de cor Amarela, Vermelha e Roxa, conforme a matriz de contagens a seguir.

Em R:

Y <- matrix(c(2, 7, 1, 1, 8, 1, 2, 5, 3,
              8, 2, 0, 8, 1, 1, 6, 3, 1,
              6, 4, 0, 8, 2, 0, 5, 4, 1,
              1, 0, 9, 0, 2, 8, 0, 3, 7), ncol = 3, byrow = TRUE)
colnames(Y) <- c("A", "V", "R")

familias <- rep(paste0("G", 1:4), each = 3)
dados <- data.frame(familias, Y)
dados

   familias A V R
1        G1 2 7 1
2        G1 1 8 1
3        G1 2 5 3
4        G2 8 2 0
5        G2 8 1 1
6        G2 6 3 1
7        G3 6 4 0
8        G3 8 2 0
9        G3 5 4 1
10       G4 1 0 9
11       G4 0 2 8
12       G4 0 3 7

Cada linha da tabela acima representa uma observação trinomial. Por exemplo, das 10 plantas da primeira repetição de G1, 2 apresentaram frutos amarelos, 7 vermelhos e 1 roxo, sendo estas categorias mutuamente excludentes, isto é, um fruto não pode pertencer a mais de uma classe de cor. Se considerarmos que as doze observações são independentes, então podemos utilizar o modelo multinomial para descrever as probabilidades de ocorrência das cores. Mais especificamente, podemos formular um modelo multinomial do tipo log-linear ou logit para estimar o efeito das famílias sobre a frequência de cores.

Mas antes, vale um exame visual das proporções de cores por família.

prop <- xtabs(Y ~ familias)
plot(prop, col = c("yellow", "red", "purple"))

Há alguma similaridade entre as famílias G2 e G3.

Formulação do modelo

Pelo modelo linear, uma observação trinomial (vetor de comprimento \(K = 3\)) seria descrita por

\[\mathsf{y}_{ij} = \mathsf{\mu}_i + \mathsf{\epsilon}_{ij}\] em que \(i = 1, 2, ..., 4\) (famílias), \(j = 1,2,3\) (repetições). O problema então está em estimar os \(\mathsf{\mu}_i\). Para isso, assumiremos que os efeitos das famílias é linear, mas não necessáriamente estão relacionados de forma direta com \(\mathsf{\mu}_i\), e sim através da função logit, que é apropriada para modelar proporções (probabilidades de ocorrência de cada cor).

\[logit(\mathsf{\mu}_i) = \log(\mathsf{\mu}_i \cdot ||\mathsf{\mu}_i||_1^{-1}) = \mathsf{x}_{ij}' \mathsf{\beta}\]

em que \(\mathsf{x}_{ij}\) representa um vetor de incidência dos efeitos das famílias sobre a observação \(\mathsf{y}_{ij}\). Isso pode ser entendido como cada linha da matriz \(\mathsf{X}\) do modelo:

X <- model.matrix( ~ familias, dados)
X

   (Intercept) familiasG2 familiasG3 familiasG4
1            1          0          0          0
2            1          0          0          0
3            1          0          0          0
4            1          1          0          0
5            1          1          0          0
6            1          1          0          0
7            1          0          1          0
8            1          0          1          0
9            1          0          1          0
10           1          0          0          1
11           1          0          0          1
12           1          0          0          1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$familias
[1] "contr.treatment"

Já \(\mathsf{\beta}\) é uma matriz de dimensão \(p \times (K-1)\) de parâmetros a serem estimados, sendo \(p = 4\) (famílias). Note que cada observação trinomial possui rank \(K-1=2\), ou dois graus de liberdade, já que são dados de proporção. Nota: nas implementações do R, os coeficientes (parâmetros) referentes à primeira classe são automaticamente zerados. Esse procedimento é similar ao que ocorre com a matriz \(\mathsf{X}\).

Dessa forma, a primeira categoria (\(k=1\), cor Amarela) é definida como categoria base, de modo que a função logit ganha sentido de log-risco relativo de ocorrência das demais categorias.

\[logit(\mathsf{\mu}_i) = \begin{bmatrix} \log \left( \frac{P(Vermelha)}{P(Amarela)} \right) \\ \log \left( \frac{P(Roxa)}{P(Amarela)} \right) \end{bmatrix}' = \begin{bmatrix} \beta_{0V} + \beta_{iV} \\ \beta_{0R} + \beta_{iR} \end{bmatrix}'\]

em que \(\beta_{0V}\) e \(\beta_{0R}\) são as constantes ou interceptos associados às ocorrências médias de cores Vermelha e Roxa, respectivamente; enquanto \(\beta_{iV}\) e \(\beta_{iR}\) representam os efeitos de uma dada família \(i\) sobre a ocorrência das cores Vermelha e Roxa.

E tem-se:

\[\mathsf{\mu}_i = \exp(\mathsf{x}_{ij}' \mathsf{\beta}) \cdot || \exp(\mathsf{x}_{ij}' \mathsf{\beta}) ||_1^{-1}\] aqui representando a probabilidade associada a cada classe.

Estimação por verossimilhança

Pelo modelo multinomial, a probabilidade de obter uma determinada observação \(\mathsf{y}_{ij}\), dados os parâmetros \(\mathsf{\mu}_i\) e \(n_{ij}\) (total de cada observação \(ij\)), é

\[P(\mathsf{y}_{ij} | \mathsf{\mu}_i, n_{ij}) = \frac{n_{ij}!}{\prod_{k=1}^K y_{ijk}!}{\prod_{k=1}^K \mu_{ik}^{y_{ijk}}} \]

Nota: \(\mu_{ik}\) é o mesmo para todas as repetições \(j\).

Agora, tomando \(P(\mathsf{y}_{ij} | \mathsf{\mu}_i, n_{ij})\) como \(L(\mathsf{\mu}_i | \mathsf{y}_{ij})\), a verossimilhança associada à observação \(\mathsf{y}_{ij}\), tem-se a seguinte log-verossimilhança:

\[\log L(\mathsf{\mu}_i | \mathsf{y}_{ij}) = \log(n_{ij}!) - \sum_{k=1}^K \log(y_{ijk}!) + \sum_{k=1}^K y_{ijk}\log(\mu_{ik}) \] mas como apenas o último termo depende de parâmetro a ser estimado, então a função de log-verossimilhança do modelo fica:

\[\log L(\mathsf{\mu} | \mathsf{Y}) \propto \sum_{i=1}^I \sum_{j=1}^J \sum_{k=1}^K y_{ijk}\log(\mu_{ik})\] Uma vez que \(\mathsf{\mu}\) é função de \(\mathsf{\beta}\), vejamos como construir a log-verossimilhança no R:

nlogL <- function(betas) {
   b <- matrix(betas, nrow = ncol(X))
   p_num <- exp(X %*% cbind(0, b))
   mat_p <- sweep(p_num, 1L, apply(p_num, 1L, sum), FUN = "/")
   ll <- Y * log(mat_p)
   -sum(ll)
}

Agora podemos utilizar um otimizador geral para encontrar numericamente as estimativas da matriz \(\mathsf{\beta}\). Mas para isso, é necessário inicializar tal matriz.

init <- matrix(0, nrow = ncol(X), ncol = ncol(Y) - 1,
    dimnames = list(colnames(X), colnames(Y)[-1]))
init

            V R
(Intercept) 0 0
familiasG2  0 0
familiasG3  0 0
familiasG4  0 0

nlogL(init)

[1] 131.8335

Perceba que o valor obtido refere-se ao negativo da função de log-verossimilhança, de forma que essa possa ser minimizada, que é o padrão da função optim(). Assim, as estimativas de máxima verossimilhança (resultado em $par) de \(\mathsf{\beta}\) são obtidas por:

mle <- optim(init, nlogL, method = "BFGS")
mle

$par
                     V             R
(Intercept)  1.3860658 -0.0002492433
familiasG2  -2.6853721 -2.3975151827
familiasG3  -2.0279414 -2.9438624442
familiasG4   0.2232172  3.1781057873

$value
[1] 88.7042

$counts
function gradient 
      25       14 

$convergence
[1] 0

$message
NULL

Risco relativo

Os riscos relativos são obtidos exponenciando os coeficientes \(\hat{\mathsf{\beta}}\), e são interpretados tomando como referência a categoria base, isto é, a cor Amarela. Então, faz sentido definir como ponto de corte, para interpretação, o valor 1.

Por exemplo, o risco da família G1 produzir frutos Vermelhos em vez de Amarelos é de \(\exp(1.386) \approx 4\) vezes. Já o risco da família G2 produzir frutos Vermelhos em vez de Amarelos é baixo, de \(\exp(1.386 - 2.685) \approx 0.27\) vezes. O risco da família G3 produzir frutos Roxos em vez de Amarelos é de \(\exp(-0.00025 + 3.178) \approx 24\) vezes, muito alto.

Probabilidades esperadas

Pelo modelo ajustado, as probabilidades esperadas em cada observação são:

p_num <- exp(X %*% cbind(A = 0, mle$par))
p_esp <- sweep(p_num, 1L, apply(p_num, 1L, sum), FUN = "/")
round(p_esp, 4)

        A      V      R
1  0.1667 0.6666 0.1667
2  0.1667 0.6666 0.1667
3  0.1667 0.6666 0.1667
4  0.7333 0.2000 0.0667
5  0.7333 0.2000 0.0667
6  0.7333 0.2000 0.0667
7  0.6333 0.3333 0.0333
8  0.6333 0.3333 0.0333
9  0.6333 0.3333 0.0333
10 0.0333 0.1667 0.8000
11 0.0333 0.1667 0.8000
12 0.0333 0.1667 0.8000

Por exemplo, a probabilidade esperada de frutos vermelhos na família G2 (linhas 4 a 6) vem de

\[logit(\hat{\mu}_{2V}) = \hat{\beta}_{0V} + \hat{\beta}_{2V}\] e então:

\[\begin{align} \hat{\mu}_{2V} &= \frac{\exp(\hat{\beta}_{0V} + \hat{\beta}_{2V})}{\exp(\hat{\beta}_{0A} + \hat{\beta}_{2A}) + \exp(\hat{\beta}_{0V} + \hat{\beta}_{2V}) + \exp(\hat{\beta}_{0R} + \hat{\beta}_{2R})} \\ &= \frac{\exp(1.386 - 2.685)}{ \exp(0 + 0) + \exp(1.386 - 2.685) + \exp(-0.00024 - 2.397)} \\ &\approx 0.20 \end{align}\]

Teste da razão de verossimilhanças

O principal objetivo do estudo exemplificado é avaliar a hipótese de nulidade de que as famílias apresentam as mesmas proporções de classes de cor dos frutos, ou seja,

\[H_0: \mathsf{\mu}_1 = \mathsf{\mu}_2 = \mathsf{\mu}_3 = \mathsf{\mu}_4\]

E isso pode ser respondido através de um TRV, cuja estatística é dada pela diferença da log-verossimilhança máxima (com \(\mathsf{\mu}(\hat{\beta})\)) e sob \(H_0\).

\[L_R = 2[\log L(\hat{\mathsf{\beta}}) - \log L(\mathsf{\beta}_{H_0})]\] Precisamos então definir \(\mathsf{\beta}_{H_0}\). Faremos isso zerando os coeficientes relativos às famílias, deixando apenas as constantes (interceptos).

betas_H0 <- mle$par
betas_H0[2:4, ] <- 0
betas_H0

                   V             R
(Intercept) 1.386066 -0.0002492433
familiasG2  0.000000  0.0000000000
familiasG3  0.000000  0.0000000000
familiasG4  0.000000  0.0000000000

Assim,

LR <- 2 * (-nlogL(mle$par) + nlogL(betas_H0))
LR

[1] 138.9259

Que é esperado ter distribuição qui-quadrado com \(\nu = 8 - 2 = 6\) graus de liberdade. Note: \(\nu\) é a diferença de número de parâmetros em \(\mathsf{\beta}\) e \(\mathsf{\beta}_{H_0}\). Com isso, calcula-se o valor-\(p\):

1 - pchisq(LR, df = 6)

[1] 0

Concluíndo com isso (\(p < 0.01\)) que a proporção não é a mesma nas quatro famílias.

Sem rodeios

Sem rodeios, usando a implementação do pacote nnet (Venables & Ripley, 2002), ajusta-se o mesmo modelo com:

library(nnet)
modelo <- multinom(Y ~ familias, dados)

# weights:  15 (8 variable)
initial  value 131.833475 
iter  10 value 88.706145
final  value 88.704203 
converged

modelo

Call:
multinom(formula = Y ~ familias, data = dados)

Coefficients:
    (Intercept) familiasG2 familiasG3 familiasG4
V  1.3860658306  -2.685372  -2.027941  0.2232172
R -0.0002491309  -2.397515  -2.943862  3.1781056

Residual Deviance: 177.4084 
AIC: 193.4084 

Já o TRV pode ser obtido por:

library(car)
Anova(modelo)

# weights:  6 (2 variable)
initial  value 131.833475 
final  value 130.382068 
converged

Analysis of Deviance Table (Type II tests)

Response: Y
         LR Chisq Df Pr(>Chisq)    
familias   83.356  6   7.23e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Alternativamente, pode-se compreender o que é realizado pela Anova() realizando ‘manualmente’ o teste de modelos encaixados, em que o modelo ajustado é comparado com o modelo nulo, sem os efeitos das famílias, contendo apenas os interceptos \(\beta_{0V}\) e \(\beta_{0R}\).

modelo_nulo <- update(modelo, formula = ~ 1)

# weights:  6 (2 variable)
initial  value 131.833475 
final  value 130.382068 
converged

anova(modelo_nulo, modelo)

     Model Resid. df Resid. Dev   Test    Df LR stat.      Pr(Chi)
1        1        22   260.7641           NA       NA           NA
2 familias        16   177.4084 1 vs 2     6 83.35573 7.771561e-16

Há aqui uma ressalva em relação aos dois TRVs, pois o feito na seção anterior utiliza dos mesmos coeficientes para os interceptos (constantes) sob \(H_0\); já a abordagem feita por último, reajusta um novo modelo apenas com os interceptos, o que promove estimativas diferentes e, consequentemente, log-verossimilhanças diferentes. Isso é esperado, já que as estimativas dos parâmetros possuem correlação. Mas ambas as abordagens são válidas, mudando as interpretações com as estimativas utilizadas.

Por último, fitted(modelo) devolve uma matriz com as probabilidades esperadas para cada observação.

Intervalos de confiança

Uma forma simples de avaliar diferenças específicas entre os grupos experimentais é por meio de intervalos de 95% de confiança para as estimativas \(\hat{\mathsf{\beta}}\) de cada categoria. Observe, por exemplo, como estão sobrepostos os ICs das famílias G2 e G3, sugerindo não haver diferenças (\(p < 0.05\)) entre elas.

confint(modelo)

, , V

                 2.5 %     97.5 %
(Intercept)  0.4061564  2.3659752
familiasG2  -4.0176991 -1.3530449
familiasG3  -3.2715493 -0.7843333
familiasG4  -2.1367076  2.5831421

, , R

                 2.5 %     97.5 %
(Intercept) -1.2397962  1.2392980
familiasG2  -4.3031844 -0.4918459
familiasG3  -5.3058266 -0.5818981
familiasG4   0.8249602  5.5312511

Grau de ajuste

O coeficiente de determinação (\(R^2\)) é uma medida amplamente utilizada para medir o grau de ajuste de modelos lineares clássicos. Mas não é apropriado para o modelo multinomial. Não obstante, uma alternativa é calcular uma medida análoga proposta por McFadden (1974):

\[\rho^2 = 1 - \frac{\log L(\hat{\mathsf{\beta}})}{\log L(\mathsf{\beta}_{H_0})}\]

r2 = c(1 - logLik(modelo)/logLik(modelo_nulo))
r2

[1] 0.3196595

Que também é interpretado como ‘proporção [0, 1] da variação explicada’.

Um modelo para classificação

No exemplo trabalhado, não é tarefa difícil classificar as famílias quanto à cor dos frutos. Isso porque a variável preditora em \(\mathsf{X}\) é do tipo categórica (famílias). Observe como podemos, a partir do modelo ajustado, predizer as cores de cada observação.

predict(modelo, type = "class")

 [1] V V V A A A A A A R R R
Levels: A V R

Esse resultado condiz com a matriz de probabilidades esperadas, calculada anteriormente.

No entanto, suponha que a variável preditora em \(\mathsf{X}\) seja quantitativa, digamos doses de um bioestimulante de antocianina, interferindo assim na coloração dos frutos. Ora, já não parace mais tão trivial classificar/predizer a cor do fruto empiricamente. Mas através de um modelo de regressão multinomial, ajustado de forma muito similar ao apresentado aqui, a tarefa poderia ser executada. Mais do que isso, o modelo poderia ainda ser utilizado para realizar predições com novos valores de dose do bioestimulante que não foram testadas (desde que estejam dentro do intervalo testado).

Miscelânea

A forma de declaração do modelo multinomial da seção anterior não é a única. Outra opção é declarar a variável resposta na foma de factor() (no exemplo, com \(K = 3\) níveis/categorias de cor), em vez da matriz de contagens em que cada coluna representa uma categoria. Essa abordagem permite realizar comparações múltiplas com implementações como a do pacote emmeans.

Podemos fazer este exercício. Mas antes será preciso reorganizar a tabulação dos dados.

library(reshape2)
md <- melt(dados)
head(md)

  familias variable value
1       G1        A     2
2       G1        A     1
3       G1        A     2
4       G2        A     8
5       G2        A     8
6       G2        A     6

Em seguida, mais um passo é necessário - utilizar a terceira coluna para repetir os valores das categorias das outras duas colunas.

g <- as.factor(unlist(apply(md, 1L, function(x) rep(x[1], x[3]))))
y <- as.factor(unlist(apply(md, 1L, function(x) rep(x[2], x[3]))))
length(y)

[1] 120

E agora, o ajuste do modelo. Observe como é idêntico ao ajustado da forma anterior.

modelo2 <- multinom(y ~ g)

# weights:  15 (8 variable)
initial  value 131.833475 
iter  10 value 88.706145
final  value 88.704203 
converged

modelo2

Call:
multinom(formula = y ~ g)

Coefficients:
    (Intercept)       gG2       gG3       gG4
R -0.0002491309 -2.397515 -2.943862 3.1781056
V  1.3860658306 -2.685372 -2.027941 0.2232172

Residual Deviance: 177.4084 
AIC: 193.4084 

Veja as probabilidades esperadas.

library(emmeans)
medias <- emmeans(modelo2, ~ g|y)
medias

y = A:
 g    prob     SE df lower.CL upper.CL
 G1 0.1667 0.0680  8  0.00978    0.324
 G2 0.7333 0.0807  8  0.54715    0.920
 G3 0.6333 0.0880  8  0.43045    0.836
 G4 0.0333 0.0328  8 -0.04224    0.109

y = R:
 g    prob     SE df lower.CL upper.CL
 G1 0.1667 0.0680  8  0.00976    0.324
 G2 0.0667 0.0455  8 -0.03835    0.172
 G3 0.0333 0.0328  8 -0.04224    0.109
 G4 0.8000 0.0730  8  0.63158    0.968

y = V:
 g    prob     SE df lower.CL upper.CL
 G1 0.6666 0.0861  8  0.46817    0.865
 G2 0.2000 0.0730  8  0.03159    0.368
 G3 0.3333 0.0861  8  0.13486    0.532
 G4 0.1667 0.0680  8  0.00977    0.324

Confidence level used: 0.95 

E os respectivos intervalos de 95% de confiança.

plot(medias)

Por fim, comparações aos pares podem ser também executadas.

pairs(medias)

y = A:
 contrast estimate     SE df t.ratio p.value
 G1 - G2   -0.5666 0.1056  8  -5.366  0.0030
 G1 - G3   -0.4666 0.1112  8  -4.195  0.0128
 G1 - G4    0.1334 0.0755  8   1.766  0.3542
 G2 - G3    0.1000 0.1194  8   0.837  0.8355
 G2 - G4    0.7000 0.0871  8   8.033  0.0002
 G3 - G4    0.6000 0.0939  8   6.390  0.0010

y = R:
 contrast estimate     SE df t.ratio p.value
 G1 - G2    0.1000 0.0819  8   1.221  0.6318
 G1 - G3    0.1333 0.0755  8   1.765  0.3544
 G1 - G4   -0.6333 0.0998  8  -6.345  0.0010
 G2 - G3    0.0333 0.0561  8   0.594  0.9311
 G2 - G4   -0.7333 0.0861  8  -8.520  0.0001
 G3 - G4   -0.7666 0.0800  8  -9.577  0.0001

y = V:
 contrast estimate     SE df t.ratio p.value
 G1 - G2    0.4666 0.1129  8   4.134  0.0139
 G1 - G3    0.3333 0.1217  8   2.738  0.0962
 G1 - G4    0.5000 0.1097  8   4.557  0.0080
 G2 - G3   -0.1333 0.1129  8  -1.181  0.6542
 G2 - G4    0.0333 0.0998  8   0.334  0.9862
 G3 - G4    0.1667 0.1097  8   1.519  0.4703

P value adjustment: tukey method for comparing a family of 4 estimates 

Referências

McFadden, D. (1973) Conditional logit analysis of qualitative choice behavior. Frontiers of Econometrics, ed. by P. Zarembka. New York: Academic Press.

Venables, W. N., Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0

---