Modelos Estadísticos. Grado Biotecnología

Introducción

Un tabla de contingencia se usa para estudiar una tabla de clasificación de un conjunto de factores donde se recogen el número de casos de cada una de las combinaciones de dichos factores.

En este tipo de modelos las variables respuestas y predictoras son todas categóricas, medidas en una escala nominal u ordinal, aunque como veremos las construcción de este tipo de modelos se centran en estudiar las frecuencias de ocurrencia de cada combinación de los factores. El objetivo básico en el análisis de tablas de contingencia es estudiar si existe alguna relación/asociación entre los factores de clasificación considerados. posteriormente, si dicha relación existe, habrá que describir cómo es para tratar de predecirla. Para cuantificar el grado de asociación entre los factores de clasificación, los modelos log-lineales proporcionarán predicciones de las frecuencias observadas en las celdas y las probabilidades asociadas a la clasificación, en función de dichos factores y la interacción entre ellos.

Ejemplo Melanoma. Los datos siguientes provienen de un estudio de pacientes con una froma de cáncer de piel llamado melanoma maligno. En una muestra de 400 pacientes se recogió información sobre la localización del rumor y su tipo histológico. Los datos son el número de pacientes (fecuency) en cada combinación de tipo de tumor (type) y localización (site). El objelivo básico del análisis es investigar si determinados tipos de tumor están asociados a ciertas localizaciones. Realiza un análisis prelimianr de la información recogida y ajusta el modelo que correspnda al objetivo del estudio. El objelivo básico del análisis es investigar si determinados tipos de tumor están asociados a ciertas localizaciones. Ambas variables, tipo de tumor y localización, tienen pues, el mismo rol en el análisis, y sólo nos interesa determinar si están asociadas o no.

Tipo / Localización	Cabeza-Cuello	Tronco	Extremidades
Melanoma de Hutchinson	22	2	10
Melanoma Superficial	16	54	115
Nodular	19	33	73
Indeterminado	11	17	28

Ejemplo Vacuna de la gripe. En un estudio prospectivo sobre una nueva vacuna para la gripe, los pacientes fueron asignados aleatoriamente a dos grupos. A los pacientes de uno de los grupos se les trató con la nueva vacuna y a los otros se les dió un placebo salino. Las respuestas fueron los niveles de anticuerpos inhibidores de hemoglutinina (HIA) encontrados en la sangre seis semanas después de la vacunación. Los datos se encuentran en el banco de datos siguiente. El objelivo del estudio es investigar el efecto de la nueva vacuna, esto es, comprobar si el hecho de dar placebo o vacuna provoca diferente respuesta HIA. Así, la variable HIA es la variable a explicar en función del tipo de tratamiento que ha recibido el paciente y las frecuencias de las celdas.

Tratamiento / Respuesta	Pequeño	Moderado	Grande
Placebo	25	8	5
Vacuna	6	18	11

Ejemplo Aspirina. En un estudio retrospectivo de casos-controles, se reunió a un gupo de pacientes con úlcera y se buscó a otro grupo de individuos de similares características a los anteriores en edad, sexo y estatus socio económico, sobre los que no se sabía que tuvieran úlcera péptica. Los pacientes con úlcera fueron clasificados de acuerdo a la localización de la úlcera: gástrica o duodenal. Se les preguntó a todos si consumían apirina. El objetivo del análisis es averiguar si existe alguna relación entre el consumo de aspirina y la existencia de algún tipo de úlcera. Si existe tal asociación, querríamos comprobar si el consumo de aspirina es más hahitual en los pacientes con úlcera gástrica que en los de úlcera duodenal. Interesa predecir el consumo de aspirina (variable a explicar) en función de si un individuo tiene úlcera o no. y si la tiene, de qué tipo es.

Tipo Úlcera	Úlcera / Consumo de aspirina	No	Si
Gástrica	Si	39	25
	No	62	6
Duodenal	Si	49	8
	No	53	8

Ejemplo Aspiraciones. En este banco de datos se muestran una parte de los resultados de un estudio en los Estados Unidos con el que se pretendía investigar el grado de asociación entre las aspiraciones de los estudiantes de bachillerato por proseguir con estudios universitarios y su entorno social, medido en términos del estatus socio-económico de su familia y el hecho de que recibieran apoyo en su familia para continuar estudiando. El objetivo es predecir la proporción de individuos que pretenden continuar en la universidad en función del estatus económico de la familia y el apoyo que reciben en ella.

Estatus	Motivación familiar / Aspiraciones Universitarias	No	Si
Baja	Baja	749	35
	Alta	233	133
Media-Baja	Baja	627	38
	Alta	330	303
Media-Alta	Baja	420	37
	Alta	374	467
Alta	Baja	153	26
	Alta	266	800

Ejemplo Contraceptivos. En la tabla aparecen los datos de un estudio sobre distintos usos contraceptivos en diferentes grupos generacionales. Entendiendo la variable método contraceptivo como variable respuesta y la variable edad como predictiva, una cuestión de interés puede ser predecir el uso anticonceptivo de una mujer en función de su edad. Esto lo podemos hacer comparando la proporción de mujeres que usa cada uno de los usos anticonceptivos con la de aquellas que no usan ninguno frente a los que si usan o son estériles para los diferentes grupos de edad.

Edad / Método contraceptivo	Esterilización	Otros	Ninguno
15-19	3	61	232
20-24	80	137	400
25-29	216	131	301
30-34	268	76	203
35-39	197	50	188
40-44	150	24	164
45-49	91	10	183

Ejemplo Satisfacción laboral. Los datos que se presentan corresponden a un estudio en el que se pretendía concluir sobre la relación entre el grado de satisfacción en el trabajo y los ingresos percibidos. El objetivo es predecir el grado de satisfacción en el trabajo en función de los ingresos percibidos (en dólares).

Ingresos / Satisfacción	Muy Insatisfecho	Poco Insatisfecho	Moderadamente Satisfecho	Muy Satisfecho
< 6000	20	24	80	82
6000 - 15000	22	38	104	125
15000 - 25000	13	28	81	113
> 25000	7	18	54	92

Modelos Log-lineales

Las cuestiones básicas de interés en el análisis de tablas de contingencia son detectar asociación entre los factores de clasificación y predecir las frecuencias medias observadas. Veamos cómo se pueden resolver estas cuestiones a partir de la modelización log-lineal de las frecuencias esperadas. Consideramos, para presentar estos modelos, únicamente dos vías (factores) de clasificación, F1 y F2. En esta situación el número de individuos que esperamos clasificar en la celda (i, j) se calcula a través del producto de los sucesos independientes “número esperado de individuos clasificados en el nivel i de F1”, \(\mu_{i+}\), y “número esperado de individuos clasificados en el nivel j de F2”, \(\mu_{+j}\) esto es,

\[E(Y_{ij}) = \mu_{ij} = \mu_{i+}*\mu_{+j}\] El efecto multiplicativo de los factores para predecir la frecuencia observada da lugar, de modo natural, a modelos logarítmicos aditivos (modelos log-lineales) de la forma:

\[log(E(Y_{ij})) = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij}\] donde \(\alpha_i\) representa el efecto principal del factor fila F1 , \(\beta_j\) el correspondiente al factor columna F2 y el \((\alpha\beta)_{ij}\) efecto de interacción entre ambos. Este modelo se denomina modelo saturado, ya que tiene tantos parámetros como celdas estamos considerando en la combinación de los factores, y tiene poco interés desde el punto de vista estadístico ya que cada parámetro del modelo se estima con un único valor en la muestra de datos. Contrastar independencia versus asociación será equivalente a contrastar interacción nula versus no nula en el modelo log-lineal, es decir debemos estudiar los modelos:

El modelo saturado dado por \(Frecuencias \sim F1 * F2\)
El modelo sin interacción dado por \(Frecuencias \sim F1 + F2\)

dados por las expresiones

\[\begin{array}{ll} M_0: log(E(Y_{ij})) &= \mu + \alpha_i + \beta_j\\ M_1: log(E(Y_{ij})) &= \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij}\\ \end{array}\]

En la práctica para contrastar independencia ajustamos un modelo con únicamente los efectos principales y valoramos la bondad del ajuste conseguido. En conclusión, reconocer asociación entre los factores de clasificación de una tabla de contingencia es equivalente a rechazar interacción nula entre ellos en el modelo log-lineal con el que se predicen las frecuencias observadas en las celdas.

Cuando consideramos la tabla de contingencia asociada con tres factores el modelo saturado vendrá dado por:

\[Frecuencias \sim F1 + F2 + F3 + F1:F2 + F1:F3 + F2:F3 + F1:F2:F3\] donde el efecto de interacción triple \(F1:F2:F3\) nos permite estudiar si existe asociación entre los tres factores conjuntamente. En caso de que el contraste de bondad de ajuste del modelo reducido \[M_0: Frecuencias \sim F1 + F2 + F3 + F1:F2 + F1:F3 + F2:F3\] no resulte significativo, podremos abordar el estudio de las interacciones dobles de la siguiente forma:

la interacción \(F1:F2\) nos permite estudiar la asociación entre ambos factores dado F3
la interacción \(F1:F3\) nos permite estudiar la asociación entre ambos factores dado F2
la interacción \(F2:F3\) nos permite estudiar la asociación entre ambos factores dado F1

Por tanto, debemos estudiar la bondad de ajuste de los modelos \[\begin{array}{ll} M_1:& Frecuencias \sim F1 + F2 + F3 + F1:F3 + F2:F3\\ M_2:& Frecuencias \sim F1 + F2 + F3 + F1:F2 + F2:F3\\ M_3:& Frecuencias \sim F1 + F2 + F3 + F1:F2 + F1:F3\\ \end{array}\]

para detectar las posibles asociaciones. Si el factor F1 Actúa como variable respuesta y los otros dos factores como predictoras, le modelo M3 no tiene sentido práctico, ya que en realidad nuestro interés radica en la asociación entre respuesta y predictoras.

Así podemos seguir aumenta el número de factores en el modelo y estudiando las posibles asociaciones entre ellos. Se trata pues de un procedimiento de comparación de efectos o modelos de tipo manual y no automático como habíamos hecho en temas anteriores.

Para escribir de forma reducida lo modelos anteriores usamos la siguiente notación (que también nos vale para el ajuste de los modelos en R). \[\begin{array}{ll} M_0:& Frecuencias \sim (F1+F2+F3)^2\\ M_1:& Frecuencias \sim F2*F3 + F1 + F1:F3\\ M_2:& Frecuencias \sim F2*F3 + F1 + F1:F2\\ M_3:& Frecuencias \sim F1*F3 + F2 + F1:F2\\ \end{array}\]

A continuación se presenta el análisis de algunos de los ejemplos, ajustando los modelos de asociación-independencia que correspondan según la pregunta u objetivo de interés en cada caso. Antes de comenzar cargamos las librerías de trabajo.

Melanoma

¿Existe asociación entre el tipo de tumor y su localización?

Esta pregunta equivale a estudiar si el efecto de interacción entre tipo de tumor y su localización resulta significativo, es decir, estudiar la bondad de ajuste del modelo sin interacción entre ambos. Cargamos los datos y representamos la información proporcionada:

# Leemos los datos
glm_multi_01 = read_csv("https://goo.gl/yeCXsv", col_types = "cci")
# Gráfico
ggplot(glm_multi_01, aes(x = site, y=frequency, fill = type)) + 
  stat_identity(geom = "bar", position = "dodge") +
  theme_bw()

Podemos ver claramente que hay ciertos tumores que tienen una mayor predisposición a aparecer en una parte del cuerpo que en otra. Pasamos a ajustar el modelo sin interacción:

# Ajuste del modelo sin interacción
modelo <- glm(frequency ~ type + site, family = poisson, data = glm_multi_01)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 2.050453e-09

Dado que el valor es inferior a 0.05 podemos concluir que hay evidencias estadísticas a favor de la asociación entre los factores tipo de tumor y localización, es decir, determinados tipos de tumor se presentan en ciertas localizaciones con más frecuencia que en otras.

Aspirina

En el estudio se controlaron las variables úlcera y localización. A los pacientes clasificados en cada una de las combinaciones de estas dos variables se les preguntó sobre el consumo de aspirina, que es la variable aleatoria y por tanto la variable a explicar en función de las otras dos. Las cuestiones que nos planteamos son:

En el grupo gástrico ¿la úlcera está asociada con el consumo de aspirina? ¿y en el grupo duodenal?. Para ambas cuestiones estamos planteando si, dado un grupo especifico (locali), el consumo de aspirina (aspirina) está relacionado con la existencia de úlcera (ulcera). Ambas preguntas se pueden contestar estudiando el efecto de interacción entre el consumo de aspirina y el desarrollo de una úlcera (aspirina:ulcera).
Si dicha asociación existe, cabe preguntarse si dicha relación es similar en los dos grupos de pacientes

# Cargamos los datos
locali <- gl(2,2,8,labels = c("Gástrica","Duodenal"))
ulcera <- gl(2,1,8,labels = c("Si","No"))
aspirina <- gl(2,4,8,labels = c("No","Si"))
frecuencia <- c(39,62,49,53,25,6,8,8)
glm_multi_04 <- data.frame(locali,ulcera,aspirina,frecuencia)
# Gráfico
ggplot(glm_multi_04, aes(x = aspirina, y = frecuencia, fill = locali)) + 
  stat_identity(geom = "bar", position = "dodge") +
  facet_wrap(~ ulcera) +
  labs(xlab ="Conusmo de aspirina", title = "Desarrollo de una úlcera", fill = "Localización") +
  theme_bw()

De nuevo se observan comportamientos distintos para la combinación de los factores considerados. En este caso el estudio de la independencia total nos llevaría a analizar un modelo donde eliminamos la interacción de orden tres. Dicho modelo se puede obtener y analizar con:

# Ajuste del modelo sin interacción
modelo <- glm(frecuencia ~ (locali + ulcera + aspirina)^2, family = poisson, data = glm_multi_04)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0.01219027

Dado que el pvalor resulta significativo podemos concluir que existe asociación entre los factores considerados. Pasamos ahora a estudiar cada una de las cuestiones planteadas ahora que sabemos que los factores están relacionados. Para saber si la existencia del desarrollo de una úlcera está asociado con el consumo de aspirina para cada tipo de úlcera, ajustamos un modelo sin la interacción aspirina:ulcera

# Ajuste del modelo sin interacción
modelo <- glm(frecuencia ~ locali*ulcera + aspirina + aspirina:locali, family = poisson, data = glm_multi_04)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0.000143616

Se rechaza claramente independencia entre el consumo de aspirina y la existencia de úlcera para ambos tipos de úlcera.

Pasamos ahora a la segunda cuestión, es decir, tratamos de averiguar si existe asociación entre el consumo de aspirina y el tipo de úlcera desarrollada. Para ello debemos analizar la interacción entre aspirina:locali cuando sabemos que se ha desarrollado una úlcera

# Ajuste del modelo sin interacción
modelo <- glm(frecuencia ~ locali*ulcera + aspirina + aspirina:ulcera, family = poisson, data = glm_multi_04)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0.005147627

Rechazamos la hipótesis de independencia a favor de que el consumo de aspirina es diferente en pacientes con distinto tipo de úlcera, pero la pregunta es ¿Cómo podemos cuantificar esas asociaciones?

Modelos Logit-Multinomial

Si en una tabla de contingencia identificamos uno de los factores de clasificación como variable respuesta y el resto como predictoras, sin duda, más que simplemente concluir que existe asociación entre ellos, interesará modelizar los datos para predecir la probabilidad de que se dé una u otra respuesta en la variable respuesta en función de las distintas combinaciones en las variables explicativas. Existe un conjunto de modelos log-lineales, útiles para investigar asociación, que son equivalentes a modelos logit, útiles para predecir una variable a través de otras; son los modelos logit-multinomial.

Consideremos una tabla de contingencia en la que están involucradas como variables de clasificación, una variable respuesta R con I niveles de clasificación, y una variable respuesta A con J niveles de clasificación. Denotamos por \(\pi_{i|j}\) las probabilidades de respuesta \(R = i\) dada cierta clasificación en la variable predictora, \((A = j)\), esto es, \[\pi_{i|j} = P(R = i | A = j).\] Los modelos logit-multinomial vienen definidos en términos de los log-odds (logaritmo del cociente) de las probabilidades condicionadas de dos niveles de respuesta R, dadas las diferentes combinaciones de niveles de las variables predictoras. Al predecir cocientes de probabilidades, se suele elegir uno de los niveles de la variable a explicar como nivel de referencia, generalmente el primero o el último. Si se toma como categoría de referencia la última de la variable a explicar de R (en este caso I), los log-odds a predecir serán: \[log\left(\frac{\pi_{i|j}}{\pi_{I|j}}\right), \text{para cada nivel i de R}\]

El ajuste de estos modelos logit-multinomiales se realiza a través de modelos log-lineales equivalentes con los que se predicen las frecuencias observadas \(Y_{ij}\) en la tabla de contingencia. Es preciso puntualizar que la equivalencia entre estos modelos logit y los log-lineales es cierta asumiendo que los totales marginales correspondientes a las variables predictivas son fijos. Básicamente es decir que las variables predictoras han sido controladas en el estudio con el fin principal de predecir la variable respuesta R. Al asumir esto, obligamos a que los efectos asociados a dichas variables predictoras hayan de aparecer siempre (junto con sus interacciones) en cualquier modelo log-lineal que ajustemos.

Dicha relación viene expresada como: \[log\left(\frac{\pi_{i|j}}{\pi_{I|j}}\right) = log\left(\frac{E(Y_{ij})}{E(Y_{Ij})}\right) = log(E(Y_{ij})) - log(E(Y_{Ij}))\] donde los \(E(Y_{ij})\) para ca \(i\) son los valores predichos del modelo log-lineal para cada combinación de los niveles de todas las variables involucradas. En realidad, esto significa que a la hora de estimar el log-odds sólo debemos tener en cuenta los coeficientes del modelo que hagan referencia a los niveles considerados de la variable respuesta, y a la interacción (si está presente en el modelo) entre respuesta y predictora, es decir: \[log\left(\frac{\pi_{i|j}}{\pi_{I|j}}\right) = (\alpha_i -\alpha_I) + (\alpha\beta_{ij} - \alpha\beta_{Ij})\] donde los \(\alpha\) son los efectos asociados con los niveles de las categorías de la respuesta y los \(\alpha\beta\) son los efectos asociados con la interacción entre respuesta y predictora. Esta relación tiene implicaciones importantes ya que para estimar los logits debemos tener en cuenta únicamente los coeficientes del modelo que están relacionados con la respuesta, a través de sus efectos principales o de sus interacciones con las predictoras. Este principio nos sirve para todas las situaciones de análisis que planteamos a continuación en función de las características de la variable respuesta y la predictora o predictoras.

Respuesta binaria nominal y predictora nominal

Imaginemos que la variable respuesta esta compuesta por dos categorías de carácter nominal (\(i=1, 2\)=. Con la estructura de estimación definida en el punto anterior si la respuesta tiene únicamente dos categorías tendríamos que \(\pi_{1|j} = 1 - \pi_{2|j}\), de forma que la ecuación anterior quedaría:

\[log\left(\frac{\pi_{2|j}}{\pi_{1|j}}\right) = log\left(\frac{\pi_{2|j}}{1-\pi_{2|j}}\right) = (\alpha_2 -\alpha_1) + (\alpha\beta_{2j} - \alpha\beta_{1j})\]

Si denotamos por \(\eta_{2|j}\) al \(log\left(\pi_{2|j}/\pi_{1|j}\right)\), a partir de la expresión anterior resulta posible obtener la probabilidad de cada nivel de la respuesta mediante: \[\pi_{2|j} = \frac{exp(\eta_{2|j})}{1+exp(\eta_{2|j})}\]

Respuesta multinomial nominal y predictora nominal

Cuando tenemos más de dos categorías el procedimiento es bastante similar. Si denotamos por \[log\left(\frac{\pi_{i|j}}{\pi_{I|j}}\right) = \eta_{i|j} = (\alpha_i -\alpha_I) + (\alpha\beta_{ij} - \alpha\beta_{Ij})\] al logit de las categorías \(i\) e \(I\), donde esta última es la de referencia, podemos obtener la probabilidad condicionada de la categoría \(i\) como: \[\pi_{i|j} = \frac{exp(\eta_{i|j})}{\sum_{l=1}^I exp(\eta_{l|j})}\] para cualquier valor de \(i\). Puesto que \(exp(\eta_{I|j}) = 1\) tendríamos: \[\pi_{i|j} = \frac{exp(\eta_{i|j})}{1 + \sum_{l \neq I} exp(\eta_{l|j})}\] Esta expresión es similar a la de dos categorías salvo por la modificación del denominador, donde se incluyen los predictores asociados a todas las categorías de la respuesta.

Respuesta nominal, predictora ordinal

Esta situación es miliar a la anterior salvo por el carácter ordinal de la predictora. Para estimar dichos modelos debemos incluir una variable ficticia (\(S\)) que contendrá los scores asociados con las categorías de la predictora. Habitualmente se usan scores consecutivos de 1 hasta n (categorías de la respuesta), donde el 1 refleja el valor más bajo y el n el más alto. Tendríamos entonces la variable \(v_1=1, v_2=2,...,v_J=J\). El logit asociado a esta situación viene dado por:

\[log\left(\frac{\pi_{i|j}}{\pi_{I|j}}\right) = \eta_{i|j} = (\alpha_i -\alpha_I) + (\alpha\beta_{ij} - \alpha\beta_{Ij})v_j\]

de forma que

\[\frac{\pi_{i|j}}{\pi_{I|j}} = exp(\eta_{i|j}) = exp((\alpha_i -\alpha_I) + (\alpha\beta_{ij} - \alpha\beta_{Ij})v_j)\] y podemos obtener la probabilidad de una categoría como: \[\pi_{i|j} = \frac{exp(\eta_{i|j})}{1 + \sum_{l \neq I} exp(\eta_{l|j})}\]

Respuesta ordinal, predictora nominal

Cuando la variable respuesta tiene carácter ordinal el proceso de estimación debe variar para tener en cuenta dicho carácter ordinal. En esta situación no obtenemos los logit de una categoría con respecto a la de referencia sino que comparamos dos categorías consecutivas \[log\left(\frac{\pi_{i+1|j}}{\pi_{i|j}}\right)\] para valorar el incremento que sufre la probabilidad de un categoría con respecto a su categoría superior.

Para estimar dichos modelos debemos incluir una variable ficticia (\(S\)) que contendrá los scores asociados con las categorías de la respuesta. Habitualmente se usan scores consecutivos de 1 hasta n (categorías de la respuesta), donde el 1 refleja el valor más bajo y el n el más alto. Tendríamos entonces la variable \(u_1=1, u_2=2,...,u_n=n\). Si tenemos únicamente una variable predictora el logit anterior se puede obtener fácilmente ya que:

\[log(E(Y_{ij})) = \mu +\alpha_i + \beta_j + \theta_j S_i\] donde \(\mu\) es la interceptación del modelo, \(\alpha_i\) es el efecto asociado con el nivel \(i\) de la respuesta, \(\beta_j\) es el efecto asociado con el nivel \(j\) de la predictora, y \(\theta\) es el efecto asociado con la variable de scores. En la práctica el logit para dos categorías consecutivas de la respuesta dependen únicamente de los coeficientes del modelo que afectaban a la respuesta y a la interacción entre respuesta y predictora, es decir:

\[log\left(\frac{\pi_{i+1|j}}{\pi_{i|j}}\right) = \eta_{i+1|j} =(\alpha_{i+1} - \alpha_i)+ \theta_j (u_{i+1} - u_i)\ \] de forma que la relación entre ambas probabilidades viene dada por: \[\frac{\pi_{i+1|j}}{\pi_{i|j}} = exp(\eta_{i+1|j}) = exp((\alpha_{i+1} - \alpha_i) + \theta_j (u_{i+1} - u_i))\]

Si queremos obtener la relación entre las probabilidades de las categorías \(i+2\) e \(i\) basta con considerar que: \[\frac{\pi_{i+2|j}}{\pi_{i|j}} = \frac{\pi_{i+2|j}}{\pi_{i+1|j}}\frac{\pi_{i+1|j}}{\pi_{i|j}}\] que se puede calcular como:

\[\frac{\pi_{i+2|j}}{\pi_{i|j}} = exp(\eta_{i+2|j}) = exp((\alpha_{i+2} - \alpha_i) + \theta_j (u_{i+2} - u_i))\]

De esta forma podemos obtener la relación entre cualquier para de categorías de la respuesta.

Respuesta ordinal, predictora ordinal

Imaginemos que tenemos una variable respuesta y otra predictora (ambas de carácter ordinal). La modelización en esta situación pasa por definir variables ficticias de scores tanto para la predictora como la respuesta. Denotamos por \(S\) y \(U\) a las variables de scores para la respuesta y la predictora respectivamente. Calculamos la variable producto \(P\) en ambos scores que será introducida en el modelo para reflejar el posible efecto de interacción entre respuesta y predictora. En esta situación tenemos que: \[log(E(Y_{ij})) = \mu +\alpha_i + \beta_j + \theta P_{ij} = \mu +\alpha_i + \beta_j + \theta S_iU_j\] donde \(\mu\) es la interceptación del modelo, \(\alpha_i\) es el efecto asociado con el nivel \(i\) de la respuesta, \(\beta_j\) es el efecto asociado con el nivel \(j\) de la predictora, y \(\theta\) es el efecto asociado con la variable \(P_{ij} = S_iU_j\). En la práctica se puede demostrar que que el logit para dos categorías consecutivas de la respuesta dependen únicamente de los coeficientes del modelo que afectaban a la respuesta y a la interacción entre respuesta y predictora, es decir:

\[log\left(\frac{\pi_{i+1|j}}{\pi_{i|j}}\right) = \eta_{i+1|j} = (\alpha_{i+1} - \alpha_i) + \theta (u_{i+1} - u_i)v_j\] de forma que la relación entre ambas probabilidades viene dada por: \[\frac{\pi_{i+1|j}}{\pi_{i|j}} = exp(\eta_{i+1|j}) = exp((\alpha_{i+1} - \alpha_i) + \theta (u_{i+1} - u_i)v_j)\]

\[\frac{\pi_{i+2|j}}{\pi_{i|j}} = exp(\eta_{i+2|j}) = exp((\alpha_{i+2} - \alpha_i) + \theta (u_{i+2} - u_i)v_j)\]

De esta forma podemos obtener la relación entre cualquier par de categorías de la respuesta. Actuado de esta forma resulta posible calcular la probabilidad de cada categoría de la respuesta sin más que fijar una categoría de referencia, y obtener todas las probabilidades asociadas con la de referencia. Procederíamos entonces como en el caso multinacional.

Todas estas ecuaciones de estimación se pueden generalizar sin problemas cuando existe más de una variable predictora, sin más que considerar las posibles interacciones entre la respuesta y las predictoras. En caso de tratarse de variables de tipo ordinal construiremos los scores correspondientes para analizar dichos efectos de interacción.

Análisis preliminar

A continuación se presentan los logits empíricos (estimados a partir de las proporciones observadas en cada una de las casillas de la tabla de contingencia) para algunos de los ejemplos más sencillos con los que venimos trabajando en este tema.

Vacuna de la gripe

Si tomamos como variable respuesta la obtenida tras el tratamiento, y usamos como referencia la categoría small podemos obtener los log-odds y representaros mediante:

# Cargamos los datos del ejemplo
glm_multi_03 <- read_csv("https://goo.gl/jPcz7x", col_types = "cci")
# Creamos variable con valores de referencia y obtenemos los log-odds asociados
referencia.val <- unlist(dplyr::select(filter(glm_multi_03,response == "small"),frequency))
glm_multi_03 <- glm_multi_03 %>% 
  mutate(referencia = rep(referencia.val,c(3,3)),
         lodds = log(frequency/referencia))
# Gráfico
# Seleccionamos las categorías a representar eliminando la de referencia
datos <- filter(glm_multi_03,response != "small")
ggplot(datos,aes(x = treatment, y = lodds, group = response, color = response)) +
  geom_point() + 
  geom_line() +
  theme_bw()

Se observa un comportamiento lineal y paralelo para los log-odds correspondientes a la comparación de cada categoría de la respuesta con respecto a la de referencia.

Aspiraciones Universitarias

Si tomamos como variable respuesta las aspiraciones universitarias en función del estatus económico y de la motivación familiar, y usamos como referencia la categoría No podemos obtener el log-odd de la categoría Si mediante:

# Cargamos los datos del ejemplo
estatus <- gl(4,2,16,labels = c("Baja","Media-Baja","Media-Alta","Alta"))
motivacion <- gl(2,1,16,labels=c("Baja","Alta"))
aspiraciones <- gl(2,8,16,labels = c("No","Si"))
frecuencia <- c(749,233,627,330,420,374,153,266,35,133,38,303,37,467,26,800)
glm_multi_05 <- data.frame(estatus,motivacion,aspiraciones,frecuencia)
# Creamos variable con valores de referencia y obtenemos los log-odds asociados
referencia.val <- unlist(dplyr::select(filter(glm_multi_05,aspiraciones == "No"),frecuencia))
glm_multi_05 <- glm_multi_05 %>% 
  mutate(referencia = rep(referencia.val,2),
         lodds = log(frecuencia/referencia))
# Gráfico
# Seleccionamos las categorías a representar eliminando la de referencia
# Como la respuesta sólo tiene una categoría no hace falta utilizarla en el gráfico
datos <- filter(glm_multi_05,aspiraciones != "No")
ggplot(datos,aes(x = estatus, y = lodds, group = motivacion,color = motivacion)) +
  geom_point() + 
  geom_line() +
  theme_bw()

Se observa una tendencia creciente (puede incluso que no lineal) conforme aumenta el estatus económico, que va relacionada directamente con una motivación más alta en el entorno familiar (comportamiento casi paralelo para ambas motivaciones).

Contraceptivos

Si tomamos como variable respuesta el método anticonceptivo utilizado, y usamos como referencia la categroria Ninguno podemos obtener los log-odds de la otras dos categorías como:

# Cargamos los datos del ejemplo
edad <- gl(7,1,21,labels = c("15-19","20-24","25-29","30-34","35-39","40-44","45-49"))
metodo <- gl(3,7,21,labels=c("Esterilación","Otros","Ninguno"))
frecuencia <- c(3,80,216,268,197,150,91,61,137,131,76,50,24,10,232,400,301,203,188,164,183)
glm_multi_06 <- data.frame(edad,metodo,frecuencia)
# Creamos variable con valores de referencia y obtenemos los log-odds asociados
referencia.val <- unlist(dplyr::select(filter(glm_multi_06,metodo == "Ninguno"),frecuencia))
glm_multi_06 <- glm_multi_06 %>% 
  mutate(referencia = rep(referencia.val,3),
         lodds = log(frecuencia/referencia))
# Gráfico
# Seleccionamos las categorías a representar eliminando la de referencia
# Como la respuesta sólo tiene una categoría no hace falta utilizarla en el gráfico
datos <- filter(glm_multi_06,metodo != "Ninguno")
ggplot(datos,aes(x = edad, y = lodds, group = metodo,color = metodo)) +
  geom_point() + 
  geom_line() +
  theme_bw()

Se observan claramente tendencias cuadráticas distintas para cada método asociadas con el grupo de edad considerado. El problema en este topo de situaciones es que como la variable edad es un factor no podemos incluir un efecto cuadrático en el modelo, Sin embargo, aprovecharemos el carácter ordinal de dicha variable para incluir una nueva variable que contenga la marca de clase de cada grupo de edad, y que nos permita ajustar modelos cuadráticos sobre dicha variable.

Conclusión

La utilización de variables ficticias numéricas para sustituir una variable categórica de tipo ordinal es una forma de proceder es muy habitual en los modelos log-lineales, donde se siempre se consideran los denominados scores o valores numéricos asociados a variables de carácter ordinal. En el caso anterior los scores son: 17, 22, 27, 32, 37, 42, y 47. Cuando la variable ordinal viene asociada con categorías y no con números, se suele asociar un código numérico con el valor más bajo de 1 indicando la categoría más baja, y vamos aumentando de 1 en 1 hasta acabar con todos los niveles del factor. En el ejemplo de las aspiraciones universitarias podríamos considerar los scores 1, 2, 3 y 4 para representar las categorías de la variable status económico con 1 igual a “Baja”, y 4 igual a “Alta”.

Por tanto, el carácter ordinal de la respuesta o de las predictoras resulta muy relevante a la hora de obtener el modelo de asociación buscado. En el punto siguiente utilizamos todo el análisis realizado hasta el momento para obtener el modelo más adecuado en lo ejemplos presentados.

Ajuste y estimación

Aspirina

En este ejemplo tenemos la variable consumo de aspirina (\(A\)) de tipo nominal con dos posibles categorías que actúa como respuesta, y dos posibles variables predictoras de tipo nominal localización de la úlcera (\(L\)) y si el sujeto tiene o no úlcera (\(E\)). Ajustamos un modelo con todas las interacciones dobles, ya que en este caso resulta imposible plantear el modelo con las interacción triple, ya que constituye el modelo saturado. La bondad del ajuste planteado viene dado por:

# Ajuste del modelo sin interacción
modelo <- glm(frecuencia ~ (locali + ulcera + aspirina)^2, family = poisson, data = glm_multi_04)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0.01219027

Rechazamos que el efecto de localización y el tipo de úlcera sea distinto para los sujetos que tomaban aspirina de los que no, es decir, existe paralelismo en el comportamiento de los sujetos que tomaban aspirina y los que no. Este comportamiento ya lo pudimos apreciar al representar os logits empíricos para estos datos. Para obtener los logits y probabilidades asociadas a cada categoría de la respuesta extraemos los coeficientes del modelo ajustado.:

summary(modelo)

## 
## Call:
## glm(formula = frecuencia ~ (locali + ulcera + aspirina)^2, family = poisson, 
##     data = glm_multi_04)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -0.5393   0.4486   0.5073  -0.4661   0.7280  -1.2085  -1.0829   1.4673  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.74870    0.15007  24.979  < 2e-16 ***
## localiDuodenal             0.06977    0.20415   0.342  0.73254    
## ulceraNo                   0.32091    0.19384   1.656  0.09782 .  
## aspirinaSi                -0.67905    0.24682  -2.751  0.00594 ** 
## localiDuodenal:ulceraNo   -0.10574    0.26147  -0.404  0.68590    
## localiDuodenal:aspirinaSi -0.70005    0.34603  -2.023  0.04306 *  
## ulceraNo:aspirinaSi       -1.14288    0.35207  -3.246  0.00117 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 127.749  on 7  degrees of freedom
## Residual deviance:   6.283  on 1  degrees of freedom
## AIC: 59.898
## 
## Number of Fisher Scoring iterations: 4

Si queremos estudiar la probabilidad de consumo de aspirina para cualquier combinación de las predictoras debemos trabajar con el logit: \[log\left(\frac{\pi_{A = SI|jk}}{\pi_{A = no|jk}}\right) = (\alpha_{A=Si} - \alpha_{A=No}) + (\theta_{A=Si,jk} - \theta_{A=No,jk})\] donde \(\alpha\) hace referencia al efecto principal asociado con la respuesta, y \(\theta\) hace referencia al coeficiente de interacción entre la respuesta y los niveles \(j\) y \(k\) de \(L\) y \(U\) respectivamente. Por ejemplo, para obtener el logit para un sujeto con \(U = No\) y \(L = Duodenal\) tendríamos que calcular:

\[\eta_\text{A = Si| U = NO; L = Duodenal} = (\alpha_{A=Si} - \alpha_{A=No}) + (\theta_{A=Si,No,Duodenal} - \theta_{A=No,No,Duodenal})\] donde \(\alpha_{A=Si} = -0.67905379\), \(\alpha_{A=No} = 0\), y \(\theta_{A=Si,No,Duodenal}\), \(\theta_{A=No,No,Duodenal}\) representan las interacciones de la respuesta y predictoras para dicha combinación que vienen dadas por:

\[\begin{array}{ll} \theta_{A=Si,No,Duodenal} &= \theta_{A=Si,No} + \theta_{A=Si,Duodenal} = - 1.14287716 - 0.70004571 = -1.842923\\ \theta_{A=No,No,Duodenal} &= \theta_{A=No,No} + \theta_{A=No,Duodenal} = 0 + 0 = 0\\ \end{array}\] de forma que: \[\eta_\text{[A = Si| U = NO; L = Duodenal]} = - 0.67905379 - 1.842923 = -2.521977\]

El resto de logits se estiman de forma similar teniendo en cuenta que \(\alpha_{A=No} = 0\) y \(\theta_{A=No,jk = 0}\) para cualquier combinación de \(j\) y \(k\), de forma que: \[\begin{array}{lll} \eta_\text{[A = Si| U = NO; L = Gástrica]} & = \alpha_{A=Si} + \theta_{A=Si,No} + \theta_{A=Si,Gástrica} &= - 0.67905379 - 1.14287716 + 0 = -1.821931\\ \eta_\text{[A = Si| U = Si; L = Duodenal]} & = \alpha_{A=Si} + \theta_{A=Si,Si} + \theta_{A=Si,Duodenal}&= - 0.67905379 + 0 - 0.70004571 = -1.3791\\ \eta_\text{[A = Si| U = Si; L = Gástrica]} & = \alpha_{A=Si} + \theta_{A=Si,Si} + \theta_{A=Si,Gástrica}&= - 0.67905379 + 0 + 0 = - 0.67905379\\ \end{array}\]

Podemos obtener ahora las probabilidades de consumo de aspirina para cada combinación de las predictoras:

\[\begin{array}{lll} \pi_\text{[A = Si| U = NO; L = Duodenal]} & = exp(-2.521977)/(1+exp(-2.521977)) &= 0.0743318\\ \pi_\text{[A = Si| U = NO; L = Gástrica]} & = exp(-1.821931)/(1+exp(-1.821931)) &= 0.1392023\\ \pi_\text{[A = Si| U = Si; L = Duodenal]} & = exp(-1.3791)/(1+exp(-1.3791)) &= 0.2011536\\ \pi_\text{[A = Si| U = Si; L = Gástrica]} & = exp(- 0.67905379)/(1+exp(- 0.67905379)) &= 0.3364725\\ \end{array}\]

A la vista de estos resultados podemos concluir que la probabilidad de consumir aspirina dado un enfermo de úlcera duodenal es de 0.20, mientras que para un enfermo de úlcera gástrica dicha probabilidad es de 0.33. Podemos concluir que la probabilidad del consumo de aspirina es superior en los enfermos con úlcera gástrica que en la duodenal, y mayor en los sujetos enfermos que en los controles.

A continuación vemos como obtener de forma directa las probabilidades de cada categoría de la respuesta utilizando la función multinom de la librería nnet. Dicha función se escribe como:

modelo <- multinom(respuesta ~ predictoras, weights = frecuencia_var, data = data_set)

donde predictoras indica el modelo de las relaciones entre las predictoras, frecuencia_var identifica la variable que contiene las frecuencias de cada combinación de repuesta y predictoras, y data_set es el conjunto de datos de trabajo. Para obtener las probabilidades asociadas con este modelo usamos la función:

predict(modelo,type = "probs")

Estas funciones solo permiten la obtención de probabilidades en modelos donde la respuesta es de tipo nominal. Veamos su funcionamiento sobre este ejemplo:

# cargamos la librería
library(nnet)
# Ajustamos el modelo correspondiente al glm asociado
mmulti <- multinom(aspirina ~ locali+ulcera, weights = frecuencia,data = glm_multi_04)

## # weights:  4 (3 variable)
## initial  value 173.286795 
## iter  10 value 113.074342
## iter  10 value 113.074342
## iter  10 value 113.074342
## final  value 113.074342 
## converged

# Añadimos las probabilidades del modleo ajustado
glm_multi_04$probabilidad <- predict(mmulti,type = "probs") 
# Seleccionamos los datos para Aspirina = Si
datos.prob <- dplyr::select(filter(glm_multi_04,aspirina == "Si"),locali,ulcera,aspirina,probabilidad)
datos.prob

# Representamos los datos obtenidos
# probabilidad máxima inferior a 0.4
ggplot(datos.prob,aes(x = locali, y = probabilidad, group = ulcera, color = ulcera)) +
  geom_point() + 
  geom_line() +
  labs(title = "Probabilidad del consumo de aspirina dado:", 
       y = "Probabilidad", 
       x = "Localización de la úlcera") +
  scale_y_continuous(breaks = seq(0,0.4,by = 0.05),limits = c(0,0.4)) +
  theme_bw()

El gráfico refleja claramente el comportamiento de las probabilidades de consumo de aspirina de acuerdo a las dos variables predictoras consideradas.

Aspiraciones Universitarias

En primer lugar establecemos el modelo con interacciones dobles para estudiar la bondad del ajuste, es decir, conocer si podemos prescindir de la interacción de orden triple. Esto implicaría que no rechazamos la hipótesis de que el efecto de la motivación es la proporción de individuos que aspiran a la universidad sea el mismo para cualquier estatus económico, tal y como pudimos ver en el comportamiento de los log-odds (tendencia paralelas). Dado que las variables predictoras pueden considerarse como ordinales, podríamos considerar sus scores y utilizarlos para construir el modelo de predicción teniendo en cuenta una posible tendencia cuadrática asociada con el estatus económico, tal y como vimos en los logits empíricos. En primer lugar consideramos ambas variables como nominales.

# Ajuste del modelo sin interacción triple
modelo <- glm(frecuencia ~ (estatus + motivacion + aspiraciones)^2, family = poisson, data = glm_multi_05)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0.664965

EL modelo ajustado es adecuado y no rechazamos la hipótesis asociada con la interacción triple. Vamos a obtener los logit asociados con este modelo tomando como referencia la categoría “No” de la variable respuesta “Aspiraciones”. El resumen del modelo ajustado viene dado por:

# Ajuste del modelo sin interacción triple
summary(modelo)

## 
## Call:
## glm(formula = frecuencia ~ (estatus + motivacion + aspiraciones)^2, 
##     family = poisson, data = glm_multi_05)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7  
## -0.15119   0.27320   0.04135  -0.05691  -0.04446   0.04719   0.32807  
##        8         9        10        11        12        13        14  
## -0.24539   0.73044  -0.35578  -0.16639   0.05952   0.15116  -0.04217  
##       15        16  
## -0.75147   0.14245  
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       6.62426    0.03602 183.918  < 2e-16 ***
## estatusMedia-Baja                -0.18496    0.05304  -3.487 0.000489 ***
## estatusMedia-Alta                -0.58183    0.05931  -9.810  < 2e-16 ***
## estatusAlta                      -1.62046    0.08450 -19.178  < 2e-16 ***
## motivacionAlta                   -1.19117    0.07166 -16.622  < 2e-16 ***
## aspiracionesSi                   -3.19497    0.11850 -26.962  < 2e-16 ***
## estatusMedia-Baja:motivacionAlta  0.55410    0.09469   5.852 4.87e-09 ***
## estatusMedia-Alta:motivacionAlta  1.07056    0.09649  11.095  < 2e-16 ***
## estatusAlta:motivacionAlta        1.78588    0.11444  15.606  < 2e-16 ***
## estatusMedia-Baja:aspiracionesSi  0.42013    0.11768   3.570 0.000357 ***
## estatusMedia-Alta:aspiracionesSi  0.73851    0.11382   6.488 8.69e-11 ***
## estatusAlta:aspiracionesSi        1.59311    0.11527  13.820  < 2e-16 ***
## motivacionAlta:aspiracionesSi     2.68292    0.09867  27.191  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3211.0014  on 15  degrees of freedom
## Residual deviance:    1.5755  on  3  degrees of freedom
## AIC: 141.39
## 
## Number of Fisher Scoring iterations: 3

Dado que la estructura de estimación es similar a la del ejemplo anterior procedemos directamente con la estimación de las probabilidades de cada categoría de la respuesta.

# Ajustamos el modelo correspondiente al glm asociado
mmulti <- multinom(aspiraciones ~ estatus + motivacion, weights = frecuencia,data = glm_multi_05)

## # weights:  6 (5 variable)
## initial  value 3459.497578 
## iter  10 value 2346.838055
## final  value 2346.837322 
## converged

# Añadimos las probabilidades del modleo ajustado
glm_multi_05$probabilidad <- predict(mmulti,type = "probs") 
# Seleccionamos los datos para aspiraciones = Si
datos.prob <- dplyr::select(filter(glm_multi_05,aspiraciones == "Si"),estatus,motivacion,aspiraciones,probabilidad)
datos.prob

# Gráfico de resultados
# probabilidad máxima inferior a 0.8
ggplot(datos.prob,aes(x = estatus, y = probabilidad, group = motivacion, color = motivacion)) +
  geom_point() + 
  geom_line() +
  labs(title = "Probabilidad de aspiraciones universitarias dado:", 
       y = "Probabilidad", 
       x = "Estatus económico") +
  scale_y_continuous(breaks = seq(0,0.8,by = 0.05),limits = c(0,0.8)) +
  theme_bw()

La probabilidad de tener aspiraciones universitarias aumenta con el estatus económico de la familia y con la motivación que esta aporta. Por ejemplo, la probabilidad de aspiraciones universitarias para una familia con estatus económico bajo y motivación baja es de 0.04, mientras que cuando el estatus económico es alto y la motivación familiar es alta, dicha probabilidad se sitúa en 0.75.

Contraceptivos

Dado que sólo tenemos una predictora vamos a ajustar un modelo sólo con efectos principales. Estudiamos la bondad del ajuste de dicho modelo para conocer el patrón de asociación entre las variables. Si desechamos la interacción doble concluiríamos que no existe asociación entre edad y método contraceptivo. En caso contrario deberíamos establecer el modelo de asociación entre ambas variables. El gráfico de log-odds ya mostraba cierto grado de asociación de tipo cuadrático con la edad. Dado que dicha variable es de tipo ordinal, si ajustamos un modelo se interacción con dicha variable deberíamos considerar la introducción de scores para obtener la predicción de la probabilidad de cada categoría de la respuesta.

# Ajuste del modelo sin interacción triple
modelo <- glm(frecuencia ~ edad + metodo, family = poisson, data = glm_multi_06)
# Bondad del asjute del modelo
1-pchisq(modelo$deviance,modelo$df.residual)

## [1] 0

Rechazamos la hipótesis de independencia. Creamos la variable numérica asociada con la edad y ajustamos un nuevo modelo con tendencia cuadrática. En este modelo incluimos la edad como factor principal, y los scores de la edad en interacción con la variable método anticonceptivo.

# Creamos los scores numéricos tomando el punto medio de cada intervalo
glm_multi_06$edad.num <- rep(c(17, 22, 27, 32, 37, 42, 47),3)
# Ajuste del modelo con efecto cuadrático en edad
# Debemos incluir el efecto edad como factor principal para mantener la estructura del modelo
modelo2 <- glm(frecuencia ~ edad + metodo*(edad.num + I(edad.num^2)), family = poisson, data = glm_multi_06)

Comparamos los dos modelos obtenidos para validar la inclusión de la variable de scores de edad

anova(modelo,modelo2,test = "Chisq")

El contraste resulta significativo indicando que el modelo de scores es mejor que el modelo de independencia. Veamos los coeficientes del modelo:

summary(modelo2)

## 
## Call:
## glm(formula = frecuencia ~ edad + metodo * (edad.num + I(edad.num^2)), 
##     family = poisson, data = glm_multi_06)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.47791  -0.59984   0.01931   0.48106   1.69573  
## 
## Coefficients: (2 not defined because of singularities)
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  2.2555965  0.1771772  12.731  < 2e-16 ***
## edad20-24                    2.1425538  0.1228478  17.441  < 2e-16 ***
## edad25-29                    3.0427215  0.1733427  17.553  < 2e-16 ***
## edad30-34                    3.2791741  0.1931540  16.977  < 2e-16 ***
## edad35-39                    3.1641959  0.1968141  16.077  < 2e-16 ***
## edad40-44                    2.7848473  0.1928054  14.444  < 2e-16 ***
## edad45-49                    2.1318918  0.1816355  11.737  < 2e-16 ***
## metodoOtros                  7.8467822  0.9194645   8.534  < 2e-16 ***
## metodoNinguno               12.2657333  0.7350314  16.687  < 2e-16 ***
## edad.num                            NA         NA      NA       NA    
## I(edad.num^2)                       NA         NA      NA       NA    
## metodoOtros:edad.num        -0.4406721  0.0592899  -7.433 1.07e-13 ***
## metodoNinguno:edad.num      -0.6999859  0.0448533 -15.606  < 2e-16 ***
## metodoOtros:I(edad.num^2)    0.0049747  0.0009263   5.371 7.85e-08 ***
## metodoNinguno:I(edad.num^2)  0.0097327  0.0006588  14.773  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1521.567  on 20  degrees of freedom
## Residual deviance:   20.475  on  8  degrees of freedom
## AIC: 182.06
## 
## Number of Fisher Scoring iterations: 4

Con una variable respuesta nominal multinomial y una predictora ordinal las ecuaciones de los logits (tomando la categoría “Ninguno” como referencia) vienen dados por:

\[log\left(\frac{\pi_{\text{Esterilizacion,j}}}{\pi_{\text{Ninguno,j}}}\right) = (\alpha_{\text{Esterilizacion}} - \alpha_{\text{Ninguno}}) + (\theta_{Esterilizacion,j} - \theta_{Ninguno,j}) v_j + (\theta^2_{Esterilizacion,j} - \theta^2_{Ninguno,j}) v^2_j\]

\[log\left(\frac{\pi_{\text{Otros,j}}}{\pi_{\text{Ninguno,j}}}\right) = (\alpha_{\text{Otros}} - \alpha_{\text{Ninguno}}) + (\theta_{Otros,j} - \theta_{Ninguno,j}) v_j + (\theta^2_{Otros,j} - \theta^2_{Ninguno,j}) v^2_j\]

donde los \(\alpha\) representan los efectos principales, \(\theta\) los efectos de las interacciones con el efecto lineal del score, \(\theta^2\) los efectos de las interacciones con el efecto cuadrático del score, y los \(v_j\) son los scores de la edad. Sustituyendo tendríamos:

\[log\left(\frac{\pi_{\text{Esterilizacion,j}}}{\pi_{\text{Ninguno,j}}}\right) = (0 - 12.2657333) + (0 + 0.6999859) v_j + (0 - 0.0097327) v^2_j = - 12.2657333 + 0.6999859 v_j - 0.0097327 v^2_j\]

\[log\left(\frac{\pi_{\text{Otros,j}}}{\pi_{\text{Ninguno,j}}}\right) = (7.8467822 - 12.2657333) + (-0.4406721 + 0.6999859) v_j + (0.0049747 - 0.0097327) v^2_j = - 4.418951 + 0.2593138 v_j - 0.004758 v^2_j\]

Con estas expresiones resulta posible obtener la probabilidad de cada categoría de la respuesta (salvo la de referencia) en función de la edad. Procedemos de forma directa con la formulación presentada para los modelos logit con más de dos categorías en la variable respuesta. Veamos el gráfico resultante:

# Secuencia de edad
edad <- seq(15,49,1)
# Valores eta para cada categoría comparada con la de referencia
EvsN <- exp(- 12.2657333 + 0.6999859*edad - 0.0097327*edad^2)
OvsN <- exp(- 4.418951 + 0.2593138*edad - 0.004758*edad^2)
# Probabilidad condicionada de cada categoria
Esteril <- EvsN/(1 + EvsN + OvsN)
Otros <- OvsN/(1 + EvsN + OvsN)
datos <- data.frame(edad,Esteril,Otros)
# Reorganizamos los datos para representarlos
datos.comb <- datos %>% 
  gather(Metodo,probabilidad,Esteril:Otros) 
# Gráfico
ggplot(datos.comb,aes(x = edad, y = probabilidad, group = Metodo, color = Metodo)) +
  geom_line() +
  labs(title = "Probabilidad de uso del método anticonceptivo:", 
       y = "Probabilidad", 
       x = "Edad") +
  scale_y_continuous(breaks = seq(0,0.6,by = 0.1),limits = c(0,0.6)) +
  scale_x_continuous(breaks = seq(15,50,by = 5),limits = c(15,50)) +
  theme_bw()

¿Qué conclusiones podemos extraer de este gráfico? ¿qué método predomina a los 20 años? ¿y a los 40?

Vacuna de la gripe

En este caso queremos estudiar un posible modelo de asociación entre la respuesta a la vacuna y el tipo de tratamiento seguido. La variable respuesta es de tipo ordinal y asignamos la variable de scores de 1 a 3, donde 1 refleja efecto pequeño y 3 refleja un efecto grande. En este caso modelizaremos los logits de categorías consecutivas: Moderado vs Pequeño, y Grande vs Moderado. Planteamos el análisis de independencia (sin interacción entre el tratamiento y la variable de scores) frente al de asociación (interacción entre tratamiento y scores):

# Creamos variables de scores
glm_multi_03$efecto.s <-c() # Creamos un vector vacio que rellenamos con los scores
glm_multi_03$efecto.s[glm_multi_03$response == "small"] <- 1
glm_multi_03$efecto.s[glm_multi_03$response == "moderate"] <- 2
glm_multi_03$efecto.s[glm_multi_03$response == "large"] <- 3
# Ajuste del modelo de independencia
modelo <- glm(frequency ~ treatment + response + efecto.s, family = poisson, data = glm_multi_03)
# Ajuste del modelo de asociación
modelo2 <- glm(frequency ~ treatment + response + treatment:efecto.s, family = poisson, data = glm_multi_03)
# Comparamos ambos modelos
anova(modelo,modelo2,test = "Chisq")

Rechazamos el modelo de independencia a favor del modelo de asociación. Veamos las estimaciones del modelo

summary(modelo2)

## 
## Call:
## glm(formula = frequency ~ treatment + response + treatment:efecto.s, 
##     family = poisson, data = glm_multi_03)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6  
##  0.3932  -1.1834   0.9992  -0.7112   0.9759  -0.5468  
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                4.87993    0.79782   6.117 9.56e-10 ***
## treatmentvaccine          -2.32157    0.69261  -3.352 0.000803 ***
## responsemoderate           0.09282    0.32850   0.283 0.777519    
## responsesmall             -0.48965    0.44919  -1.090 0.275683    
## treatmentplacebo:efecto.s -1.25109    0.36266  -3.450 0.000561 ***
## treatmentvaccine:efecto.s       NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 23.8072  on 5  degrees of freedom
## Residual deviance:  4.3106  on 1  degrees of freedom
## AIC: 39.439
## 
## Number of Fisher Scoring iterations: 4

Para este modelo Los logits asociados vienen dados por:

\[log\left(\frac{\pi_{\text{moderate,j}}}{\pi_{\text{small,j}}}\right) = (\alpha_{\text{moderate}} - \alpha_{\text{small}}) + \theta_j(u_{moderate} - u_{small})\] \[log\left(\frac{\pi_{\text{large,j}}}{\pi_{\text{moderate,j}}}\right) = (\alpha_{\text{large}} - \alpha_{\text{moderate}}) + \theta_j(u_{large} - u_{moderate})\]

que al sustituir por las correspondientes estimaciones y scores proporciona:

\[log\left(\frac{\pi_{\text{moderate,j}}}{\pi_{\text{small,j}}}\right) = (0.09282 - (-0.48965)) + \theta_j (2 -1) = 0.58247 + \theta_j \] \[log\left(\frac{\pi_{\text{large,j}}}{\pi_{\text{moderate,j}}}\right) = (0 - 0.09282) + \theta_j (3 - 2) = -0.09282 + \theta_j\]

con \(\theta_1 = -1.25109\) para el placebo y \(\theta_2 = 0\) para el tratamiento. De esta forma las relaciones entre las probabilidades vienen dadas por:

\[\frac{\pi_{\text{moderate,j}}}{\pi_{\text{small,j}}} = exp(0.58247 + \theta_j)\] \[\frac{\pi_{\text{large,j}}}{\pi_{\text{moderate,j}}} = exp(-0.09282 + \theta_j)\]

Además resulta posible obtener: \[\frac{\pi_{\text{large,j}}}{\pi_{\text{small,j}}} = exp((0 - (-0.48965)) + \theta_j (3 - 1)) = exp(0.48965 + 2\theta_j)\]

Evaluando dichas expresiones tenemos que para los sujetos tratados: \[\pi_{\text{moderate}} = 1.8 * \pi_{\text{small}} \text{; } \pi_{\text{large}} = 0.9 * \pi_{\text{moderate}} \text{; y } \pi_{\text{large}} = 1.6 * \pi_{\text{small}}\] mientras que en los controles las relaciones son: \[\pi_{\text{moderate}} = 0.5 * \pi_{\text{small}} \text{; } \pi_{\text{large}} = 0.3 * \pi_{\text{moderate}} \text{; y } \pi_{\text{large}} = 0.2 * \pi_{\text{small}}\]

Podemos concluir que hay un mayor cambio en los sujetos tratados que en los controles, además de que antes ya habíamos indicado que existía una asociación entre el tratamiento y la respuesta. Si tomamos como categoría de referencia “small” podemos obtener las probabilidades de cada categoría partir de las expresiones de \(\frac{\pi_{\text{large,j}}}{\pi_{\text{small,j}}}\) y \(\frac{\pi_{\text{moderate,j}}}{\pi_{\text{small,j}}}\).

tratamiento <- c("Control","Tratado")
eta.LM <- c(exp(0.48965 + 2*(-1.25109)),exp(0.48965))
eta.MS <- c( exp(0.58247 - 1.25109),exp(0.58247))
# Obtenemos a partir de los anteriores
eta.LS <- eta.LM*eta.MS
# Datos
datos <-data.frame(tratamiento,eta.LM,eta.MS,eta.LS)
datos <- dplyr::select(datos,tratamiento,eta.MS,eta.LS)
datos

# Probabilidades
orden<-c("small","moderate","large")
probabilidades <- datos %>% 
  mutate(large = eta.LS/(1+eta.LS+eta.MS),
         moderate = eta.MS/(1+eta.LS+eta.MS),
         small = 1- large - moderate) %>%
  dplyr::select(tratamiento,large,moderate,small) %>%
  gather(`large`,`moderate`,`small`,key = efecto, value = probabilidad)
probabilidades

# Gráfico 
ggplot(probabilidades,aes(x = efecto, y = probabilidad, group = tratamiento, color = tratamiento)) +
  geom_line() +
  labs(y = "Probabilidad", 
       x = "Efecto") +
  scale_y_continuous(breaks = seq(0,0.7,by = 0.1),limits = c(0,0.7)) +
  scale_x_discrete(limits=orden) +
  theme_bw()

¿Cómo interpretamos el gráfico?

Satisfacción laboral

En este caso queremos estudiar un posible modelo de asociación entre la satisfacción laboral en función de los ingresos percibidos. Tanto la variable respuesta como la predictora se encuentran en escala ordinal. Debemos introducir los scores asociados a ambas variables para poder considerarlos en el modelo de asociación. Creamos scores en escala continua de 1 a 4 para la satisfacción (\(u_i\)), y de 1 a 4 para los ingresos (\(v_j\)) (ya que los intervalos considerados no tienen la misma amplitud en este caso). La interacción entre ambas variables se reflejará mediante el producto de las variables de scores creadas.

glm_multi_07 = read_csv("https://goo.gl/5T0nh0", col_types = "ciiii")
# Construimos la tabla donde cada fila recoje la frecuenia observada para combinación de las variables
glm_multi_07 = glm_multi_07 %>% 
  gather(`Muy insatisfecho`,`Poco insatisfecho`,`Moderadamente satisfecho`,`Muy satisfecho`,key = "Estado", value = frecuencia)
# Generamos las variables de scores y su interacción
# En lugar de ir selccionando cada  valor lo hacemos de golpe
glm_multi_07 <- glm_multi_07 %>%
  mutate(ingresos.s = as.numeric(gl(4,1,16)),
         estado.s = as.numeric(gl(4,4,16)),
         asociacion = ingresos.s*estado.s)
# Ajuste del modelo de independencia
modelo <- glm(frecuencia ~ Ingresos + Estado, family = poisson, data = glm_multi_07)
# Ajuste del modelo de asociación
modelo2 <- glm(frecuencia ~ Ingresos + Estado + asociacion, family = poisson, data = glm_multi_07)
# Comparamos ambos modelos
anova(modelo,modelo2,test = "Chisq")

Rechazamos el modelo de independencia a favor del modelo de asociación que obtenemos con los scores. Veamos el resumen del modelo:

summary(modelo2)

## 
## Call:
## glm(formula = frecuencia ~ Ingresos + Estado + asociacion, family = poisson, 
##     data = glm_multi_07)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.03435  -0.17539   0.00509   0.27032   0.57984  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              3.98067    0.13123  30.334  < 2e-16 ***
## Ingresos> 25000         -1.26389    0.36669  -3.447 0.000567 ***
## Ingresos15000 - 25000   -0.57692    0.24745  -2.331 0.019729 *  
## Ingresos6000 - 15000    -0.01053    0.14344  -0.073 0.941479    
## EstadoMuy insatisfecho  -1.13006    0.20826  -5.426 5.76e-08 ***
## EstadoMuy satisfecho    -0.01809    0.11720  -0.154 0.877350    
## EstadoPoco insatisfecho -0.82254    0.13813  -5.955 2.60e-09 ***
## asociacion               0.11194    0.03641   3.075 0.002108 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 445.7627  on 15  degrees of freedom
## Residual deviance:   2.3859  on  8  degrees of freedom
## AIC: 107.42
## 
## Number of Fisher Scoring iterations: 4

Las ecuaciones de los logit para este modelo vienen dadas por \[\frac{\pi_{\text{Muy Satisfecho|j}}}{\pi_{\text{Moderadamente Satisfecho|j}}} = (\alpha_{\text{Muy Satisfecho|j}} - \alpha_{\text{Moderadamente Satisfecho|j}} ) + \theta (u_{\text{Muy Satisfecho|j}} - u_{\text{Moderadamente Satisfecho|j}})v_j\]

\[\frac{\pi_{\text{Moderadamente Satisfecho|j}}}{\pi_{\text{Poco Insatisfecho|j}}} = (\alpha_{\text{Moderadamente Satisfecho|j}} - \alpha_{\text{Poco Insatisfecho|j}} ) + \theta (u_{\text{Moderadamente Satisfecho|j}} - u_{\text{Poco Insatisfecho|j}})v_j\]

\[\frac{\pi_{\text{Poco Insatisfecho|j}}}{\pi_{\text{Muy Insatisfecho|j}}} = (\alpha_{\text{Poco Insatisfecho|j}} - \alpha_{\text{Muy Insatisfecho|j}} ) + \theta (u_{\text{Poco Insatisfecho|j}} - u_{\text{Muy Insatisfecho|j}})v_j\]

Sustituyendo por los coeficientes estimados del modelo tenemos:

\[log\left(\frac{\pi_{\text{Muy Satisfecho|j}}}{\pi_{\text{Moderadamente Satisfecho|j}}}\right) = (-0.0180 - 0) + 0.11194*(4-3)* v_j = -0.0180 + 0.11194 * v_j\]

\[log\left(\frac{\pi_{\text{Moderadamente Satisfecho|j}}}{\pi_{\text{Poco Insatisfecho|j}}}\right) = (0 - (-0.82254) + 0.11194*(3-2)* v_j = 0.82254 + 0.11194 * v_j\]

\[log\left(\frac{\pi_{\text{Poco Insatisfecho|j}}}{\pi_{\text{Muy Insatisfecho|j}}}\right) = (-0.82254 - (-1.13006)) + 0.11194*(2-1)* v_j = 0.30752 + 0.11194 * v_j\]

El coeficiente \(\theta\) refleja el crecimiento que sufre el logit por cada aumento del score (cambio entre los niveles) de la variable predictora. A continuación se calculan los ratios (y se representan gráficamente) de las probabilidades para los diferentes niveles de la predictora, identificando los niveles de la respuesta como “Muy Insatisfecho = S1”, “Poco Insatisfecho = S2”, “Moderadamente Satisfecho = S3”, y “Muy Satisfecho = S4”:

# scores de ingresos
score.ingresos <- 1:4
ingresos <- c("<6000","6000-15000","15000-25000",">25000")
# Relaciones entre probabilidades
S4vsS3 <- exp(-0.0180 + 0.11194 * score.ingresos)
S3vsS2 <- exp(0.82254 + 0.11194 * score.ingresos)
S2vsS1 <- exp(0.30752 + 0.11194 * score.ingresos)
# Utilizando estas
S4vsS2 <- S4vsS3*S3vsS2
S4vsS1 <- S4vsS3*S3vsS2*S2vsS1
# Banco de datos
datos <- data.frame(score.ingresos,ingresos,S4vsS3,S3vsS2,S2vsS1,S4vsS2,S4vsS1)
datos.comb <- dplyr::select(datos,score.ingresos,ingresos,S4vsS3,S4vsS2,S4vsS1)
datos.comb

Podemos ver que el ratio entre la probabilidad de muy satisfecho y muy insatisfechos para los que ganan menos de 6000 dólares es de 4.25 mientras que se va hasta los 11.65 para los que ganan más de 25000 dólares. Se aprecia un crecimiento en los ratios en función de los ingresos

Vamos a obtener las probabilidades de cada categoría a partir de los ratios obtenidos tomando como referencia la categoría “Muy Insatisfecho”.

probabilidades <- datos%>% 
  mutate(eta21 = S2vsS1,
         eta31 = S3vsS2*eta21,
         eta41 = S4vsS3*eta31,
         Muy.satisfecho = eta41/(1+eta41+eta31+eta21),
         Moderadamente.satisfecho = eta31/(1+eta41+eta31+eta21),
         Poco.insatisfecho = eta21/(1+eta41+eta31+eta21),
         Muy.insatisfecho = 1 - Muy.satisfecho - Moderadamente.satisfecho - Poco.insatisfecho) %>%
  dplyr::select(score.ingresos,ingresos,Muy.satisfecho,Moderadamente.satisfecho,Poco.insatisfecho,Muy.insatisfecho) %>%
  gather(`Muy.satisfecho`, `Moderadamente.satisfecho`, `Poco.insatisfecho`, `Muy.insatisfecho`, key = Estado,value = Probabilidad)
# Gráfico 
ggplot(probabilidades,aes(x = score.ingresos, y = Probabilidad, group = Estado, color = Estado)) +
  geom_line() +
  labs(y = "Probabilidad", 
       x = "Ingresos") +
  scale_y_continuous(breaks = seq(0,0.6,by = 0.1),limits = c(0,0.6)) +
  scale_x_discrete(limits = ingresos) +
  theme_bw()

¿Qué conclusiones podemos extraer de este gráfico?

Diagnóstico

El diagnóstico de este tipo de modelos se parece al diagnóstico de los modelos de poisson. Dado que sólo tenemos un valor predicho y un residuo por cada combinación de niveles el diagnóstico nos limitaremos a:

Realizar el gráfico de residuos vs ajustados del predictor lineal.
Realizar el gráfico de valores predichos (frecuencias predichas) vs observados (frecuenias observadas).

Aspirina

Obtnemos los valores ajustados (frecuencias predichas) y residuos para este modelo.

# Ajuste del modelo sin interacción
modelo <- glm(frecuencia ~ (locali + ulcera + aspirina)^2, family = poisson, data = glm_multi_04)
# Datos para el análisis
datos.comb <- glm_multi_04 %>% 
  dplyr::select(frecuencia,locali,ulcera,aspirina) %>%
  mutate(residuos = resid(modelo), # residuos del predictor lineal
         aju.frec = predict(modelo,type ="response"), # predicción de frecuencias
         aju.lineal = predict(modelo))  # predicción del predictor lineal
# Gráfico Residuos vs ajustados
ggplot(datos.comb,aes(aju.lineal,residuos)) + 
  geom_point() +
  labs(x = "Ajustados", y = "Residuos") + 
  geom_hline(yintercept = 0) +
  theme_bw()

# Gráfico Observados vs ajustados
ggplot(datos.comb,aes(aju.frec,frecuencia)) + 
  geom_point() +
  labs(x = "Ajustados", y = "Observados") + 
  geom_abline(intercept = 0, slope = 1) +
  theme_bw()

El modelo ajustado muestra un buen ajuste.

Obtén el diagnótico para el resto de modelos tratados en este tema.

Modelos log-lineales para tablas de contingencia

Introducción

Modelos Log-lineales

Melanoma

Aspirina

Modelos Logit-Multinomial

Respuesta binaria nominal y predictora nominal

Respuesta multinomial nominal y predictora nominal

Respuesta nominal, predictora ordinal

Respuesta ordinal, predictora nominal

Respuesta ordinal, predictora ordinal

Análisis preliminar

Vacuna de la gripe

Aspiraciones Universitarias

Contraceptivos

Conclusión

Ajuste y estimación

Aspirina

Aspiraciones Universitarias

Contraceptivos

Vacuna de la gripe

Satisfacción laboral

Diagnóstico

Aspirina

Bibliografía