Los modelos que hemos visto hasta ahora suponian que la variable respuesta fuese normal, o binomial en el caso de la regresión logítica.En este caso vamos a ir un paso más alla y estudiar que pasa cuando las distribuciones no son ni la normal ni la binomial.
¿Qúe vamos a estudiar?
En el último caso la variable no es ni normal ni binomial, sino que es una especie de combinación de ambas.Nos interesa he hecho de que el evento haya sucedido o no y al mismo tiempo queremos saber el tiempo que ha pasado hasta que sucede ese evento.
Me voy a centrar mucho en cómo escribimos el modelo, cómo lo formulamos matemáticamente, cómo lo interpretámos, y en la significación estadísitica tanto de los parámetros, como de las variables en si.
No es lo mismo que un parametro sea significativo, y que la variable asociada a ese modelo sea significativo
En un modelo multinomial, la variable de interés tiene mas de dos categorías.
Cuando teníamos 2 categorías -en el caso de la reresión logística- planteabamos un modelo en base al logit de la probabilidad del evento. \[log(\frac{P(Y=1|X)}{P(Y=0|X)})\] Teniamos dos posibles respuestas, nuestra evento de interes lo denotabamos como Y=1, la alternativa como Y=0. La probabilidad de exito sería \(P(Y=1|X)\) y lo que modelizamos es el logaritmo de esa probabilidad del evento, dividido por la probabilidad del no-evento. Y a todo eso le llamábamos el logit, de ahi el nombre de la regresión logística.
El tema ahora es que nuestra varaible no toma solo los valores 0 y 1. Sino que puede tomar ahsta k valores, siendo \(k>1\). Para visualizarlo, vamos a asumir que nuestra VR toma 3 posibles valores \(k = 0,1, 2\).
Si consideramos que nuestra v.r. toma los valores 0, 1 y 2, lo que tenemos ahora son 2 ecuaciones logit. \[log(\frac{P(Y=1|X)}{P(Y=0|X)}) = \beta_0^1 + \beta_1^1X1 + ... + \beta_q^1X_q )\] \[log(\frac{P(Y=2|X)}{P(Y=0|X)} = \beta_0^2 + \beta_1^2X1 + ... + \beta_q^2X_q)\] Los betas de cada modelo son distintos
Siempre vamos a tener el nuemro de valores -1 de ecuaciones logit.
Siempre vamos a tener que definir cual es la categoría de referencia de nuestra variable respuesta. Lo normal es que cojamos en cero cómo referencia. Y la idea aqui es que yo genero submodelos de regresión logistica donde todos mis posibles valores de k los comparo con el \(k=0\), y para cada uno de esos genero una ecuación.
Las betas de cada modelo son distintos. Yo tengo un modelo donde las covariables que van a explicar la varibale respuesta son hasta q. Y cada una lleva asociado un parametro. Ese parámetro toma valores diferentes en la primera ecuacion y en la segunda ecuacion. En el caso de la regresión multinomial estamos asumiendo que esa relacion que tiene mi variable explcativa con el logit en cada una de las comparaciones que hacemos es distinta. Lo cual nos lleva a la posibilidad de que una variable no explique una comparativa -el cero contra el 1- pero que no me explique el cero contra el 2.
Sea \(Y\) una variable respuest con distribución multinomial que toma valores \(k=1,...,K\), con probabilidades \(p_1,...,p_k\) y \(X = (X_1,...,X_q)\) es el vector de q variables explicativas. Sea \((x_i,y_i)_{i=1}^n\) ( Una muestra de n individuos donde he observado datos para mis q variables explicativas y para mi variable respuesta que toma valores de 1 hasta k )
El siguiente paso es estimar esos parámetros. ¿Cómo lo hacemos? Mediante la método de maxima verosimilitud. ¿Eso que quiere decir? Que tenemos que definir cual es la función de verosimilitud.
Lo que me dice la funcion de verosimilitud lo que me dice es cómo de probable es que en base a ese modelo yo volviese a obtener exactamente esos valores que ya he observado. Cómo de probable es que yo vuelva a observar lo que ya he obseravdo.Yo quiero que esa probabilidad sea máxima.
Lo que tenemos que hacer es que esa funcion de verosimilitud sea maima respecto a nuestros parámetros y para eso derivamos con respecto a todos esos parámetros y generamos una serie de ecuacione que no son lineales y necesitámos métodos iterativos para resolverlos. En la practica lo que se hace es aplicar el logaritmo sobre esa función por que siempre es mucho más facil la funcion maximo del logaritmo que de la función verosimil sin el logaritmo.
En el caso de una variable respuesta multinomial la expresión para el logaritmo de la verosiilitud es la siguiente: \[ l(\beta) = \sum_{i=1}^{n}\sum_{k=1}^{K-1}y_{ik}x_i^T\beta^k-log(1+\sum_{k=1}^{K-1}x_i\beta^k)\] Dónde las \(x\) son datos que hemos observado, las i son datos que hemos observado, y lo unico que no conocemos ahi son los parametros beta. Entonces una vez tengo esa ecuacion, la puedo deribar con respecto a cada uo de esos parámetros beta, igualar a cero, y entocnes obengo los parámetros que me maximizan esa logverosimilitud.
¿Qué significación estadítica tienen esos parámetros?
Es decir, esos parámetros los hemos estimado. Ahora yo quiero saber si ese valor que yo he estimado es significativo. O dicho de otra forma: ¿realmente hay una relación entre esa variable y mi variable respuesta?
Lo mas intuitivo sería plantear un contraste de hipotesis sencillo.Si mi variable explicativa es continua, puedo utilizar un test de wall, donde la hipótesis nula es \(\beta=0\).
Pero si mi variable explicativa es categorica, no me vale con ese test.En un test donde tengo dos betas asociadas a una variable, necesito un test que me haga ese contraste en su conjunto. Y ese test es el test de la razón de verosimilitud.
En este test lo que compara son las deviance del modelo que tiene esa variable y el modelo que no la tiene. La deviance es un parámetro en el que tenemos en cuenta el valor de la verosimilitud de ese modelo. Lo que quiero saber es que si al meter esa variable x2 en el modelo, la regreson que obtengo aumenta la verosimilitud del modelo. Si lo que gana es estadisticamente significativo, me quedo con la varaible. Pero si le meto una variable al modelo para no ganar nada, me quedo sin esa varaible.
Lo que hace el test es comparar los parametros deviance del modelo con esa variable y sin esa varaible. Y lo que estamos asumiendo por detrás es que esa diferencia de deviance -a la que llamamos estadístico G- sigue una distribucion chi cuadrado. Esa función nos permite calcular un p-valor asocaido a ese test.
El odds ratio es el parámetro que vamos a utilizar para interpretar el modelo. Lo que hace es comparar la probabilidad de éxito frente a la probabilidad de fracaso para los diferentes valores que toman las variables explicativas.
\[OR_k(x_0,x_1) = \frac{P(Y=k|X=x_1)/P(Y=0|X=x_1)}{ P(Y=k|X=x_0)/P(Y=0|X=x_0)}, k=1,2.\] Donde \(k\) indica el nivel de la variable respuesta con el que estámos comparándo y \(x_o\) y \(x_1\) son los valores de las variables explicativas. Si solo tenémos una varibale explicativa dicotómoica que toma valores 0 y 1, entonces tendríamos \(OR_k(0,1)\)
De todas las variables que yo tengo, me tengo que quedar con aquellas que mejor me expliquen el modelo. Y no quiero 10 si me valen 5. Queremos modelos que sean lo mejor posible con el menor número de variables posible.
ANALIZAR LAS VARIABLES DE FORMA UNIVARIANTE. De una en una. cada una, si solo metemos esa variable en el modelo, si son significativas o no lo son.
CONSTRUIR UN MODELO MULTIVARIANTE ¿Qué varaibles meto? Claramente todas aquellas que sean significativas en el univariante. Pero es que además, aquellas que en el univariante no sean significativas pero tenga un p-valor por debajo de 0.25 las voy a meter también. Por que puede pasar que una variable por si misma no sea significativa, pero junto con otra si.
COMPARAR MODELOS Tenemos que comparar modelos,tenemos que ver cual es mejor.Uno de los criterios va a ser el AIC.
A PARTIR DE AQUI SEGUIMOS CON EL EJEMLPLO DE REGRESIÓN MULTINOMIAL.
los datos que vamosa utilizar provienen de un estudio cuyo objetivo era estudiar que variables influian en el conocimiento y la aptitud y el comportamiento de las mujeres hacia las mamografías son datos de 412 mujeres
setwd("/home/alex/Escritorio/Máster Matemáticas/MODELIZACIÓN ESTADÍSTICA/IRANZU BARRIO")
mamexp <- read.table("mamexp.txt",header=TRUE)
names(mamexp) # para saber que variables tengo, cuales son los nombres que tengo en mi base de datos, y es una manera sencilla de verificar que hemos cargado bien los datos## [1] "me" "symp" "pb" "hist" "bse" "dect"MOELO MULTINOMIAL (3 CATEGORIAS–> 2 ECUACIONES LOGIT)
lo primero que tenemos que hacer es plantear el modelo y ver si cada una de esas variables son o no son significativas
todas aquellas variables que tomen valores pero que cada uno de esos valores identifique una categoría, R piensa por defecto que es una variable continua que toma valores 0 y 1, no sabe que es una variable categórica. YO le tengo que decir a R que esa variable no es numérica, que es categórica, con la función factor.
# Convertimos variables categóricas en factor
mamexp$symp <- factor(mamexp$symp)
mamexp$hist <- factor(mamexp$hist)
mamexp$bse <- factor(mamexp$bse)
mamexp$dect <- factor(mamexp$dect)
# de esta manera R ya sabe que 0 es una categoría y 1 otra
attach(mamexp)
table(me)## me
## 0 1 2
## 234 104 74Ajustamos modelos multinomiales univariantes( con 2 de las variables con hist y dect) La funcion multinom nos da las funciones logit de las variables que ajustamos, referenciadas al cero.
# Cargamos la librería para usar la función multinom()
library("nnet")
# función de R para ajustar modelos de regresión multinomial:
# la variable hist me dice si hay historial de cancer de mama en la familia 0=No, 1=Si
# eso quiere decir que necesito un único beta.
multi1_hist <- multinom(me~hist)## # weights: 9 (4 variable)
## initial value 452.628263
## iter 10 value 396.169969
## iter 10 value 396.169969
## final value 396.169969
## convergedmulti1_hist## Call:
## multinom(formula = me ~ hist)
##
## Coefficients:
## (Intercept) hist1
## 1 -0.9509794 1.256531
## 2 -1.2504803 1.009384
##
## Residual Deviance: 792.3399
## AIC: 800.3399La variable me tiene 3 valores 0=nunca, 1=último año y 2= más de un año. Por tanto, los logits creados a la salida de la funcion multinom para la variable hist son último año vs nunca (1/0) y más de un año vs nunca (2/0). R elige automaticamente el cero como categoría de referencia.
multi1_dect <- multinom(me~dect)## # weights: 12 (6 variable)
## initial value 452.628263
## iter 10 value 389.200621
## final value 389.200537
## convergedmulti1_dect## Call:
## multinom(formula = me ~ dect)
##
## Coefficients:
## (Intercept) dect1 dect2
## 1 -2.565201 0.7062589 2.1062453
## 2 -1.178617 -0.3926042 0.1977986
##
## Residual Deviance: 778.4011
## AIC: 790.4011# si no tuviese hecho el attach la función sería asi: multinom(me~hist, data = namexp)La variable dect tiene 3 valores (0,1,2). Cómo tengo 3 categoría necesito 2 betas.Pero cómo estoy en un modelo de regresión multinomial, mi v.r. toma 3 valores y por lo tanto tengo 2 ecuaciones, voy a necesitar un total de 4 betas. 2 por el numero de ecuaciones que tenga.
Siguiente pregunta: he estimado, tengo los valores correspondiente para mis betas. La siguiente pregunta es: ¿Son esta disticamente significativas? De verdad esa variable hist o esa variable det en esos modelos univariantes me explican mi varaible respuesta?–> necesitamos un test de la razon de verosimilitud, why? Por que yo tengo 2 ecuaciones, y a la hora de saber si esa variable es estadisticamente significativa, yo lo que quiero saber es si explica en su conjunto toda mi variable respuesta. Luego ya analizare esos betas uno a uno para saber en que casos en concreto esas difrencias son significativas. Puede haber que en unas comparaciones no lo sea, y en otras que si lo sea, pero quiero saber si en su conjunto la variable es significativa.
¿Cómo la hacemos? ANOVA. Hay que tener cuidado en utilizarla bien ya que es una funcion que esta desarrollada de tal forma que sirve para comparar muchos modelos en situaciones diferentes.
Para saber si esa variable es significativa a o no, la tengo que comparar con el modelo nulo, donde no hago nada, donde no le introduzo ninguna variable, me quedo solo con mis probabilidades asignadas originales en base a mis datos y no las intento modelizar en base a otras variables. Por eso tenemos que ajustar el modelo nulo (multi0), donde mi variable respuesta sigue siendo la misma pero no tengo covariables ¿Y como le idgo a r que no tengo covariables? gusanito 1. Despues comparo mi modelo nulo con el modelo que si tiene esa variable ¿y que hago? comparar la deviance, las verosimilitudes de los 2 modelos.
La hipótesis nula siempre va a ser que los betas son cero.En este caso, lo que estoy asumiendo en esa hipotesis nula es que todos los betas son cero. y la alternativa es que yo tengo algo distinto de cero. y en el momento que yo tengo algo distinto de 0 ya consigo una significación estadística.
# Test de la razón de verosimilitud. Significaciónn de variables
multi0 <- multinom(me~1)## # weights: 6 (2 variable)
## initial value 452.628263
## final value 402.599009
## convergedanova(multi1_hist,multi0)# comparo modelo nulo con el modelo que si tiene esa variable## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 822 805.1980 NA NA NA
## 2 hist 820 792.3399 1 vs 2 2 12.85808 0.001614001#Pr(chi)=0.001614001<0.05-->TRUE. Es decir, el p-valor es menor al nivel de significación alpha=0.05. Por tanto
# rechazamos h0, osea, que rechazamos que beta0 sea nula, y por tanto puedo decir que la variable hist es estadisticamente significativa# otra cosa que tenemos que tener en cuenta es la diferencia entre grados de libertad de los dos modelos.- la distribucion shi cuadrado utiliza esos grados de libertas-
anova(multi1_dect,multi0)## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 822 805.1980 NA NA NA
## 2 dect 818 778.4011 1 vs 2 4 26.79694 2.184912e-052.184912e-05<0.05, rechazamos h0, variable Si ES signficativa
# Estimaciónn de OR (ODS RATIO)
exp(coef(multi1_hist))# este comando calcula las ods ratio del coef hist. ## (Intercept) hist1
## 1 0.3863624 3.513213
## 2 0.2863672 2.743911Interpretacion: las mujeres con historia familiar de cancer de pecho tienen un ods de haberse hecho una mamo en el ultimo año que es 3.51 veces mayor a las mujeres sin historia familiar. Recordad que cuando hacemos el odds ratio yo estoy comparando la probabilidad de exito con la probabilidad de fracaso en aquellas que tienen historia familiar y en aquellas que no. Un odds ratio de 1 quiere decir que esos ratios son iguales. Eso querria decir que mi probabilidad de 0 a 1, de pasar de la categoria 0 a 1 en mi variable respuesta no cambiaria entre las muejeres con historia y sin historia.Entonces, todo lo que sea de un odds ratio por encima de 1 quiere decir más riesgo.Más probable que se de esa categoría con la que estámos calculando el logit. En este caso sería “último año”
Las mujeres con historia familiar tienen un odds ratio 2.74 mayor de haberse hecho una mamo hace mas de un año que las que no tienen historial familiar.
exp(coef(multi1_dect))## (Intercept) dect1 dect2
## 1 0.07690374 2.0263961 8.217330
## 2 0.30770391 0.6752959 1.218717exp(coef(multi1_dect))[1,]# aqui me quedo solo con la primera parte, los parametros del logit del ultimo año vs nunca## (Intercept) dect1 dect2
## 0.07690374 2.02639611 8.21732991A mi me interesa interpretar aquellos betas que son distintos de cero.Yo aqui estoy estudiando la probabilidad de hacerse una mamografía el último año vs no hacersela nunca. Y estoy metiendo mi variable dect que tiene 3 categorías en el modelo. Cojo la categoría cero como referencia y por lo tanto tengo: un beta para dect=1 y un beta para dect=2. Si me quedo con esta primera ecuación, el primer beta lo que me esta comparando es: en la situacion en la que comparo ultimo año vs nunca, que efecto tengo cuando paso del valor de dect=0 al valor de dect=1, y el segundo beta me dice que efecto tengo cuando paso del valor dect=0 a dect=2 en esta primera ecuación.
Tengo 4 betas distintos. Y puede pasar que alguno de ellos sea significativo pero otros no.(test de wall) El test de wall es equivalente al intervalo de confianza que podemos calcular para esos parámetros. Si yo calculo un intervalo de confianza al 95% para esos parametros,¿qué tendría que mirar para saber si son estadisticamente significativo?- Que el cero esté en el intervalo.
Si el cero esta en el intervalo significa que ese beta puede ser cero, indicando que no es estadisticamente significativo. Porque cuando el cero no esté en el intervalo de confianza de ese beta, querra decir que ese beta es significativo. Porque yo ya he visto que la variable en su conjunto es significativa. Ahora lo que quiero saber es exactamente en que comparaciones en concreto tengo esas diferencias significativas.
# Estimación intervalo de confianza para los parámetros
IC_dect <- confint(multi1_dect, level=0.95)
IC_dect[,,1] # estimaciones de los parámetros del logit Ultimo año vs Nunca (0/1)## 2.5 % 97.5 %
## (Intercept) -4.5993833 -0.5310183
## dect1 -1.4169276 2.8294454
## dect2 0.0551994 4.1572913Si el cero esta en intervalo de confianza del beta, pero tambien lo puede mirar si el 1 esta en el intervalo de confianza del Odds ratio, por que recordad: que el ods ratio es la exponencial del beta. si beta es cero, el ods ratio es 1
He sacado los intervalos de confianza para los betas de esa primera ecuación. Cómo tengo el cero en el intervalo de confianza de dect1 eso quiere decir que ese beta podría ser cero. Lo cual quiere decir que cuando estoy comparando ultimo año vs nunca, de dect=0 a dect=1 no hay evidencias estadisticamente significativas. Sin embargo si que las hay cuando comparo dect=0 con dect=2.EN ESTE LOGIT.
¿Cómo lo interpreto? el odds de hacerse una mamo dentro del ultimo año entre las mujeres que piensan que es muy probable que la mamo detecte un nuevo caso de cancer(dect=2), es 8.22 veces mayor que el ods entre las mujeres que piensan que no es probale que la mamo lo detecte (dect=1).
Tendríamos que hacer lo mismo si queremos comparar el logit de más de un año vs nuca
IC_dect[,,2] # estimaciones de los parámetros del logit Más de un año vs Nunca## 2.5 % 97.5 %
## (Intercept) -2.2992503 -0.0579843
## dect1 -1.6359055 0.8506970
## dect2 -0.9656586 1.3612557En este caso, los intervalos de confianza de los betas, ambos incluyen al cero. Por lo tanto ninguno de esos 2 betas es estadisticamente significativo. Lo cual quiere decir que de los 4 betas que hemos estimado solo tenemos 1 que es estadisticamente significativo. Pero en su conjunto, obtenemos que esa variable ES estadisticamente significativa.
# En base al test de la razón de la verosimilitud, todos los modelos univariantes son significativos:
multi1_symp <- multinom(me~symp)## # weights: 15 (8 variable)
## initial value 452.628263
## iter 10 value 370.716115
## final value 370.706395
## convergedmulti1_bse <- multinom(me~bse)## # weights: 9 (4 variable)
## initial value 452.628263
## final value 394.032437
## convergedmulti1_pb <- multinom(me~pb)## # weights: 9 (4 variable)
## initial value 452.628263
## final value 384.972360
## convergedanova(multi1_symp,multi0)## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 822 805.1980 NA NA NA
## 2 symp 816 741.4128 1 vs 2 6 63.78523 7.63456e-12anova(multi1_bse,multi0)## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 822 805.1980 NA NA NA
## 2 bse 820 788.0649 1 vs 2 2 17.13314 0.0001903642anova(multi1_pb,multi0)## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 822 805.1980 NA NA NA
## 2 pb 820 769.9447 1 vs 2 2 35.2533 2.212299e-08# después tendríamos que montar el mejor modelo de entre todas las variablesDespués tendríamos que montar el mejor modelo de entre todas las variables.
Entonces vamos a generar un modelo multivariante en el que vamos a meter todas las variables que queramos meter.
Para meter variables, directamente las puedo sumar:
# Ajuste del modelo múltiple con todas las variables explicativas
multi2 <- multinom(me ~ hist+dect+symp+bse+pb)## # weights: 30 (18 variable)
## initial value 452.628263
## iter 10 value 348.174100
## iter 20 value 346.951267
## final value 346.950964
## convergedsummary(multi2)## Call:
## multinom(formula = me ~ hist + dect + symp + bse + pb)
##
## Coefficients:
## (Intercept) hist1 dect1 dect2 symp2 symp3 symp4
## 1 -2.9990607 1.366274 0.01700706 0.9041717 0.1101408 1.9249158 2.457230
## 2 -0.9860258 1.065446 -0.92448094 -0.6906153 -0.2901346 0.8172916 1.132224
## bse1 pb
## 1 1.291772 -0.2194421
## 2 1.052183 -0.1482079
##
## Std. Errors:
## (Intercept) hist1 dect1 dect2 symp2 symp3 symp4
## 1 1.539293 0.4375239 1.1619394 1.1268643 0.9228324 0.7776651 0.775400
## 2 1.111829 0.4593986 0.7137332 0.6871025 0.6440622 0.5397872 0.547666
## bse1 pb
## 1 0.5299080 0.07551477
## 2 0.5149934 0.07636886
##
## Residual Deviance: 693.9019
## AIC: 729.9019La salida ya se nos complica un poco. Sigo teniendo dos filas porque tengo 2 logits, pero se me van acumulando las columnas.
confint(multi2)## , , 1
##
## 2.5 % 97.5 %
## (Intercept) -6.0160192 0.01789780
## hist1 0.5087428 2.22380478
## dect1 -2.2603524 2.29436653
## dect2 -1.3044418 3.11278514
## symp2 -1.6985774 1.91885907
## symp3 0.4007202 3.44911153
## symp4 0.9374738 3.97698584
## bse1 0.2531716 2.33037280
## pb -0.3674483 -0.07143587
##
## , , 2
##
## 2.5 % 97.5 %
## (Intercept) -3.16516983 1.193118231
## hist1 0.16504089 1.965850477
## dect1 -2.32337238 0.474410491
## dect2 -2.03731145 0.656080761
## symp2 -1.55247327 0.972204160
## symp3 -0.24067187 1.875255164
## symp4 0.05881781 2.205629242
## bse1 0.04281458 2.061551754
## pb -0.29788816 0.001472263Ahora yo quiero saber con qué variables me tengo que quedar. Lo que puedo hacer es estudiar los intervalos de confianza para todos esos betas en cada una de las ecuaciones, y me fijo en cuales el cero no está incluido. Pero eso lo que me dice es si esos betas son significativos, pero no estamos hablando de la varible puesto que me falta la segunda parte de la ecuación. O sea que tengo que hacer lo mismo con la siguiente ecuación.
Si hay una variable la cual todas sus betas son cero, va a ser una candidata a no entrar en el modelo. Peroque tenemos que hacer? o estudiar el test de la razon de verosimilitud para todas las posibles combinaciones, o utilizar una función,que se llama step, y lo que hace es analizar diferentes modelos, partiendo de ese modelo grande que nosotros hemos introducido, para ver cual seria mejor. y el criterio que utiliza esa funcion step para ver cual es el mejor, es el valor del AIC.
step(multi2)## Start: AIC=729.9
## me ~ hist + dect + symp + bse + pb
##
## trying - hist
## # weights: 27 (16 variable)
## initial value 452.628263
## iter 10 value 353.888326
## iter 20 value 352.579296
## final value 352.579290
## converged
## trying - dect
## # weights: 24 (14 variable)
## initial value 452.628263
## iter 10 value 351.526518
## final value 351.173900
## converged
## trying - symp
## # weights: 21 (12 variable)
## initial value 452.628263
## iter 10 value 366.301013
## final value 366.161513
## converged
## trying - bse
## # weights: 27 (16 variable)
## initial value 452.628263
## iter 10 value 352.379941
## iter 20 value 351.786439
## final value 351.786432
## converged
## trying - pb
## # weights: 27 (16 variable)
## initial value 452.628263
## iter 10 value 353.100061
## iter 20 value 352.019446
## iter 20 value 352.019443
## iter 20 value 352.019443
## final value 352.019443
## converged
## Df AIC
## <none> 18 729.9019
## - dect 14 730.3478
## - bse 16 735.5729
## - pb 16 736.0389
## - hist 16 737.1586
## - symp 12 756.3230## Call:
## multinom(formula = me ~ hist + dect + symp + bse + pb)
##
## Coefficients:
## (Intercept) hist1 dect1 dect2 symp2 symp3 symp4
## 1 -2.9990607 1.366274 0.01700706 0.9041717 0.1101408 1.9249158 2.457230
## 2 -0.9860258 1.065446 -0.92448094 -0.6906153 -0.2901346 0.8172916 1.132224
## bse1 pb
## 1 1.291772 -0.2194421
## 2 1.052183 -0.1482079
##
## Residual Deviance: 693.9019
## AIC: 729.9019Si no le quito ningun parámetro a mi modelo, tengo un AIC de 729, que resutla ser el más pequeño de todos.
multi3 <- multinom(me ~hist+symp+bse+pb)## # weights: 24 (14 variable)
## initial value 452.628263
## iter 10 value 351.526518
## final value 351.173900
## convergedstep(multi3)## Start: AIC=730.35
## me ~ hist + symp + bse + pb
##
## trying - hist
## # weights: 21 (12 variable)
## initial value 452.628263
## iter 10 value 357.587784
## final value 357.041154
## converged
## trying - symp
## # weights: 15 (8 variable)
## initial value 452.628263
## iter 10 value 373.863782
## final value 373.834205
## converged
## trying - bse
## # weights: 21 (12 variable)
## initial value 452.628263
## iter 10 value 356.025814
## final value 355.756885
## converged
## trying - pb
## # weights: 21 (12 variable)
## initial value 452.628263
## iter 10 value 358.173731
## final value 357.670364
## converged
## Df AIC
## <none> 14 730.3478
## - bse 12 735.5138
## - hist 12 738.0823
## - pb 12 739.3407
## - symp 8 763.6684## Call:
## multinom(formula = me ~ hist + symp + bse + pb)
##
## Coefficients:
## (Intercept) hist1 symp2 symp3 symp4 bse1 pb
## 1 -2.181670 1.373809 0.1530116 2.0855801 2.618928 1.2666132 -0.2502638
## 2 -1.733803 1.115519 -0.3074208 0.8418555 1.142751 0.9816765 -0.1380179
##
## Residual Deviance: 702.3478
## AIC: 730.3478Pero si le quito la variable dect, que en modelo univariante si que era significativa (aunque solo una de sus betas lo era), el AIC no se reduce pero es muy parecido.
¿En este caso que hacemos?
Pues ajusto el modelo con todas las variables, y ajusto mi modelo con todas las variables menos dect.Y los comparo. Y cuando los comparo lo que estoy haciendo es el test de la razon de verosimilitud.
# Test de la razón de verosimilitud (comparación de modelos)
anova(multi2,multi3)## Model Resid. df Resid. Dev Test Df LR stat.
## 1 hist + symp + bse + pb 810 702.3478 NA NA
## 2 hist + dect + symp + bse + pb 806 693.9019 1 vs 2 4 8.445871
## Pr(Chi)
## 1 NA
## 2 0.07654504Y obtengo un \(p-valor = 0.07654504\). ¿Meto o no meto dect?
\(H_0: vector de betas =0\) por tanto cuando comparo dos modelos, mi hipotesis nula siempre va a ser el modelo pequeño. Alternativa: $H_1: $ modelo grade, vector de betas distinto de cero.
cómo el p valor del anova sale superior al nivel de significacion de 0.05, aceptamos h0 Si el p-valor es mayor al nivel de significación y acepto H0, quiere decir que el vector de betas es igual a cero, por tanto no entra la variable dect en el modelo al no ser significativa.