Hola soy Alexander y este es mi proyecto final del curso de R. Yo soy estudiante de ingeniería Ambiental de la EPN, escogi esta carrera por el deseo inefable de cuidar el medio ambiente, aprender a explotar los recursos naturales de forma sostenible, colaborar con el cumplimiento de la legislación vigente y por la necesidad imperante de solucionar los graves problemas ambientales.
En el presente trabajo se va a analizar una base de datos que contienen información relacionada con las variables (Edad, Colesterol, Índice de Masa Corporal, Tensión Arterial Diastólica y Género) de 70 pacientes. A continuación se definen las principales variables encontradas en la base de datos seleccionada, con el fin de tener un conocimiento previo del tema y asi estimar cual es la variable independiente y dependiente para el análisis.
| Nombre | Descripcion |
|---|---|
| Colesterol | Es una sustancia cerosa y parecida a la grasa que se encuentra en todas las células de su cuerpo. Su hígado produce colesterol. También se encuentra en algunos alimentos, como la carne y los productos lácteos. Su cuerpo necesita algo de colesterol para funcionar bien. Pero si tiene demasiado colesterol en la sangre, tiene un mayor riesgo de enfermedad arterial coronaria. |
| IMC | Es un número que se calcula con base en el peso y la estatura de la persona. El IMC es un indicador de la gordura bastante confiable para la mayoría de las personas. El IMC se usa como una herramienta de detección para identificar posibles problemas de salud de los adultos |
| Tensión Arterial Diastólica | La tensión arterial diastólica o (la baja), es la presión que la sangre ejerce cuando el corazón se relaja para volver a llenarse de sangre. La tensión arterial diastólica suele aumentar hasta los 50 años y a partir de esa edad tiende a disminuir. |
Lectura de la base de datos
library(readxl)
tension <- read_excel("~/Curso_R_muy_principiante/Proyecto final/pacientes.xlsx",
skip = 2)
tension$PACIENTE <- as.character(tension$PACIENTE)
row.names(tension) <- tension$PACIENTE
tension$EDAD <- as.integer(tension$EDAD)
tension$COLESTEROL <- as.integer(tension$COLESTEROL)
tension$IMC <- as.double(tension$IMC)
tension$TAD <- as.integer(tension$TAD)
tension$GENERO <- as.factor(tension$GENERO)
kable(head(tension), align = "c") %>%
kable_styling(bootstrap_options = c("striped","bordered"), full_width = F, position = "center")
| PACIENTE | EDAD | COLESTEROL | IMC | TAD | GENERO |
|---|---|---|---|---|---|
| 1 | 42 | 292 | 31.64 | 97 | Hombre |
| 2 | 64 | 235 | 30.80 | 90 | Hombre |
| 3 | 47 | 200 | 25.61 | 80 | Hombre |
| 4 | 56 | 200 | 26.17 | 75 | Mujer |
| 5 | 54 | 300 | 31.96 | 100 | Hombre |
| 6 | 48 | 215 | 23.18 | 67 | Hombre |
La información copilada en la tabla anterior muestra la siguiente estructura.
str(tension)
## Classes 'tbl_df', 'tbl' and 'data.frame': 70 obs. of 6 variables:
## $ PACIENTE : chr "1" "2" "3" "4" ...
## $ EDAD : int 42 64 47 56 54 48 57 52 67 46 ...
## $ COLESTEROL: int 292 235 200 200 300 215 216 254 310 237 ...
## $ IMC : num 31.6 30.8 25.6 26.2 32 ...
## $ TAD : int 97 90 80 75 100 67 NA 70 105 70 ...
## $ GENERO : Factor w/ 2 levels "Hombre","Mujer": 1 1 1 2 1 1 2 1 1 2 ...
summary(tension)
## PACIENTE EDAD COLESTEROL IMC
## Length:70 Min. :42.00 Min. :175.0 Min. :19.10
## Class :character 1st Qu.:49.00 1st Qu.:214.2 1st Qu.:22.36
## Mode :character Median :56.00 Median :230.0 Median :25.38
## Mean :55.24 Mean :236.8 Mean :25.47
## 3rd Qu.:60.00 3rd Qu.:254.0 3rd Qu.:27.81
## Max. :68.00 Max. :315.0 Max. :33.91
## NA's :2
## TAD GENERO
## Min. : 65.00 Hombre:41
## 1st Qu.: 75.00 Mujer :29
## Median : 80.00
## Mean : 81.65
## 3rd Qu.: 90.00
## Max. :105.00
## NA's :1
Gracias a la función summary observamos un resumen de cada varible de la base de datos. Se tiene 4 variables numericas y una variable categórica. Tambien nos indica cuantos valores NA existen en cada varible.
edad_est <- c(paste("Minimo: ", min(tension$EDAD)),
paste("Media: ", round(mean(tension$EDAD),2)),
paste("Máximo: ", max(tension$EDAD)),
paste("Desviación Estandar: ", round(sd(tension$EDAD),2)),
paste("Primer Cuartil :", quantile(tension$EDAD,prob=0.25)))
edad_est
## [1] "Minimo: 42" "Media: 55.24"
## [3] "Máximo: 68" "Desviación Estandar: 7.09"
## [5] "Primer Cuartil : 49"
La variable Edad indica que los pacientes estan entre los 42 y 68 años.
colesterol_est <- c(paste("Minimo: ", min(tension$COLESTEROL)),
paste("Media: ", round(mean(tension$COLESTEROL),2)),
paste("Máximo: ", max(tension$COLESTEROL)),
paste("Desviación Estandar: ", round(sd(tension$COLESTEROL),2)),
paste("Primer Cuartil :", quantile(tension$COLESTEROL,prob=0.25)))
colesterol_est
## [1] "Minimo: 175" "Media: 236.77"
## [3] "Máximo: 315" "Desviación Estandar: 34.6"
## [5] "Primer Cuartil : 214.25"
La variable Colesterol, indica que el nivel de colesterol en los pacientes es en promedio 236.27 mg/dl que sobrepasa el nivel normal de colesterol normal que es menos de 200 mg/dl (Leer más)
imc_est <-c(paste("Minimo: ", min(tension$IMC, na.rm = TRUE)),
paste("Media: ", round(mean(tension$IMC, na.rm = TRUE),2)),
paste("Máximo: ", max(tension$IMC, na.rm = TRUE)),
paste("Desviación Estandar: ", round(sd(tension$IMC, na.rm = TRUE),2)),
paste("Primer Cuartil :", round(quantile(tension$IMC, prob=0.25, na.rm = TRUE),2))
)
imc_est
## [1] "Minimo: 19.1" "Media: 25.47"
## [3] "Máximo: 33.91" "Desviación Estandar: 3.96"
## [5] "Primer Cuartil : 22.36"
La media del IMC de los 70 pacientes es de 25.54 que está dentro de los valores correspondientes a “sobrepeso”. (Leer más)
tad_est <- c(
paste("Minimo: ", min(tension$TAD, na.rm = TRUE)),
paste("Media: ", round(mean(tension$TAD, na.rm = TRUE),2)),
paste("Máximo: ", max(tension$TAD, na.rm = TRUE)),
paste("Desviación Estandar: ", round(sd(tension$TAD, na.rm = TRUE),2)),
paste("Primer Cuartil :", quantile(tension$TAD, prob=0.25, na.rm = TRUE))
)
tad_est
## [1] "Minimo: 65" "Media: 81.65"
## [3] "Máximo: 105" "Desviación Estandar: 11.32"
## [5] "Primer Cuartil : 75"
La presión arterial diástolica normal es menor de 80mmHg, para este caso se encontro pacientes que superan en mas de 20 mmHg los valores normales. (Leer más)
table(tension$GENERO)
##
## Hombre Mujer
## 41 29
El género al tratarse de una variable finita de dos levels: Hombre y Mujeres. Solo se calculo su frecuencia, encontrandose que hay mas hombres en la base de datos.
boxplot(tension$IMC~tension$GENERO, col = c("yellow", "blue"), ylab = "IMC" )
Los diagramas de caja indica que el promedio de IMC en hombre es ligeramnete mayor que en las mujeres, en cuanto a la dispersion y desviacion de datos, los pacientes hombres tienen mejor distribucion que las mujeres es mas cercana a una distribución normal. La distancia entre cuantiles, media, minimos y maximos en hombres es mu similar mientras que en mujeres no lo es.
ggplot(data=tension, aes(x=tension$PACIENTE, y=tension$EDAD)) +
geom_bar(stat="identity") +
labs(title="EDAD") +
xlab("Paciente") +
ylab("EDAD") +
scale_fill_manual(values=c('green','yellow'))
La edad es una variable discreta, por lo que se grafico la edad de cada paciente en un diagrama de barras, observandose que la mayor parte de pacientes tienen edades inferiores a los 50 años.
ggplot(data=tension, aes(x=tension$PACIENTE, y=tension$COLESTEROL)) +
geom_bar(stat="identity") +
labs(title="COLESTEROL") +
xlab("Paciente") +
ylab("COLESTEROL") +
scale_fill_manual(values=c('black','lightgray'))
El nivel de colesterol de la mayoria de pacientes esta entre el primer y tercer cuantil. Se observa que hay un porcentaje similar entre pacientes que superan y tienenn un nivel normal de colesterol.
ggplot(data=tension, aes(x=tension$PACIENTE , y=tension$TAD)) +
geom_bar(stat="identity") +
labs(title="TENSION ARTERIAL DIASTOLICA") +
xlab("Paciente") +
ylab("TAD") +
scale_fill_manual(values=c('black','lightgray'))
De acuerdo al diagrama de barras del TAD son pocos los caso en pacientes que superan el valor normal de presion arterial diastólica (80 mmHg)
Para limpiar la información. Se utilizó missmap () del paquete de Amelia. Como lo visto en clases, previamente con summary se demostró que hay varios NA en algunas variables tal como se comprobó con missmap, luego de ello se eliminó los NA usando la funcion na.omit()
Compruebe si hay NA en el marco de datos.
missmap(tension,col=c('yellow','black'),y.at=1,y.labels='',legend=TRUE)
tension <- na.omit(tension) #ELIMINACIÓN DE LOS NA
La correlación es cualquiera de una amplia clase de relaciones estadísticas que involucran dependencia, aunque en el uso común a menudo se refiere a la medida en que dos variables tienen una relación lineal entre sí.
Los diagramas de correlación son una gran forma de explorar datos y ver si hay algún término de interacción.
corrplot(cor(select(tension,-GENERO,-PACIENTE)))
Interpretación “TAD es la variale dependiente”. Existe correlacion lineal directa entre TAD y las variables edad, colesterol y TAD. El TAD aumenta con el aumento en edad (bajo), colesterol (medio-alto), IMC (medio).
tension %>%
ggplot(aes(TAD)) +
stat_density() +
theme_bw()
Las visualizaciones anteriores revelan que las densidades máximas de TAD están entre 70 y 80.
tension %>%
select(c(EDAD, COLESTEROL, IMC, TAD)) %>%
melt(id.vars = "TAD") %>%
ggplot(aes(x = value, y = TAD, colour = variable)) +
geom_point(alpha = 0.7) +
stat_smooth(aes(colour = "black")) +
facet_wrap(~variable, scales = "free", ncol = 2) +
labs(x = "Variable Value", y = "Nivel") +
theme_minimal()
Los resultados del gráfico anterior están en correlación con el corrplot. Sin embargo se observa que hay una alta dispersion de los puntos respecto a la correlaciòn con TAD.
Para el analisis ANOVA se escogio la variable dependiente TAD y una variable categórica que en este caso corresponde a la variable GENERO.
Llamar la funcion ANOVA
fm = aov( lm(tension$TAD ~ tension$GENERO) )
fm
## Call:
## aov(formula = lm(tension$TAD ~ tension$GENERO))
##
## Terms:
## tension$GENERO Residuals
## Sum of Squares 120.320 8405.799
## Deg. of Freedom 1 65
##
## Residual standard error: 11.37189
## Estimated effects may be unbalanced
Resumen:
summary(fm)
## Df Sum Sq Mean Sq F value Pr(>F)
## tension$GENERO 1 120 120.3 0.93 0.338
## Residuals 65 8406 129.3
Delimitación de la región de aceptación y rechazo.
qf(0.05, 0, 70-1, lower.tail = F)
## [1] NaN
Valores del estadístico qf es NaN (“Not a Number”) por lo que se desconoce la region de rechazo. En nuetro caso 0.93 es mucho menor que el valor crítico obtenido. Sin embargo, se puede deducir que no existe buena correlacion entre el genero y el nivel de TAD.
media <- mean(tension$TAD[tension$GENERO == "Hombre"])
valor_t <- pt(0.05/2, 70 - 1)
sp <- sqrt(120) #desviación típica de la varianza muestral común
ee <- valor_t * (sp/ sqrt(2)) #error de estimación
media
## [1] 82.53846
Límite superior e inferior del intervalo de confianza de la media de TAD en los HOMBRES.
La raíz cuadrada de la media de los cuadrados del error, además de proporcionarnos una estimación de la varianza muestral de todos los datos, se utiliza en la obtención de los intervalos de confianza de las medias en cada uno de los grupos de interés.
Por ejemplo, este sería el intervalo de confianza de la media del TAD en hombres, con un nivel de confianza del 95%
Límite superior e inferior del intervalo de confianza de la media de TAD en los HOMBRES.
Limite_Superior <- media + ee
Limite_Inferior <- media - ee
limite <- data.frame(Limite_Superior,Limite_Inferior)
kable(limite, align = "c") %>%
kable_styling(bootstrap_options = c("striped","bordered"), full_width = F, position = "center")
| Limite_Superior | Limite_Inferior |
|---|---|
| 86.48841 | 78.58851 |
Si hemos detectado diferencias significativas entre las medias de las poblaciones. ¿Sería posible saber cuáles son los grupos que generan estas diferencias?
intervals = TukeyHSD(fm)
intervals
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm(tension$TAD ~ tension$GENERO))
##
## $`tension$GENERO`
## diff lwr upr p adj
## Mujer-Hombre -2.717033 -8.342609 2.908543 0.3383328
plot(intervals)
Como el grafico contienen a cero son estadisticamente iguales (TAD Y Genero)
plot(fm$residuals)
summary(fm$residuals)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -17.538 -8.680 -2.538 0.000 10.179 22.462
boxplot(fm$residuals)
No existen datos dispersos.
hist(fm$residuals)
qqnorm(fm$residuals)
qqline(fm$residuals)
Los datos no tienden a una distribución normal, debido a que p-value en menor 0.05.
shapiro.test(fm$residuals)
##
## Shapiro-Wilk normality test
##
## data: fm$residuals
## W = 0.94407, p-value = 0.004598
Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:
boxplot(fm$residuals~tension$GENERO, col = c("yellow", "blue"))
Tienen diferente ancho de cajas por lo que no hay igualdad de varinzas entre las categorias de la variable género.
desviaciones <- tapply(fm$residuals, tension$GENERO, sd)
max(desviaciones) / min(desviaciones)
## [1] 1.201014
Comparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente) y si es menor o igual a 2 si hay homocedasticidad. En este caso se tiene un valor de 1.201014 lo que significa que si hay homocedasticidad.
bartlett.test(fm$residuals ~ tension$GENERO)
##
## Bartlett test of homogeneity of variances
##
## data: fm$residuals by tension$GENERO
## Bartlett's K-squared = 1.0159, df = 1, p-value = 0.3135
De acuerdo, al indice de Bartlett no hay homocedasticidad ya que p-value es mayor a 0.05.
A partir de los residuos del modelo s comprobo que el modelo ANOVA no es adecuado. Por que no se cuumplierón los tres supuestos: independencia, homocedasticidad y normalidad.
set.seed(123)
split <- sample.split(tension,SplitRatio =0.75)
train <- subset(tension, split==TRUE)
test <- subset(tension,split==FALSE)
Retire las variables EDAD e intercepto por no fueron significativas en la primera corrida del modelo.
model <- lm(TAD ~ -1 + COLESTEROL + IMC, data = train)
summary(model)
##
## Call:
## lm(formula = TAD ~ -1 + COLESTEROL + IMC, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.3547 -6.4538 0.1365 6.2741 19.6198
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## COLESTEROL 0.20959 0.03746 5.595 1.42e-06 ***
## IMC 1.30311 0.34981 3.725 0.000564 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.555 on 43 degrees of freedom
## Multiple R-squared: 0.9902, Adjusted R-squared: 0.9898
## F-statistic: 2182 on 2 and 43 DF, p-value: < 2.2e-16
Las variables COLESTEROL y IMC tienen la menor porbabilidad no tener buena correlación.
Permite visualizar nuestro modelo de regresión lineal trazando los residuos. La diferencia entre el valor observado de la variable dependiente (y) y el valor predicho (y) se denomina residual (e).
res <- residuals(model)
# Convertir residuos en un DataFrame
#
res <- as.data.frame(res)
Un modelo es aceptable si: * Las variables son estadisticamente significativas * Si el R^2>=0.5 * Residuos que se ajusten a una N(0,1) distribucion normal de 0 a 1, (Varianza=1) * Residuos no correlacionados, deben estar diversos.
ggplot(res,aes(res)) + geom_histogram(fill='blue',alpha=0.5)
No tiene se asemaeja a una distribucion normal.
plot(model)
Probemos nuestro modelo prediciendo en nuestro conjunto de datos de prueba.
test$predicted.TAD <- predict(model,test)
pl1 <-test %>%
ggplot(aes(TAD,predicted.TAD)) +
geom_point(alpha=0.5) +
stat_smooth(aes(colour='black')) +
xlab('VALOR ACTUAL DE TAD') +
ylab('VALOR ESTIMADO DE TAD')+
theme_bw()
ggplotly(pl1)
Usando Root Mean Square Error, una medida estandarizada de cuán lejos estábamos con nuestros valores predichos.
ee_modelo <- test$TAD-test$predicted.TAD #Calculamos el error entre si es cercano a uno es un mododelo bueno
rmse <- sqrt(mean(ee_modelo)^2)
Conclusión El Root Mean Square Error (RMSE) para nuestro modelo es 5.100186 y los resultados pueden mejorarse aún más utilizando la extracción de variables y entrenando el modelo.