Trabajo 1: Modelos RLS
Prueba
Estudio de caso: Modelos RLS
Los datos proporcionan las estimaciones del número medio de dientes cariados, perdidos y empastados (DMFT) a la edad de 12 años, y el consumo medio anual de azúcar en los cinco años anteriores para 90 países. El estudio se plantea diferentes objetivos pero el principal es conocer si el econusmo medio de azúcar influye en el DMFT. Las variables consideradas son:
- Country: País.
- Indus: si el país se considera un país industrializado; un factor con niveles Ind (industrializado) o NonInd (no industrializado).
- Sugar: el consumo medio anual de azúcar en kilogramos por persona y año, calculado a lo largo de los cinco años (o tanto como esté disponible) anteriores a la encuesta.
- DMFT: estimaciones del número medio de dientes cariados, perdidos y empastados a los 12 años.
Fuente: M. Woodward and A. R. P. Walker (1994) Sugar consumption and dental caries: evidence from 90 countries. British Dental Journal, 176, 297–302.
Cargamos los datos correspondientes:
Para este banco de datos se deben responder las cuestiones siguientes:
- ¿Cuál es el objetivo principal del estudio? ¿Se pueden establecer objetivos secundarios que puedan resultar relevantes? ¿Qué variable está asociada con el objetivo principal del estudio?
Parte 1.
- Representa gráficamente la posible relación entre el consumo de azúcar y DMFT, y extrae las conclusiones correspondientes.
- Ajusta un modelo adecuado (atendiendo al gráfico anterior) y comenta el modelo obtenido y su capacidad explicativa.
- Comprueba si el modelo considerado verifica las hipótesis del modelo de RLS.
- ¿Cuál será la predicción de DMFT para valores de Sugar igual a 15, 30 y 45? ¿Cómo podemos interpretar esas predicciones?
Parte 2.
- Representa gráficamente la posible relación entre el consumo de azúcar y DMFT identificando por la industrialización del país, y extrae las concluiones correspondientes.
- Ajusta un modelo adecuado (atendiendo al gráfico anterior) para cada grupo de países y comenta los modelos obtenidos y su capacidades explicativas.
- Comprueba si los modelos considerados verifican las hipótesis del modelo de RLS.
- ¿Cuál será la predicción de DMFT para valores de Sugar igual a 15, 30 y 45 en cada uno de los modelos? ¿Cómo podemos intereptar esas predicciones?
Objetivos
El objetivo principal del estudio es determinar si la estimación del número medio de dientes cariados, perdidos y empastados a los 12 años (DMFT) se puede relacionar con el consumo medio anuka de azúacar (Sugar) en un conjunto de países.
Como objetivo secundario nos podríamos plantear si la posible relación entre DMFT y Sugar puede verse influenciada porque el país sea industrializado o no, lo que está intimamente relacionado con los hábitos higiénicos y el nivel de salud de cada uno de los países considerados.
La variable que está vinculada con el objetivo principal del estudio es DMFT ya que es la que tomamos como variable respuesta.
A continuación se muestran los resultados obtenidos con los diferentes análisis realziados para este banco de datos, mientras que en el apéndice se muestra el código empleado.
Análisis gráfico del modelo completo
Para el análisis de la posible relación entre DMFT y Sugar realizamos en primer lugar el correspondiente gráfico de dispersión (al tratarse ambas variables de tipo numérico).
A la vista del gráfico podemos concluir que parece existir cierta dependencia entre DMFT y Sugar, aunque como se puede apreciar claramente la dispersión del número estimado de piezas dentales con problemas es más variable conforme aumenta el consumo medio anual de azúcar. Además, la posible tendencia entre DMFT y Sugar parece más la correspondiente a un modelo exponencial que a un modelo lineal. Nos planteamos por tanto dos posibles modelos iniciales en el tratamaiento de estos datos:
\[M1: DMFT = \beta_0 + \beta_1* Sugar + \epsilon\] \[M2: log(DMFT) = \beta_0 + \beta_1* Sugar + \epsilon\]
Procedemos con el ajuste de ambos modelos y comparamos la capacidad explicativa de cada uno de ellos. Comenzamos con los resultados obtneidos apra el modelo 1:
| Parameter | Coefficient | Pr(>|t|) |
|---|---|---|
| (Intercept) | 1.3 (0.69, 1.91) | < 0.001 |
| Sugar | 0.05 (0.03, 0.06) | < 0.001 |
El modelo de crecimiento lineal entre DMFT y Sugar obtenido viene dado por:
\[\widehat{DMFT} = 1.3 + 0.05*Sugar\] indicando que el número estimado de piezas dentales con problemas es aproximadamante 1 , intervalo de confianza entre 0 y 2, caundo el consumo promedio de azúcar es de 0 kg, mientras que dicho número crece en un unidad cuando aumentamos en 20 kg el consumo medio de azúcar al año. Por tanto, la diferencia estimada del número de piezas dentales entre dos países con una diferencia de consumo medio anual de azúcar de 60 kilos es de 4 puezas dentales. Con respecto a la capacidad explicativa del modelo tenemos que:
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.225 | 0.217 | 1.4 | 25.6 | 2.28e-06 | 1 | -157 | 319 | 327 | 171 | 88 | 90 |
donde podemos ver que el modelo resulta significativo (test F de la regresión con p-valor \(< 0.05\)), indicando que hay asociación entre DMFT y el consumo de azúcar, pero sin embargo el \(R^2\) asociado es del 22% lo que implica que aunque el modelo es significativo la capacidad explicativa es muy baja, debido principalmente a la gran variabilidad observada entre los diferentes países.
Pasamos al análisis del modelo con crecimiento exponencial.
| Parameter | Coefficient | Pr(>|t|) |
|---|---|---|
| (Intercept) | 0.15 (-0.09, 0.4) | 0.222 |
| Sugar | 0.02 (0.01, 0.03) | < 0.001 |
En este caso el modelo obtenido viene dado por:
\[log(DMFT) = 0.15 + 0.02*Sugar\]
donde podemos ver que el número de piezas dentales con problemas con un consumo de azúcar de 0 kg es igual a \(exp(0.15)\), mientras que el número se incrementa a razón de \(exp(0.02)\) por cada kilo que aumentamos nuestro consumo de azúcar. En cuanto a la calidad del modelo tenemos que:
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.285 | 0.276 | 0.56 | 35 | 6.21e-08 | 1 | -74.5 | 155 | 163 | 27.6 | 88 | 90 |
indicando que el modelo es significativo (p-valor \(< 0.05\)) y con una capacidad explicativa de 28%, es decir, tenemos una ganancia de un 6% con respecto al modelo de crecimiento lineal. En esta situación parece preferible el modelo de crecimiento exponencial pero nos faltaría verificar si se cumplen las hipótesis del modelo, con repecto a los residuos otenidos en cada uno de esos modelos: linealidad, normalidad, varianza constante, e independencia. Para verificar las hipótesis utlizamos un gráfico de dispersión de los residuos versus valores ajustados (linealidad), el test de Shapiro-Wilks (normalidad), test de Breush-Pagan (varianza constante), y el test de Durbin-Watson (independencia). En la tabla siguiente se muestran los p-valores asociados a cada uno de los tests considerados para ambos modelos.
| Modelo/Hipótesis | Varianza constante | Normalidad | Independencia |
|---|---|---|---|
| Crecimiento lineal | 0.1358 | 3.038e-05 | 0.2703 |
| Crecimiento exponencial | 0.07945 | 0.2317 | 0.1748 |
Como se puede ver el modelo de crecimiento lineal no verifica la hipótesis de Normalidad (p-valor \(< 0.05\)) indicando el incumplimiento de dicha hipótesis e invalidando el modelo obtenido. Por otro lado el modelo de crecimiento exponencial si verifica todas las hipótesis, incluyendo la de linealidad como se ve en el gráfico siguiente:
Podriamos tratar de utilizar Box-Cox sobre el modelo de crecimiento lienal para tratar de estabilizar la varianza pero dado que el modelo de crecimiento exponencial proporciona mejores resultados, vamos a ver únicamente cual es la transformación que obtendríamos.
Como se puede ver la trasnformación de raíz cuadrada (\(lambda = 0.5\)) se encuentra dentro del intervalo de confianza de las posibles transformaciones, aunque como se puede comprobar su calidad es similar al modelo de crecimeinto exponencial. Nos quedamos con el modelo de crecimeinto exponencial.
Analizamos ahora la predicción del número de piezas dentales con problemas para conusmos de azúcar de 10, 30 y 45. Obtenemos la predicción del modelo y finalmente expresamos en la escala de la varaible original. En términos del modelo obtenemos:
| Sugar | fit | lwr | upr |
|---|---|---|---|
| 10 | 0.36 | 0.18 | 0.55 |
| 30 | 0.79 | 0.67 | 0.9 |
| 45 | 1.1 | 0.95 | 1.26 |
y en la escala original
| Sugar | fit | lwr | upr |
|---|---|---|---|
| 10 | 1.44 | 1.19 | 1.73 |
| 30 | 2.19 | 1.95 | 2.47 |
| 45 | 3.01 | 2.57 | 3.53 |
lo que indica que el número de piezas dentales con defectos para el conjunto de países con un consumo de medio anual de azúcar de 10 kg se sítuan entre 1 y 2, mientras que los que consumen 30 kg están entre 2 y 3, y los de 45 kg entre 3 y 4. Se puede ver que a pesar de usar un modelo de crecimiento exponencial los cambios observados no son muy grandes debido en gran medida a que el crecimiento observado es muy suave (como se puede ver en el gráfico de predicción a contiuación).
Para el segundo grupo de análisis vamos a tener en cuenta la varaible que indica si el país está industrializado o no. Para ello realizamos en primer lugar un gráfico de dispersión donde identificamos cada tipo de país para verificar si tienen comportamientos similares o diferentes.
Como se aprecía claramente los países menos industrializados con consumos de azúcar en el promedio anual más bajos tienen sin embargo más problemas dentales, debido segurmente a las condiciones y hábitos higiénicos en dichos paises. Por otro lado, los países más industrializados tienen unos hábitos alimenticios que les hacen consumir mucho más azúcar, y que acarrean un mayor número de problemas dentales a pesar de unos hábitos higiénicos más adecuados.
Los paises no industrilizados parecen seguir el mismo patrón de comportamiento que observabamos cuando trabajamos con todos los datos, mientras que los industrializados parecen tener un comportmaiento más aleatorio, no mostrando ningún patrón claro.
En esta situación optamos por mantener el modelo de crecimiento exponencial para ambos grupos de paises ya que parece que ese comportamiento se podría estar reproduciendo. Dividimos el conjunto de datos en los dos grupos y ajustamos los modelos de crecimiento en cada caso.
A continuación se presentan los resultados para los datos de los países industrializados:
| Parameter | Coefficient | Pr(>|t|) |
|---|---|---|
| (Intercept) | 1.39 (0.64, 2.13) | < 0.001 |
| Sugar | -0.01 (-0.02, 0.01) | 0.495 |
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.0174 | -0.019 | 0.388 | 0.479 | 0.495 | 1 | -12.6 | 31.3 | 35.4 | 4.06 | 27 | 29 |
El modelo obtenido viene dado por:
\[\widehat{log(DMFT)} = 1.39 - 0.01 *Sugar\]
pero con una capacidad explicativa del 1% y donde podemos ver que el coeficiente asociado con Sugar es considerado igual a cero (test F de la regresión no significativo, pvalor \(> 0.05\)). Esto implica que para este conjunto de datos no se puede establecer un modelo entre consumo de azúcar y número estimado de piezas dentales con problemas. Parece pues que los defectos observados con estos países se debe atribuir a otro tipo de factores o efectos que no han sido considerados en el estudio.
Pasamos ahora al modelo para los países no industrializados.
| Parameter | Coefficient | Pr(>|t|) |
|---|---|---|
| (Intercept) | 0.1 (-0.18, 0.37) | 0.494 |
| Sugar | 0.02 (0.01, 0.03) | < 0.001 |
| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.261 | 0.248 | 0.604 | 20.8 | 2.64e-05 | 1 | -54.8 | 116 | 122 | 21.5 | 59 | 61 |
En este caso el modelo obtenido es
\[\widehat{log(DMFT)} = 0.1 + 0.02 *Sugar\]
con una capacidad explicativa del 26% y donde podemos ver que el coeficiente asociado con Sugar es considerado distinto de cero (test F de la regresión no significativo, pvalor \(< 0.05\)). Por tanto, para los países no industrializados si podemos establecer un modelo de crecimiento exponencial muy similar al del modelo para todos los datos.
Los tests de diagnóstico asociados con este modelo resultan no significativos indicando que e modelo verifica todas las hipótesis necesarias, y que por tanto es válido para el proceso de predicción, a pesar de su baja capacidad explicativa.
Los valroes de predicción (que se pueden obtener con un código simialr al del modleo completo) proporcionan resultados muy simialares a los del modelo completo.
Apéndice
Gráfico del modelo completo
ggplot(datos, aes(x = Sugar, y = DMFT)) +
geom_point() +
labs(x = "Consumo medio de azúcar (en Kg)",
y = "Número medio de dientes con problemas a los 12 años")Ajustes iniciales
# Modelo 1
fit1 <- lm(DMFT ~ Sugar, data = datos)
glm_coef(fit1)
glance(fit1)
# Modelo 2
datos$lDMFT <- log(datos$DMFT)
fit2 <- lm(lDMFT ~ Sugar, data = datos)
glm_coef(fit2)
glance(fit2)# Diagnóstico 1
dg1 <-fortify(fit1)
# Test de diagnóstico
bptest(fit1)
# Test de diagnóstico
shapiro.test(dg1$.stdresid)
# Test de diagnóstico
dwtest(fit1, alternative = "two.sided")
# Diagnóstico 2
dg2 <-fortify(fit2)
# Test de diagnóstico
bptest(fit2)
# Test de diagnóstico
shapiro.test(dg2$.stdresid)
# Test de diagnóstico
dwtest(fit2, alternative = "two.sided")dg2 <-fortify(fit2)
ggplot(dg2, aes(x = .fitted, y = .stdresid)) +
geom_point() +
stat_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula 'y ~ x'newpred <- data.frame(Sugar = c(10, 30, 45))
# Predicción para la media de la respuesta
# Opción interval = "confidence"
newdata <- data.frame(newpred,
predict(fit2,
newpred,
interval = "confidence"))
round(newdata, 2)# Gráfico del ajuste
plot_model(fit2, "pred", terms = ~Sugar,
ci.lvl = 0.95,
show.data = TRUE,
axis.title = c("Consumo medio de azúcar (en Kg)", "Número medio de dientes con problemas a los 12 años"),
title = " ")# Predicción para la media de la respuesta
newdata <- data.frame(Sugar = seq(min(datos$Sugar),
max(datos$Sugar),
length = 50))
newdata <- data.frame(newdata, predict(fit2,
newdata,
interval = "confidence"))
# Deshacemos la transformación para volver a la escala de viscosidad
newdata[,2:4] <- exp(newdata[,2:4])
# Gráfico del ajuste
ggplot(newdata, aes(x = Sugar, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
geom_point(data = datos, aes(x = Sugar, y = DMFT)) +
labs(x = "Consumo medio de azúcar (en Kg)",
y = "Número medio de dientes con problemas a los 12 años")Modelo identificando tipo de país
ggplot(datos, aes(x = Sugar, y = DMFT, color = Indus)) +
geom_point() +
labs(x = "Consumo medio de azúcar (en Kg)",
y = "Número medio de dientes con problemas a los 12 años",
col = "Tipo de País")datos1 <- datos[datos$Indus == "Ind",]
datos2 <- datos[datos$Indus == "NonInd",]
fit1 <- lm(lDMFT ~ Sugar, data = datos1)
glm_coef(fit1)
glance(fit1)
fit2 <- lm(lDMFT ~ Sugar, data = datos2)
glm_coef(fit2)
glance(fit2)