25 Modelo multinomial

En este capítulo se muestran como ajustar un modelo multinomial (nominal y ordinal) con efectos aleatorios.

Ejemplo: modelo multinomial nominal

Este ejemplo está basado en este tutorial. El objetivo es construir un modelo multinomial nominal para explicar la variable Response que tiene cuatro niveles: BareNoun, NumeralNoun, QuantAdjNoun y QuantNoun. En esta situación los niveles de Response no tienen un orden, por esa razón es que este ejemplo encaja en un modelo multinomial nominal.

Los datos que se van a utilizar en este ejemplo corresponden a un experimento en el cual a un grupo de personas con cierta habilidad de lenguaje (English_Mono, German_Mono, L2_Advanced, L2_Intermediate), se les presentaron diez figuras (Item) y ellos debían identificar si la figura se refería a un sustantivo Bare, Numeral, QuantAdj o Quant. Las diez figuras son una muestra de todas las figuras que se podrían tener para el experimento, por eso se pueden considerar que la figura (Item) puede tener una influencia en la respuesta.

El código para acceder a la base de datos se muestra a continuación.

# Para cargar los datos
pict <- base::readRDS(url("https://slcladal.github.io/data/pict.rda", "rb"))
# Primeras observaciones de la base de datos
head(pict, n=6)

##   Id Participant       Group Item     Response Gender Age
## 1  1        G001 German_Mono    1  NumeralNoun   Male  18
## 2  2        G002 German_Mono    3  NumeralNoun   Male  18
## 3  3        G003 German_Mono    4  NumeralNoun   Male  18
## 4  4        G004 German_Mono    6 QuantAdjNoun   Male  18
## 5  5        G005 German_Mono    8  NumeralNoun   Male  18
## 6  6        G006 German_Mono    9 QuantAdjNoun   Male  18

Solución.

Antes de ajustar el modelo vamos a explorar un poco los datos. Vamos a crear tablas de frecuencia para la variable respuesta, las covariables y la variable de agrupación Item.

# Tabla para la respuesta
table(pict$Response)

## 
##     BareNoun  NumeralNoun QuantAdjNoun    QuantNoun 
##          394          510          152           34

# Tabla para las covariables
table(pict$Gender)

## 
## Female   Male 
##    740    350

table(pict$Group)

## 
##    English_Mono     German_Mono     L2_Advanced L2_Intermediate 
##             230             310             340             210

# Tabla para la agrupación
table(pict$Item)

## 
##   1   3   4   6   8   9  11  12  13  15 
## 109 109 109 109 109 109 109 109 109 109

Vamos a retirar la palabra Noun de las respuestas, de esta manera nos quedará más compacta una tabla que haremos al final.

library(stringr)
pict$Response <- str_replace(string=pict$Response, 
                             pattern="Noun.*", 
                             replacement="")
pict$Response <- factor(pict$Response)

Vamos ahora a crear el modelo multinomial nominal con la ayuda de la función mblogit del paquete mclogit de Elff (2022). Las covariables que vamos a utilizar para predecir a Response son Gender, Group y Age. La variable Item se refiere a la figura o fotografía mostrada a los participantes y esta variable se va a utilizar como fuente del intercepto aleatorio \(b0\) que se va a incluir en el modelo. El intercepto aleatorio se le indica a mblogit por medio de random = ~ 1 | Item.

library(mclogit)
mod1 <- mblogit(formula = Response ~ Gender + Group + Age, 
                random = ~ 1 | Item, 
                data = pict)

## 
## Iteration 1 - deviance = 1648.607 - criterion = 0.8003659
## Iteration 2 - deviance = 1564.408 - criterion = 0.08391308
## Iteration 3 - deviance = 1544.629 - criterion = 0.02825377
## Iteration 4 - deviance = 1539.952 - criterion = 0.005299121
## Iteration 5 - deviance = 1538.684 - criterion = 0.000183948
## Iteration 6 - deviance = 1538.481 - criterion = 1.807221e-07
## Iteration 7 - deviance = 1538.475 - criterion = 2.534822e-13
## converged

Con el modelo ajustado podemos explorar qué tan bien logra clasificar las observaciones de los datos de entrenamiento. A continuación el código para calcular las probabilidades estimadas \(\hat{P}\), las respuestas estimadas \(\widehat{Resp}\) por el modelo, una tabla de confusión tradicional y la tasa de precisión del modelo con los datos de entrenamiento.

# Probabilidades estimadas para cada nivel de Response
probs <- predict(object=mod1, type="response")
# Para identificar el nivel de Response con mayor probabilidad
aux_fun <- function(x) names(x)[which.max(x)]
Response_hat <- apply(X=probs, MARGIN=1, FUN=aux_fun)
# Matriz de confusion
Response_hat <- factor(Response_hat, levels=levels(pict$Response))
tabla <- table(True_Response=pict$Response, Response_hat)
tabla

##              Response_hat
## True_Response Bare Numeral Quant QuantAdj
##      Bare      308      55     0       31
##      Numeral    78     429     0        3
##      Quant       6      26     0        2
##      QuantAdj   45      37     0       70

# Accuracy
sum(diag(tabla)) / sum(tabla)

## [1] 0.740367

Es posible usar el modelo ajustado y la misma función predict para hacer predicciones de la respuesta a items que estén en la base de datos pero con algunos cambios en las covariables. Supongamos que queremos predecir la respuesta nominal para los items 1, 3, 9 y 11 cuando los valores de las covariables son las que se muestran a continuación en el objeto newdata.

newdata <- data.frame(Item=factor(c(1, 3, 9, 11)), 
                      Gender=factor(c("Male", "Female", "Female", "Male")),
                      Group=factor(c("English_Mono", "Geman_Mono", "L2_Advanced", "L2_Intermediate")),
                      Age=c(18, 20, 23, 17))
newdata

##   Item Gender           Group Age
## 1    1   Male    English_Mono  18
## 2    3 Female      Geman_Mono  20
## 3    9 Female     L2_Advanced  23
## 4   11   Male L2_Intermediate  17

Para hacer la predicción usamos la función predict de la siguiente manera.

probs <- predict(object=mod1, newdata=newdata, type="response")
probs

##            Bare     Numeral       Quant    QuantAdj
## [1,] 0.09863091 0.782680376 0.005847403 0.112841308
## [2,] 0.54910245 0.433187674 0.008833600 0.008876273
## [3,] 0.40685403 0.007197128 0.030881057 0.555067783
## [4,] 0.80135086 0.130020315 0.010356409 0.058272417

En el objeto anterior están las probabilidades de elegir Bare, Numeral, Quant o QuantAdj para los cuatro items. La probabilidad con mayor valor indica la posible elección de cada item. Para extraer la respuesta asociada con la mayor probabilidad podemos hacer lo siguiente.

aux_fun <- function(x) names(x)[which.max(x)]
apply(X=probs, MARGIN=1, FUN=aux_fun)

## [1] "Numeral"  "Bare"     "QuantAdj" "Bare"

De la salida anterior se obtienen las posibles respuestas de cada item cuando los valores de sus covariables son las del objeto newdata.

Cuando se use predict con un nuevo conjunto de datos, debemos asegurarnos de que la estructura del nuevo conjunto de datos coincida con la estructura de los datos de entrenamiento. Es decir, que si usa covariable usad en el modelo es un factor, en el nuevo conjunto de datos debe ser un factor también.

Ejemplo: modelo multinomial ordinal

Este ejemplo fue tomado de la sección 10.3 del libro de Agresti (2018). En este ejemplo se analizan los resultados de un ensayo clínico que comparó un fármaco hipnótico activo con un placebo en dos ocasiones en el tratamiento de pacientes con insomnio. La respuesta (hora de conciliar el sueño) cae en una de cuatro categorías ordenadas que son 1 (< 20 min), 2 (20–30 min), 3 (30-60 min) y 4 (> 60 min). Como la respuesta de cada paciente está en una de las cuatro categorías ordenadas, es que se utiliza un modelo multinomial ordinal.

El código para acceder a la base de datos se muestra a continuación.

Insomnia <- read.table("http://www.stat.ufl.edu/~aa/cat/data/Insomnia.dat",
                       header=TRUE)
Insomnia$response <- factor(Insomnia$response) # response var for clmm must be factor
head(Insomnia)

##   case treat occasion response
## 1    1     1        0        1
## 2    1     1        1        1
## 3    2     1        0        1
## 4    2     1        1        1
## 5    3     1        0        1
## 6    3     1        1        1

Solución.

El modelo de interés en este ejemplo es

\[ \text{logit}[P(Y_{it} \leq j) ] = b0_i + \beta_{0j} + \beta_1 t + \beta_2 x + \beta_3 (t \times x), \]

siendo \(i\) el subíndice que identifica a los pacientes \(i=1, 2, \ldots, 239\) y \(j\) el subíndice que identifica la respuesta ordenada 1, 2, 3 o 4. La variable \(t\) se refiere a la ocasión y \(x\) al tratamiento. El intercepto aleatorio se denota como \(b0\) y el intercepto fijo para cada nivel de \(Y\) se denota por \(\beta_{0j}\).

Vamos ahora a crear el modelo multinomial ordinal con la ayuda de la función clmm del paquete ordinal de Christensen (2023). Las covariables que vamos a utilizar para predecir a response son occasion, treat y la interacción entre ambas. La variable case se refiere al paciente y esta variable se va a utilizar como fuente del intercepto aleatorio \(b0\) que se va a incluir en el modelo. El intercepto aleatorio se le indica a clmm por medio del término (1|case) en la fórmula del modelo.

library(ordinal)
fit <- clmm(response ~ (1|case) + occasion + treat + occasion:treat,
            nAGQ=20, data=Insomnia) 
summary(fit)

## Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite 
## quadrature approximation with 20 quadrature points 
## 
## formula: response ~ (1 | case) + occasion + treat + occasion:treat
## data:    Insomnia
## 
##  link  threshold nobs logLik  AIC     niter     max.grad cond.H 
##  logit flexible  478  -592.97 1199.94 401(1549) 1.03e-04 6.1e+01
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  case   (Intercept) 3.628    1.905   
## Number of groups:  case 239 
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## occasion       -1.60158    0.28336  -5.652 1.59e-08 ***
## treat          -0.05785    0.36630  -0.158  0.87450    
## occasion:treat -1.08130    0.38046  -2.842  0.00448 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2  -3.4896     0.3588  -9.727
## 2|3  -1.4846     0.2903  -5.114
## 3|4   0.5613     0.2701   2.078

Los efectos fijos se puden obtener usando la función coef.

coef(fit)

##            1|2            2|3            3|4       occasion          treat 
##    -3.48959549    -1.48463671     0.56127844    -1.60158084    -0.05785368 
## occasion:treat 
##    -1.08129512

Los efectos aleatorios se puden obtener usando la función ranef.

b0 <- ranef(fit)$case

Usando la información anterior se puede escribir el modelo ajustado de la siguiente manera:

\[ \text{logit}[P(Y_{it} \leq 1) ] = b0_i -3.49 -1.60 t -0.06 x -1.08 (t \times x), \] \[ \text{logit}[P(Y_{it} \leq 2) ] = b0_i -1.49 -1.60 t -0.06 x -1.08 (t \times x), \]

\[ \text{logit}[P(Y_{it} \leq 3) ] = b0_i + 0.56 -1.60 t -0.06 x -1.08 (t \times x), \] El intercepto aleatorio \(\tilde{b}0\) predicho está almacenado en el objeto b0 y se puede utilizar para completar las expresiones anteriores y así calcular las probabilidades.

Ejemplo: simulando observaciones de un modelo multinomial nominal

En este ejemplo vamos a simular observaciones \(n_i=50\) observaciones para \(G=5\) grupos que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar la forma para estimar los parámetros del modelo.

\[\begin{align*} y_{ij} | b_0 &\sim Multinomial(\left\lbrace 1, 2, 3 \right\rbrace, \left\lbrace \pi_1, \pi_2, \pi_3 \right\rbrace) \\ \pi_{1i} &= \frac{e^{1.62-0.11x_{ij}+b0}}{1+e^{1.62-0.11x_{ij}+b0}+e^{5.70.2.47x_{ij}+b0}} \\ \pi_{2i} &= \frac{e^{5.70-2.47x_{ij}+b0}}{1+e^{1.62-0.11x_{ij}+b0}+e^{5.70-2.47x_{ij}+b0}} \\ \pi_{3i} &= \frac{1}{1+e^{1.62-0.11x_{ij}+b0}+e^{5.70-0.2.47x_{ij}+b0}} \\ b_{0} &\sim N(0, \sigma^2_{b0}=4) \\ x_{ij} &\sim U(1, 4) \end{align*}\]

El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(1.62, -0.11, 5.70, -2.47, \sigma_{b0}=2)^\top\).

Solución.

El código para simular las 500 observaciones se muestra a continuación. Observe que se fijó la semilla en dos ocasiones para que el lector pueda replicar el ejemplo y obtener los mismos resultados.

ni <- 10
G <- 200
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
x <- runif(n=nobs, min=1.2, max=3.9)
b0 <- rnorm(n=G, mean=0, sd=sqrt(0)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni)             # El mismo intercepto aleatorio pero repetido
num1 <- exp(1.62-0.11*x+b0)
num2 <- exp(5.70-2.47*x+b0)
num3 <- 1
den <- num1 + num2 + num3
pi1 <- num1 / den
pi2 <- num2 / den
pi3 <- num3 / den
probs <- cbind(pi1, pi2, pi3)
y <- apply(X=probs, MARGIN=1, which.max)
y <- factor(y, levels=3:1)
datos <- data.frame(obs, grupo, x, y)
table(datos$y)

## 
##    3    2    1 
##    0  417 1583

References

Agresti, Alan. 2018. An Introduction to Categorical Data Analysis. John Wiley & Sons.

Christensen, Rune Haubo Bojesen. 2023. Ordinal: Regression Models for Ordinal Data. https://github.com/runehaubo/ordinal.

Elff, Martin. 2022. Mclogit: Multinomial Logit Models, with or Without Random Effects or Overdispersion. http://mclogit.elff.eu.