---

Modelos Estadísticos. Grado Biotecnología

---

Resolución Ejercicios Modelos Lineales Generalizados para respuesta binomial

Abstract

Guia para la resolución de ejercicios sobre modelos lineales. Se muestra el código necesario para resolver cada uno de los ejercicios planteados. Se recomienda ir respondiendo a las preguntas que se van planteando a lo largo de cada uno de los ejercicios.

---

Librerías y funciones

---

# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)
library(reshape2)
library(lmtest)
library(mgcv)
library(MASS)
library(modelr)
library(survival)
library(survminer)

---

Caso 1

---

En un experimento se sometió a cierto número de cucarachas a cinco horas de exposición a disulfato de carbono gaseoso a varias concentraciones. Se pretendía investigar la relación existente entre la dosis de disulfato administrada y la resistencia de los insectos; si existe tal relación, determinar la dosis a la cual es posible garantizar el exterminio del 50% de los insectos.

caso1 = read_csv("https://goo.gl/E2MlSZ", col_types = "dii")
# Calculamos los vivos para el ajuste de modelos
caso1 = caso1 %>% 
  mutate(alive = number - dead, probabilidad = dead / number)

Representamos la probabilidad de muerte en función de la dosis utilizada.

ggplot(caso1, aes(x = dose, y = probabilidad)) +
  geom_point() +
  labs(x = "Dosis", y = "Probabilidad de muerte") +
  geom_smooth(method = loess, se = FALSE) +
  theme_bw()

Estimación del modelo

¿Qué tipo de modelo consideramos para esta situación experimental? ¿El modelo obtenido es adecuado para explicar la probabildiad de muerte en función de la concentración de disulfato de carbono gaseoso? ¿Cuál es la ecuación del modelo ajustado? ¿Cómo interpretamos la interceptación del modelo? ¿Y la pendiente?

# Variable repuesta agrupada
yres <- cbind(caso1$dead,caso1$alive)
M1 <- glm(yres  ~ dose, family = binomial(link = logit), caso1)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

Obtención de la dosis letal 50.

# Dosis letal 0.5.  Despejamos de la ecuación estimada
prb <- 0.5
LD50 <- (log(prb /(1-prb)) + 60.717) / 34.270
LD50

Diagnóstico del modelo

Calculamos los valores para el diagnóstico del modelo y realizamos los gráficos de diagnóstico. ¿Qué opinas sobre los gráficos de diagnóstico obtenidos? ¿Necesitamos modificar nuestro modelo de partida? Justifica tu respuesta.

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","dose",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) + 
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Proponemos un modelo polinómico de grado 2 en la dosis. ¿El modelo obtenido es adecuado para explicar la probabilidad de muerte en función de la concentración de disulfato de carbono gaseoso? ¿Cuál es la ecuación del modelo ajustado? ¿Cómo interpretamos la interceptación del modelo? ¿Y las pendientes? ¿Cómo obtenemos la dosis letal 50 en esta situación? ¿Qué valor porporciona?

# Variable repuesta agrupada
yres <- cbind(caso1$dead,caso1$alive)
M1 <- glm(yres  ~ dose + I(dose^2), family = binomial(link = logit), caso1)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

Realizamos el diagnóstico el nuevo modelo. ¿Qué opinas de los gráficos obtenidos? ¿podemos considerar como válido el modelo obtenido o debemos plantear una nueva modificación?

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","dose",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Predicción del modelo

Para obtener la predicción en este caso obtenemos las cantidades de interés, valoramos la coincidencia entre éxitos estimados y observados, y representamos la probabilidad de éxito asociada con la dosis utilizada.

¿Qué podemos comentar sobre los valores predichos por el modelo? ¿Cuál es la probabilidad de morir a concentraciones de dosis iagual a 1.7, 1.8 y 1.9? ¿Qué dosis propones usar para conseguir una mortalidad superior al 90%?

# Predicción
prediccion <- predict.glm(M1,type = "response", se.fit = TRUE)
# éxitos estimados
exitos <- round(prediccion$fit*caso1$number,0)
exitos
# comparativa
comparativa <- data.frame(dose = caso1$dose, original = caso1$dead, predicho = exitos)
comparativa

Curva de probabilidad de morir

# Valores paa predicción
newdata <- data.frame(dose = rep(seq(min(caso1$dose), max(caso1$dose), .005)))
# Obtenemos la predicción para el modelo ajustado
newdata <- data.frame(newdata, 
                      predict(M1, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza y añadimos la dosis 
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit)
)
# Gráfico de la predicción con log dosis
ggplot(newdata, aes(x = dose, y = fit)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
  scale_x_continuous(breaks = seq(1.5,2,0.025)) +
  scale_y_continuous(breaks = seq(0,1,0.1)) +
  labs(x = "Dosis", 
       y = "Probabilidad de morir", 
       title = "Bandas de predicción") +
  theme_bw()

---

Caso 2

---

Se realiza un experimento in vitro para estimar el número de anteras embriogénicas de las especies de plantas Datura innoxia Mill bajo dos condiciones experimentales. El primer tratamiento consiste en almacenar a 3° C durante 48 horas, y el segundo consiste en un control donde no se aplica ningún tratamiento. Además se considera una variable que representa los tres valores de fuerza de centrifugación. Las variables registradas son total (número de anteras totales), embryogenic (numero de anteras embriogénicas), storage (tipo de almacenamiento), centrifuge (velocidad centrifugación). Es de interés en el análisis investigar si efectivamente se demostraba un mayor número de anteras embriogénicas para las diferentes condiciones experimentales.

caso2 = read_csv("https://goo.gl/6P3zRr", col_types = "ddcd")
# Recodificación del factor y variable de no embryogenic
caso2 = caso2 %>% 
  mutate(storage=fct_recode(storage,"Control" = "1","treatment" = "2"), 
         nembrig = total - embryogenic)

Estimación del modelo

¿Qué tipo de modelo de partida consideramos para esta situación experimental? ¿Cuál es el modelo obtenido tras la selección de efectos? ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de los efectos presentes en é? ¿El modelo obtenido es adecuado para explicar la probabilidad antera embriogénica en función del tipo de almacenamiento y la velocidad de centrifugación?

# Variable repuesta agrupada
yres <- cbind(caso2$embryogenic,caso2$nembrig)
M1 <- glm(yres  ~ centrifuge*storage, family = binomial(link = logit), caso2)
M1 <- step(M1, test = "Chisq")

# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

Diagnóstico del modelo

Calculamos los valores para el diagnóstico del modelo y realizamos los gráficos de diagnóstico. ¿Qué opinas sobre los gráficos de diagnóstico obtenidos? ¿Necesitamos modificar nuestro modelo de partida? Justifica tu respuesta.

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","centrifuge", "storage",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) + 
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Predicción del modelo

Para obtener la predicción en este caso obtenemos las cantidades de interés, valoramos la coincidencia entre éxitos estimados y observados, y representamos la probabilidad de éxito asociada con la dosis utilizada.

¿Qué podemos comentar sobre los valores predichos por el modelo?

# Predicción
prediccion <- predict.glm(M1,type = "response", se.fit = TRUE)
# éxitos estimados
exitos <- round(prediccion$fit*caso2$total,0)
exitos
# comparativa
comparativa <- data.frame(centrifuge = caso2$centrifuge, 
                          storage = caso2$storage, 
                          original = caso2$embryogenic, 
                          predicho = exitos)
comparativa

Estimamos la curva de probabilidad de antera embriogénica. ¿Cuál es la probabilidad de una antera embriogénica para los dos tipos de almacanmiento y velocidades de centrifugación de 100, 200, o 300? ¿Que información podemos extarer de dichas probabilidades? ¿Qué opción de combinaciones de las predictoras resultará más óptimo para obtener anteras embriogénicas?

# Valores paa predicción
newdata <- data.frame(centrifuge = rep(seq(40,350,1),each=2), 
                      storage = factor(c("Control", "treatment")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, 
                      predict(M1, newdata, type = "response", se.fit=TRUE))
# Gráfico de la predicción con log dosis
ggplot(newdata, aes(x = centrifuge, y = fit, color = storage)) +
  geom_line() +
  scale_x_continuous(breaks = seq(40,360,20)) +
  scale_y_continuous(breaks = seq(0,1,0.01)) +
  labs(x = "Velocidad centrifugación", 
       y = "Probabilidad de antera embriogénica", 
       title = "Curva de probabilidad") +
  theme_bw()

---

Caso 3

---

Se realiza un ensayo clínico para determinar en un grupo de personas mayores su estado psiquiátrico. Para cada sujeto se realiza un análisis completo y se clasifica cada uno en función de si muestra rasgos de senilidad (calificados como 1) o no (calificados como 0). Por otro lado se les pasa el test de escala de inteligencia de adultos para saber si la puntuación obtenida puede ser un indicador de si la persona tiene rasgos de senilidad o no.

caso3 = read_csv("https://goo.gl/6E8fhd", col_types = "dc")
# Convertimos la respuesta en variable numérica para el ajuste de modelos 
# y como factor para representar los datos
caso3 = caso3 %>% 
  mutate(senility = 1*(senility==1),
         senilfct = fct_recode(as.factor(senility),"Senil" = "1","No senil" = "0"))

Realizamos un gráfico inicial de los datos registrados. ¿Qué opinas sobre la relación entre la puntuación del test y la calsificacón del sujeto?

ggplot(caso3) +
  geom_boxplot(aes(senilfct, score)) + 
  labs(x = "", y = "Puntuación del test") +
  theme_bw()

Estimación del modelo

¿Qué tipo de modelo consideramos para esta situación experimental? ¿El modelo obtenido es adecuado para explicar la probabilidad de caracterizar a un sujeto como senil en función del apuntuación del test? ¿Cuál es la ecuación del modelo estimada y cómo se interpretan los coeficientes de dicho modelo?

# Modelo para respuesta individualizada
M1 <- glm(senility  ~ score, family = binomial(link = logit), caso3)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

Diagnóstico del modelo

Calculamos los valores para el diagnóstico del modelo y realizamos los gráficos de diagnóstico. ¿Qué opinas sobre los gráficos de diagnóstico obtenidos? ¿Necesitamos modificar nuestro modelo de partida? Justifica tu respuesta.

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","score",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) + 
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Predicción

Dado que tenemos datos individualizados obtenemos el porcentaje de clasificación correcta. ¿Qué opinas sobre la predicción del modelo? ¿Resultinteresante utilizar la puntuación del test para clsificar a un sujeto como senil? Justifica tus respuestas.

# Obtención de predicción de la respuesta
prediccion <- predict(M1,type = "response")
# Clasificacmos a cada sujeto como éxito o fracaso
clasificado <- 1*(prediccion>=0.5)
tabla <- table(caso3$senilfct,clasificado)
tabla
# Porcentaje de calsificación correcta
round(100*sum(diag(tabla))/sum(tabla),2)

---

Caso 4

---

Se realiza un estudio para conocer que el impacto de la bomba de hiroshima en al aparición de caso de leucemia. Para ello se registro para todos los sujetos que presentaron algún tipo de cáncer el grado de radiación al que fue sometida la persona. Las variables que aparece son los conteos del número de casos de leucemia y de otros tipos de cáncer para los diferentes niveles de radiación. Es de interés en el análisis investigar la influencia del grado de radiación en la aparición de un mayor número de casos de leucemia.

caso4 = read_csv("https://goo.gl/ZDIWVC", col_types = "cdii")
# Calculamos la probabilidad de aparición de leucemia 
caso4 <- caso4 %>% mutate(probabilidad = leukemia / total)

Representamos las probabilidades obtenidas

# Construimos el orden de las categorias
orden <- c("0","1 to 9","10 to 49","50 to 99","100 to 199","200 +")
# Representamos los datos
ggplot(caso4, aes(x =radiation, y = probabilidad)) + 
  geom_point() + 
  scale_x_discrete(limits = orden) +
  labs(x = "Nivel de radiación",
       y = "Proababilidad de caso de leucemia") +
  theme_bw()

Estimación del modelo

¿Qué tipo de modelo de partida consideramos para esta situación experimental? ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de los efectos presentes en é? ¿El modelo obtenido es adecuado para explicar la probabilidad de observar un caso de leucemia en función de la radiación a la que se ha visto sometida el sujeto?

# Variable repuesta agrupada
yres <- cbind(caso4$leukemia,caso4$other)
M1 <- glm(yres  ~ radiation, family = binomial(link = logit), caso4)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)
# Tabla ANOVA asociada al factor
anova(M1,test = "Chisq")

Dado que sólo tenemos una observación por nivel del factor resulta imposible hacer un análisis adecuado para este tipo de datos ya que tenemos un completamente saturado (tantas observaciones como parámetros). Consideramos la variable numérica de marcas de clase asociada con el factor. ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de sus coeficientes? ¿El modelo obtenido es adecuado para explicar la probabilidad de observar un caso de leucemia en función de la radiación a la que se ha visto sometida el sujeto?

# Variable repuesta agrupada
caso4$radionum <- c(0,5,29.5,74.5,149.5,200)
M1 <- glm(yres  ~ radionum, family = binomial(link = logit), caso4)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

Diagnóstico del modelo

Calculamos los valores para el diagnóstico del modelo y realizamos los gráficos de diagnóstico. ¿Qué opinas sobre los gráficos de diagnóstico obtenidos? ¿Necesitamos modificar nuestro modelo de partida? Justifica tu respuesta.

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","radionum",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) + 
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Predicción del modelo

Estimamos la curva de probabilidad de consumo de drogas en función de los valores de las predictoras. ¿Cuál es la probabilidad de cosumode observar un caso de leucemia para los niveles de radiación de 50, 100, 150, y 200? ¿Que información podemos extraer de dichas probabilidades?

# Valores paa predicción
newdata <- data.frame(radionum = seq(0,200,1))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, 
                      predict(M1, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza y añadimos la dosis 
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit)
)
# Gráfico de la predicción con log dosis
ggplot(newdata, aes(x = radionum, y = fit)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
  scale_x_continuous(breaks = seq(0,200,10)) +
  scale_y_continuous(breaks = seq(0,1,0.05)) +
  labs(x = "Nivel de radiación", 
       y = "Probabilidad de observar un caso de leucemia", 
       title = "Curva de probabilidad") +
  theme_bw()

---

Caso 5

---

Los datos siguientes describen los patrones de comportamiento en el consumo de drogas psicotrópicas en una muestra de individuos del Oeste de Londres. Los investigadores se plantean las preguntas siguientes:

• ¿Hay diferencias por ’sexo’ en el consumo de drogas?
• ¿Cómo influye la ’edad’ para explicar el consumo de psicotrópicos?
• ¿La edad influye igual en hombres y en mujeres?

Obtén la expresión y el valor de las predicciones sobre el consumo de psicotrópicos en hombres y en mujeres de 17 y 52 años con el modelo ajustado.

sexo <- c(rep("H",4),rep("M",4))
edad <- c("16-29","30-44","45-64","65-74","16-29","30-44","45-64","65-74")
usa <- c(21,32,70,43,46,89,169,51)
nousa <- c(683,596,705,295,738,700,847,196)
caso5 = data.frame(sexo, edad, usa, nousa)
# Calculamos la probabilidad de consumo de drogas para cada combinación
caso5 <- caso5 %>% mutate(total = usa + nousa, 
                          probabilidad = usa /total)

Representamos los datos observados. A la vista del gráfico ¿puedes responder a alguna de las preguntas planteadas en la investigación?

ggplot(caso5,aes(x = edad, y = probabilidad, color = sexo)) + 
  geom_point() +
  labs(x = "Edad",
       y = "Probabilidad de uso de drogas") +
  theme_bw()

Estimación del modelo

Para la estimación del modelo convertimos la variable ordinal (edad) en numérica para ajustar el modelo, dado que si no tendríamos un modelo saturado (mismas observaciones que combinaciones de los factores).

caso5$edadnum <- c(22.5,37,54.5,69.5)

Ajustamos el modelo con la nueva variable. ¿Qué tipo de modelo de partida consideramos para esta situación experimental? ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de los efectos presentes en él? ¿El modelo obtenido es adecuado para explicar la probabilidad de conusmo de drgogas? ¿Cómo respondemos a las preguntas planteadas en la investigación?

# Variable repuesta agrupada
yres <- cbind(caso5$usa,caso5$nousa)
M1 <- glm(yres  ~ edadnum * sexo, family = binomial(link = logit), caso5)
# Selección del modelo
M1 <- step(M1, test = "Chisq")
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)
# Tabla ANOVA asociada al factor
anova(M1, test = "Chisq")

Diagnóstico del modelo

Calculamos los valores para el diagnóstico del modelo y realizamos los gráficos de diagnóstico. ¿Qué opinas sobre los gráficos de diagnóstico obtenidos? ¿Necesitamos modificar nuestro modelo de partida? Justifica tu respuesta.

# Valores para el diagnóstico
diagnostico <- fortify(M1)
# Valores predichos y residuos en la escala original 
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico
# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","edadnum", ".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable)) + 
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Predicción del modelo

Estimamos la curva de probabilidad del consumo de psicotrópicos en hombres y en mujeres de 17 y 52 años con el modelo ajustado

# Valores paa predicción
newdata <- data.frame(edadnum = rep(seq(16,74,1),each=2), 
                      sexo = factor(c("H", "M")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, 
                      predict(M1, newdata, type = "response", se.fit=TRUE))
# Gráfico de la predicción con log dosis
ggplot(newdata, aes(x = edadnum, y = fit, color = sexo)) +
  geom_line() +
  scale_x_continuous(breaks = seq(16,74,2)) +
  scale_y_continuous(breaks = seq(0,1,0.01)) +
  labs(x = "Edad", 
       y = "Probabilidad conusmo de drogas", 
       title = "Curva de probabilidad") +
  theme_bw()

---

Caso 6

---

Se desea estudiar la sensibilidad de un test basado en el diagnóstico de la tuberculosis a través de una prueba basada en rayos X . Los datos registrados aparecen a continuación.

tuberculosis <- c(22,8)
total <- c(51,1739)
rayosx <- c("Positivo","Negativo")
caso6 = data.frame(tuberculosis, total, rayosx)
caso6 <- caso6 %>% mutate(notuberculosis = total - tuberculosis)

Estimación del modelo

Ajustamos el modelo correspondiente a la información experimental. ¿Cuál es la sensibilidad del test a partir del modelo estimado?¿y la especificidad?

# Variable repuesta agrupada
yres <- cbind(caso6$tuberculosis, caso6$notuberculosis)
M1 <- glm(yres  ~ rayosx, family = binomial(link = logit), caso6)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

---

Caso 7

---

En el Hospital de Yale-New Heaven, en Connecticut, se llevó a cabo un estudio para investigar la relación entre los nacimientos prematuros (el niño nazca antes de 37 semanas de gestación o su peso sea inferior a 2500 g.) y la edad de la madre. La población de estudio consistió en 175 madres de niños nacidos únicos y prematuros, y 303 madres de niños no prematuros. Los datos se presentan a continuación. ¿Hay alguna relación entre la edad de la madre y el hecho de que un niño nazca prematuro? ¿Cuál es el grupo de edad con mayor riesgo?

edad <- c("14-17","18-19","20-24","25-29","+30")
casos <- c(15,22,47,56,35)
controles <- c(16,25,62,122,18)
caso7 = data.frame(edad, casos, controles)
# Calculamos varaible numérica asociada a edad
caso7 <-caso7 %>% mutate(edadnum = c(15.5, 18.5, 22, 27, 30))

Estimación del modelo

Ajustamos el modelo correspondiente a la información experimental. ¿Qué tipo de modelo de partida consideramos para esta situación experimental? ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de los efectos presentes en él? ¿El modelo obtenido es adecuado para explicar la probabilidad de conusmo de drgogas? ¿Cómo respondemos a las preguntas planteadas en la investigación?

# Variable repuesta agrupada
yres <- cbind(caso7$casos, caso7$controles)
M1 <- glm(yres  ~ edadnum, family = binomial(link = logit), caso7)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)

---

Caso 8

---

En el banco de datos siguiente se presentan los resultados de una encuesta realizada en 1998. A cada sujeto de una muestra de 300 adultos se le pidió que opinara sobre qué política consideraba adecuada implantar respecto al uso de tabaco en lugares públicos. Las opciones planteadas son

• Opción 1: Sin restricciones
• Opción 2: Fumar sólo en áreas exclusivas
• Opción 3: No fumar nunca
• Opción 4: No opina

¿Hay alguna relación entre la actitud frente al tabaco y el nivel de estudios?

nivel <- c("Est. Superiores", "Secundaria", "Primaria")
opt1 <- c(5,15,15)
opt2 <- c(44,100,40)
opt3 <- c(23,30,10)
opt4 <- c(3,5,10)
total <- c(75,150,75)
caso8 = data.frame(nivel, opt1, opt2, opt3, opt4, total)

Nos interesamos únicamente por la opción 3 de no fumar nunca. Calculamos la respuesta de acuerdo a esta opción. Más tarde estudiaremos los modelos logit para trabajar con todas las opciones de respuesta a la vez.

caso8 <- caso8 %>% mutate(nofumar = opt4, fumar = opt1 + opt2)

Estimación del modelo

Ajustamos el modelo correspondiente a la información experimental. ¿Qué tipo de modelo de partida consideramos para esta situación experimental? ¿Cuál es la ecuación del modelo estimada y cómo se interpreta cada uno de los efectos presentes en él? ¿Cómo respondemos a las preguntas planteadas en la investigación?

# Variable repuesta agrupada
yres <- cbind(caso8$nofumar, caso8$fumar)
M1 <- glm(yres  ~ nivel, family = binomial(link = logit), caso8)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)
# Resumen del modelo ajustado
summary(M1)
# Tabla ANOVA de efecto
anova(M1, test="Chisq")

Realiza el diagnóstico de este modelo. En caso de detectar algún problema, modifica el modelo y interpreta los resultados nuevamente. ¿Qué ocurre si creamos una variable ordinal para representar el nivel de estudios y ajustamos ese nuveo modelo? ¿Cómo respondemos a las preguntas planteadas en la investigación?

# Creamos vector numérico
caso8 <- caso8 %>% mutate(nivelnum = c(3,2,1))
# Variable repuesta agrupada
yres <- cbind(caso8$nofumar, caso8$fumar)
M1 <- glm(yres  ~ nivelnum, family = binomial(link = logit), caso8)
# Bondad de ajuste
1-pchisq(M1$deviance,M1$df.residual)

## [1] 0.1146155

# Resumen del modelo ajustado
summary(M1)

## 
## Call:
## glm(formula = yres ~ nivelnum, family = binomial(link = logit), 
##     data = caso8)
## 
## Deviance Residuals: 
##      1       2       3  
##  1.007  -1.106   0.502  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -1.0959     0.6721  -1.631   0.1030  
## nivelnum     -0.7850     0.3816  -2.057   0.0397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7.0732  on 2  degrees of freedom
## Residual deviance: 2.4894  on 1  degrees of freedom
## AIC: 16.851
## 
## Number of Fisher Scoring iterations: 4

# Tabla ANOVA de efecto
anova(M1, test="Chisq")

Realiza e interpreta el diagnóstico para este modelo

---

Bibliografía

---

---

Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.