Regresión ordinal y sus aplicaciones

Regresi´on ordinal y sus aplicaciones FACULTAD DE MATEM´ ATICAS DEPARTAMENTO DE ESTAD´ ISTICA E INVESTIGACI´ ON OPERATIVA Dirigido por: Juan Manuel Mu˜noz Pichardo Miguel Arias Ben´ıtez Junio 2018 “It only ends once. Anything that happens before that is just progress.” Jacob, Lost. ´ Indice general Resumen ....................................... v Summary ....................................... vii 1. Introducci´on a la Regresi´on Categ´orica 1 1.1. Distribuci´on Multinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1. Familia exponencial . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. Modelo lineal generalizado . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1. Estimaci´on de los par´ametros . . . . . . . . . . . . . . . . . . . . . 5 1.2.2. Bondaddeajuste............................ 7 1.2.3. Regiones de confianza para β..................... 8 1.2.4. Residuos................................. 9 1.3. Modelos logit para respuestas nominales . . . . . . . . . . . . . . . . . . . 11 2. Modelos de regresi´on ordinal 13 2.1. Modelo de Odds Proporcionales . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.1. Presentaci´on del modelo . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.2. Odds Ratios e Intervalos de Confianza . . . . . . . . . . . . . . . . 16 2.1.3. Extensi´on del Modelo Ordinal a kvariables . . . . . . . . . . . . . 17 2.1.4. Funci´on de Probabilidad para el Modelo Ordinal . . . . . . . . . . 18 2.1.5. Estimaci´on de par´ametros . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.6. Residuos................................. 19 2.2. Otros modelos ordinales de inter´es . . . . . . . . . . . . . . . . . . . . . . 20 2.2.1. Modelo de Ratios Continuados . . . . . . . . . . . . . . . . . . . . 20 2.2.2. Modelo Logit de Categor´ıas Adyacentes . . . . . . . . . . . . . . . 20 3. Modelado en R 23 3.1. Librer´ıa“ordinal”................................ 23 3.2. EjemploenR.................................. 25 3.2.1. Datos wine ............................... 25 3.2.2. Predicciones sobre wine ........................ 30 4. Ejemplos Reales 31 4.1. Regresi´on Log´ıstica Ordinal Aplicada a la Identificaci´on de Factores de Riesgo para C´ancer de Cuello Uterino . . . . . . . . . . . . . . . . . . . . 31 iii iv 4.2. Predicci´on del rendimiento en una asignatura empleando la regresi´on log´ısticaordinal................................. 31 4.3. El impacto de las relaciones interpersonales en la satisfacci´on laboral general 32 4.4. Las percepciones de riesgo de los consumidores en alimentos l´acteos: aplicaci´on de una regresi´on log´ıstica ordinal . . . . . . . . . . . . . . . . . . . 32 Bibliograf´ıa 33 Resumen El objetivo de este trabajo es la construcci´on de una base te´orica con la que formular y desarrollar la forma general de algunos modelos de regresi´on ordinal, a su vez que motivar el estudio con una serie de ejemplos reales donde se emplearon los m´etodos descritos. Para ello, a modo de introducci´on en el primer cap´ıtulo, se exponen una serie de conceptos te´oricos sobre la idea estad´ıstica de la regresi´on; desde la definici´on y algunos conceptos de la distribuci´on multinomial, hasta una descripci´on breve del modelo que tomaremos como base, el Modelo Lineal Generalizado. En el cap´ıtulo segundo se presenta el Modelo de Odds Proporcionales, tambi´en conocido como Modelo Logit Acumulado. Primero se estudia e ilustra el caso bivariante para posteriormente su extensi´on a k-variables, donde se presenta la funci´on de probabilidad y sus residuos, a su vez que se estiman los par´ametros y se construyen los intervalos de confianza. En el tercer cap´ıtulo se trabajar´a sobre un conjunto de datos en R para ilustrar la construcci´on y los resultados al aplicar las t´ecnicas de regresi´on descritas en estos. Los datos empleados representan un experimento sobre ciertos factores que determinan la acidez en el vino; en nuestro caso, la temperatura y el contacto entre el zumo y las pieles de las uvas cuando se extrae de ellas. Por ´ultimo, en el cap´ıtulo cuarto, se describen una serie de estudios reales de diversas ´ındoles sobre los que se aplicaron m´etodos de regresi´on ordinales con el objetivo de extraer una serie de conclusiones, por ejemplo, para la identificaci´on de factores de riesgo en el c´ancer de cuello uterino. v Summary The main objective of this work is the construction of a theoretical basis with which to formulate and develop the general form of some models of ordinal regression, at the same time as motivating the study with some real examples where the described methods were used. For this, as an introduction, in the first chapter some theoretical concepts on the statistical idea of regression are exposed; from the definition and some concepts of the multinomial distribution, to a brief description of the model that we will take as a base, the Generalized Linear Model. In the second chapter we present the Proportional Odds Model, also known as the Accumulated Logit Model. First, the bivariate case is studied and illustrated for its extension to kvariables, where the probability function and its residuals are presented, while the parameters are estimated and confidence intervals are constructed. In the third chapter we will work on a set of data in R to illustrate the construction and the results when applying the regression techniques described. The data used represent an experiment on certain factors that determine the bitterness in the wine; in our case, the temperature and the contact between the juice and the skins of the grapes when it is extracted from them. Finally, in the fourth chapter, we describe some real studies of various types on which ordinal regression methods were applied in order to dig out a series of conclusions, for example, for the identification of risk factors in cancer of the cervix. vii 6Introducci´on a la Regresi´on Categ´orica donde            µi(β) = h(xt iβ) σ2 i(β) = a(φ)v(h(xt iβ)) V(µ) = ∂2b(θ)/∂σ2 Di(β) = ∂h(xt iβ)/∂η con ηi=xt iβ Otros conceptos importantes a tener en cuenta en la estimaci´on m´aximo-veros´ımil del vector de par´ametros son: Matriz de informaci´on de Fisher esperada: F(β) = Cov s(β) = X i Fi(β) Fi(β) = xixt iwi(β)wi(β) = D2 i(β)σ−2 i(β) Matriz de Fisher observada: Fobs(β) = −∂2l(β) ∂β∂βt verific´andose que F(β) = E(Fobs(β)) Para las funciones de enlace naturales σ(µi) = xt iβ, las matrices se simplifican de la forma: s(β) = 1 a(φ)X i xi[yi−µi(β)] F(β) = 1 a(φ)X i V(µi(β))xixt iF(β) = Fobs(β) La obtenci´on de la estimaci´on de m´axima-verosmilitud se plantea generalmente como las soluciones de la ecuaci´on de verosimiltud s(ˆ β) = 0 lo que corresponde a un m´aximo local, es decir, con la matriz de segundas derivadas Fobs(ˆ β) definida positiva. Las ecuaciones resultantes no suelen ser lineales y para resolverlas necesitaremos de m´etodos num´ericos iterativos como el de Fisher Scoring o el de M´ınimos Cuadrados Ponderados Iterativos, cuyas iteraciones se definen a partir de un estimador inicial ˆ β0por: ˆ β(k+1) =ˆ β(k)+F−1(ˆ β(k))s(ˆ β(k))con :k= 0,1,2, ... Obteniendo a trav´es de estos m´etodos las estimaciones de los par´ametros del modelo ˆ β, estas estimaciones m´aximo-veros´ımiles tienen las propiedades de consistencia, eficiencia asint´otica y distribuci´on normal asint´otica. Introducci´on a la Regresi´on Categ´orica 7 En el caso en el que el par´ametro de dispersi´on sea desconocido, se puede considerar el siguiente estimador consistente: ˆ φ=1 n−pX i [yi−µi(ˆ β)]2 v(µi(ˆ β)) 1.2.2. Bondad de ajuste Una vez estimados los par´ametros debemos valorar cuan bueno es nuestro modelo, es decir valorar la discrepancia entre los datos observados y los datos esperados. De esta manera, determinar cuantos t´erminos son necesarios en la estructura lineal para una descripci´on ´optima de los datos intentando no saturarlo de variables explicativas que har´an un modelo bien ajustado pero de dificil comprensi´on ni, en caso contrario, un defecto de variables que har´an un modelo de f´acil interpretaci´on pero de pobre ajuste es otro de los problemas que nos plantearemos. Trataremos de constuir un modelo intermedio entre el modelo saturado y el modelo nulo, donde el primero se refiere al modelo en el que el n´umero de par´ametros es igual al n´umero de observaciones (ninguna simplificaci´on) y el segundo que es el modelo m´as simple en el que solo se usa el par´ametro µ, el valor esperado para todas las observaciones (simplificaci´on total, asume efecto nulo de las variables explicativas) En el modelo lineal generalizado, la bondad de ajuste se puede evaluar de distintas formas, entre ellas destacan: La funci´on o estad´ıstico desviaci´on D(y;µ)=2{l(y;y)−l(ˆµ;y)} Es la distancia entre el logaritmo de la funci´on verosimilitud del modelo saturado y el modelo con el que se est´a trabajando. Un valor peque˜no de la desviaci´on indica que para un n´umero menor de par´ametros, se obtiene un ajuste tan bueno como cuando se ajusta el modelo saturado. Si el modelo es correcto el estad´ıstico se distribuye asint´oticamente segun una χ2 n−p con n−pgrados de libertad [15]. D(y, ˆµ)∼χ2 n−p 8Introducci´on a la Regresi´on Categ´orica Coeficiente de determinaci´on R2: La medida R2se define como la proporci´on de la varianza total de la variable explicada por la regresi´on. El R2, tambi´en llamado coeficiente de determinaci´on, refleja la bondad del ajuste de un modelo a la variable que pretender explicar. Este coeficiente viene dado por: R2= 1 −D(y, ˆµ) D(y, ˆµ0) donde D(y, ˆµ) y D(y, ˆµ0) son las funciones de desviaci´on del modelo ajustado y nulo respectivamente. Se verifica que 0 ≤R2≤1 Estad´ıstico Chi-cuadrado de Pearson: χ2= n X i=1 (yi−ˆµi)2 V(ˆµi donde V(ˆµ) es la funci´on varianza estimada para la distribuci´on de la variable objetivo. En cuanto a lo que a procesos de seleci´on de modelos se refiere mencionaremos el Criterio de informaci´on Akaike (AIC), en el caso general el AIC es: AIC =k−2ln(ˆ L) donde kes el n´umero de par´ametros y ˆ Les el m´aximo valor de la funci´on de verosmilitud para el modelo estimado. El modelo se optimiza minimizando el valor de AIC, este modelo recompensa la bondad de ajuste y penaliza el aumento de la cantidad de par´ametros estimados. 1.2.3. Regiones de confianza para β A continuaci´on construiremos un intervalo de confianza realizando inferencias sobre el vector de par´ametros desconocidos βde dimensi´on p, la mayoria de cuestiones consiguen reformularse a traves de una hip´otesis lineal de la forma Cβ, siendo Cuna matriz de rango total s≤pyξun vector de constantes conocido de dimensi´on s Ho:Cβ =ξ H1:Cβ 6=ξ Para este procedimiento se puede usar el Estad´ıstico de Wald entre otros. ´ Este se basa en la distribuci´on asint´otica del vector ˆ βy est´a definido por: ξW= [Cˆ β−ξ]T[CF−1ˆ β(C)0][Cˆ β−ξ] Introducci´on a la Regresi´on Categ´orica 9 determina la distancia ponderada entre el estimador Cˆ βy su valor determinado por la hip´otesis nula. F−1(ˆ β) denota la estimaci´on de la matriz de informaci´on de Fisher de ˆ β Asint´oticamente y bajo hip´otesis nula este estad´ıstico se distribuye como una distribuci´on Chi-cuadrado con sgrados de libertad χ2 s Usando este estad´ıstico, una regi´on de confianza para βcon un nivel de confianza del 100(1 −α) % viene dada por: {β∈Rp|(ˆ β−β)T[V ar(ˆ β)]−1(ˆ β−β)< χ2 p,1−α} 1.2.4. Residuos En la pr´actica podemos encontrar el problema de que a´un habiendo escogido cuidadosamente un modelo, al ajustarlo a un conjunto de datos el resultado sea insatisfactorio. Las desviaciones sistem´aticas se originan por haber escogido inadecuadamente la funci´on de enlace o las variables explicativas. Las diferencias aisladas pueden darse por puntos extremos o porque estos realmente sean err´oneos, ´estos se conocen como outliers. La comprobaci´on de la adecuaci´on del modelo es un requisito fundamental que se realiza sobre el conjunto de datos para encontrar posibles fallos en las suposiciones hechas por el modelo, as´ı como los outliers que puedan interferir desproporcionadamente en los resultados del ajuste. Como en la regresi´on lineal, los residuos son los utilizados para verificar esta adecuaci´on del modelo. Expresan la discrepancia entre una observacci´on y su valor ajustado y tambi´en pueden indicar la presencia de valores an´omalos que puedan requerir un estudio m´as concreto. Entre otros residuos los m´as destacados son: El residuo b´asico: Definido como la diferecnia entre le valor observado, yi, de la variable respuesta y el valor ajustado, ˆyi, por el modelo. rb i=yi−ˆyicon i = 1, ..., n 10 Introducci´on a la Regresi´on Categ´orica El residuo de Pearson: Es la contribuci´on individual al estad´ıstico χ2de Pearson, se define como: rp i=yi−ˆµi qˆ φ V ar( ˆµi) con i = 1, ..., n siendo ˆ φun estimador consistente del par´ametro escala φ. Y su versi´on studentizada viene dada por: rsp i=yi−ˆµi qˆ φ V ar( ˆµi)(1 −hi) con i = 1, ..., n siendo hiel elemento diagonal de la matriz H, donde: H=W1/2X(XTWX)−1XTW1/2 con Wuna matriz diagonal cuyos elementos de la diagonal principal son: wi=1 V ar(µi)∂µi ∂η 2 La ventaja de este residuo studentizado frente al anterior reside en que la captaci´on de la variabilidad de los datos es mejor debido a que usa el valor hi, este es ´util para medir la influencia de cada observaci´on. El residuo desviaci´on: Se define como: rD i=sign(yi−ˆµipdii= 1, ..., n dies el llamado componente desviaci´on,di= 2(l(yi, yi)−l( ˆµi, yi) Y su versi´on studentizada: rsD0 i=rD i qˆ φ(1 −hi) donde hies el i-´esimo elemento de la diagonal de la matriz Hyˆ φes la estimaci´on del par´ametro de escala φ. Introducci´on a la Regresi´on Categ´orica 11 1.3. Modelos logit para respuestas nominales Al igual que se describe en Mar´ın [14], se denota Ccomo el n´umero de categor´ıas de la variable Yy{π1, ..., πj}las probabilidades de cada respuesta, satisfaciendo X j πj= 1. Se parte de nobservaciones independientes extraidas. La distribuci´on de probabilidad del n´umero de observaciones de las C categor´ıas sigue una distribuci´on multinomial tal y como nombramos anteriormente. Esta modeliza la probabilidad de cada una de las posibles maneras en que nobservaciones pueden repartirse entre C categor´ıas. Al ser la respuesta nominal, el orden entre las categor´ıas es irrelevante; este es el problema central que trataremos m´as adelante. Se toma una categor´ıa como respuesta base, por ejemplo la ´ultima categor´ıa (C), y se define un modelo logit con respecto a ella: logπc πC=αc+βcx donde c= 1, ..., C −1. El modelo tiene C−1 ecuaciones con sus propios par´ametros, y los efectos var´ıan con respecto la categor´ıa que se ha tomado como base. Cuando C= 2, el modelo equivale a una ´unica ecuaci´on log(π1/π2) = logit(π1) y se obtiene el modelo de regresi´on log´ıstica est´andar. La ecuaci´on general logit con respecto a la categor´ıa base C determina tambi´en los logits para cualquier pareja de categor´ıas. Considerando c1yc2dos categor´ıas cualesquiera tenemos: logπc1 πc2=logπc1/πC πc2/πC=logπc1 πC−logπc2 πC = (αc1−βc1x)−(αc2−βc2x) = (αc1−αc2)+(βc1−βc2)x. De este modo, la ecuaci´on para las categor´ıas c1yc2tiene tambi´en la forma α+βx donde α= (α= (αc1−αc2) y β= (βc1−βc2)x Cap´ıtulo 2 Modelos de regresi´on ordinal En este cap´ıtulo, la regresi´on log´ıstica est´andar ser´a extendida para soportar variables respuesta que tengan m´as de dos categor´ıas ordenadas. Cuando las categor´ıas de la variable respuesta tengan un orden natural la regresi´on log´ıstica ordinal ser´a la elegida como forma m´as ´optima de estudiar estos datos. Se estudiar´a la forma matem´atica general del modelo de regresi´on log´ıstica ordinal, as´ı como se desarrollar´a su interpretaci´on, las f´ormulas para la odds ratio, intervalos de confianza, t´ecnicas para tests de hip´otesis y para el estudio de la significaci´on de la variable objetivo. Para el desarrollo de este cap´ıtulo se ha utilizado diversa bibliograf´ıa recogida al final de la memoria, especialmente Harrell [9], Kleinbaum y otros [12] y Kleinbaum y Klein[11] 2.1. Modelo de Odds Proporcionales El modelo log´ıstico ordinal que se va a desarrollar es el llamado de Odds Proporcionales omodelo de odds proporcionales, tambi´en conocido como el Modelo Logit Acumulado. Para ilustrar la idea del modelo odds proporcional asumimos que tenemos una variable respuesta con cinco categor´ıas y consideramos las cuatro posibles formas de dividirlas en s´olo dos categor´ıas respetando el orden natural. Por ejemplo, todas estas divisiones ser´ıan posibles respetando el orden: 0 1 2 3 4 01 234 0 1 2 3 4 0 1 2 3 4 13 14 Modelos de regresi´on ordinal Generalmente, si una variable respuesta ordinal Dtiene Gcategor´ıas (D= 0,1,2, ..., G− 1), entonces hay G−1 formas de dicotomizar la respuesta: (D≥1 ´o D < 1; D≥2 ´o D < 2, ..., D ≥G−1 ´o D < G −1). Para un suceso aleatorio S, se define su “odds” o “ventaja” como la raz´on entre la probabilidad de ocurrencia y la probabilidad de no ocurrencia. Con la categorizaci´on de D, se puede definir la “odds” o “ventaja” de que D≥g dividida por la probabilidad de que D < g, i.e. odds (D≥g) = P(D≥g) P(D < g)donde g = 1,2,3, ..., G −1 El modelo odds proporcional hace una importante suposici´on. Bajo este modelo, el odds ratio que eval´ua el efecto de una variable explicativa para cualquiera de las divisiones o categorizaciones anteriores ser´a el mismo independientemente de donde se realice el punto de corte sobre las categorias. Suponemos que tenemos una variable respuesta con cinco niveles y una variable explicativa dicot´omica (E= 0, E = 1). Entonces, bajo la suposici´on de odds proporcionales, el odds ratio que compara categor´ıas iguales o mayores que 1 y categor´ıas menores que 1 es el mismo que el que compara categor´ıas mayores o iguales a 4 con categor´ıas menores que 4. Formalmente: OR(D≥1) = odds [(D≥1)|E= 1] odds [(D≥1)|E= 0] =odds [(D≥4)|E= 1] odds [(D≥4)|E= 0] =OR(D≥4) En otras palabras, el odds ratio es invariante al punto utilizado para la dicotomizaci´on. Esto implica que si hay Gcategor´ıas en la respuesta, solo hay un par´ametro (β) para cada una de las variables predictoras o explicativas. Sin embargo sigue habiendo constantes separadas (αg) para cada una de las G−1 comparaciones. Esto contrasta con la regresi´on log´ıstica polit´omica1, donde hay G−1 par´ametros para cada variable predictora, as´ı como constantes separadas para cada una de las G−1 comparaciones. En resumen: Variable Par´ametro Constante α1, α2, ..., αG−1 X1β1 Cuadro 2.1: Ordinal Variable Par´ametro Constante α1, α2, ..., αG−1 X1β11, β21, ..., βG−1 Cuadro 2.2: Polit´omica La hip´otesis de invarianza del odds ratio en cuanto a los puntos de corte no es la misma que suponer que el odds dado para un patr´on de exposici´on es invariante. Usando 1No solo dicot´omica, multi-categ´orica Modelos de regresi´on ordinal 15 el ejemplo anterior, para una realizaci´on dada de E(e.j. , E = 0), el odds que compara categor´ıas mayores o iguales a 1 con las menores no es igual al odds que compara categor´ıas mayores o iguales a 4 con las menores. odds (D≥1) 6=odds (D≥4) donde, para E= 0 odds(D≥1) =P(D≥1)|E= 0) P(D≥1)|E= 0) 6=P(D≥4)|E= 0) P(D≥4)|E= 0) =odds(D≥4) pero OR (D≥1) = OR (D≥4) 2.1.1. Presentaci´on del modelo Procedemos ahora a presentar la forma del modelo odds proporcional con una respuesta Dde Gniveles (D= 0,1,2, ..., G −1) y una variable explicativa X1.El modelo expresa la probabilidad de que la variable respuesta est´e es una categor´ıa igual o superior agen funci´on de la variable explicativa X1como sigue: P(D≥g|X1) = 1 1 + exp[−(αg+β1X1)], g = 1,2, ..., G −1 Por tanto, la probabilidad de que la variable respuesta est´e en una categor´ıa inferior ages: P(D < g |X1) = exp[−(αg+β1X1)] 1 + exp[−(αg+β1X1)] El modelo puede ser definido equivalentemente en t´erminos del odds de una desigualdad. Si sustituimos la formula P(D≥ |X1) por la expresi´on para el odds entonces: odds(D≥g|X1) = P(D≥g|X1) 1−P(D≥g|X1)=P(D≥g|X1) P(D < g |X1)= =exp(αg+β1X1) = eαg·eβ1X1 Cap´ıtulo 3 Modelado en R En este cap´ıtulo se aborda la aplicaci´on del modelo objeto del presente trabajo desde dos aspectos distintos: La aplicaci´on de los m´etodos de inferencia a trav´es de la librer´ıa “ordinal” de R. Una ilustraci´on de su aplicaci´on e interpretaci´on de resultados. 3.1. Librer´ıa “ordinal” Para el modelaje se usar´a la librer´ıa “ordinal” [4], su objetivo consiste en la implementaci´on del modelo de odds proporcionales [2.1] y otros modelos ordinales. Las funciones m´ass relevantes para la regresi´on ordinal son: clm Ajusta mediante modelos acumulados como el modelo de odds proporcional. El modelo permite varias funciones de enlace y umbrales estructurados que restringen los puntos de corte equidistante o simetricamente dispuestos alrededor de los umbrales centrales. Se usa una modificaci´on del algoritmo de Newton para optimizar la funci´on de m´axima verosimilitud. Se pueden a˜nadir estructuras determinadas para las distribuci´on de los puntos de corte. A trav´es de esta funci´on podemos obtener los vectores de los coeficientes de regre- si´on, los vectores de las constantes α, las probabilidades ya ajustadas, etc. anova.clm Comparaci´on de modelos acumulados a trav´es de contrastes de raz´on de verosimilitudes. 23 24 Modelado en R confint Calcula intervalos de confianza a partir de la funci´on de m´axima verosimilitud de uno o m´as par´ametros. Obtenemos una matriz donde sus columnas proporcionan los intervalos para cada par´ametro. Tambi´en tiene la opci´on de crear la gr´afica en el caso de que queramos estudiarla en busca de problemas de linealidad. convergence Verifica la precisi´on de las estimaciones de los par´ametros de los modelos acumulados. El n´umero correcto decimales y n´umero de d´ıgitos significativos se da para las estimaciones de m´axima verosimilitud de los par´ametros en un modelo de enlace acumulado creado con la funci´on clm. Se obtiene informaci´on sobre la convergencia, errores de estimaci´on de los par´ametros estimados. predict.clm Se obtienen los valores esperados para un modelo previamente creado con la orden clm. Requiere del modelo “clm” y de un “data frame” donde buscar las variables con las que predecir. Devuelve una lista con las predicciones o valores ajustados en el caso de que no se le introduzca un “data frame” con las variables explicativas sobre las que apoyarse. Ofrece la opci´on de que se muestren los intervalos de confianza y los errores. Modelado en R 25 3.2. Ejemplo en R 3.2.1. Datos wine Vamos a considerar los datos de acidez de ciertos vinos y los procesos seguidos en Thompson [18], Fahrmeir y otros [7], Christensen [3] y basandonos en el ejemplo expuesto en Christensen [5] se ilustrar´a en R el modelaje. Los datos de Randall (1989), disponibles en el paquete “ordinal” de R, est´an presentados en el cuadro 3.1, disponible como el conjunto de datos wine en el paquete ordinal. library("ordinal") data(wine) head(wine) ## response rating temp contact bottle judge ## 1 36 2 cold no 1 1 ## 2 48 3 cold no 2 1 ## 3 47 3 cold yes 3 1 ## 4 67 4 cold yes 4 1 ## 5 77 4 warm no 5 1 ## 6 60 4 warm no 6 1 str(wine) ## 'data.frame': 72 obs. of 6 variables: ## $ response: num 36 48 47 67 77 60 83 90 17 22 ... ## $ rating : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 3 4 4 4 5 5 1 2 ... ## $ temp : Factor w/ 2 levels "cold","warm": 1 1 1 1 2 2 2 2 1 1 ... ## $ contact : Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 2 2 1 1 ... ## $ bottle : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 2 ... ## $ judge : Factor w/ 9 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ... El conjunto de datos representa un experimento sobre ciertos factores que determinan la amargura del vino; donde 1 = “lo menos amargo” y 5 = “lo m´as amargo”, dos factores de tratamiento (temperatura y contacto) cada una con dos niveles, la temperatura y el contacto entre el zumo y las pieles de las uvas cuando se extrae de ellas. Nueve jueces evaluaron cada vino de dos botellas para cada una de las cuatro condiciones de tratamiento, por lo tanto, hay 72 observaciones en total. La variable objetivo Ya estudiar ser´a la variable “rating” ∈ {1, ..., 5}que es una categorizaci´on de la variable “response” la cual califica la acidez de los vinos. Vamos a ajustar el siguiente modelo acumulado para los datos wine: logit(P(Yi≤j)) = θj−β1(temperaturai)−β2(contactoi) i= 1, ..., n j = 1, ..., J −1(3.1) 26 Modelado en R Juez Temperatura Contacto Botella 1 2 3 4 5 6 7 8 9 fr´ıo no 1 2 1 2 3 2 3 1 2 1 fr´ıo no 2 3 2 3 2 3 2 1 2 2 fr´ıo s´ı 3 3 1 3 3 4 3 2 2 3 fr´ıo s´ı 4 4 3 2 2 3 2 2 3 2 templado no 5 4 2 5 3 3 2 2 3 3 templado no 6 4 3 5 2 3 4 3 3 2 templado s´ı 7 5 5 4 5 3 5 2 3 4 templado s´ı 8 5 4 4 3 3 4 3 4 4 Cuadro 3.1: Calificaciones de la amargura de algunos vinos blancos. Los datos han sido tomados de Randall (1989). Este es un modelo para la probabilidad acumulada de que la calificaci´on i-´esima caiga sobre la categor´ıa j-´esima o superior, donde iindica cada observaci´on (n= 72) y los ´ındices j= 1, ..., J reflejan la categor´ıa respuesta (J= 5). El par´ametro θjes el punto de corte para el j-´esimo modelo acumulado, logit(P(Yi≤j)). Este modelo es el modelo de odds proporcional descrito en la secci´on 2.1 de este trabajo. Con el comando clm modelizamos a partir del modelo de odds proporcionales previamente desarrollado, ajustado mediante el m´etodo de m´axima verosmilitud fm1<-clm(rating ~temp +contact, data=wine) fm1 ## formula: rating ~ temp + contact ## data: wine ## ## link threshold nobs logLik AIC niter max.grad cond.H ## logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 ## ## Coefficients: ## tempwarm contactyes ## 2.503 1.528 ## ## Threshold coefficients: ## 1|2 2|3 3|4 4|5 ## -1.344 1.251 3.467 5.006 Podemos obtener informaci´on adicional mediante summary Modelado en R 27 summary(fm1) ## formula: rating ~ temp + contact ## data: wine ## ## link threshold nobs logLik AIC niter max.grad cond.H ## logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## tempwarm 2.5031 0.5287 4.735 2.19e-06 *** ## contactyes 1.5278 0.4766 3.205 0.00135 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Threshold coefficients: ## Estimate Std. Error z value ## 1|2 -1.3444 0.5171 -2.600 ## 2|3 1.2508 0.4379 2.857 ## 3|4 3.4669 0.5978 5.800 ## 4|5 5.0064 0.7309 6.850 El primer resultado es la tabla de coeficientes con estimaciones de par´ametros, errores est´andar y los p-valores basados en el m´etodo de Wald. Las estimaciones mediante el m´etodo de m´axima verosmilitud para los par´ametros son: ˆ β1= 2,50,ˆ β2= 1,53,{ˆ θj}={−1,34,1,25,3,47,5,01} El n´umero de iteraciones Newton-Raphson se da a continuaci´on de niter. Tenemos max.grad que es el gradiente absoluto m´aximo de la funci´on de log-verosimilitud con respecto a los par´ametros. Un gradiente absoluto peque˜no es una condici´on necesaria para la convergencia del modelo. El procedimiento iterativo indicar´a convergencia siempre que el gradiente absoluto m´aximo est´e por debajo de clm.control()$gradTol ## [1] 1e-06 Los coeficientes para la temperatura y el contacto son positivos, lo que indica que una temperatura m´as alta y m´as contacto aumenta la amargura del vino, es decir, la calificaci´on en categor´ıas superiores es m´as probable. La odds ratio del suceso Y≥jes exp(βtratamiento), por lo que la odds ratio de acidez que clasifica en la categor´ıa jo superior a temperaturas templadas frente a las fr´ıas es 28 Modelado en R exp(coef(fm1)[5]) ## tempwarm ## 12.22034 Los p-valores para los coeficientes de ubicaci´on dados por el summary se basan en el estad´ıstico de Wald. Las pruebas de raz´on de verosimilitud proporcionan pruebas m´as precisas. Estas se puede obtener con el m´etodo anova, por ejemplo, la prueba de ratios de m´axima verosimilitud del contacto es fm2 <- clm(rating ~temp, data=wine) anova(fm2,fm1) ## Likelihood ratio tests of cumulative link models: ## ## formula: link: threshold: ## fm2 rating ~ temp logit flexible ## fm1 rating ~ temp + contact logit flexible ## ## no.par AIC logLik LR.stat df Pr(>Chisq) ## fm2 5 194.03 -92.013 ## fm1 6 184.98 -86.492 11.043 1 0.0008902 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 que en este caso produce un p-valor ligeramente menor. De manera equivalente, podemos usar Drop1 para obtener pruebas de raz´on de verosimilitud de las variables explicativas mientras se controlan el resto variables: drop1(fm1,test="Chi") ## Single term deletions ## ## Model: ## rating ~ temp + contact ## Df AIC LRT Pr(>Chi) ## <none> 184.98 ## temp 1 209.91 26.928 2.112e-07 *** ## contact 1 194.03 11.043 0.0008902 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Las pruebas de raz´on de verosimilitud de las variables explicativas ignorando el resto vienen dadas por el m´etodo add1: Modelado en R 29 fm0 <- clm(rating ~1,data=wine) add1(fm0, scope =~temp +contact, test ="Chi") ## Single term additions ## ## Model: ## rating ~ 1 ## Df AIC LRT Pr(>Chi) ## <none> 215.44 ## temp 1 194.03 23.4113 1.308e-06 *** ## contact 1 209.91 7.5263 0.00608 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 En este caso, estas ´ultimas pruebas no son tan fuertes como las pruebas que controlan la otra variable. Los intervalos de confianza se obtienen a partir del comando confint como sigue confint(fm1) ## 2.5 % 97.5 % ## tempwarm 1.5097627 3.595225 ## contactyes 0.6157925 2.492404 Estos se basan en los perfiles de la funcion de verosimilitud y generalmente son bastante precisos. Menos preciso, pero m´as simple y con intervalos de confianza sim´etricos basados en los errores est´andar de los par´ametros (tambi´en llamados, intervalos de confianza de Wald) se obtienen con: confint(fm1, type="Wald") ## 2.5 % 97.5 % ## 1|2 -2.3578848 -0.330882 ## 2|3 0.3925794 2.109038 ## 3|4 2.2952980 4.638476 ## 4|5 3.5738541 6.438954 ## tempwarm 1.4669081 3.539296 ## contactyes 0.5936345 2.461961 Adem´as del enlace logit, que es el predeterminado, el probit, el log-log, el log-log complementario y el de cauchit tambi´en se pusen usar para el ajuste. Por ejemplo, un modelo de odds proporcionales como el anterior pero con el enlace log-log ser´ıa: 30 Modelado en R fm.cll <- clm(rating ~contact +temp, data=wine, link="cloglog") En resumen, el primer modelo que creamos supone que los {βj}son constantes para todos los valores de las variables explicativas restantes, en nuestro ejemplo temperatura ycontacto. Esto se conoce como la suposici´on de odds proporcionales osuposici´on de pendientes iguales. 3.2.2. Predicciones sobre wine Los valores ajustados se obtienen con la orden fitted(fm1) y producen probabilidades ajustadas, es decir, la i-´esima probabilidad ajustada ser´ıa la probabilidad de que la i- ´esima observaci´on caiga sobre la categor´ıa prevista. Las predicciones sobre qu´e categor´ıa respuesta tiene mas probabilidad sobre las que caer la i-´esima observaci´on son: pfm1=predict(fm1,type="class") pfm1$fit[1:15] ## [1]223333442233334 ## Levels: 1 2 3 4 5 Digamos que solo queremos las predicciones para las cuatro posibles combinaciones entre temperatura y contacto, entonces ser´ıa: combinaciones <- expand.grid(temp=levels(wine$temp), contact=levels(wine$contact)) cbind(combinaciones, predict(fm1, newdata=combinaciones)$fit) ## temp contact 1 2 3 4 5 ## 1 cold no 0.206790132 0.57064970 0.1922909 0.02361882 0.00665041 ## 2 warm no 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701 ## 3 cold yes 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711 ## 4 warm yes 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469 Los errores est´andar y los intervalos de confianza para las predicciones tambi´en se pueden calcular. Por ejemplo, para las primeras cuatro observaciones; las predicciones, los errores est´andar y los intervalos de confianza al 95 % ser´ıan: h=head(do.call("cbind",predict(fm1, se.fit=TRUE,interval=TRUE))) h[1:4,] ## fit se.fit lwr upr ## [1,] 0.57064970 0.08683884 0.39887109 0.7269447 ## [2,] 0.19229094 0.06388672 0.09609419 0.3477399 ## [3,] 0.44305990 0.07939754 0.29746543 0.5991420 ## [4,] 0.09582084 0.04257593 0.03887676 0.2173139 Cap´ıtulo 4 Ejemplos Reales En esta secci´on se nombrar´an y describir´an brevemente algunas situaciones reales que fueron estudiadas con modelos ordinales. 4.1. Regresi´on Log´ıstica Ordinal Aplicada a la Identificaci´on de Factores de Riesgo para C´ancer de Cuello Uterino Este estudio fue realizado por Evaristo Navarro,An´ıbal Verbel,Delia Robles yKennedy Hurtado, en Barranquilla, Colombia y publicado el 25 de Agosto de 2014. (Navarro y otros [16]) Seg´un los autores, la identificaci´on de factores de riesgo para c´ancer de cuello uterino es determinante a la hora de establecer diagn´osticos efectivos que, en un momento dado, pueden ser determinantes para salvar vidas. Desde esta perspectiva se realiz´o este estudio sobre una muestra constituida por 105 pacientes. En el estudio fue considerada como variable objetivo el C´ancer de cuello uterino (CCU) y como variables explicativas los factores relacionados con la paridad (Edad (ED), N´umero de Hijos Nacidos Vivos (NHV), N´umero de Hijos Nacidos Muertos (NHM), tipo de parto (TP) y tipo de embarazo (TE)). Tambi´en se incluyeron las caracter´ısticas de la conducta sexual (Enfermedades ven´ereas (EV)). De manera general se observa que el riesgo de tener c´ancer de cuello uterino es mayor cuando aumenta el n´umero de hijos en partos por ces´area y se ha perdido un hijo. 4.2. Predicci´on del rendimiento en una asignatura empleando la regresi´on log´ıstica ordinal Este estudio fue realizado por Jobany J. Heredia,Aida G. Rodr´ıguez yJos´e A. Vilalta, llevado por el departamento de Ingenier´ıa Industrial, Facultad de Ingenier´ıa Industrial del Instituto Superior Polit´ecnico “Jos´e Antonio Echeverr´ıa” en La Habana, Cuba. A˜no de publicaci´on, 2014. (Heredia y otros [10]) 31