Proporcional de probabilidades de asunción en el ordinal regresión logística en R con los paquetes VGAM y rms

Question

La suposición de que el ordinal de la regresión logística es la parte proporcional de probabilidades de asunción. El uso de R y el 2 paquetes mencionados tengo 2 formas de comprobar eso, pero tengo preguntas en cada uno de ellos.

1) Utilizando el paquete rms

Dada la siguiente comandos

library(rms)
ddist <- datadist(Ki67,Cyclin_E)
options(datadist='ddist')
f <- lrm(grade ~Ki67+Cyclin_E);f
sf <- function(y)
c('Y>=1'=qlogis(mean(y >= 1)),'Y>=2'=qlogis(mean(y >= 2)),'Y>=3'=qlogis(mean(y >= 3)))
s <- summary(grade ~Ki67+Cyclin_E, fun=sf)
plot(s,which=1:3,pch=1:3,xlab='logit',main='',xlim=c(-2.5,2.5))

Tengo

lrm(formula = grade ~ Ki67 + Cyclin_E)

Frequencies of Missing Values Due to Each Variable
   grade     Ki67 Cyclin_E 
       0        0        3 


                     Model Likelihood     Discrimination    Rank Discrim.    
                        Ratio Test            Indexes          Indexes       

Obs            42    LR chi2     11.38    R2       0.268    C       0.728    
 1             11    d.f.            2    g        1.279    Dxy     0.456    
 2             15    Pr(> chi2) 0.0034    gr       3.592    gamma   0.458    
 3             16                         gp       0.192    tau-a   0.308    
max |deriv| 1e-07                         Brier    0.166                     


         Coef    S.E.   Wald Z Pr(>|Z|)
y>=2     -0.1895 0.8427 -0.22  0.8221  
y>=3     -2.0690 0.9109 -2.27  0.0231  
Ki67      0.0971 0.0330  2.94  0.0033  
Cyclin_E -0.0076 0.0227 -0.33  0.7387 

El s tabla de: (lamentablemente no sé cómo subir un gráfico en I)

grade    N=45

+--------+-------+--+----+---------+----------+
|        |       |N |Y>=1|Y>=2     |Y>=3      |
+--------+-------+--+----+---------+----------+
|Ki67    |[ 2, 9)|12|Inf |0.6931472|-1.0986123|
|        |[ 9,16)|12|Inf |0.3364722|-2.3978953|
|        |[16,24)|10|Inf |2.1972246| 0.0000000|
|        |[24,44]|11|Inf |2.3025851| 1.5040774|
+--------+-------+--+----+---------+----------+
|Cyclin_E|[ 3,16)|15|Inf |1.0116009|-0.1335314|
|        |[16,22)| 7|Inf |1.7917595|-0.9162907|
|        |[22,33)|10|Inf |1.3862944|-0.8472979|
|        |[33,80]|10|Inf |0.4054651|-0.4054651|
|        |Missing| 3|Inf |      Inf| 0.6931472|
+--------+-------+--+----+---------+----------+
|Overall |       |45|Inf |1.1284653|-0.4054651|
+--------+-------+--+----+---------+----------+

Donde para el Ki67 veo que 3 de las 4 diferencias logit(P[Y> = 2])-logit(P[Y> = 3]) son cerca de 2. Sólo el último es muy baja (alrededor de 0.8). Pero aquí Ki67 es continua y no categórica, así que no sé si los resultados de la tabla son correctos y no hay ningún valor de p para decidir. Por mi forma de correr lo anterior en el programa SPSS y yo no se rechaza la hipótesis.

2) Utilizando el paquete VGAM

Aquí el uso de los siguientes comandos tengo el modelo bajo el supuesto de probabilidades proporcionales

library(VGAM)
fit1 <- vglm(grade ~Ki67+Cyclin_E,family=cumulative(parallel=T))
summary(fit1)

Y los resultados

Coefficients:
                Estimate Std. Error  z value
(Intercept):1  0.1894723   0.820442  0.23094
(Intercept):2  2.0690395   0.886732  2.33333
Ki67          -0.0970972   0.032423 -2.99467
Cyclin_E       0.0075887   0.021521  0.35261

Number of linear predictors:  2 

Names of linear predictors: logit(P[Y< = 1]), logit(P[Y< = 2])

Dispersion Parameter for cumulative family:   1

Residual deviance: 79.86801 on 80 degrees of freedom

Log-likelihood: -39.93401 on 80 degrees of freedom

Number of iterations: 5 

Mientras que el uso de los siguientes comandos tengo el modelo sin la asunción de probabilidades proporcionales

fit2 <- vglm(grade ~Ki67+Cyclin_E,family=cumulative(parallel=F))

donde por desgracia yo receice el siguiente mensaje

Mensaje de advertencia: En vglm.en forma(x = x, y = y), w = w, offset = offset, Xm2 = Xm2, : convergencia no se obtiene en 30 iteraciones

Sin embargo, si el tipo I summary(fit2) obtengo los resultados, pero de nuevo no sé si son correctos. Mi intención era utilizar los siguientes comandos y obtener la respuesta, pero sé que yo a dudar de si esto es correcto (por cierto, si lo hago llego p-value=0.6.

pchisq(deviance(fit1)-deviance(fit2),
df=df.residual(fit1)-df.residual(fit2),lower.tail=FALSE)

Así que, respecto de los métodos mencionados anteriormente, ¿alguien sabe si los resultados que obtengo son válidos o, en el caso de la VGAM paquete es hay alguna manera de aumentar el número de itterations?¿Hay alguna otra manera de comprobarlo?

Answer 1

Recomiendo parcial residual de las parcelas con el rms del paquete lrm y residuals.lrm funciones. Usted también puede ajustar una serie de modelos binarios utilizando diferentes puntos de corte para $Y$ y trazar el registro de los odds ratios frente de corte.

Answer 2

Hay errores estándar robustos disponibles en el paquete? o el sándwich de paquete de calcular para usted?

El sucio secreto acerca de la comprobación del modelo es casi nunca hay buena potencia estadística para detectar significativas violaciones de su modelo. Usted todavía puede comprobar que ellos, pero más como post-hoc de análisis de sensibilidad.

El mejor es el uso de métodos que son agnósticos en cuanto a si los supuestos bajo los cuales el modelo es 'óptimo' son correctos. E. g. de riesgos proporcionales casi siempre son violados en cox modelos lineales y modelos más probable es que nunca se corresponden con la verdadera base lineal de la estructura. Sin embargo, con errores estándar robustos, no importa, y que todavía puede hacer declaraciones acerca de las tendencias.