Trazado de los resultados del análisis de regresión logística ordenada

Question

He intentado representar los resultados de un análisis de regresión logística ordenada calculando las probabilidades de aprobar cada categoría de respuesta de la variable dependiente (escala Likert de 6 puntos, que va de "1" a "6"). Sin embargo, he recibido probabilidades extrañas cuando he calculado las probabilidades basándome en esta fórmula: $\rm{Pr}(y_i \le k|X_i) = \rm{logit}^{-1}(X_i\beta)$ .

A continuación puede ver cómo he intentado calcular exactamente las probabilidades y trazar los resultados del modelo de regresión logística ordenada ( m2 ) que instalé utilizando el polr función ( MASS paquete). Las probabilidades ( probLALR ) que calculé y utilicé para trazar una "puntuación media esperada" son desconcertantes ya que la puntuación media esperada en el gráfico aumenta a lo largo del continuo RIV.st mientras que el coeficiente para RIV.st es negativo (-0,1636). Habría esperado que la puntuación media esperada disminuyera debido al efecto principal negativo de RIV.st y la irrelevancia de los términos de interacción para la condición de baja admiración y baja rivalidad (LALR) del actual diseño 2 por 2 (primer factor = f.adm ; segundo factor = f.riv ; codificación ficticia 0 y 1).

¿Alguna idea de cómo dar sentido al patrón encontrado? ¿Es ésta la forma correcta de calcular las probabilidades? La forma en que utilicé los interceptos en la fórmula para calcular las probabilidades podría ser problemática (cf, Coeficiente negativo en la regresión logística ordenada ).

m2 <- polr(short.f ~ 1 + f.adm*f.riv + f.adm*RIV.st + f.riv*RIV.st, data=sampleNS)

# f.adm  = dummy (first factor of 2 by 2 design);
# f.riv  = dummy (second factor of 2 by 2 design);
# RIV.st = continuous predictor (standardized)
summary(m2)
Coefficients:
                Value Std. Error t value
f.adm1         1.0203    0.14959  6.8203
f.riv1        -0.8611    0.14535 -5.9240
RIV.st        -0.1636    0.09398 -1.7403
f.adm1:f.riv1 -1.2793    0.20759 -6.1625
f.adm1:RIV.st  0.0390    0.10584  0.3685
f.riv1:RIV.st  0.6989    0.10759  6.4953

Intercepts:
    Value    Std. Error t value 
1|2  -2.6563   0.1389   -19.1278
2|3  -1.2139   0.1136   -10.6898
3|4  -0.3598   0.1069    -3.3660
4|5   0.9861   0.1121     8.7967
5|6   3.1997   0.1720    18.6008

Aquí puedes ver cómo he intentado calcular las probabilidades ( probLALR ) para 1 de las 4 condiciones del diseño 2 por 2:

inv.logit  <- function(x){ return(exp(x)/(1+exp(x))) }
Pred       <- seq(-3, 3, by=0.01)
b = c(-2.6563,-1.2139,-0.3598,0.9861,3.1997) # intercepts of model m2
a = c(1.0203,-0.8611,-0.1636,-1.2793,0.0390,0.6989) # coefficients of m2
probLALR   <- data.frame(matrix(NA,601,5))
for (k in 1:5){ 
    probLALR[,k] <- inv.logit(b[k] + a[1]*0 + a[2]*0 + 
                               a[3]*Pred  + a[4]*0*0 + 
                               a[5]*Pred*0 + a[6]*Pred*0)
}

plot(Pred,probLALR[,1],type="l",ylim=c(0,1))   # p(1)
lines(Pred,probLALR[,2],col="red")             # p(1 or 2)
lines(Pred,probLALR[,3],col="green")           # P(1 or 2 or 3)
lines(Pred,probLALR[,4],col="orange")          # P(1 or 2 or 3 or 4)
lines(Pred,probLALR[,5],col="orange")          # P(1 or 2 or 3 or 4 or 5)

# option response functions:

orc = matrix(NA,601,6)
orc[,6] = 1-probLALR[,5]        # prob of 6
orc[,5]= probLALR[,5]-probLALR[,4]  # prob of 5
orc[,4]= probLALR[,4]-probLALR[,3]  # prob of 4
orc[,3]= probLALR[,3]-probLALR[,2]  # prob of 3
orc[,2]= probLALR[,2]-probLALR[,1]  # prob of 2
orc[,1]= probLALR[,1]           # prob of 1

plot(Pred,orc[,1],type="l",ylim=c(0,1))   # p(1)
lines(Pred,orc[,2],col="red")             # p(2)
lines(Pred,orc[,3],col="green")           # P(3)
lines(Pred,orc[,4],col="orange")          # P(4)
lines(Pred,orc[,5],col="purple")          # P(5)
lines(Pred,orc[,6],col="purple")          # P(6)

# mean score

mean = orc[,1]*1+orc[,2]*2+orc[,3]*3+orc[,4]*4+orc[,5]*5+orc[,6]*6
plot(Pred,mean,type="l",xlab="RIV.st",ylab="expected mean score",ylim=c(1,6))

Answer 1

Gracias a Wilco Emons por la siguiente solución al problema:

En polr el modelo de enlace acumulativo está parametrizado de forma un poco diferente que en el libro de Agresti al que se hace referencia. El problema se puede resolver fácilmente cambiando la línea de código:

probLALR[,k] <- inv.logit(b[k] +  a[1]*0 + a[2]*0 + a[3]*Pred + a[4]*0*0 + 
                                  a[5]*Pred*0 + a[6]*Pred*0                ) 

en

probLALR[,k] <- inv.logit(b[k] - (a[1]*0 + a[2]*0 + a[3]*Pred + a[4]*0*0 + 
                                  a[5]*Pred*0 + a[6]*Pred*0)               ) 

Gracias también a Achim Zeileis por su sugerencia de utilizar predict(m2, type="prob") ¡! A continuación encontrará una forma de calcular las probabilidades mediante el predict() función:

Pred             <- seq(-3, 3, by=0.01)
PRED.LALR        <- data.frame(rep(NA,601))
PRED.LALR$f.adm  <- as.factor(rep(0,601))
PRED.LALR$f.riv  <- as.factor(rep(0,601))
PRED.LALR$RIV.st <- Pred
prob.LALR        <- predict(m2,PRED.LALR,type="prob")

scoreLALR        <- prob.LALR[,1]*1 + prob.LALR[,2]*2 + prob.LALR[,3]*3 + 
                    prob.LALR[,4]*4 + prob.LALR[,5]*5 + prob.LALR[,6]*6
plot(Pred, scoreLALR, col="green", ylim=c(1,6))