Definiendo una regresión multinomial previa. Estudio de caso con `MCMCglmm`

Question

Me resulta difícil entender cómo definir la distribución previa para una regresión multinomial.

• ¿En qué unidad se debe establecer la probabilidad previa dado que la respuesta realmente no tiene unidades (solo categorías)? ¿Debería ser una probabilidad, o tal vez un logaritmo de probabilidades?
• ¿Depende las dimensiones de la prior del número de niveles en la variable de respuesta categórica?
• ¿Cuál es el significado de la matriz de varianza-covarianza para una variable de respuesta categórica?

---

Aquí hay un subconjunto de mis datos (donde los nombres de las variables y resultados han sido cambiados). ¡Sería genial un ejemplo de trabajo con MCMCglmm en estos datos!

df=read.table(text="y     x1       x2
    1  amarillo 106.00 6.190476
    2  amarillo 120.00 5.254762
    3  amarillo  57.00 6.202381
    4  amarillo 115.33 5.652381
    5  amarillo 175.00 6.154762
    6  amarillo  74.00 8.285714
    7  amarillo 104.67 3.766667
    8  amarillo  95.50 7.976190
    9  amarillo 108.00 8.792857
    10 amarillo 121.33 7.935714
    11 amarillo  66.67 6.969048
    12 amarillo  30.00 7.333333
    13 amarillo  45.00 6.811905
    14 amarillo  70.00 7.550000
    15 amarillo  48.00 7.316667
    16 amarillo 211.00 4.650000
    17 amarillo  69.00 8.369048
    18 amarillo 110.50 6.621429
    19 amarillo 203.00 6.095238
    20 amarillo  75.33 8.211905
    21 amarillo 207.33 6.211905
    22 amarillo  54.00 7.961905
    23 amarillo  74.00 7.019048
    24 amarillo 113.00 4.221429
    25 amarillo  23.00 7.942857
    26 amarillo  80.00 7.511905
    27 amarillo 257.00 7.878571
    28 amarillo 211.00 7.754762
    29 amarillo  99.00 8.016667
    30 amarillo 120.00 7.728571
    31 amarillo 222.50 5.840476
    32 amarillo  44.00 4.209524
    33 amarillo  63.00 6.614286
    34 amarillo  57.00 8.669048
    35 amarillo 223.33 7.033333
    36 amarillo 128.00 6.754762
    37 amarillo 128.00 5.561905
    38 amarillo 121.00 7.471429
    39 amarillo  70.00 7.445238
    40 amarillo  85.67 5.261905
    41 amarillo 113.33 8.509524
    42 amarillo  82.00 6.697619
    43      rojo 207.33 4.180952
    44      rojo 167.67 5.302381
    45      rojo 366.50 7.102381
    46      rojo 230.00 4.942857
    47      rojo 201.00 5.754762
    48      rojo 226.00 9.076190
    49      rojo 193.33 7.066667
    50      rojo 170.00 7.314286
    51      rojo 361.33 7.502381
    52      azul 154.00 4.342857
    53      rojo 199.33 6.361905
    54      azul  97.00 7.750000
    55      azul  82.33 6.209524
    56      azul  55.67 5.321429
    57      azul  47.50 5.911905
    58      azul  15.67 7.185714
    59      azul  96.50 6.452381
    60      azul 202.33 8.576190
    61      azul 157.00 6.669048
    62      azul 117.33 5.828571
    63      azul 105.67 8.485714
    64      azul 108.67 5.714286
    65      azul 296.67 5.852381
    66      azul 206.50 6.826190
    67      azul  88.50 6.178571
    68      azul 163.00 7.833333
    69      azul 151.50 8.983333")

y aquí hay una llamada a MCMCglmm para la cual las priors predeterminadas conducen a un mensaje de error

set.seed(12)
m = MCMCglmm(y ~  -1  + trait:(x1) + trait:(x2)  , rcov = ~ us(trait):units,
 data = df, family = "categorical", verbose = TRUE, burnin = 8000,
 nitt = 40000, thin = 50)

Estructura G/R mal condicionada (CN = 24007848728601288.000000):
 Use priors adecuadas si no las tiene o reescale los datos si los tiene

Answer 1

¡Buena pregunta! En realidad he intentado comprender lo mismo durante bastante tiempo, así que simplemente compartiré mis experiencias. Todavía soy un novato en esta área, así que espero no haber cometido errores en la notación.

siempre, y quiero decir siempre, estandariza tus variables continuas para que tengan una media=0 y una sd=1 (o incluso sd=2). Consulta algunos de los posts o artículos de Andrew Gelman sobre esto. Solo busca Andrew Gelman estandarización, hay muchos buenos papers y posts.

Considera tus coeficientes como log(odds-ratios) (en referencia a una categoría de referencia, la explicación sigue). Para una discusión más detallada, consulta esta respuesta. Andrew Gelman también tiene algunas recomendaciones sobre priors, como la cauchy, o normal(0,1). Sus papers son sobre regresión logística, pero encuentro que estas recomendaciones también se extienden a la regresión multi-resultado.

La dimensión del prior realmente depende del número de resultados. Si tienes tres resultados, tienes estas tres dependencias lineales:

$ y_1 = \beta_{0,1} + \sum_i\beta_{i,1}*x_{i,1} $ y $ y_2 = \beta_{0,2} + \sum_i\beta_{i,2}*x_{i,2} $ y $ y_3 = \beta_{0,3} + \sum_i\beta_{i,3}*x_{i,3} $

El segundo subíndice denota el resultado. Normalmente pondrías un prior en cada uno de estos coeficientes. Ilustraré con la intersección, $\beta_{0,k}$,

$\beta_{0,1} \sim \mathcal{N}(0,1)$ y $\beta_{0,2} \sim \mathcal{N}(0,1)$ y $\beta_{0,3} \sim \mathcal{N}(0,1)$

Solo para completar las matemáticas, la probabilidad del resultado $k$ es,

$p_k = \frac{y_k}{\sum_k y_k}$ y la verosimilitud es $y \sim \text{Multinomial}(p)$

Por favor ten en cuenta que a menudo obligarías a que uno de los resultados sea el resultado "de referencia" debido a problemas de identificabilidad. Wiki tiene una descripción detallada de por qué. Prácticamente, esto significa que obligas a que los coeficientes de la categoría de referencia sean 0, y por lo tanto $p_{referencia} = 1$, ya que $exp(0) = 1$.

Nunca he usado MCMCglmm así que no puedo responder nada específico sobre eso, lo siento.