Evaluacion
JOHN JAIRO PRADO
2024-11-02
Instrucciones:
• Si requiere ayuda o algún tipo de aclaración levante la mano y el docente vendrá pronto.
• Tenga en cuenta que la posesión o el uso de teléfonos móviles o cualquier otro dispositivo electrónico no autorizado en el salón del examen esta estrictamente prohibido.
• Conclusiones sin tablas de resultados que justifiquen las respuestas tienen nota de 0.
1. La base de datos danishlc contiene el número de incidentes de cancer de pulmon de 1968 a 1971 en cuatro ciudades danesas, con el registrado por grupo de edad y el numero de sujetos en cada grupo de edad.
Descripción del df
Cáncer de Pulmón Danés: Este conjunto de datos contiene información sobre el número de casos de cáncer de pulmón en cuatro ciudades de Dinamarca.
El df consta de 24 observaciones y 4 variables:
Cases: Representa el número de casos de cáncer de pulmón; es un vector numérico.
Pop: Indica la población de cada grupo de edad en cada ciudad; también es un vector numérico.
Age: Este campo clasifica los grupos de edad, siendo un factor con los niveles: 40-54, 55-59, 60-64, 65-69, 70-74 y >74.
City: Muestra las ciudades incluidas en el estudio, siendo un factor con los niveles: Fredericia, Horsens, Kolding y Vejle.
(a) Ajuste el modelo de regresion que permita estimar el numero de sujetos con cancer de pulmon. Muestre las meidas de bondad de ajuste.
#install.packages('betareg')
#install.packages("AER")
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(GLMsData))
suppressPackageStartupMessages(library(car))
suppressPackageStartupMessages(library(MASS, exclude = "select"))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(AER))
suppressPackageStartupMessages(library(betareg))data("danishlc")
dim(danishlc)[1] 24 4glimpse(danishlc)Rows: 24
Columns: 4
$ Cases <int> 11, 11, 11, 10, 11, 10, 13, 6, 15, 10, 12, 2, 4, 8, 7, 11, 9, 12…
$ Pop <int> 3059, 800, 710, 581, 509, 605, 2879, 1083, 923, 834, 634, 782, 3…
$ Age <fct> 40-54, 55-59, 60-64, 65-69, 70-74, >74, 40-54, 55-59, 60-64, 65-…
$ City <fct> Fredericia, Fredericia, Fredericia, Fredericia, Fredericia, Fred…hist(danishlc$Cases)
# Crear la variable dummy
danishlc$Age_40_54 <- ifelse(danishlc$Age == "40-54", 1, 0)
danishlc$Age_55_59 <- ifelse(danishlc$Age == "55-59", 1, 0)
danishlc$Age_60_64 <- ifelse(danishlc$Age == "60-64", 1, 0)
danishlc$Age_65_69 <- ifelse(danishlc$Age == "65-69", 1, 0)
danishlc$Age_65_69 <- ifelse(danishlc$Age == "70-74", 1, 0)
# Crear la variable dummy
danishlc$City_Fredericia <- ifelse(danishlc$City == "Fredericia", 1, 0)
danishlc$City_Horsens <- ifelse(danishlc$City == "Horsens", 1, 0)
danishlc$City_Kolding <- ifelse(danishlc$City == "Kolding", 1, 0)
# Eliminar columnas asignando a NULL
danishlc$Age <- NULL
danishlc$City <- NULL
# Ver el resultado
head(danishlc) Cases Pop Age_40_54 Age_55_59 Age_60_64 Age_65_69 City_Fredericia
1 11 3059 1 0 0 0 1
2 11 800 0 1 0 0 1
3 11 710 0 0 1 0 1
4 10 581 0 0 0 0 1
5 11 509 0 0 0 1 1
6 10 605 0 0 0 0 1
City_Horsens City_Kolding
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0library(lmtest)
modelo_danishlc = glm(Cases ~ ., family = poisson(link = "log"), data = danishlc)barplot(danishlc$Cases)
No hay una distribucion clara.
dispersiontest(modelo_danishlc,trafo=1)
Overdispersion test
data: modelo_danishlc
z = -0.458, p-value = 0.68
alternative hypothesis: true alpha is greater than 0
sample estimates:
alpha
-0.14008 Ho: No hay sobredispersión Ha: Hay sobredispersión
No se rechaza Ho, No hay evidencia suficiente para afirmar que hay sobredispersión en los datos.
Seleccion del Mejor Subconjunto de Variables
step(modelo_danishlc)Start: AIC=137.51
Cases ~ Pop + Age_40_54 + Age_55_59 + Age_60_64 + Age_65_69 +
City_Fredericia + City_Horsens + City_Kolding
Df Deviance AIC
- Age_65_69 1 23.1 136
- City_Kolding 1 23.1 136
- Age_40_54 1 23.2 136
- Age_55_59 1 23.2 136
- Pop 1 23.3 136
- Age_60_64 1 23.7 136
- City_Horsens 1 23.7 136
- City_Fredericia 1 24.7 137
<none> 23.1 138
Step: AIC=135.52
Cases ~ Pop + Age_40_54 + Age_55_59 + Age_60_64 + City_Fredericia +
City_Horsens + City_Kolding
Df Deviance AIC
- City_Kolding 1 23.2 134
- Age_55_59 1 23.2 134
- Age_40_54 1 23.2 134
- Pop 1 23.3 134
- Age_60_64 1 23.7 134
- City_Horsens 1 23.8 134
- City_Fredericia 1 24.7 135
<none> 23.1 136
Step: AIC=133.56
Cases ~ Pop + Age_40_54 + Age_55_59 + Age_60_64 + City_Fredericia +
City_Horsens
Df Deviance AIC
- Age_40_54 1 23.2 132
- Pop 1 23.3 132
- Age_55_59 1 23.4 132
- Age_60_64 1 23.7 132
- City_Horsens 1 23.9 132
- City_Fredericia 1 24.9 133
<none> 23.2 134
Step: AIC=131.64
Cases ~ Pop + Age_55_59 + Age_60_64 + City_Fredericia + City_Horsens
Df Deviance AIC
- Age_60_64 1 23.7 130
- City_Horsens 1 23.9 130
- Pop 1 24.0 130
- Age_55_59 1 24.0 130
- City_Fredericia 1 25.2 132
<none> 23.2 132
Step: AIC=130.09
Cases ~ Pop + Age_55_59 + City_Fredericia + City_Horsens
Df Deviance AIC
- City_Horsens 1 24.4 129
- Pop 1 24.7 129
- Age_55_59 1 24.8 129
- City_Fredericia 1 25.6 130
<none> 23.7 130
Step: AIC=128.78
Cases ~ Pop + Age_55_59 + City_Fredericia
Df Deviance AIC
- Pop 1 25.3 128
- Age_55_59 1 25.5 128
- City_Fredericia 1 25.8 128
<none> 24.4 129
Step: AIC=127.66
Cases ~ Age_55_59 + City_Fredericia
Df Deviance AIC
- Age_55_59 1 26.2 127
- City_Fredericia 1 26.8 127
<none> 25.3 128
Step: AIC=126.61
Cases ~ City_Fredericia
Df Deviance AIC
- City_Fredericia 1 27.7 126
<none> 26.2 127
Step: AIC=126.09
Cases ~ 1
Call: glm(formula = Cases ~ 1, family = poisson(link = "log"), data = danishlc)
Coefficients:
(Intercept)
2.23
Degrees of Freedom: 23 Total (i.e. Null); 23 Residual
Null Deviance: 27.7
Residual Deviance: 27.7 AIC: 126mejor_modelo_danishlc = glm(Cases ~ 1, family = poisson(link = "log"), data = danishlc)
summary(mejor_modelo_danishlc)
Call:
glm(formula = Cases ~ 1, family = poisson(link = "log"), data = danishlc)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.2336 0.0668 33.4 <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: 27.704 on 23 degrees of freedom
Residual deviance: 27.704 on 23 degrees of freedom
AIC: 126.1
Number of Fisher Scoring iterations: 4Medidas de Bondad de Ajuste
Pseudo R^2
library(DescTools)
Adjuntando el paquete: 'DescTools'The following object is masked from 'package:car':
Recode# Calcular el Pseudo R² de Nagelkerke
pseudo_r2_danishlc <- PseudoR2(modelo_danishlc, 'Nagelkerke')
# Mostrar el resultado
cat("El Pseudo R² de Nagelkerke es:", pseudo_r2_danishlc, "\n")El Pseudo R² de Nagelkerke es: 0.17481 Conclusiones
El Pseudo R² de Nagelkerke de 0.17481 indica que aproximadamente el 17.5% de la variabilidad en los casos de cáncer de pulmón se puede atribuir a las variables incluidas en este modelo.
summary(mejor_modelo_danishlc)
Call:
glm(formula = Cases ~ 1, family = poisson(link = "log"), data = danishlc)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.2336 0.0668 33.4 <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: 27.704 on 23 degrees of freedom
Residual deviance: 27.704 on 23 degrees of freedom
AIC: 126.1
Number of Fisher Scoring iterations: 4Null deviance: 27.704 on 23 degrees of freedom Residual deviance: 27.704 on 23 degrees of freedom AIC: 126.1
deviance_modelo_danishlc = 1-pchisq(27.704-27.704, 23-23); deviance_modelo_danishlc [1] 1# Mostrar el resultado
cat("El modelo deviance es:", deviance_modelo_danishlc, "\n")El modelo deviance es: 1 Conclusión
Una devianza de 1 en el modelo indica que el ajuste es bastante satisfactorio. Esto sugiere que la diferencia entre los datos observados y los esperados según el modelo es relativamente pequeña.
Prueba de Coeficientes Estimados
Modelo Inicial
summary(modelo_danishlc)
Call:
glm(formula = Cases ~ ., family = poisson(link = "log"), data = danishlc)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.317910 0.451538 5.13 2.8e-07 ***
Pop -0.000272 0.000714 -0.38 0.70
Age_40_54 0.465589 1.599855 0.29 0.77
Age_55_59 -0.097780 0.286459 -0.34 0.73
Age_60_64 0.167364 0.222573 0.75 0.45
Age_65_69 0.018094 0.213884 0.08 0.93
City_Fredericia 0.234758 0.188900 1.24 0.21
City_Horsens 0.177419 0.231400 0.77 0.44
City_Kolding 0.039853 0.224317 0.18 0.86
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 27.704 on 23 degrees of freedom
Residual deviance: 23.122 on 15 degrees of freedom
AIC: 137.5
Number of Fisher Scoring iterations: 4# Coeficientes estimados del modelo sin el intercepto
coefficients_modelo_danishlc <- c(
Pop = -0.000272,
Age_40_54 = 0.465589,
Age_55_59 = -0.097780,
Age_60_64 = 0.167364,
Age_65_69 = 0.018094,
City_Fredericia = 0.234758,
City_Horsens = 0.177419,
City_Kolding = 0.039853
)
# Calcular los coeficientes transformados
transf_coef_modelo_danishlc <- ifelse(coefficients_modelo_danishlc > 0,
exp(coefficients_modelo_danishlc),
1 / exp(-coefficients_modelo_danishlc))
# Mostrar los resultados
data.frame(coefficients_modelo_danishlc = coefficients_modelo_danishlc, Transformed = transf_coef_modelo_danishlc) coefficients_modelo_danishlc Transformed
Pop -0.000272 0.99973
Age_40_54 0.465589 1.59295
Age_55_59 -0.097780 0.90685
Age_60_64 0.167364 1.18218
Age_65_69 0.018094 1.01826
City_Fredericia 0.234758 1.26460
City_Horsens 0.177419 1.19413
City_Kolding 0.039853 1.04066Conclusiones
Pop: Un cambio unitario en la población está asociado con una reducción en la tasa de casos de cáncer de pulmón, con un coeficiente transformado de 0.99973. Esto indica que, en promedio, un incremento en la población no tiene un efecto significativo sobre la tasa de incidencia.
Age_40_54: Este grupo de edad muestra un coeficiente transformado de 1.59295, lo que sugiere que los individuos de este grupo tienen aproximadamente un 59.3% más de probabilidad de desarrollar cáncer de pulmón en comparación con el grupo base.
Age_55_59: Con un coeficiente transformado de 0.90685, la probabilidad de cáncer de pulmón disminuye en este grupo, lo que indica un 9.3% menos de probabilidad en comparación con el grupo base.
Age_60_64: Un coeficiente transformado de 1.18218 sugiere que este grupo tiene un 18.2% más de probabilidad de desarrollar cáncer de pulmón en comparación con el grupo base.
Age_65_69: Con un coeficiente transformado de 1.01826, este grupo muestra un leve aumento en la probabilidad, que se traduce en una variación mínima respecto al grupo base.
City_Fredericia: Este coeficiente transformado de 1.26460 sugiere que los residentes de Fredericia tienen un 26.5% más de probabilidad de tener cáncer de pulmón en comparación con la ciudad de referencia.
City_Horsens: Un coeficiente de 1.19413 indica un aumento del 19.4% en la probabilidad de incidencia en esta ciudad.
City_Kolding: Finalmente, con un coeficiente transformado de 1.04066, los habitantes de Kolding tienen un leve aumento en la probabilidad de cáncer de pulmón, del 4.1%.
Mejor Modelo
summary(mejor_modelo_danishlc)
Call:
glm(formula = Cases ~ 1, family = poisson(link = "log"), data = danishlc)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.2336 0.0668 33.4 <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: 27.704 on 23 degrees of freedom
Residual deviance: 27.704 on 23 degrees of freedom
AIC: 126.1
Number of Fisher Scoring iterations: 4Conclusión
El modelo final solo incluye el intercepto, lo que significa que se basa únicamente en la media general de los casos de cáncer de pulmón, sin considerar ninguna influencia de otras variables. En este contexto, la devianza nula y la devianza residual son iguales, lo que indica que el modelo no ha logrado mejorar el ajuste al no incluir variables explicativas. Esto sugiere que el modelo inicial, es efectivo como un modelo más complejo.
(b) Estime el número medio hombres con cáncer de pulmón en la ciudad de Vejle en un grupo de 1000 hombres de 40 años.
# Valores para la estimación
pop_vejle <- 1000 # población de 1000 hombres
age_40_54 <- 1 # grupo de edad 40-54
age_55_59 <- 0 # no incluye este grupo
age_60_64 <- 0 # no incluye este grupo
age_65_69 <- 0 # no incluye este grupo
city_fredericia <- 0 # no es Fredericia
city_horsens <- 0 # no es Horsens
city_kolding <- 0 # no es Kolding
# Calcular la predicción
log_mu <- coef(modelo_danishlc)[1] +
coef(modelo_danishlc)[2] * pop_vejle +
coef(modelo_danishlc)[3] * age_40_54 +
coef(modelo_danishlc)[4] * age_55_59 +
coef(modelo_danishlc)[5] * age_60_64 +
coef(modelo_danishlc)[6] * age_65_69 +
coef(modelo_danishlc)[7] * city_fredericia +
coef(modelo_danishlc)[8] * city_horsens +
coef(modelo_danishlc)[9] * city_kolding
# Exponenciar para obtener el número esperado de casos
num_esperado <- exp(log_mu)
# Imprimir la conclusión
cat("Se estima que en la ciudad de Vejle, para un grupo de 1000 hombres de 40 años,",
"el número medio de casos de cáncer de pulmón es:", round(num_esperado, 2), "\n")Se estima que en la ciudad de Vejle, para un grupo de 1000 hombres de 40 años, el número medio de casos de cáncer de pulmón es: 12.32 (c) Muestre los residuos de pearson y la deviance. ¿Existe algun dato atıpico?
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelo_danishlc)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_modelo_danishlc,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 1 valor atipico en el modelo inicial (Rango entre 10 y 15) y en el modelo final (Rango entre 10 y 15), sin embargo no parecen afectar al modelo.
# Supongamos que resid_modelo_danishlc es un data frame que contiene los residuos
# Asegúrate de que tienes ambos tipos de residuos calculados
resid_modelo_danishlc <- data.frame(
deviance = residuals(modelo_danishlc, type = "deviance"),
pearson = residuals(modelo_danishlc, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modelo_danishlc = resid_modelo_danishlc[resid_modelo_danishlc$deviance > 2 & resid_modelo_danishlc$pearson > 2, ]
resul_filt_modelo_danishlc[1] deviance pearson
<0 rows> (o 0- extensión row.names)La filtración de residuos en resultados_filt_inv_blocks no devolvió filas, lo que indica que no hay residuos que cumplan con las condiciones especificadas (deviance > 2 y pearson > 2). Esto sugiere que, en tu modelo, los residuos están dentro de un rango que no se considera problemático según esos umbrales.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_modelo_danishlc)
StudRes Hat CookD
1 0.54247 0.041667 0.013503
2 0.54247 0.041667 0.013503
9 1.74642 0.041667 0.156090
12 -2.95895 0.041667 0.261410
13 -2.00519 0.041667 0.138266Conclusiones
La observación 12 es la única que cumple con el criterio de ser atípica basada en los residuos estudiantizados.
Todas las observaciones tienen un leverage bajo, lo que es positivo.
Todas las observaciones tienen valores de Cook’s D que están muy por debajo de 1, por lo que no parecen tener una influencia desproporcionada en el modelo.
Por lo tanto, solo la observación 12 sería considerada como un punto a revisar más a fondo, mientras que las demás parecen comportarse normalmente.
2. En la base de datos trout se busca determinar el numero de huevos de trucha muertos por una toxina aplicada en los animales.
Descripcion de la BD: El DF describe efecto del cianato de potasio en los huevos de trucha y describe la supervivencia de los huevos de trucha expuestos al cianato de potasio.
Variables
Conc: Representa la concentración de cianato de potasio, medida en miligramos por litro (mg/litro). Variable numérica.
When: Indica el momento en que se aplica el tóxico. Se clasifica en dos niveles: Now (Ahora): Cuando se aplica el tóxico de inmediato. Later (Más tarde): Después de que los huevos han sido sometidos a un proceso de endurecimiento en agua. Variable categórica (factor).
Number: Número total de huevos utilizados en la prueba. Variable numérica.
Dead: Descripción: Número de huevos que resultaron muertos después de la exposición al cianato de potasio, medido a los 19 días. Variable numérica.
Los datos recopilados muestran el número de huevos de trucha que murieron a los 19 días tras ser expuestos al cianato de potasio (kscn). En cada frasco, se separaron los huevos en dos grupos: la mitad fue sometida a un proceso de endurecimiento en agua antes de que se aplicara el tóxico, mientras que la otra mitad fue expuesta de inmediato. Este diseño experimental permite investigar de manera efectiva cómo la concentración del cianato y el momento de su aplicación influyen en la tasa de supervivencia de los huevos de trucha. Así, el análisis de estos datos es crucial para comprender el impacto del cianato de potasio en el desarrollo de esta especie y para evaluar posibles riesgos ambientales asociados a su uso.
(a) Ajuste el “mejor” modelo de regresion con para estimar el numero de huevos dañados como variable de interes. Muestre las medidas de bondad de ajuste.
data("trout")
dim(trout)[1] 48 4glimpse(trout)Rows: 48
Columns: 4
$ Conc <int> 90, 90, 90, 90, 180, 180, 180, 180, 360, 360, 360, 360, 720, 72…
$ When <fct> Later, Later, Later, Later, Later, Later, Later, Later, Later, …
$ Number <int> 111, 97, 108, 122, 68, 109, 109, 118, 98, 110, 129, 103, 83, 87…
$ Dead <int> 8, 10, 10, 9, 4, 6, 11, 6, 6, 5, 9, 17, 2, 3, 16, 9, 60, 47, 49…trout$When_Now <- ifelse(trout$When == "Now", 1, 0)
trout$When <- NULL
head(trout) Conc Number Dead When_Now
1 90 111 8 0
2 90 97 10 0
3 90 108 10 0
4 90 122 9 0
5 180 68 4 0
6 180 109 6 0library(lmtest)
modelo_huevos = glm(Dead ~ ., family = poisson(link = "log"), data = trout)barplot(trout$Dead)
No hay una distribucion clara.
dispersiontest(modelo_huevos,trafo=1)
Overdispersion test
data: modelo_huevos
z = 2.92, p-value = 0.0018
alternative hypothesis: true alpha is greater than 0
sample estimates:
alpha
7.8078 Ho: No hay sobredispersión Ha: Hay sobredispersión
Se rechaza Ho, hay evidencia suficiente para afirmar que hay sobredispersión en los datos.
Con base a esta Conclusión vamos a tomar el modelo de regresion binomial negativa
modelo_nb_huevos = glm.nb(Dead ~ ., data = trout, link = log)Seleccion del Mejor Subconjunto de Variables
step(modelo_nb_huevos)Start: AIC=368.18
Dead ~ Conc + Number + When_Now
Df Deviance AIC
- When_Now 1 54.1 366
<none> 54.0 368
- Number 1 59.6 372
- Conc 1 120.0 432
Step: AIC=366.3
Dead ~ Conc + Number
Df Deviance AIC
<none> 54.0 366
- Number 1 59.5 370
- Conc 1 119.8 430
Call: glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012342402,
link = log)
Coefficients:
(Intercept) Conc Number
0.593366 0.000848 0.012387
Degrees of Freedom: 47 Total (i.e. Null); 45 Residual
Null Deviance: 135
Residual Deviance: 54 AIC: 368mejor_modelo_nb_huevos = glm.nb(Dead ~ Conc + Number, data = trout, init.theta = 2.012342402, link = log)
summary(mejor_modelo_nb_huevos)
Call:
glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.593340 0.597246 0.99 0.320
Conc 0.000848 0.000110 7.70 1.4e-14 ***
Number 0.012387 0.005187 2.39 0.017 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.0123) family taken to be 1)
Null deviance: 135.249 on 47 degrees of freedom
Residual deviance: 54.031 on 45 degrees of freedom
AIC: 368.3
Number of Fisher Scoring iterations: 1
Theta: 2.012
Std. Err.: 0.480
2 x log-likelihood: -360.298 Medidas de Bondad de Ajuste
Pseudo R^2
library(DescTools)
# Calcular el Pseudo R² de Nagelkerke
pseudo_r2_danishlc <- PseudoR2(modelo_nb_huevos, 'Nagelkerke')
# Mostrar el resultado
cat("El Pseudo R² de Nagelkerke es:", pseudo_r2_danishlc, "\n")El Pseudo R² de Nagelkerke es: 0.63181 Conclusiones
El Pseudo R² de Nagelkerke de 0.17481 indica que aproximadamente el 17.5% de la variabilidad en los casos de cáncer de pulmón se puede atribuir a las variables incluidas en este modelo.
summary(mejor_modelo_nb_huevos)
Call:
glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.593340 0.597246 0.99 0.320
Conc 0.000848 0.000110 7.70 1.4e-14 ***
Number 0.012387 0.005187 2.39 0.017 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.0123) family taken to be 1)
Null deviance: 135.249 on 47 degrees of freedom
Residual deviance: 54.031 on 45 degrees of freedom
AIC: 368.3
Number of Fisher Scoring iterations: 1
Theta: 2.012
Std. Err.: 0.480
2 x log-likelihood: -360.298 Null deviance: 27.704 on 23 degrees of freedom Residual deviance: 27.704 on 23 degrees of freedom AIC: 126.1
deviance_modelo_nb_huevos = 1-pchisq(27.704-27.704, 23-23); deviance_modelo_nb_huevos [1] 1# Mostrar el resultado
cat("El modelo deviance es:", deviance_modelo_nb_huevos, "\n")El modelo deviance es: 1 Conclusión
Una devianza de 1 en el modelo indica que el ajuste es bastante satisfactorio. Esto sugiere que la diferencia entre los datos observados y los esperados según el modelo es relativamente pequeña.
Prueba de Coeficientes Estimados
Modelo Inicial
summary(modelo_nb_huevos)
Call:
glm.nb(formula = Dead ~ ., data = trout, link = log, init.theta = 2.017737155)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.599994 0.597221 1.00 0.315
Conc 0.000847 0.000110 7.70 1.4e-14 ***
Number 0.012667 0.005282 2.40 0.016 *
When_Now -0.076648 0.222508 -0.34 0.730
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.0177) family taken to be 1)
Null deviance: 135.56 on 47 degrees of freedom
Residual deviance: 54.02 on 44 degrees of freedom
AIC: 370.2
Number of Fisher Scoring iterations: 1
Theta: 2.018
Std. Err.: 0.482
2 x log-likelihood: -360.177 # Coeficientes estimados del modelo sin el intercepto
coefficients_modelo_nb_huevos <- c(
Conc = 0.000847,
Number = 0.012667,
When_Now = -0.076648
)
# Calcular los coeficientes transformados
transf_coef_modelo_nb_huevos <- ifelse(coefficients_modelo_nb_huevos > 0,
exp(coefficients_modelo_nb_huevos),
1 / exp(-coefficients_modelo_nb_huevos))
# Mostrar los resultados
data.frame(coefficients_modelo_nb_huevos = coefficients_modelo_nb_huevos, Transformed = transf_coef_modelo_nb_huevos) coefficients_modelo_nb_huevos Transformed
Conc 0.000847 1.00085
Number 0.012667 1.01275
When_Now -0.076648 0.92622Conclusiones
Concentración de cianato (Conc): Indica que un incremento en la concentración de cianato de potasio tiene un efecto positivo en el número de huevos muertos, aunque el factor transformado 1.00085 sugiere que este efecto es bastante pequeño y cercano a 1, lo que implica que el cambio es casi despreciable en términos prácticos.
Número de huevos utilizados (Number): El coeficiente de 0.012667 sugiere que, a medida que aumenta el número de huevos utilizados, también se espera un ligero incremento en el número de huevos muertos, con un factor transformado de 1.01275. Esto significa que por cada huevo adicional utilizado, el número de huevos muertos aumenta en aproximadamente un 1.3%, lo que es una Alteración en el contexto de un gran número de huevos.
Momento de aplicación (When_Now): Un coeficiente de -0.076648 implica que aplicar el tóxico “ahora” (en lugar de “más tarde”) reduce el número de huevos muertos, con un factor transformado de 0.92622. Esto sugiere que el riesgo de mortalidad es un 7.4% menor cuando el tóxico se aplica de inmediato, lo que podría indicar que la exposición temprana tiene un efecto protector.
Mejor Modelo
summary(mejor_modelo_nb_huevos)
Call:
glm.nb(formula = Dead ~ Conc + Number, data = trout, init.theta = 2.012337563,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.593340 0.597246 0.99 0.320
Conc 0.000848 0.000110 7.70 1.4e-14 ***
Number 0.012387 0.005187 2.39 0.017 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.0123) family taken to be 1)
Null deviance: 135.249 on 47 degrees of freedom
Residual deviance: 54.031 on 45 degrees of freedom
AIC: 368.3
Number of Fisher Scoring iterations: 1
Theta: 2.012
Std. Err.: 0.480
2 x log-likelihood: -360.298 # Coeficientes estimados del modelo sin el intercepto
coef_mejor_modelo_nb_huevos <- c(
Conc = 0.000848,
Number = 0.012387
)
# Calcular los coeficientes transformados
transf_coef_mejor_modelo_nb_huevos <- ifelse(coef_mejor_modelo_nb_huevos > 0,
exp(coef_mejor_modelo_nb_huevos),
1 / exp(-coef_mejor_modelo_nb_huevos))
# Mostrar los resultados
data.frame(coef_mejor_modelo_nb_huevos = coef_mejor_modelo_nb_huevos, Transformed = transf_coef_mejor_modelo_nb_huevos) coef_mejor_modelo_nb_huevos Transformed
Conc 0.000848 1.0008
Number 0.012387 1.0125Conclusión
Concentración de cianato (Conc): El coeficiente transformado de 0.000848 indica que un aumento en la concentración de cianato de potasio se asocia con un leve incremento en el número de huevos muertos. El factor transformado de 1.0008 sugiere que, por cada incremento unitario en la concentración, el número de huevos muertos aumenta en aproximadamente un 0.08%. Esto señala que, aunque hay una relación positiva, el impacto es muy sutil en términos prácticos.
Número de huevos utilizados (Number): El coeficiente de 0.012387 implica que un incremento en el número de huevos utilizados se relaciona con un aumento en el número de huevos muertos. El factor transformado de 1.0125 indica que, por cada huevo adicional utilizado, se espera un incremento del 1.25% en el número de huevos muertos. Este resultado sugiere que el efecto del tamaño de la muestra es más relevante y tiene un impacto notable en la mortalidad de los huevos.
Los coeficientes entre el modelo inicial y el mejor modelo son bastante similares para las variables de concentración y número de huevos. Sin embargo, la eliminación de la variable “When_Now” en el mejor modelo sugiere que esta variable no aporta valor predictivo adicional y puede haber sido un factor de confusión en el modelo inicial.
(b) Realice una estimacion con Con = 90, When = later y Number = 234.
# Crear un nuevo dataframe con los valores de interés
nueva_data <- data.frame(Conc = 90, Number = 234, When_Now = factor("Later", levels = c("Now", "Later")))
# Realizar la predicción
prediccion <- predict(mejor_modelo_nb_huevos, newdata = nueva_data, type = "response")
# Mostrar el resultado
cat("La estimación del número de huevos de trucha dañados, bajo las condiciones de una concentración de cianato de potasio de 90 mg/litro, aplicado más tarde, y con un total de 234 huevos, es de aproximadamente:", round(prediccion), "\n")La estimación del número de huevos de trucha dañados, bajo las condiciones de una concentración de cianato de potasio de 90 mg/litro, aplicado más tarde, y con un total de 234 huevos, es de aproximadamente: 35 (c) Existe datos atıpicos o influyentes?. Muestre los residuales y las medidas de influencia.
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelo_nb_huevos)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_modelo_nb_huevos,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 4 valores atipicos en el modelo inicial (Rango entre 30 y 40) y en el modelo final se observa 2 valores atípicos (Rango entre 35 y 40).
# Supongamos que resid_modelo_nb_huevos es un data frame que contiene los residuos
# Asegúrate de que tienes ambos tipos de residuos calculados
resid_modelo_nb_huevos <- data.frame(
deviance = residuals(modelo_nb_huevos, type = "deviance"),
pearson = residuals(modelo_nb_huevos, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modelo_nb_huevos = resid_modelo_nb_huevos[resid_modelo_nb_huevos$deviance > 2 & resid_modelo_nb_huevos$pearson > 2, ]
resul_filt_modelo_nb_huevos deviance pearson
38 2.5809 4.3342
39 2.2148 3.4911Los valores de devianza (2.5809 y 2.2148) indican la discrepancia entre el modelo ajustado y los datos observados. Un valor más bajo de devianza sugiere un mejor ajuste del modelo a los datos.
os residuos de Pearson obtenidos (4.3342 y 3.4911) indican que el modelo puede no estar ajustándose adecuadamente a los datos. Dado que estos valores son relativamente altos, sugieren que hay discrepancias significativas entre los valores observados y los valores predichos por el modelo.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_modelo_nb_huevos)
StudRes Hat CookD
26 0.51619 0.262014 0.0449878
38 2.46973 0.023612 0.1400845
39 2.08432 0.032006 0.1237827
40 -2.49186 0.102878 0.0683893
45 -0.13786 0.136187 0.0011561Conclusiones
Residuos Estudiantizados (StudRes): Un valor absoluto mayor que 2 o 3 podría indicar un posible problema en el ajuste del modelo. En este caso las observaciones (38, 39, 40) son preocupantes
Hat (Leverage): Un valor de leverage (Hat) mayor a 0.2 puede indicar que el punto tiene un alto potencial de influir en el modelo. En este caso la observación 26 tiene un leverage elevado, lo que podría ser problemático.
Cook’s Distance (CookD): Un valor mayor a 1 puede sugerir que el punto tiene un impacto desproporcionado en la estimación del modelo. En este caso todos los valores de CookD están por debajo de 1, lo que indica que no hay puntos que influyan desproporcionadamente en el modelo.
3. Un grupo de 24 aficionados al ciclismo de nivel principiante se dieron cita en un parque para una carrera de ciclismo que consiste en andar en bicicleta durante 30 minutos a lo largo de un circuito de 1.3 millas. El ganador es el que pedalea la distancia mas larga. La informacion registrada para cada participante es el genero (M/F), si tuvo una experiencia previa en carreras ası (sı/no), autoevaluacion de habilidades para terminar la carrera y hacerlo bien (en una escala de 10 puntos, siendo 10 la confianza mas alta) y la distancia recorrida en bicicleta (en millas). [?] Los datos se presentan a continuacion:
| gender | prior expr | self eval | distance |
|---|---|---|---|
| F | no | 2 | 1.90 |
| F | no | 2 | 2.10 |
| F | yes | 8 | 3.80 |
| F | yes | 4 | 3.00 |
| M | no | 5 | 4.20 |
| F | yes | 10 | 8.20 |
| F | no | 3 | 3.10 |
| F | no | 4 | 2.40 |
| F | no | 5 | 4.60 |
| M | yes | 6 | 8.70 |
| F | no | 6 | 4.70 |
| M | yes | 7 | 4.20 |
| F | no | 7 | 4.40 |
| F | yes | 3 | 3.10 |
| M | yes | 10 | 6.40 |
| F | yes | 4 | 3.20 |
| F | no | 6 | 5.10 |
| M | no | 10 | 5.90 |
| F | no | 6 | 5.00 |
| M | yes | 3 | 3.60 |
| F | no | 7 | 4.40 |
| M | yes | 10 | 11.20 |
| F | yes | 3 | 3.00 |
| M | yes | 7 | 4.30 |
(a) Ajuste el “mejor” modelo de regresion que estime la distancia recorrida. Muestre las medidas de bondad de ajuste. Interprete los coeficientes estimados del modelo.
datos_ciclismo <- data.frame(
gender = c("F", "F", "F", "F", "M", "F", "F", "F", "F", "M", "F", "M", "F", "F", "M", "F", "F", "M", "F", "M", "F", "M", "F", "M"),
prior_expr = c("no", "no", "yes", "yes", "no", "yes", "no", "no", "no", "yes", "no", "yes", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "yes", "yes", "yes"),
self_eval = c(2, 2, 8, 4, 5, 10, 3, 4, 5, 6, 6, 7, 7, 3, 10, 4, 6, 10, 6, 3, 7, 10, 3, 7),
dist = c(1.90, 2.10, 3.80, 3.00, 4.20, 8.20, 3.10, 2.40, 4.60, 8.70, 4.70, 4.20, 4.40, 3.10, 6.40, 3.20, 5.10, 5.90, 5.00, 3.60, 4.40, 11.20, 3.00, 4.30)
)
# Mostrar el data frame
head(datos_ciclismo) gender prior_expr self_eval dist
1 F no 2 1.9
2 F no 2 2.1
3 F yes 8 3.8
4 F yes 4 3.0
5 M no 5 4.2
6 F yes 10 8.2# Convertir a variables dummy
datos_ciclismo <- datos_ciclismo %>%
mutate(
gender_M = ifelse(gender == "M", 1, 0),
pr_expr_yes = ifelse(prior_expr == "yes", 1, 0)
) %>%
select(-gender, -prior_expr) # Eliminar las columnas originales
head(datos_ciclismo) self_eval dist gender_M pr_expr_yes
1 2 1.9 0 0
2 2 2.1 0 0
3 8 3.8 0 1
4 4 3.0 0 1
5 5 4.2 1 0
6 10 8.2 0 1modeloinv=glm(dist~.,family=inverse.gaussian(link="log"),data=datos_ciclismo)
modeloga=glm(dist~.,family=Gamma(link="log"),data=datos_ciclismo)# Definición de la función
evaluar_modelo <- function(L_cicl) {
if (L_cicl > 0) {
return(paste("El modelo preferido es el inversamente gaussiano. (L_cicl =", round(L_cicl, 4), ")"))
} else if (L_cicl < 0) {
return(paste("El modelo preferido es el gamma. (L_cicl =", round(L_cicl, 4), ")"))
} else {
return("Los modelos son equivalentes (L_cicl es 0).")
}
}
# Supongamos que L_cicl ya ha sido calculado
L_cicl <- log(logLik(modeloinv)/logLik(modeloga))
resultado <- evaluar_modelo(L_cicl)
# Imprimir el resultado
cat(resultado)El modelo preferido es el gamma. (L_cicl = -0.0584 )Seleccion del Mejor Subconjunto de Variables
step(modeloga)Start: AIC=76.8
dist ~ self_eval + gender_M + pr_expr_yes
Df Deviance AIC
- pr_expr_yes 1 1.30 75.0
- gender_M 1 1.39 76.4
<none> 1.28 76.8
- self_eval 1 3.27 103.9
Step: AIC=75.11
dist ~ self_eval + gender_M
Df Deviance AIC
<none> 1.30 75.1
- gender_M 1 1.46 75.6
- self_eval 1 3.34 104.0
Call: glm(formula = dist ~ self_eval + gender_M, family = Gamma(link = "log"),
data = datos_ciclismo)
Coefficients:
(Intercept) self_eval gender_M
0.667 0.127 0.188
Degrees of Freedom: 23 Total (i.e. Null); 21 Residual
Null Deviance: 4.46
Residual Deviance: 1.3 AIC: 75.1mejor_modeloga=glm(dist ~ self_eval + gender_M + pr_expr_yes, family=Gamma(link="log"),data=datos_ciclismo)
summary(mejor_modeloga)
Call:
glm(formula = dist ~ self_eval + gender_M + pr_expr_yes, family = Gamma(link = "log"),
data = datos_ciclismo)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6494 0.1386 4.69 0.00014 ***
self_eval 0.1266 0.0234 5.41 2.7e-05 ***
gender_M 0.1652 0.1310 1.26 0.22205
pr_expr_yes 0.0571 0.1143 0.50 0.62282
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.068413)
Null deviance: 4.4596 on 23 degrees of freedom
Residual deviance: 1.2797 on 20 degrees of freedom
AIC: 76.8
Number of Fisher Scoring iterations: 5Medidas de Bondad de Ajuste
Pseudo R^2
library(DescTools)
# Calcular el Pseudo R² de Nagelkerke
pseudo_r2_cicl <- PseudoR2(modeloga, 'Nagelkerke')
# Mostrar el resultado
cat("El Pseudo R² de Nagelkerke es:", pseudo_r2_cicl, "\n")El Pseudo R² de Nagelkerke es: 0.73201 Conclusiones
El Pseudo R² de Nagelkerke es el 73.2% de la variabilidad en la distancia recorrida por los participantes en tu estudio. lo que indica que las variables independientes incluidas en el modelo (Autoevaluacion, genero Experiencia Previa) tienen un impacto significativo en la distancia que recorren.
summary(mejor_modeloga)
Call:
glm(formula = dist ~ self_eval + gender_M + pr_expr_yes, family = Gamma(link = "log"),
data = datos_ciclismo)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6494 0.1386 4.69 0.00014 ***
self_eval 0.1266 0.0234 5.41 2.7e-05 ***
gender_M 0.1652 0.1310 1.26 0.22205
pr_expr_yes 0.0571 0.1143 0.50 0.62282
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.068413)
Null deviance: 4.4596 on 23 degrees of freedom
Residual deviance: 1.2797 on 20 degrees of freedom
AIC: 76.8
Number of Fisher Scoring iterations: 5Prueba de la Deviance
Null deviance: 4.4596 on 23 degrees of freedom Residual deviance: 1.2797 on 20 degrees of freedom AIC: 76.8
deviance_modeloga = 1-pchisq(4.4596-1.2797, 23-20)
# Mostrar el resultado
cat("La deviance tiene un valor de:", deviance_modeloga, "\n")La deviance tiene un valor de: 0.36471 Conclusión
La deviance indica un buen ajuste del modelo a los datos.
Prueba de Coeficientes Estimados
Modelo Inicial
summary(modeloga)
Call:
glm(formula = dist ~ ., family = Gamma(link = "log"), data = datos_ciclismo)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6494 0.1386 4.69 0.00014 ***
self_eval 0.1266 0.0234 5.41 2.7e-05 ***
gender_M 0.1652 0.1310 1.26 0.22205
pr_expr_yes 0.0571 0.1143 0.50 0.62282
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.068413)
Null deviance: 4.4596 on 23 degrees of freedom
Residual deviance: 1.2797 on 20 degrees of freedom
AIC: 76.8
Number of Fisher Scoring iterations: 5# Coeficientes estimados del modelo sin el intercepto
coefficients_modeloga <- c(
self_eval = 0.1266,
gender_M = 0.1652,
pr_expr_yes = 0.0571
)
# Calcular los coeficientes transformados
transf_coef_modeloga <- ifelse(coefficients_modeloga > 0,
exp(coefficients_modeloga),
1 / exp(-coefficients_modeloga))
# Mostrar los resultados
data.frame(coefficients_modeloga = coefficients_modeloga, Transformed = transf_coef_modeloga) coefficients_modeloga Transformed
self_eval 0.1266 1.1350
gender_M 0.1652 1.1796
pr_expr_yes 0.0571 1.0588Mejor Modelo
summary(mejor_modeloga)
Call:
glm(formula = dist ~ self_eval + gender_M + pr_expr_yes, family = Gamma(link = "log"),
data = datos_ciclismo)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6494 0.1386 4.69 0.00014 ***
self_eval 0.1266 0.0234 5.41 2.7e-05 ***
gender_M 0.1652 0.1310 1.26 0.22205
pr_expr_yes 0.0571 0.1143 0.50 0.62282
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.068413)
Null deviance: 4.4596 on 23 degrees of freedom
Residual deviance: 1.2797 on 20 degrees of freedom
AIC: 76.8
Number of Fisher Scoring iterations: 5# Coeficientes estimados del modelo sin el intercepto
coefficients_mejor_modeloga <- c(
self_eval = 0.1266,
gender_M = 0.1652,
pr_expr_yes = 0.0571
)
# Calcular los coeficientes transformados
transf_coef_mejor_modeloga <- ifelse(coefficients_mejor_modeloga > 0,
exp(coefficients_mejor_modeloga),
1 / exp(-coefficients_mejor_modeloga))
# Mostrar los resultados
data.frame(coefficients_mejor_modeloga = coefficients_mejor_modeloga, Transformed = transf_coef_mejor_modeloga) coefficients_mejor_modeloga Transformed
self_eval 0.1266 1.1350
gender_M 0.1652 1.1796
pr_expr_yes 0.0571 1.0588Conclusiones
la coincidencia de resultados entre el modelo inicial y el mejor modelo subraya la estabilidad de las relaciones observadas entre las variables y la distancia recorrida.
(b) Estime la distancia que puede recorrer un caballero, sin experiencia en carrera y con una autoevaluacion de 6.
# Crear un nuevo dataframe con los valores para predecir
nuevo_dato <- data.frame(self_eval = 6, gender_M = 1, pr_expr_yes = 0)
# Realizar la predicción
prediccion <- predict(modeloga, nuevo_dato, type = "response")
# Imprimir la predicción
cat("La distancia estimada es:", prediccion, "unidades.\n")La distancia estimada es: 4.8254 unidades.(c) Existe datos atıpicos o influyentes?. Muestre los residuales y las medidas de influencia.
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modeloga)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_modeloga,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 1 valor atipico en el modelo inicial (aprox 10) y en el modelo final (aprox 10), sin embargo no parecen afectar al modelo.
# Supongamos que resid_modeloga es un data frame que contiene los residuos
# Asegúrate de que tienes ambos tipos de residuos calculados
resid_modeloga <- data.frame(
deviance = residuals(modeloga, type = "deviance"),
pearson = residuals(modeloga, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modeloga = resid_modeloga[resid_modeloga$deviance > 2 & resid_modeloga$pearson > 2, ]
resul_filt_modeloga[1] deviance pearson
<0 rows> (o 0- extensión row.names)La filtración de residuos en resultados_filt_inv_blocks no devolvió filas, lo que indica que no hay residuos que cumplan con las condiciones especificadas (deviance > 2 y pearson > 2). Esto sugiere que, en tu modelo, los residuos están dentro de un rango que no se considera problemático según esos umbrales.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_modeloga)
StudRes Hat CookD
3 -1.62792 0.20055 0.116588
6 0.64833 0.32277 0.051222
10 3.04644 0.15087 0.377718
18 -1.36921 0.30240 0.157088Conclusiones
Dado que todos los puntos cumplen con las tres reglas de diagnóstico, Se puede concluir que estas observaciones no presentan problemas de influencia y son consistentes con el ajuste del modelo.
4. Los datos siguientes son de un estudio de 60 sujetos: 30 diabéticos y 30 no diabéticos que tenía como objetivo determinar el efecto de la obesidad, los antecedentes familiares y practicar algún tipo de actividad física en el desarrollo de la diabetes. IMC hace referencia al índice de masa corporal, ATF hace referencia a realizar (1=si) algún tipo de actividad física y HF hace referencia a los antecedentes familiares de la enfermedad.
Sujetos Diabéticos (SD)
Sujetos No Diabéticos (SND)
| IMC | SD_HF | SD_ATF | IMC | SND_HF | SND_ATF |
|---|---|---|---|---|---|
| 22.10 | 1 | 1 | 26.70 | 0 | 1 |
| 31.30 | 0 | 0 | 24.40 | 0 | 1 |
| 33.80 | 1 | 0 | 29.40 | 0 | 0 |
| 33.70 | 1 | 1 | 26.00 | 0 | 0 |
| 23.10 | 1 | 1 | 24.20 | 1 | 0 |
| 26.80 | 1 | 0 | 29.70 | 0 | 0 |
| 32.30 | 1 | 0 | 30.20 | 0 | 1 |
| 31.40 | 1 | 0 | 23.40 | 0 | 1 |
| 37.60 | 1 | 0 | 42.40 | 0 | 0 |
| 32.40 | 1 | 0 | 25.80 | 0 | 0 |
| 29.10 | 0 | 1 | 39.80 | 0 | 1 |
| 28.60 | 0 | 1 | 31.60 | 0 | 0 |
| 35.90 | 0 | 0 | 21.80 | 1 | 1 |
| 30.40 | 0 | 0 | 24.20 | 0 | 1 |
| 39.80 | 0 | 0 | 27.80 | 1 | 1 |
| 43.30 | 1 | 0 | 37.50 | 1 | 1 |
| 32.50 | 0 | 0 | 27.90 | 1 | 1 |
| 28.70 | 0 | 1 | 25.30 | 1 | 0 |
| 30.30 | 0 | 0 | 31.30 | 0 | 1 |
| 32.50 | 1 | 0 | 34.50 | 1 | 1 |
| 32.50 | 1 | 0 | 25.40 | 0 | 1 |
| 21.60 | 1 | 1 | 27.00 | 1 | 1 |
| 24.40 | 0 | 1 | 31.10 | 0 | 0 |
| 46.70 | 1 | 0 | 27.30 | 0 | 1 |
| 28.60 | 1 | 1 | 24.00 | 0 | 0 |
| 29.70 | 0 | 0 | 33.50 | 0 | 0 |
| 29.60 | 0 | 1 | 20.70 | 0 | 0 |
| 22.80 | 0 | 0 | 29.20 | 1 | 1 |
| 34.80 | 1 | 0 | 30.00 | 0 | 1 |
| 37.30 | 1 | 0 | 26.50 | 0 | 0 |
(a) Ajuste el “mejor” modelo de regresión que explique el tener diabetes. Muestre las medidas de bondad de ajuste.
obesidad <- data.frame(
SD_IMC_1 = c(22.10, 31.30, 33.80, 33.70, 23.10, 26.80, 32.30, 31.40, 37.60, 32.40, 29.10, 28.60, 35.90, 30.40, 39.80, 43.30, 32.50, 28.70, 30.30, 32.50, 32.50, 21.60, 24.40, 46.70, 28.60, 29.70, 29.60, 22.80, 34.80, 37.30),
SD_HF = c(1, 0, 1, 1, 1, 1, 1, 1, 1, 1,0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1),
SD_ATF = c(1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0),
SND_IMC_2 = c(26.70, 24.40, 29.40, 26.00, 24.20, 29.70, 30.20, 23.40, 42.40, 25.80, 39.80, 31.60, 21.80, 24.20, 27.80, 37.50, 27.90, 25.30, 31.30, 34.50, 25.40, 27.00, 31.10, 27.30, 24.00, 33.50, 20.70, 29.20, 30.00, 26.50),
SND_HF_2 = c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),
SND_ATF_2 = c(1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0)
)
# Crear la variable 'diabetico' basándote en antecedentes familiares y criterios adicionales
obesidad$diabetico <- ifelse(obesidad$SD_HF == 1 | obesidad$SD_IMC_1 >= 30, 1, 0)
obesidad$diabetico <- ifelse(obesidad$SD_HF == 1 | obesidad$SD_IMC_1 >= 30, 1, 0)mod_diabetes <- glm(diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 + SND_ATF_2, data = obesidad, family = binomial)
# Resumen del modelo
summary(mod_diabetes)
Call:
glm(formula = diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 +
SND_ATF_2, family = binomial, data = obesidad)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.183 5.323 0.03 0.973
SD_IMC_1 0.236 0.163 1.45 0.148
SD_ATF -1.341 1.284 -1.04 0.296
SND_IMC_2 -0.208 0.121 -1.71 0.086 .
SND_HF_2 -0.207 1.359 -0.15 0.879
SND_ATF_2 1.953 1.371 1.42 0.154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32.596 on 29 degrees of freedom
Residual deviance: 20.812 on 24 degrees of freedom
AIC: 32.81
Number of Fisher Scoring iterations: 6Seleccion del Mejor Subconjunto de Variables
step(mod_diabetes)Start: AIC=32.81
diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 + SND_ATF_2
Df Deviance AIC
- SND_HF_2 1 20.8 30.8
- SD_ATF 1 22.0 32.0
<none> 20.8 32.8
- SND_ATF_2 1 23.2 33.2
- SD_IMC_1 1 23.5 33.5
- SND_IMC_2 1 24.3 34.3
Step: AIC=30.84
diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_ATF_2
Df Deviance AIC
- SD_ATF 1 22.0 30.0
<none> 20.8 30.8
- SND_ATF_2 1 23.2 31.2
- SD_IMC_1 1 23.7 31.7
- SND_IMC_2 1 24.3 32.3
Step: AIC=29.97
diabetico ~ SD_IMC_1 + SND_IMC_2 + SND_ATF_2
Df Deviance AIC
<none> 22.0 30.0
- SND_IMC_2 1 24.9 30.9
- SND_ATF_2 1 25.5 31.5
- SD_IMC_1 1 28.6 34.6
Call: glm(formula = diabetico ~ SD_IMC_1 + SND_IMC_2 + SND_ATF_2, family = binomial,
data = obesidad)
Coefficients:
(Intercept) SD_IMC_1 SND_IMC_2 SND_ATF_2
-3.200 0.302 -0.186 2.079
Degrees of Freedom: 29 Total (i.e. Null); 26 Residual
Null Deviance: 32.6
Residual Deviance: 22 AIC: 30mejor_mod_diabetes = glm(diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 + SND_ATF_2, family = binomial(link = "logit"), data = obesidad)Medidas de Bondad de Ajuste
Pseudo R^2
library(DescTools)
# Calcular el Pseudo R² de Nagelkerke
pseudo_r2_danishlc <- PseudoR2(mod_diabetes, 'Nagelkerke')
# Mostrar el resultado
cat("El Pseudo R² de Nagelkerke es:", pseudo_r2_danishlc, "\n")El Pseudo R² de Nagelkerke es: 0.49024 Conclusiones
Pseudo R² de Nagelkerke de 0.49024 indica que aproximadamente el 49% de la variabilidad en la variable diabetes se puede atribuir a las variables incluidas en este modelo
summary(mejor_mod_diabetes)
Call:
glm(formula = diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 +
SND_ATF_2, family = binomial(link = "logit"), data = obesidad)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.183 5.323 0.03 0.973
SD_IMC_1 0.236 0.163 1.45 0.148
SD_ATF -1.341 1.284 -1.04 0.296
SND_IMC_2 -0.208 0.121 -1.71 0.086 .
SND_HF_2 -0.207 1.359 -0.15 0.879
SND_ATF_2 1.953 1.371 1.42 0.154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32.596 on 29 degrees of freedom
Residual deviance: 20.812 on 24 degrees of freedom
AIC: 32.81
Number of Fisher Scoring iterations: 6Prueba de la Deviance
Null deviance: 32.596 on 29 degrees of freedom Residual deviance: 20.812 on 24 degrees of freedom AIC: 32.81
deviance_mod_diabetes = 1-pchisq(32.596-20.812, 29-24)
# Mostrar el resultado
cat("El modelo deviance es:", deviance_mod_diabetes, "\n")El modelo deviance es: 0.03787 Conclusión
El valor de deviance de 0.03787 es bastante bajo y sugiere que el modelo se ajusta bien a los datos.
Prueba de Coeficientes Estimados
Modelo Inicial
summary(mod_diabetes)
Call:
glm(formula = diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 +
SND_ATF_2, family = binomial, data = obesidad)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.183 5.323 0.03 0.973
SD_IMC_1 0.236 0.163 1.45 0.148
SD_ATF -1.341 1.284 -1.04 0.296
SND_IMC_2 -0.208 0.121 -1.71 0.086 .
SND_HF_2 -0.207 1.359 -0.15 0.879
SND_ATF_2 1.953 1.371 1.42 0.154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32.596 on 29 degrees of freedom
Residual deviance: 20.812 on 24 degrees of freedom
AIC: 32.81
Number of Fisher Scoring iterations: 6# Coeficientes estimados del modelo sin el intercepto
coefficients_mod_diabetes <- c(
SD_IMC_1 = 0.236,
SD_ATF = -1.341,
SND_IMC_2 = -0.208,
SND_HF_2 = -0.207,
SND_ATF_2 = 1.953
)
# Calcular los coeficientes transformados
transf_coef_mod_diabetes <- ifelse(coefficients_mod_diabetes > 0,
exp(coefficients_mod_diabetes),
1 / exp(-coefficients_mod_diabetes))
# Mostrar los resultados
data.frame(coefficients_mod_diabetes = coefficients_mod_diabetes, Transformed = transf_coef_mod_diabetes) coefficients_mod_diabetes Transformed
SD_IMC_1 0.236 1.26617
SD_ATF -1.341 0.26158
SND_IMC_2 -0.208 0.81221
SND_HF_2 -0.207 0.81302
SND_ATF_2 1.953 7.04981Conclusiones
El aumento de una unidad en SD_IMC_1 está asociado con un aumento del 26.6% en las probabilidades de ser diabético, manteniendo constantes las demás variables. La transformación indica que el efecto es multiplicativo en términos de odds (odds ratio).
Un aumento de una unidad en SD_ATF reduce las probabilidades de ser diabético en aproximadamente un 26.1%. Esto refleja una reducción en la probabilidad de diabetes.
Un incremento en SND_IMC_2 también reduce las probabilidades de diabetes en un 18.8%. Aunque no es tan significativo como SD_ATF, sigue indicando una relación inversa.
Un aumento en SND_HF_2 está asociado con una reducción del 18.6% en las probabilidades de diabetes. Refuerza la tendencia de que mayores valores en esta variable se asocian con menores probabilidades.
Un aumento en SND_ATF_2 incrementa las probabilidades de ser diabético en aproximadamente un 195% (7.05 en términos de odds), indicando que esta variable tiene un efecto muy fuerte en la probabilidad de diabetes. Esto implica que es un factor de riesgo considerable.
Nota: Para obtener el cambio porcentual, puedes usar la fórmula Cambio% =(e^coef−1)×100.
Mejor Modelo
summary(mejor_mod_diabetes)
Call:
glm(formula = diabetico ~ SD_IMC_1 + SD_ATF + SND_IMC_2 + SND_HF_2 +
SND_ATF_2, family = binomial(link = "logit"), data = obesidad)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.183 5.323 0.03 0.973
SD_IMC_1 0.236 0.163 1.45 0.148
SD_ATF -1.341 1.284 -1.04 0.296
SND_IMC_2 -0.208 0.121 -1.71 0.086 .
SND_HF_2 -0.207 1.359 -0.15 0.879
SND_ATF_2 1.953 1.371 1.42 0.154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 32.596 on 29 degrees of freedom
Residual deviance: 20.812 on 24 degrees of freedom
AIC: 32.81
Number of Fisher Scoring iterations: 6Conclusión
Indica de que el modelo ha capturado adecuadamente las relaciones en los datos.
El conjunto de variables elegidas es el adecuado para el modelo.
(b) Estime la probabilidad de que un sujeto con antecedentes familiares de la enfermedad, que realiza actividad física y con un IMC = 26 desarrolle la enfermedad.
# Crear un nuevo dataframe con los valores de interés
nueva_data <- data.frame(SD_IMC_1 = 26, SD_HF = 1, SD_ATF = 1,
SND_IMC_2 = 0, SND_HF_2 = 0, SND_ATF_2 = 0)
# Realizar la predicción
prediccion <- predict(mejor_mod_diabetes, newdata = nueva_data, type = "response")
# Mostrar el resultado
cat("La probabilidad de que un sujeto con IMC de 26, antecedentes familiares de diabetes y que realiza actividad física desarrolle diabetes es de aproximadamente:", round(prediccion * 100, 2), "%\n")La probabilidad de que un sujeto con IMC de 26, antecedentes familiares de diabetes y que realiza actividad física desarrolle diabetes es de aproximadamente: 99.32 %(c) Estime los residuos de Pearson y Deviance, ¿existe algún dato atípico?
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(mod_diabetes)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_mod_diabetes,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 1 valor atipico en el mejor modelo (Rango entre 25 y 30), sin embargo no parecen afectar al modelo.
Calculo de los Residuos
resid_mod_diabetes <- data.frame(
deviance = residuals(mod_diabetes, type = "deviance"),
pearson = residuals(mod_diabetes, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_mod_diabetes = resid_mod_diabetes[resid_mod_diabetes$deviance > 2 & resid_mod_diabetes$pearson > 2, ]
resul_filt_mod_diabetes[1] deviance pearson
<0 rows> (o 0- extensión row.names)La filtración de residuos en resultados_filt_inv_blocks no devolvió filas, lo que indica que no hay residuos que cumplan con las condiciones especificadas (deviance > 2 y pearson > 2). Esto sugiere que, en el modelo, los residuos están dentro de un rango que no se considera problemático según esos umbrales.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_mod_diabetes)
StudRes Hat CookD
9 1.4565 0.55946 0.37025
11 -1.2706 0.57239 0.28728
27 -2.2296 0.24517 0.33322
28 -2.2957 0.39324 0.62202Conclusiones
Todos los puntos analizados (9, 11, 27 y 28) cumplen las reglas de diagnóstico para los residuos estudentizados, el apalancamiento y la distancia de Cook, indicando que no son puntos influyentes en el modelo.
5. El departamento de Administración de la Salud está analizando la utilización de los recursos hospitalarios. Se recuperaron datos sobre 30 hospitales elegidos al azar en EEUU y se está analizando el porcentaje de pacientes de la sala de emergencias (ER) que fueron hospitalizados, la ubicación del hospital (urbano o rural), tipo de hospital (público o privado) y el tamaño del hospital por número de camas. Los datos son:
| percent | hosp | hosp loc | hosp type | N beds |
|---|---|---|---|---|
| 1 | 17 | rural | private | 56 |
| 2 | 39 | rural | public | 144 |
| 3 | 38 | urban | public | 61 |
| 4 | 48 | rural | public | 186 |
| 5 | 30 | rural | private | 132 |
| 6 | 25 | urban | private | 589 |
| 7 | 5 | urban | public | 53 |
| 8 | 4 | rural | private | 73 |
| 9 | 48 | rural | private | 154 |
| 10 | 4 | urban | public | 38 |
| 11 | 26 | rural | private | 318 |
| 12 | 15 | urban | public | 35 |
| 13 | 28 | urban | private | 184 |
| 14 | 34 | urban | private | 173 |
| 15 | 31 | urban | public | 63 |
| 16 | 4 | urban | public | 91 |
| 17 | 6 | urban | public | 77 |
| 18 | 39 | urban | private | 237 |
| 19 | 41 | urban | private | 56 |
| 20 | 45 | rural | public | 43 |
| 21 | 13 | urban | public | 64 |
| 22 | 42 | rural | public | 193 |
| 23 | 28 | urban | private | 363 |
| 24 | 31 | urban | public | 600 |
| 25 | 48 | rural | public | 468 |
| 26 | 41 | rural | public | 311 |
| 27 | 9 | urban | public | 65 |
| 28 | 13 | urban | private | 44 |
| 29 | 44 | urban | public | 479 |
| 30 | 16 | rural | public | 72 |
(a) Modele la proporción de pacientes hospitalizados como variable respuesta. Muestre las medidas de bondad de ajuste.
options(scipen = 10)
# Supongamos que tienes los datos en un dataframe llamado 'hospitales'
hospital <- data.frame(
percent_hosp = c(17, 39, 38, 48, 30, 25, 5, 4, 48, 4,
26, 15, 28, 34, 31, 4, 6, 39, 41, 45,
13, 42, 28, 31, 48, 41, 9, 13, 44, 16),
hosp_loc = c("rural", "rural", "urban", "rural", "rural",
"urban", "urban", "rural", "rural", "urban",
"rural", "urban", "urban", "urban", "urban",
"urban", "urban", "urban", "urban", "urban",
"urban", "rural", "urban", "urban", "rural",
"rural", "urban", "urban", "urban", "rural"),
hosp_type = c("private", "public", "public", "public", "private",
"private", "public", "private", "private", "public",
"private", "public", "private", "private", "public",
"public", "public", "private", "private", "private",
"public", "public", "private", "private", "public",
"public", "public", "private", "private", "public"),
N_beds = c(56, 144, 61, 186, 132,
589, 53, 73, 154, 38,
318, 35, 184, 173, 63,
91, 77, 237, 56, 43,
64, 193, 363, 600, 468,
311, 65, 44, 479, 72)
)
# Crear variable dummy para hosp_loc
hospital$hosp_loc <- ifelse(hospital$hosp_loc == "urban", 1, 0)
# Crear variable dummy para hosp_type
hospital$hosp_type <- ifelse(hospital$hosp_type == "public", 1, 0)
hist(hospital$percent_hosp)
#Crear prop
hospital$pro <- hospital$percent_hosp / 100
#Eliminar la variable percent_hosp
hospital$percent_hosp <- NULLAjustar el Modelo
# Ajustar el modelo beta
modelo_beta <- betareg(pro ~ N_beds + hosp_loc + hosp_type, data = hospital)
# Resumen del modelo
summary(modelo_beta)
Call:
betareg(formula = pro ~ N_beds + hosp_loc + hosp_type, data = hospital)
Quantile residuals:
Min 1Q Median 3Q Max
-2.398 -0.697 0.009 0.805 1.707
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.042176 0.313781 -3.32 0.0009 ***
N_beds 0.001883 0.000772 2.44 0.0147 *
hosp_loc -0.423366 0.259509 -1.63 0.1028
hosp_type -0.159905 0.267233 -0.60 0.5496
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 9.24 2.31 3.99 0.000065 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 20 on 5 Df
Pseudo R-squared: 0.259
Number of iterations: 22 (BFGS) + 2 (Fisher scoring) Calcular el mejor modelo
library(betareg)
# Suponiendo que ya tienes el dataframe 'hospital' con las variables
# Definimos las variables independientes
variables <- c("N_beds", "hosp_loc", "hosp_type")
# Crear un data frame para almacenar resultados
resultados <- data.frame(Modelo = character(), AIC = numeric(), stringsAsFactors = FALSE)
# Probar todas las combinaciones de las variables independientes
for (i in 1:length(variables)) {
combinaciones <- combn(variables, i, simplify = FALSE)
for (combinacion in combinaciones) {
# Crear la fórmula del modelo
formula_modelo <- as.formula(paste("pro ~", paste(combinacion, collapse = " + ")))
# Ajustar el modelo
modelo_beta <- tryCatch({
betareg(formula_modelo, data = hospital)
}, error = function(e) {
cat("Error al ajustar el modelo:", e$message, "\n")
return(NULL)
})
# Calcular el AIC si el modelo se ajustó correctamente
if (!is.null(modelo_beta)) {
aic <- AIC(modelo_beta)
resultados <- rbind(resultados, data.frame(Modelo = paste(combinacion, collapse = ", "), AIC = aic))
}
}
}
# Verificar si hay resultados
if (nrow(resultados) > 0) {
# Encontrar el mejor modelo (menor AIC)
mejor_modelo <- resultados[which.min(resultados$AIC), ]
# Mostrar el mejor modelo
cat("El mejor modelo es:", mejor_modelo$Modelo, "con AIC =", mejor_modelo$AIC, "\n")
} else {
cat("No se encontraron modelos válidos.\n")
}El mejor modelo es: N_beds, hosp_loc con AIC = -31.558 Conclusión
El mejor modelo es:
mejor_mod_beta <- betareg(pro ~ N_beds + hosp_loc, data = hospital)mejor_mod_beta <- betareg(pro ~ N_beds + hosp_loc, data = hospital), con un AIC = -31.558
Medidas de Bondad de Ajuste
Pseudo R^2
Pseudo R-squared: 0.239
Conclusiones
Un Pseudo R-squared de 0.239 en un modelo de regresión beta indica que aproximadamente el 23.9% de la variabilidad en la proporción de pacientes hospitalizados (pro) es explicada por las variables independientes incluidas en el modelo, en este caso, las variables N_beds y hosp_loc.
summary(mejor_mod_beta)
Call:
betareg(formula = pro ~ N_beds + hosp_loc, data = hospital)
Quantile residuals:
Min 1Q Median 3Q Max
-2.217 -0.700 0.049 0.826 1.825
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.150700 0.250012 -4.60 0.0000042 ***
N_beds 0.002048 0.000738 2.77 0.0056 **
hosp_loc -0.429505 0.258618 -1.66 0.0968 .
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 9.16 2.29 3.99 0.000065 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 19.8 on 4 Df
Pseudo R-squared: 0.239
Number of iterations: 17 (BFGS) + 1 (Fisher scoring) Prueba de Coeficientes Estimados
Modelo Inicial
summary(modelo_beta)
Call:
betareg(formula = formula_modelo, data = hospital)
Quantile residuals:
Min 1Q Median 3Q Max
-2.398 -0.697 0.009 0.805 1.707
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.042176 0.313781 -3.32 0.0009 ***
N_beds 0.001883 0.000772 2.44 0.0147 *
hosp_loc -0.423366 0.259509 -1.63 0.1028
hosp_type -0.159905 0.267233 -0.60 0.5496
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 9.24 2.31 3.99 0.000065 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 20 on 5 Df
Pseudo R-squared: 0.259
Number of iterations: 22 (BFGS) + 2 (Fisher scoring) # Coeficientes estimados del modelo sin el intercepto
coefficients_modelo_beta <- c(
N_beds = 0.001883,
hosp_loc = -0.423366,
hosp_type = -0.159905
)
# Calcular los coeficientes transformados
transf_coef_modelo_beta <- ifelse(coefficients_modelo_beta > 0,
exp(coefficients_modelo_beta),
1 / exp(-coefficients_modelo_beta))
# Mostrar los resultados
data.frame(coefficients_modelo_beta = coefficients_modelo_beta, Transformed = transf_coef_modelo_beta) coefficients_modelo_beta Transformed
N_beds 0.001883 1.00188
hosp_loc -0.423366 0.65484
hosp_type -0.159905 0.85222Conclusiones
N_beds (Número de Camas = 0.001883): Por cada cama adicional en el hospital, se espera que la proporción de pacientes hospitalizados aumente en aproximadamente 0.1883%. El factor transformado es 1.00188, lo que indica que un aumento en el número de camas se asocia con un ligero aumento en la proporción de hospitalizaciones.
hosp_loc (Ubicación del Hospital = -0.423366): Comparado con los hospitales urbanos, estar ubicado en un área rural se asocia con una disminución de aproximadamente 42.34% en la proporción de pacientes hospitalizados. El factor transformado es 0.65484, lo que sugiere que la proporción de hospitalización en hospitales rurales es aproximadamente 65.5% de la de hospitales urbanos.
hosp_type (Tipo de Hospital = -0.159905): Comparado con los hospitales privados, ser un hospital público se asocia con una disminución de aproximadamente 15.99% en la proporción de pacientes hospitalizados. El factor transformado es 0.85222, indicando que la proporción de hospitalización en hospitales públicos es aproximadamente 85.2% de la de hospitales privados.
Mejor Modelo
summary(mejor_mod_beta)
Call:
betareg(formula = pro ~ N_beds + hosp_loc, data = hospital)
Quantile residuals:
Min 1Q Median 3Q Max
-2.217 -0.700 0.049 0.826 1.825
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.150700 0.250012 -4.60 0.0000042 ***
N_beds 0.002048 0.000738 2.77 0.0056 **
hosp_loc -0.429505 0.258618 -1.66 0.0968 .
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 9.16 2.29 3.99 0.000065 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 19.8 on 4 Df
Pseudo R-squared: 0.239
Number of iterations: 17 (BFGS) + 1 (Fisher scoring) # Coeficientes estimados del modelo sin el intercepto
coefficients_mejor_mod_beta <- c(
N_beds = 0.002048,
hosp_loc = -0.429505
)
# Calcular los coeficientes transformados
transf_coef_mejor_mod_beta <- ifelse(coefficients_mejor_mod_beta > 0,
exp(coefficients_mejor_mod_beta),
1 / exp(-coefficients_mejor_mod_beta))
# Mostrar los resultados
data.frame(coefficients_mejor_mod_beta = coefficients_mejor_mod_beta, Transformed = transf_coef_mejor_mod_beta) coefficients_mejor_mod_beta Transformed
N_beds 0.002048 1.00205
hosp_loc -0.429505 0.65083Conclusiones
N_beds (Número de Camas = 0.002048): Por cada cama adicional en el hospital, se espera que la proporción de pacientes hospitalizados aumente en aproximadamente 0.2048%. El factor transformado es 1.00205, lo que indica que un aumento en el número de camas se asocia con un leve aumento en la proporción de hospitalizaciones.
hosp_loc (Ubicación del Hospital = -0.429505): Comparado con los hospitales urbanos, estar ubicado en un área rural se asocia con una disminución de aproximadamente 42.95% en la proporción de pacientes hospitalizados. El factor transformado es 0.65083, lo que sugiere que la proporción de hospitalización en hospitales rurales es aproximadamente 65.1% de la de hospitales urbanos.
(b) Estime la proporción de pacientes hospitalizados en un hospital público rural con 50 camas.
# Crear un nuevo dataframe con las condiciones especificadas
nuevo_hospital <- data.frame(
hosp_loc = 0, # rural
hosp_type = 1, # público
N_beds = 50 # número de camas
)
# Asegúrate de que 'modelo_beta' es el modelo que has ajustado previamente
proporcion_estimacion <- predict(modelo_beta, newdata = nuevo_hospital, type = "response")
# Mostrar la estimación
cat("La proporción estimada de pacientes hospitalizados en un hospital público rural con 50 camas es:", round(proporcion_estimacion, 4), "\n")La proporción estimada de pacientes hospitalizados en un hospital público rural con 50 camas es: 0.2483 Conclusión
La proporción estimada de pacientes hospitalizados en un hospital público rural con 50 camas es de aproximadamente 0.2483. Esto sugiere que, bajo las condiciones consideradas en el modelo, alrededor del 24.83% de los pacientes de la sala de emergencias serían hospitalizados en este tipo de hospital.
(c) ¿Existen datos atípicos? Muestre los residuos de deviance y pearson.
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelo_beta)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_mod_beta,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 1 valor atipico en el modelo inicial (Rango entre 5 y 10) y en el modelo final (Rango entre 20 y 25), sin embargo no parecen afectar al modelo.
# Supongamos que resid_modelo_beta es un data frame que contiene los residuos
# Asegúrate de que tienes ambos tipos de residuos calculados
resid_modelo_beta <- data.frame(
deviance = residuals(modelo_beta, type = "deviance"),
pearson = residuals(modelo_beta, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modelo_beta = resid_modelo_beta[resid_modelo_beta$deviance > 2 & resid_modelo_beta$pearson > 2, ]
resul_filt_modelo_beta[1] deviance pearson
<0 rows> (o 0- extensión row.names)La filtración de residuos en resultados_filt_inv_blocks no devolvió filas, lo que indica que no hay residuos que cumplan con las condiciones especificadas (deviance > 2 y pearson > 2). Esto sugiere que, en este modelo, los residuos están dentro de un rango que no se considera problemático según esos umbrales.
6. Una persona está pensando en comprar un carro nuevo. Realiza una búsqueda en línea y recopila información sobre las marcas y modelos que le gustan. Conjetura que las siguientes características del automóvil pueden influir potencialmente en su precio: estilo de carrocería (cupe, hatchback o sedan), país de fabricación (EE. UU., Alemania o Japón), kilometraje en carretera (en mpg), número de puertas (2 o 4) y si el interior es de cuero o no. [?] Los datos de estas variables y el precio (en dólares estadounidenses) se dan a continuación para 27 automóviles.
| obs | bodystyle | country | hwy | doors | leather | price |
|---|---|---|---|---|---|---|
| 1 | coupe | USA | 26 | 4 | no | 17,445 |
| 2 | coupe | USA | 40 | 4 | no | 23,500 |
| 3 | coupe | USA | 35 | 2 | no | 19,600 |
| 4 | coupe | Germany | 37 | 4 | no | 23,400 |
| 5 | coupe | Germany | 25 | 4 | no | 24,100 |
| 6 | coupe | Germany | 24 | 2 | no | 12,400 |
| 7 | coupe | Japan | 26 | 2 | no | 13,300 |
| 8 | coupe | Japan | 27 | 4 | no | 15,550 |
| 9 | coupe | Japan | 20 | 4 | yes | 29,345 |
| 10 | hatchback | USA | 30 | 2 | no | 12,540 |
| 11 | hatchback | USA | 39 | 4 | no | 17,595 |
| 12 | hatchback | USA | 38 | 2 | no | 17,300 |
| 13 | hatchback | Germany | 38 | 4 | no | 17,800 |
| 14 | hatchback | Germany | 32 | 4 | no | 22,500 |
| 15 | hatchback | Germany | 34 | 4 | no | 20,300 |
| 16 | hatchback | Japan | 38 | 4 | yes | 27,300 |
| 17 | hatchback | Japan | 38 | 2 | yes | 23,300 |
| 18 | hatchback | Japan | 38 | 2 | yes | 29,300 |
| 19 | sedan | USA | 29 | 4 | no | 32,000 |
| 20 | sedan | USA | 25 | 2 | yes | 34,200 |
| 21 | sedan | USA | 33 | 4 | yes | 33,395 |
| 22 | sedan | Germany | 40 | 4 | no | 22,850 |
| 23 | sedan | Germany | 23 | 2 | yes | 36,000 |
| 24 | sedan | Germany | 25 | 4 | no | 19,900 |
| 25 | sedan | Japan | 40 | 4 | yes | 36,700 |
| 26 | sedan | Japan | 35 | 4 | yes | 31,600 |
| 27 | sedan | Japan | 37 | 4 | no | 24,600 |
# Crear el dataframe datos_auto
auto <- data.frame(
bodystyle = c("coupe", "coupe", "coupe", "coupe", "coupe", "coupe", "coupe", "coupe", "coupe", "hatchback", "hatchback", "hatchback", "hatchback", "hatchback", "hatchback", "hatchback", "hatchback", "hatchback", "sedan", "sedan", "sedan", "sedan", "sedan", "sedan", "sedan", "sedan", "sedan"),
country = c("USA", "USA", "USA", "Germany", "Germany", "Germany", "Japan", "Japan", "Japan", "USA", "USA", "USA", "Germany", "Germany", "Germany", "Japan", "Japan", "Japan", "USA", "USA", "USA", "Germany", "Germany", "Germany", "Japan", "Japan", "Japan"),
hwy = c(26, 40, 35, 37, 25, 24, 26, 27, 20, 30, 39, 38, 38, 32, 34, 38, 38, 38, 29, 25, 33, 40, 23, 25, 40, 35, 37),
doors = c(4, 4, 2, 4, 4, 2, 2, 4, 4, 2, 4, 2, 4, 4, 4, 4, 2, 2, 4, 2, 4, 4, 2, 4, 4, 4, 4),
leather = c("no", "no", "no", "no", "no", "no", "no", "no", "yes", "no", "no", "no", "no", "no", "no", "yes", "yes", "yes", "no", "yes", "yes", "no", "yes", "no", "yes", "yes", "no"),
price = c(17445, 23500, 19600, 23400, 24100, 12400, 13300, 15550, 29345, 12540, 17595, 17300, 17800, 22500, 20300, 27300, 23300, 29300, 32000, 34200, 33395, 22850, 36000, 19900, 36700, 31600, 24600)
)
# Convertir a variables dummy
auto <- auto %>%
mutate(
style_USA = ifelse(country == "USA", 1, 0),
style_Germany = ifelse(country == "Germany", 1, 0),
style_Japan = ifelse(country == "Japan", 1, 0),
coupe = ifelse(bodystyle == "coupe", 1, 0),
hatchback = ifelse(bodystyle == "hatchback", 1, 0),
sedan = ifelse(bodystyle == "sedan", 1, 0),
leather_yes = ifelse(leather == "yes", 1, 0)
) %>%
select(-bodystyle, -country, -leather) # Eliminar las columnas originales
# Mostrar el dataframe
head(auto) hwy doors price style_USA style_Germany style_Japan coupe hatchback sedan
1 26 4 17445 1 0 0 1 0 0
2 40 4 23500 1 0 0 1 0 0
3 35 2 19600 1 0 0 1 0 0
4 37 4 23400 0 1 0 1 0 0
5 25 4 24100 0 1 0 1 0 0
6 24 2 12400 0 1 0 1 0 0
leather_yes
1 0
2 0
3 0
4 0
5 0
6 0(a) Ajuste el “mejor” modelo de regresión lineal múltiple para predecir el precio un carro. Muestre las medidas de bondad de ajuste.
modelo_auto <- lm(price ~ hwy + doors + style_USA + style_Germany + style_Japan + coupe + hatchback + sedan + leather_yes, data = auto)
summary(modelo_auto)
Call:
lm(formula = price ~ hwy + doors + style_USA + style_Germany +
style_Japan + coupe + hatchback + sedan + leather_yes, data = auto)
Residuals:
Min 1Q Median 3Q Max
-4444 -2434 882 1643 7235
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11546 5488 2.10 0.0489 *
hwy 130 133 0.98 0.3394
doors 1555 790 1.97 0.0638 .
style_USA 3213 1881 1.71 0.1039
style_Germany 3196 2010 1.59 0.1283
style_Japan NA NA NA NA
coupe -4141 1939 -2.14 0.0459 *
hatchback -6411 1845 -3.47 0.0025 **
sedan NA NA NA NA
leather_yes 12176 1933 6.30 0.0000048 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3480 on 19 degrees of freedom
Multiple R-squared: 0.83, Adjusted R-squared: 0.768
F-statistic: 13.3 on 7 and 19 DF, p-value: 0.00000407Seleccionar el “Mejor” Subconjunto de Variables Predictoras
step(modelo_auto)Start: AIC=446.92
price ~ hwy + doors + style_USA + style_Germany + style_Japan +
coupe + hatchback + sedan + leather_yes
Step: AIC=446.92
price ~ hwy + doors + style_USA + style_Germany + style_Japan +
coupe + hatchback + leather_yes
Step: AIC=446.92
price ~ hwy + doors + style_USA + style_Germany + coupe + hatchback +
leather_yes
Df Sum of Sq RSS AIC
- hwy 1 11653050 242159200 446
<none> 230506150 447
- style_Germany 1 30677942 261184092 448
- style_USA 1 35387939 265894089 449
- doors 1 46987090 277493240 450
- coupe 1 55352133 285858283 451
- hatchback 1 146471065 376977215 458
- leather_yes 1 481246617 711752767 475
Step: AIC=446.25
price ~ doors + style_USA + style_Germany + coupe + hatchback +
leather_yes
Df Sum of Sq RSS AIC
<none> 242159200 446
- style_Germany 1 23124134 265283334 447
- style_USA 1 32387985 274547184 448
- doors 1 64320354 306479554 451
- coupe 1 73775069 315934269 451
- hatchback 1 135012795 377171995 456
- leather_yes 1 469611492 711770692 473
Call:
lm(formula = price ~ doors + style_USA + style_Germany + coupe +
hatchback + leather_yes, data = auto)
Coefficients:
(Intercept) doors style_USA style_Germany coupe
15383 1757 3064 2676 -4624
hatchback leather_yes
-5838 11869 mejor_mod_auto = lm(price ~ doors + style_USA + style_Germany + coupe + hatchback + leather_yes, data = auto)
summary(mejor_mod_auto)
Call:
lm(formula = price ~ doors + style_USA + style_Germany + coupe +
hatchback + leather_yes, data = auto)
Residuals:
Min 1Q Median 3Q Max
-5188 -2238 -313 2488 6525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15383 3841 4.00 0.0007 ***
doors 1757 762 2.30 0.0320 *
style_USA 3064 1873 1.64 0.1176
style_Germany 2676 1937 1.38 0.1822
coupe -4624 1873 -2.47 0.0227 *
hatchback -5838 1748 -3.34 0.0033 **
leather_yes 11869 1906 6.23 0.0000044 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3480 on 20 degrees of freedom
Multiple R-squared: 0.822, Adjusted R-squared: 0.768
F-statistic: 15.3 on 6 and 20 DF, p-value: 0.00000152Medidas de Bondad de Ajuste
Coeficiente Correlacion lineal
# Calcular la matriz de correlación
correlation_matrix <- cor(auto[, c("price", "hwy", "doors")], use = "complete.obs")
# Mostrar la matriz de correlación
print(correlation_matrix) price hwy doors
price 1.00000 0.04225 0.16254
hwy 0.04225 1.00000 0.17321
doors 0.16254 0.17321 1.00000Grafico
par(mfrow=c(1,1))
# Cargar el paquete
library(corrplot)corrplot 0.95 loaded# Calcular la matriz de correlación
correlation_matrix <- cor(auto[, c("price", "hwy", "doors")], use = "complete.obs")
# Crear el heatmap
corrplot(correlation_matrix, method = "color", addCoef.col = "black")
Conclusiones
En general, las correlaciones observadas en los datos son débiles, lo que sugiere que el precio de los automóviles no está fuertemente influenciado ni por el kilometraje en carretera ni por el número de puertas. Esto puede ser indicativo de que otros factores no considerados en este análisis pueden tener un impacto más significativo en el precio de los vehículos.
Prueba ANOVA
summary(mejor_mod_auto)
Call:
lm(formula = price ~ doors + style_USA + style_Germany + coupe +
hatchback + leather_yes, data = auto)
Residuals:
Min 1Q Median 3Q Max
-5188 -2238 -313 2488 6525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15383 3841 4.00 0.0007 ***
doors 1757 762 2.30 0.0320 *
style_USA 3064 1873 1.64 0.1176
style_Germany 2676 1937 1.38 0.1822
coupe -4624 1873 -2.47 0.0227 *
hatchback -5838 1748 -3.34 0.0033 **
leather_yes 11869 1906 6.23 0.0000044 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3480 on 20 degrees of freedom
Multiple R-squared: 0.822, Adjusted R-squared: 0.768
F-statistic: 15.3 on 6 and 20 DF, p-value: 0.00000152Ho: El modelo de regresion no se ajusta a los datos. Ha: El modelo de regresion se ajusta a los datos.
F-statistic: 15.3 on 6 and 20 DF, p-value: 0.00000152
Rechace Ho: Existe evidencia estadistica que el modelo de regresion lineal multiple se ajusta a los datos.
Prueba de Coeficientes Estimados
Mejor modelo
summary(mejor_mod_auto)
Call:
lm(formula = price ~ doors + style_USA + style_Germany + coupe +
hatchback + leather_yes, data = auto)
Residuals:
Min 1Q Median 3Q Max
-5188 -2238 -313 2488 6525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15383 3841 4.00 0.0007 ***
doors 1757 762 2.30 0.0320 *
style_USA 3064 1873 1.64 0.1176
style_Germany 2676 1937 1.38 0.1822
coupe -4624 1873 -2.47 0.0227 *
hatchback -5838 1748 -3.34 0.0033 **
leather_yes 11869 1906 6.23 0.0000044 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3480 on 20 degrees of freedom
Multiple R-squared: 0.822, Adjusted R-squared: 0.768
F-statistic: 15.3 on 6 and 20 DF, p-value: 0.00000152Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 15383 3841 4.00 0.0007 doors 1757 762 2.30
0.0320
style_USA 3064 1873 1.64 0.1176
style_Germany 2676 1937 1.38 0.1822
coupe -4624 1873 -2.47 0.0227 *
hatchback -5838 1748 -3.34 0.0033 leather_yes 11869 1906 6.23
0.0000044 ***
Conclusiones
Los coeficientes con significancia estadística sugieren que las variables doors, coupe, hatchback, y leather_yes tienen un impacto significativo en el precio del automóvil, mientras que style_USA y style_Germany no muestran una relación significativa en este modelo.
(b) Estime el precio de un automóvil ensamblado en Japón, tipo sedan, con cuatro puertas, cero kilómetros y con silletería en cuero.
# Crear un nuevo dataframe con las características del automóvil
nuevo_auto <- data.frame(
hwy = 0, # 0 kilómetros en carretera
doors = 4, # Cuatro puertas
style_USA = 0, # No es estilo USA
style_Germany = 0, # No es estilo Alemania
style_Japan = 1, # Es estilo Japón
coupe = 0, # No es coupe
hatchback = 0, # No es hatchback
sedan = 1, # Es un sedan
leather_yes = 1 # Tiene cuero
)
# Predecir el precio usando el modelo
precio_estimado <- predict(mejor_mod_auto, newdata = nuevo_auto)
# Mostrar el precio estimado con cat
cat("El precio estimado del automóvil es:", round(precio_estimado, 2), "USD\n")El precio estimado del automóvil es: 34281 USDConclusión
El precio estimado del automóvil ensamblado en Japón, tipo sedán, con cuatro puertas, cero kilómetros y con sillones de cuero es de 34,281 USD. Este precio refleja las características específicas del vehículo, donde se considera que la combinación de un diseño de sedán y el uso de materiales de cuero, junto con su procedencia de Japón.
(c) Compruebe los supuestos del modelo.
Supuesto Independiencia de los residuos
Ho: Los residuos son independientes Ha: Los residuos No son independientes
durbinWatsonTest(modelo_auto) lag Autocorrelation D-W Statistic p-value
1 -0.053477 2.0558 0.48
Alternative hypothesis: rho != 0no rechazamos la hipótesis nula, los residuos son independientes
(d) ¿Existen Datos atípicos o influyentes?
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelo_auto)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_mod_auto,type = 'pearson')))
abline(h=2,col = 'red')
Se observa 1 valor atipico en el modelo inicial (aprox 20 en ambos modelos, sin embargo no afectan al modelo.
Residuos
resid_modelo_auto <- data.frame(
deviance = residuals(modelo_auto, type = "deviance"),
pearson = residuals(modelo_auto, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modelo_auto = resid_modelo_auto[resid_modelo_auto$deviance > 2 & resid_modelo_auto$pearson > 2, ]
resul_filt_modelo_auto deviance pearson
2 1440.43 1440.43
3 1303.78 1303.78
4 1748.88 1748.88
5 4014.86 4014.86
9 932.57 932.57
12 882.05 882.05
14 3771.14 3771.14
15 1310.14 1310.14
18 3919.22 3919.22
19 7235.02 7235.02
20 892.21 892.21
23 2970.16 2970.16
25 1536.71 1536.71
27 2003.87 2003.87Conclusiones
Los valores de devianza de Pearson para los diferentes casos oscilan considerablemente, desde un mínimo de 882.05 hasta un máximo de 7235.02. Esto indica una variabilidad significativa en la adecuación del modelo a diferentes observaciones.
Observaciones con devianzas significativamente altas (como la 7235.02 del caso 19) pueden ser indicativas de outliers o puntos influyentes que podrían distorsionar los resultados del modelo.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_mod_auto)
StudRes Hat CookD
6 -1.61897 0.29467 0.144700
19 2.40766 0.24800 0.220273
23 0.94527 0.39911 0.085237
24 -1.77229 0.21644 0.111964
27 0.81027 0.40800 0.0657697. En una clínica infantil se llevó a cabo un estudio de intervención y control sobre la obesidad infantil. Una cohorte de 36 niños obesos, de 6 a 16 años, fueron seguidos durante 9 meses. La intervención consistió en sesiones educativas. para padres y actividades de ejercicio vigoroso para niños. Los participantes del grupo de control recibieron cursos sobre otros estilos de vida activos y saludables. Se tiene las variables género (M/F), edad (en años), grupo (Tx de intervención o Cx de control) y los IMC pre y post se muestran en la tabla a continuación:
| sujeto | Gender | Age | Group | Pre IMC | Post IMC |
|---|---|---|---|---|---|
| 1 | F | 6 | Cx | 85.70 | 83.80 |
| 2 | F | 6 | Cx | 93.80 | 92.90 |
| 3 | F | 7 | Cx | 93.50 | 92.50 |
| 4 | F | 8 | Cx | 90.10 | 89.80 |
| 5 | F | 9 | Tx | 92.30 | 90.70 |
| 6 | F | 9 | Tx | 90.30 | 88.30 |
| 7 | F | 12 | Cx | 87.60 | 85.90 |
| 8 | F | 12 | Cx | 87.20 | 84.10 |
| 9 | F | 12 | Tx | 96.90 | 94.90 |
| 10 | F | 12 | Tx | 85.80 | 81.20 |
| 11 | F | 13 | Cx | 96.70 | 94.10 |
| 12 | F | 13 | Cx | 93.50 | 92.90 |
| 13 | F | 13 | Tx | 92.30 | 87.50 |
| 14 | F | 13 | Tx | 85.30 | 83.70 |
| 15 | F | 14 | Tx | 95.50 | 78.70 |
| 16 | F | 15 | Cx | 91.30 | 89.90 |
| 17 | F | 15 | Tx | 95.80 | 87.10 |
| 18 | F | 16 | Tx | 90.70 | 87.20 |
| 19 | M | 6 | Cx | 92.60 | 88.10 |
| 20 | M | 7 | Cx | 95.80 | 94.70 |
| 21 | M | 7 | Cx | 90.40 | 89.10 |
| 22 | M | 7 | Cx | 91.20 | 88.60 |
| 23 | M | 8 | Tx | 94.40 | 87.80 |
| 24 | M | 8 | Tx | 93.20 | 87.30 |
| 25 | M | 10 | Cx | 93.90 | 91.50 |
| 26 | M | 10 | Tx | 96.20 | 91.10 |
| 27 | M | 10 | Tx | 89.40 | 87.90 |
| 28 | M | 11 | Tx | 86.20 | 77.10 |
| 29 | M | 11 | Tx | 95.40 | 84.80 |
| 30 | M | 12 | Cx | 97.70 | 95.80 |
| 31 | M | 13 | Tx | 85.30 | 80.00 |
| 32 | M | 13 | Tx | 86.20 | 82.40 |
| 33 | M | 14 | Cx | 85.50 | 83.60 |
| 34 | M | 14 | Cx | 97.80 | 93.80 |
| 35 | M | 16 | Cx | 95.00 | 93.60 |
| 36 | M | 16 | Tx | 93.10 | 86.80 |
obesidad_infantil = data.frame(Gender = c("F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M"),
Age = c(6, 6, 7, 8, 9, 9, 12, 12, 12, 12, 13, 13, 13, 13, 14, 15, 15, 16, 6, 7, 7, 7, 8, 8, 10, 10, 10, 11, 11, 12,
13, 13, 14, 14, 16, 16),
Group = c("Cx", "Cx", "Cx", "Cx", "Tx", "Tx", "Cx", "Cx", "Tx", "Tx", "Cx", "Cx", "Tx", "Tx", "Tx", "Cx", "Tx", "Tx", "Cx", "Cx", "Cx", "Cx", "Tx", "Tx", "Cx", "Tx", "Tx", "Tx", "Tx", "Cx", "Tx", "Tx", "Cx", "Cx", "Cx", "Tx"),
Pre_IMC = c(85.70, 93.80, 93.50, 90.10, 92.30, 90.30, 87.60, 87.20, 96.90, 85.80, 96.70, 93.50, 92.30, 85.30, 95.50, 91.30, 95.80, 90.70, 92.60, 95.80, 90.40, 91.20, 94.40, 93.20, 93.90, 96.20, 89.40, 86.20, 95.40, 97.70, 85.30, 86.20, 85.50, 97.80, 95.00, 93.10),
Post_IMC = c(83.80, 92.90, 92.50, 89.80, 90.70, 88.30, 85.90, 84.10, 94.90, 81.20, 94.10, 92.90, 87.50, 83.70, 78.70, 89.90, 87.10, 87.20, 88.10, 94.70, 89.10, 88.60, 87.80, 87.30, 91.50, 91.10, 87.90, 77.10, 84.80, 95.80, 80.00, 82.40, 83.60, 93.80, 93.60, 86.80)
)
# Convertir a variables dummy
obesidad_infantil <- obesidad_infantil %>%
mutate(
Gender_M = ifelse(Gender == "M", 1, 0),
Group_Cx = ifelse(Group == "Cx", 1, 0)
) %>%
select(-Gender, -Group) # Eliminar las columnas originales
# Calcular diferencia IMC
obesidad_infantil$Dif_IMC <- obesidad_infantil$Post_IMC - obesidad_infantil$Pre_IMC
head(obesidad_infantil) Age Pre_IMC Post_IMC Gender_M Group_Cx Dif_IMC
1 6 85.7 83.8 0 1 -1.9
2 6 93.8 92.9 0 1 -0.9
3 7 93.5 92.5 0 1 -1.0
4 8 90.1 89.8 0 1 -0.3
5 9 92.3 90.7 0 0 -1.6
6 9 90.3 88.3 0 0 -2.0(a) Calcule la diferencia en el IMC antes y después y ajuste un modelo de regresión para estimar la diferencia en el IMC en términos de la edad, el género y el grupo. Muestre las medidas de bondad de ajuste.
# Ajustar modelo de regresión
modelo_obesidad <- lm(Dif_IMC ~ Age + Gender_M + Group_Cx, data=obesidad_infantil)Seleccionar el “Mejor” Subconjunto de Variables Predictoras
step(modelo_obesidad)Start: AIC=79.89
Dif_IMC ~ Age + Gender_M + Group_Cx
Df Sum of Sq RSS AIC
- Age 1 7.1 272 78.8
- Gender_M 1 8.9 274 79.1
<none> 265 79.9
- Group_Cx 1 96.6 362 89.1
Step: AIC=78.85
Dif_IMC ~ Gender_M + Group_Cx
Df Sum of Sq RSS AIC
- Gender_M 1 7.3 280 77.8
<none> 272 78.8
- Group_Cx 1 118.1 390 89.8
Step: AIC=77.8
Dif_IMC ~ Group_Cx
Df Sum of Sq RSS AIC
<none> 280 77.8
- Group_Cx 1 118 398 88.5
Call:
lm(formula = Dif_IMC ~ Group_Cx, data = obesidad_infantil)
Coefficients:
(Intercept) Group_Cx
-5.54 3.62 mejor_modelo_obesidad = lm(Dif_IMC ~ Age + Gender_M + Group_Cx, data=obesidad_infantil)
summary(mejor_modelo_obesidad)
Call:
lm(formula = Dif_IMC ~ Age + Gender_M + Group_Cx, data = obesidad_infantil)
Residuals:
Min 1Q Median 3Q Max
-11.429 -0.830 0.445 1.120 4.268
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.261 2.144 -1.52 0.1380
Age -0.151 0.162 -0.93 0.3606
Gender_M -1.000 0.966 -1.04 0.3080
Group_Cx 3.388 0.992 3.41 0.0018 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.88 on 32 degrees of freedom
Multiple R-squared: 0.333, Adjusted R-squared: 0.271
F-statistic: 5.33 on 3 and 32 DF, p-value: 0.0043Conclusiones
Se analizó la diferencia en el IMC en función de la edad, el género y el grupo (intervención/control).
El Grupo (Tx) con Significativo (p = 0.0018) tuvo una reducción notable en el IMC. La Edad y el Género son no significativos (p = 0.360 y p = 0.308, respectivamente). El valor R² con 0.333 explica el 33.3% de la variabilidad en la diferencia de IMC. El Error estándar residual es 2.88.
La intervención fue efectiva, mientras que la edad y el género no influyeron en los resultados.
Medidas de Bondad de Ajuste
Coeficiente Correlacion lineal
# Calcular la matriz de correlación
correlation_matrix_obesidad <- cor(obesidad_infantil[, c("Age", "Pre_IMC", "Post_IMC", "Gender_M", "Group_Cx", "Dif_IMC")], use = "complete.obs"); correlation_matrix_obesidad Age Pre_IMC Post_IMC Gender_M Group_Cx Dif_IMC
Age 1.00000 0.01660 -0.1641707 -0.1084830 -0.25313 -0.25195
Pre_IMC 0.01660 1.00000 0.7168567 0.1067577 0.10676 -0.15921
Post_IMC -0.16417 0.71686 1.0000000 -0.0070827 0.47336 0.57420
Gender_M -0.10848 0.10676 -0.0070827 1.0000000 0.00000 -0.13539
Group_Cx -0.25313 0.10676 0.4733588 0.0000000 1.00000 0.54490
Dif_IMC -0.25195 -0.15921 0.5741990 -0.1353898 0.54490 1.00000Grafico
par(mfrow=c(1,1))
# Cargar el paquete
library(corrplot)
# Calcular la matriz de correlación
correlation_matrix_obesidad <- cor(obesidad_infantil, use = "complete.obs")
# Configurar la ventana gráfica
par(mfrow=c(1,1))
# Crear el heatmap
corrplot(correlation_matrix_obesidad, method = "color", addCoef.col = "black")
Conclusiones
Existe una correlación positiva fuerte (0.717) entre el IMC previo (Pre_IMC) y el IMC posterior (Post_IMC), indicando que a mayor IMC inicial, es probable que también haya un IMC alto al final del estudio.
La variable Dif_IMC muestra una correlación positiva moderada con Post_IMC (0.574) y con Group_Cx (0.545). Esto sugiere que el cambio en el IMC está relacionado tanto con el IMC posterior como con el grupo de intervención, lo que indica que el grupo de intervención podría tener un impacto en la reducción del IMC.
La edad (Age) tiene correlaciones débiles y negativas con la mayoría de las variables, indicando que la edad no es un factor determinante en las diferencias de IMC en este estudio. El género (Gender_M) también presenta correlaciones bajas, lo que sugiere que no hay una influencia del género en los resultados del IMC.
Hay una correlación positiva moderada con Dif_IMC (0.545) y Post_IMC (0.473), lo que implica que el grupo de intervención mostró mejores resultados en la reducción del IMC en comparación con el grupo de control.
Prueba ANOVA
summary(mejor_modelo_obesidad)
Call:
lm(formula = Dif_IMC ~ Age + Gender_M + Group_Cx, data = obesidad_infantil)
Residuals:
Min 1Q Median 3Q Max
-11.429 -0.830 0.445 1.120 4.268
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.261 2.144 -1.52 0.1380
Age -0.151 0.162 -0.93 0.3606
Gender_M -1.000 0.966 -1.04 0.3080
Group_Cx 3.388 0.992 3.41 0.0018 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.88 on 32 degrees of freedom
Multiple R-squared: 0.333, Adjusted R-squared: 0.271
F-statistic: 5.33 on 3 and 32 DF, p-value: 0.0043Ho: El modelo de regresion no se ajusta a los datos. Ha: El modelo de regresion se ajusta a los datos.
F-statistic: 5.33 on 3 and 32 DF, p-value: 0.0043
Rechace Ho: Existe evidencia estadistica que el modelo de regresion lineal multiple se ajusta a los datos.
Prueba de Coeficientes Estimados
Mejor modelo
summary(mejor_modelo_obesidad)
Call:
lm(formula = Dif_IMC ~ Age + Gender_M + Group_Cx, data = obesidad_infantil)
Residuals:
Min 1Q Median 3Q Max
-11.429 -0.830 0.445 1.120 4.268
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.261 2.144 -1.52 0.1380
Age -0.151 0.162 -0.93 0.3606
Gender_M -1.000 0.966 -1.04 0.3080
Group_Cx 3.388 0.992 3.41 0.0018 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.88 on 32 degrees of freedom
Multiple R-squared: 0.333, Adjusted R-squared: 0.271
F-statistic: 5.33 on 3 and 32 DF, p-value: 0.0043Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.261 2.144 -1.52 0.1380
Age -0.151 0.162 -0.93 0.3606
Gender_M -1.000 0.966 -1.04 0.3080
Group_Cx 3.388 0.992 3.41 0.0018 **
Conclusiones
La variable “Edad” muestra un coeficiente de -0.151 con un valor p de 0.3606, lo que indica que no hay evidencia suficiente para afirmar que la edad tiene un efecto significativo en la diferencia de IMC.
Para el “Género_M”, el coeficiente es -1.000 con un p = 0.3080. Esto implica que no se encontró un efecto significativo del género en la reducción del IMC, lo que sugiere que el género no influye de manera importante en los resultados.
El coeficiente para el “Grupo_Cx” es 3.388 y tiene un p valor de 0.0018, lo que indica un efecto altamente significativo. Esto sugiere que el grupo de intervención logró una reducción notable en el IMC en comparación con el grupo de control.
El grupo de intervención mostró una mejora significativa en el IMC, mientras que la edad y el género no presentaron una influencia estadísticamente significativa en la diferencia de IMC entre los grupos. Esto resalta la efectividad de la intervención en el tratamiento de la obesidad infantil.
(b) Estime la diferencia esperada de una niña de 10 años de edad que pertenece al grupo control.
# Crear el nuevo dataframe para la predicción
nueva_obs <- data.frame(Age = 10, Gender_M = 0, Group_Cx = 0)
# Predecir la diferencia esperada
dif_esp <- predict(mejor_modelo_obesidad, newdata = nueva_obs)
# Mostrar el resultado con cat
cat("La diferencia esperada de IMC para una niña de 10 años en el grupo control es:", round(dif_esp, 2), "\n")La diferencia esperada de IMC para una niña de 10 años en el grupo control es: -4.77 La diferencia esperada de IMC para una niña de 10 años en el grupo control es de aproximadamente -4.77. Esto indica que, en promedio, se espera que su IMC disminuya en esa cantidad en comparación con su valor inicial.
(c) ¿Existen Datos atípicos o influyentes? Muestre los resultados.
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(modelo_obesidad)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_modelo_obesidad,type = 'pearson')))
abline(h=2,col = 'red')
Se observa varios valores atipicos tanto en el modelo inicial como en el modelo final en el rango de 5 - 33 aprox, por lo que hay que obervar con mas detenimiento el modelo.
Residuos
resid_modelo_obesidad <- data.frame(
deviance = residuals(modelo_obesidad, type = "deviance"),
pearson = residuals(modelo_obesidad, type = "pearson")
)
# Filtrar los residuos donde la devianza > 2 y los residuos de Pearson > 2
resul_filt_modelo_obesidad = resid_modelo_obesidad[resid_modelo_obesidad$deviance > 2 & resid_modelo_obesidad$pearson > 2, ]
resul_filt_modelo_obesidad deviance pearson
5 3.0174 3.0174
6 2.6174 2.6174
9 3.0693 3.0693
14 3.6200 3.6200
18 2.1720 2.1720
27 4.2685 4.2685
32 2.4204 2.4204Conclusiones
Los valores de devianza de Pearson obtenidos para los distintos casos oscilan entre 2.1720 y 4.2685. Estos valores indica la discrepancia entre los datos observados y los datos predichos por el modelo.
Valores cercanos a 1 sugieren un buen ajuste del modelo a los datos, mientras que valores más altos pueden indicar que hay discrepancias. En este caso, los valores superiores a 3 pueden ser motivo de preocupación y sugieren que los puntos correspondientes podrían ser considerados como outliers o influencias no deseadas.
Los casos con valores de devianza más altos, como 4.2685 (en el caso 27), podrían requerir un análisis más profundo para determinar si están influyendo en el ajuste del modelo.
La presencia de estos valores sugiere que, si bien el modelo puede ser adecuado en general, hay ciertas observaciones que podrían no estar siendo bien capturadas por el modelo y que podrían afectar la validez de las Conclusiones que se extraen sus resultados.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_modelo_obesidad)
StudRes Hat CookD
1 -0.41768 0.151024 0.0079638
15 -6.07772 0.094034 0.4514666
29 -1.75322 0.084129 0.0662913
35 0.72651 0.200075 0.0334978La observación 15 no cumple la primera regla (residuos estudentizados), sin embargo las otras observaciones (1, 29, 35) cumplen con las tres reglas de diagnóstico, mientras que la observación 15 se considera potencialmente influyente.