R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

#EXPLORACION DE DATOS 

library(readxl)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(splines)
library(mgcv)
## Cargando paquete requerido: nlme
## 
## Adjuntando el paquete: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
library(ggplot2)
library(gratia)
## Warning: package 'gratia' was built under R version 4.5.3
library(tree)
## Warning: package 'tree' was built under R version 4.5.3
datose <- read_xlsx("salud_cardiovascular.xlsx")
head(datose)
## # A tibble: 6 × 10
##      id  edad sexo    imc presion_sistolica glucosa_ayuno actividad_fisica
##   <dbl> <dbl> <chr> <dbl>             <dbl>         <dbl> <chr>           
## 1     1    19 F      23.3               102            73 activo          
## 2     2    24 M      17.8               138            98 moderado        
## 3     3    44 F      25                 116            84 activo          
## 4     4    74 F      33.1               149            91 moderado        
## 5     5    31 F      28.9                80            84 sedentario      
## 6     6    31 F      33.1               118            97 moderado        
## # ℹ 3 more variables: fumador <chr>, riesgo_cardiovascular <dbl>,
## #   riesgo_alto <chr>

modelacion de los datos

EDAD

## 
## Call:
## glm(formula = riesgo_alto ~ poly(edad, 2), family = binomial, 
##     data = datose)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.2681     0.5175  -4.383 1.17e-05 ***
## poly(edad, 2)1  75.3423    12.2910   6.130 8.80e-10 ***
## poly(edad, 2)2   4.5504     9.1195   0.499    0.618    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 243.96  on 497  degrees of freedom
## AIC: 249.96
## 
## Number of Fisher Scoring iterations: 8
## 
## Call:
## glm(formula = riesgo_alto ~ poly(edad, 3), family = binomial, 
##     data = datose)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      -5.112      2.221  -2.302   0.0213 *
## poly(edad, 3)1  173.045     70.010   2.472   0.0134 *
## poly(edad, 3)2  -62.092     46.816  -1.326   0.1847  
## poly(edad, 3)3   39.533     23.848   1.658   0.0974 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 239.37  on 496  degrees of freedom
## AIC: 247.37
## 
## Number of Fisher Scoring iterations: 10
## 
## Call:
## glm(formula = riesgo_alto ~ poly(edad, 4), family = binomial, 
##     data = datose)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -5.110971   5.861117  -0.872    0.383
## poly(edad, 4)1 173.001050 178.859243   0.967    0.333
## poly(edad, 4)2 -62.052118 156.844964  -0.396    0.692
## poly(edad, 4)3  39.512358  81.284008   0.486    0.627
## poly(edad, 4)4   0.009228  34.397474   0.000    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 239.37  on 495  degrees of freedom
## AIC: 249.37
## 
## Number of Fisher Scoring iterations: 12
## Analysis of Deviance Table
## 
## Model 1: riesgo_alto ~ edad
## Model 2: riesgo_alto ~ poly(edad, 2)
## Model 3: riesgo_alto ~ poly(edad, 3)
## Model 4: riesgo_alto ~ poly(edad, 4)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       498     244.19                       
## 2       497     243.96  1   0.2341  0.62854  
## 3       496     239.37  1   4.5852  0.03225 *
## 4       495     239.37  1   0.0000  0.99979  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SEXO

#basandose en el hecho de que la variable es categorica solo se pueden usar modelos parametricos en este caso se usa el modelo glm ya que la variable es categorica 

datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$sexo <- as.factor(datose$sexo)
#modelo
modelosexo <- glm(riesgo_alto ~ sexo,
                  data = datose,
                  family = binomial)

summary(modelosexo)
## 
## Call:
## glm(formula = riesgo_alto ~ sexo, family = binomial, data = datose)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0090     0.1396  -7.227 4.95e-13 ***
## sexoM         0.3095     0.1961   1.579    0.114    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 606.66  on 498  degrees of freedom
## AIC: 610.66
## 
## Number of Fisher Scoring iterations: 4
#dado a que el modelo mustra que el cambio de sexo no es estaditicamente significativo para explicar el riesgo alto se utiliza entonces un chisq test para observar si realmente no hay relacion estadistica 
tabla_sexo <- table(datose$sexo, datose$riesgo_alto)

print(tabla_sexo)
##    
##     bajo alto
##   F  192   70
##   M  159   79
prop.table(tabla_sexo, margin = 1)
##    
##          bajo      alto
##   F 0.7328244 0.2671756
##   M 0.6680672 0.3319328
chisq.test(tabla_sexo)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabla_sexo
## X-squared = 2.2, df = 1, p-value = 0.138
#el chisq test muestra que sexo relamente no es significativa al mostrar un p.value mayor a 0.05

imc

#factor
datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$riesgo_num <- ifelse(datose$riesgo_alto == "alto", 1, 0)
#ver si hay relacion evidente
plot(datose$imc, datose$riesgo_num, 
     main="ims vs riesgo ",
     pch=19, col="darkblue")

# como no hay se utilizara el mejor modelo 
modeloimc <- glm(riesgo_alto ~ poly(imc, 3), data = datose, family = binomial)

#ver cual modelo polinomial es mejor 
m1 <- glm(riesgo_alto ~ poly(imc, 1), data = datose, family = binomial)
m2 <- glm(riesgo_alto ~ poly(imc, 2), data = datose, family = binomial)
m3 <- glm(riesgo_alto ~ poly(imc, 3), data = datose,family = binomial)
m4 <- glm(riesgo_alto ~ poly(imc, 4), data = datose,family = binomial)
tabla_anova_imc <- anova(m1, m2, m3, m4)
print(tabla_anova_imc)
## Analysis of Deviance Table
## 
## Model 1: riesgo_alto ~ poly(imc, 1)
## Model 2: riesgo_alto ~ poly(imc, 2)
## Model 3: riesgo_alto ~ poly(imc, 3)
## Model 4: riesgo_alto ~ poly(imc, 4)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       498     590.83                       
## 2       497     587.13  1   3.6961  0.05454 .
## 3       496     584.96  1   2.1709  0.14064  
## 4       495     584.19  1   0.7750  0.37868  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#loess
span_values <- seq(0.2, 0.9, by = 0.05)

errores <- sapply(span_values, function(s) {modelolo <- loess(riesgo_num ~ imc, data =datose, span = s)
sum(modelolo$residuals^2)})

optimo_span <- span_values[which.min(errores)]

print(paste("El span que minimiza el error es:", optimo_span))
## [1] "El span que minimiza el error es: 0.2"
modeloimcl <- loess(riesgo_num ~ imc, data = datose, span = optimo_span)

#spline 
modelospimc <- smooth.spline(x = datose$imc, y = datose$riesgo_num,cv = TRUE)
## Warning in smooth.spline(x = datose$imc, y = datose$riesgo_num, cv = TRUE):
## cross-validation with non-unique 'x' values seems doubtful
#predicciones
x_range <- seq(min(datose$imc, na.rm = TRUE), max(datose$imc, na.rm = TRUE),length.out = 100)

prediccionesl <- predict(modeloimcl, newdata = data.frame(imc= x_range))

predlineal <- predict(modeloimc, newdata = data.frame(imc = x_range), type = "response")

predspline <- predict(modelospimc,  x = x_range)

#comparativa 

plot(datose$imc, datose$riesgo_num,
     pch = 16,
     col = "pink",
     main = "IMC vs riesgo cardiovascular")

lines(x_range, predlineal,col = "red", lwd = 3)
lines(x_range, prediccionesl, col = "steelblue", lwd = 3)
lines(predspline$x, predspline$y,col = "darkgreen",lwd = 3)

#comparacion estadistica 
data.frame(Modelo = c("Polinomial3", "LOESS", "Spline"),RMSE = c(sqrt(mean((datose$riesgo_num - predict(modeloimc))^2, na.rm = TRUE)),sqrt(mean((datose$riesgo_num - predict(modeloimcl))^2, na.rm = TRUE)),sqrt(mean((datose$riesgo_num- predict(modelospimc, x = datose$imc)$y)^2, na.rm = TRUE))))
##        Modelo      RMSE
## 1 Polinomial3 1.3543741
## 2       LOESS 0.4409212
## 3      Spline 0.4475874
#Aqui el modelo que reduce el error y por ende modela mejor la relacion es el modelo LOESS, sin embargo el spline se encuentra cerca en custion se reducir el error .Es entonces que se prefiere el spline sya que este proporciona una carva mas estable que no es tan sensible como el modelo LOESS

presion_sistolica

#factor
datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$riesgo_num <- ifelse(datose$riesgo_alto == "alto", 1, 0)
#observar la relacion
plot(datose$presion_sistolica, datose$riesgo_num, 
     main="presion vs riesgo ",
     pch=19, col="darkblue")

# como no hay se utilizara el mejor modelo 
modelopr <- glm(riesgo_alto ~ poly(presion_sistolica, 3), data = datose, family = binomial)

#ver cual modelo polinomial es mejor 
m1p <- glm(riesgo_alto ~ poly(presion_sistolica, 1), data = datose, family = binomial)
m2p <- glm(riesgo_alto ~ poly(presion_sistolica, 2), data = datose, family = binomial)
m3p <- glm(riesgo_alto ~ poly(presion_sistolica, 3), data = datose, family = binomial)
m4p <- glm(riesgo_alto ~ poly(presion_sistolica, 4), data = datose, family = binomial)
tablapr <- anova(m1p, m2p, m3p, m4p)
print(tablapr)
## Analysis of Deviance Table
## 
## Model 1: riesgo_alto ~ poly(presion_sistolica, 1)
## Model 2: riesgo_alto ~ poly(presion_sistolica, 2)
## Model 3: riesgo_alto ~ poly(presion_sistolica, 3)
## Model 4: riesgo_alto ~ poly(presion_sistolica, 4)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       498     595.69                       
## 2       497     593.73  1   1.9622  0.16127  
## 3       496     590.39  1   3.3360  0.06778 .
## 4       495     590.21  1   0.1847  0.66737  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#LOES
spanvaluespr <- seq(0.2, 0.9, by = 0.05)

errorespr <- sapply(spanvaluespr, function(s) {modelolopr <- loess(riesgo_num ~ presion_sistolica, data = datose, span = s)
sum(modelolopr$residuals^2)})

optimospanpr <- spanvaluespr[which.min(errorespr)]

print(paste("El span que minimiza el error es:", optimospanpr))
## [1] "El span que minimiza el error es: 0.2"
modeloprlo <- loess(riesgo_num ~ presion_sistolica, data = datose, span = optimospanpr)

#Spline
modelosppr <- smooth.spline( x = datose$presion_sistolica, y = datose$riesgo_num, cv = TRUE)
## Warning in smooth.spline(x = datose$presion_sistolica, y = datose$riesgo_num, :
## cross-validation with non-unique 'x' values seems doubtful
#predicciones
xrangepr <- seq(min(datose$presion_sistolica, na.rm = TRUE),max(datose$presion_sistolica, na.rm = TRUE),length.out = 100)

prediccionesprlo <- predict(modeloprlo,newdata = data.frame(presion_sistolica = xrangepr))

predpr <- predict(modelopr,newdata = data.frame(presion_sistolica = xrangepr))

predsplinepr <- predict(modelosppr,x = xrangepr)

#comparativa 
plot(datose$presion_sistolica, datose$riesgo_num,
     pch = 16,
     col = "purple",
     main = "Presión sistólica vs riesgo cardiovascular")

lines(xrangepr, predpr, col = "red", lwd = 3)
lines(xrangepr, prediccionesprlo, col = "steelblue", lwd = 3)
lines(predsplinepr$x, predsplinepr$y, col = "darkgreen", lwd = 3)

#comparativa estadistica 
data.frame(Modelo = c("Polinomial3", "LOESS", "Spline"),RMSE = c(sqrt(mean((datose$riesgo_num - predict(modelopr))^2, na.rm = TRUE)), sqrt(mean((datose$riesgo_num - predict(modeloprlo))^2, na.rm = TRUE)), sqrt(mean((datose$riesgo_num - predict(modelosppr, x = datose$presion_sistolica)$y)^2, na.rm = TRUE))))
##        Modelo      RMSE
## 1 Polinomial3 1.3050591
## 2       LOESS 0.4440639
## 3      Spline 0.4484260
#Aqui el modelo que reduce el error y por ende modela mejor la relacion es el modelo LOESS, sin embargo el spline se encuentra cerca en custion se reducir el error .Es entonces que se prefiere el spline sya que este proporciona una carva mas estable que no es tan sensible como el modelo LOESS

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

glucosa_ayuno

#factor
datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$riesgo_num <- ifelse(datose$riesgo_alto == "alto", 1, 0)
#ver la relacion
plot(datose$glucosa_ayuno, datose$riesgo_num,main = "glucosa en ayuno vs riesgo",pch = 19, col = "darkblue")

#como no se observa una ralacion simple se utilizara otro modelo 
modeloga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 4), data = datose, family = binomial)

m1ga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 1), data = datose, family = binomial)
m2ga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 2), data = datose, family = binomial)
m3ga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 3), data = datose, family = binomial)
m4ga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 4), data = datose, family = binomial)
m5ga <- glm(riesgo_alto ~ poly(glucosa_ayuno, 5), data = datose, family = binomial)

tablaga <- anova(m1ga, m2ga, m3ga, m4ga, m5ga)
print(tablaga)
## Analysis of Deviance Table
## 
## Model 1: riesgo_alto ~ poly(glucosa_ayuno, 1)
## Model 2: riesgo_alto ~ poly(glucosa_ayuno, 2)
## Model 3: riesgo_alto ~ poly(glucosa_ayuno, 3)
## Model 4: riesgo_alto ~ poly(glucosa_ayuno, 4)
## Model 5: riesgo_alto ~ poly(glucosa_ayuno, 5)
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       498     609.16                     
## 2       497     608.80  1  0.35432   0.5517
## 3       496     608.67  1  0.13232   0.7160
## 4       495     607.12  1  1.55088   0.2130
## 5       494     606.64  1  0.47867   0.4890
#LOESS
spanvaluesga <- seq(0.2, 0.9, by = 0.05)

erroresga <- sapply(spanvaluesga, function(s) {modelologa <- loess(riesgo_num ~ glucosa_ayuno, data = datose, span = s)
sum(modelologa$residuals^2)})


optimospanga <- spanvaluesga[which.min(erroresga)]

print(paste("El span que minimiza el error es:", optimospanga))
## [1] "El span que minimiza el error es: 0.2"
modelogalo <- loess(riesgo_num ~ glucosa_ayuno, data = datose, span = optimospanga)

#Spline
modelospga <- smooth.spline(x = datose$glucosa_ayuno,y = datose$riesgo_num,cv = TRUE)
## Warning in smooth.spline(x = datose$glucosa_ayuno, y = datose$riesgo_num, :
## cross-validation with non-unique 'x' values seems doubtful
#predicciones
xrangega <- seq(min(datose$glucosa_ayuno, na.rm = TRUE),max(datose$glucosa_ayuno, na.rm = TRUE),length.out = 100)

prediccionesgalo <- predict(modelogalo,newdata = data.frame(glucosa_ayuno = xrangega))

predga <- predict(modeloga,newdata = data.frame(glucosa_ayuno = xrangega))

predsplinega <- predict(modelospga,x = xrangega)

#comparativa
plot(datose$glucosa_ayuno, datose$riesgo_num,
     pch = 16,
     col = "lightpink",
     main = "Glucosa en ayuno vs riesgo cardiovascular")

lines(xrangega, predga, col = "red", lwd = 3)
lines(xrangega, prediccionesgalo, col = "steelblue", lwd = 3)
lines(predsplinega$x, predsplinega$y, col = "darkgreen", lwd = 3)

#comparativa estadistica
data.frame(Modelo = c("Polinomial4", "LOESS", "Spline"),RMSE = 
c(sqrt(mean((datose$riesgo_num - predict(modeloga))^2, na.rm = TRUE)),
sqrt(mean((datose$riesgo_num - predict(modelogalo))^2, na.rm = TRUE)),
sqrt(mean((datose$riesgo_num - predict(modelospga, x = datose$glucosa_ayuno)$y)^2, na.rm = TRUE))))
##        Modelo      RMSE
## 1 Polinomial4 1.2504715
## 2       LOESS 0.4498901
## 3      Spline 0.4573785
#Aqui el modelo que reduce el error y por ende modela mejor la relacion es el modelo LOESS, sin embargo el spline se encuentra cerca en custion se reducir el error .Es entonces que se prefiere el spline sya que este proporciona una carva mas estable que no es tan sensible como el modelo LOESS

Fumador

#basandose en el hecho de que la variable es categorica solo se pueden usar modelos parametricos en este caso se usa el modelo glm ya que la variable es categorica 

datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$fumador <- as.factor(datose$fumador)
#modelo
modelosfu <- glm(riesgo_alto ~ fumador,
                  data = datose,
                  family = binomial)

summary(modelosfu)
## 
## Call:
## glm(formula = riesgo_alto ~ fumador, family = binomial, data = datose)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9144     0.1155  -7.917 2.44e-15 ***
## fumadorsi     0.2099     0.2175   0.965    0.334    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 608.24  on 498  degrees of freedom
## AIC: 612.24
## 
## Number of Fisher Scoring iterations: 4
#dado a que el modelo mustra que el cambio de si es fumador o no no es estaditicamente significativo para explicar el riesgo alto se utiliza un chisq test para observar si realmente no hay relacion estadistica 
tablafu <- table(datose$fumador, datose$riesgo_alto)

print(tablafu)
##     
##      bajo alto
##   no  262  105
##   si   89   44
prop.table(tablafu, margin = 1)
##     
##           bajo      alto
##   no 0.7138965 0.2861035
##   si 0.6691729 0.3308271
chisq.test(tablafu)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tablafu
## X-squared = 0.73185, df = 1, p-value = 0.3923
#el chisq test muestra que ser fumador no es realmente significativo al mostrar un p.value mayor a 0.05

Actividad fisica

#basandose en el hecho de que la variable es categorica solo se pueden usar modelos parametricos en este caso se usa el modelo glm ya que la variable es categorica 

datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$actividad_fisica <- as.factor(datose$actividad_fisica)
#modelo
modelosac <- glm(riesgo_alto ~ actividad_fisica,
                  data = datose,
                  family = binomial)

summary(modelosac)
## 
## Call:
## glm(formula = riesgo_alto ~ actividad_fisica, family = binomial, 
##     data = datose)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -1.1550     0.2391  -4.830 1.36e-06 ***
## actividad_fisicamoderado     0.3460     0.2873   1.204    0.228    
## actividad_fisicasedentario   0.3784     0.2799   1.352    0.176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 609.16  on 499  degrees of freedom
## Residual deviance: 607.14  on 497  degrees of freedom
## AIC: 613.14
## 
## Number of Fisher Scoring iterations: 4
#dado a que el modelo mustra que el cambio de si hace o no actividad fisica noes estaditicamente significativo para explicar el riesgo alto se utiliza un chisq test para observar si realmente no hay relacion estadistica 
tablaac <- table(datose$actividad_fisica, datose$riesgo_alto)

print(tablaac)
##             
##              bajo alto
##   activo       73   23
##   moderado    128   57
##   sedentario  150   69
prop.table(tablaac, margin = 1)
##             
##                   bajo      alto
##   activo     0.7604167 0.2395833
##   moderado   0.6918919 0.3081081
##   sedentario 0.6849315 0.3150685
chisq.test(tablaac)
## 
##  Pearson's Chi-squared test
## 
## data:  tablaac
## X-squared = 1.9613, df = 2, p-value = 0.3751
#el chisq test muestra que hacer actividad fisica no es relamente significativo al mostrar un p.value mayor a 0.05

#GAM

#convertir a factores 
datose$riesgo_alto <- factor(datose$riesgo_alto, levels = c("bajo", "alto"))
datose$sexo <- as.factor(datose$sexo)
datose$fumador <- as.factor(datose$fumador)
datose$actividad_fisica <- as.factor(datose$actividad_fisica)

datose$riesgo_num <- ifelse(datose$riesgo_alto == "alto", 1, 0)
#exploracion visualde datos
par(mfrow = c(2, 4))

plot(datose$edad, datose$riesgo_num, pch = 16, col = "blue",
     main = "Riesgo alt0/Edad")
lines(lowess(datose$edad, datose$riesgo_num), col = "darkblue", lwd = 2)

plot(datose$imc, datose$riesgo_num, pch = 16, col = "pink",
     main = "Riesgo alto/IM")
lines(lowess(datose$imc, datose$riesgo_num), col = "lightblue", lwd = 2)

plot(datose$presion_sistolica, datose$riesgo_num, pch = 16, col = "aquamarine",
     main = "Riesgo alto/Presión")
lines(lowess(datose$presion_sistolica, datose$riesgo_num), col = "purple", lwd = 2)

plot(datose$glucosa_ayuno, datose$riesgo_num, pch = 16, col = "red",
     main = "Riesgo alto/Glucosa")
lines(lowess(datose$glucosa_ayuno, datose$riesgo_num), col = "black", lwd = 2)


barplot(prop.table(table(datose$sexo, datose$riesgo_alto), margin = 1),
        beside = TRUE, col = c("gold", "brown"),
        main = "Riesgo alto/Sexo")

barplot(prop.table(table(datose$fumador, datose$riesgo_alto), margin = 1),
        beside = TRUE, col = c("yellow", "green"),
        main = "Riesgo alto/Fumador")

barplot(prop.table(table(datose$actividad_fisica, datose$riesgo_alto), margin = 1),
        beside = TRUE, col = c("darkgreen", "darkblue"),
        main = "Riesgo alto/Actividad física")

par(mfrow = c(1, 1))

#gam
modelogam <- gam(
  riesgo_alto ~ 
    s(edad, k = 10) +
    s(imc, k = 10) +
    s(presion_sistolica, k = 10) +
    s(glucosa_ayuno, k = 10) +
    sexo +
    fumador +
    actividad_fisica,
  data = datose,
  family = binomial,
  method = "REML"
)
summary(modelogam)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## riesgo_alto ~ s(edad, k = 10) + s(imc, k = 10) + s(presion_sistolica, 
##     k = 10) + s(glucosa_ayuno, k = 10) + sexo + fumador + actividad_fisica
## 
## Parametric coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -8.1698     1.2247  -6.671 2.55e-11 ***
## sexoM                        1.0234     0.4998   2.048 0.040593 *  
## fumadorsi                    3.2808     0.7169   4.576 4.73e-06 ***
## actividad_fisicamoderado     1.5616     0.7021   2.224 0.026135 *  
## actividad_fisicasedentario   2.6142     0.7609   3.436 0.000591 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df Chi.sq  p-value    
## s(edad)              1.000  1.000 54.614  < 2e-16 ***
## s(imc)               3.231  4.042 36.038 4.37e-07 ***
## s(presion_sistolica) 3.682  4.592 15.658  0.00573 ** 
## s(glucosa_ayuno)     1.000  1.000  2.238  0.13472    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.832   Deviance explained = 81.8%
## -REML = 62.379  Scale est. = 1         n = 500
#evaluacion del modelo conjunto
par(mfrow = c(2, 2))
gam.check(modelogam)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-1.097051e-05,1.838552e-06]
## (score 62.37901 & scale 1).
## Hessian positive definite, eigenvalue range [5.504187e-06,0.6052596].
## Model rank =  41 / 41 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                        k'  edf k-index p-value  
## s(edad)              9.00 1.00    1.00    0.42  
## s(imc)               9.00 3.23    1.01    0.62  
## s(presion_sistolica) 9.00 3.68    0.95    0.09 .
## s(glucosa_ayuno)     9.00 1.00    0.97    0.22  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(1, 1))
#evaluacion de valores parciales 

par(mfrow = c(2, 2))

plot(modelogam, shade = TRUE, shade.col = "blue",
     seWithMean = TRUE, scale = 0, residuals = TRUE,
     pch = 16, cex = 0.4, col = "darkblue",
     select = 1, main = "s(edad) — efecto parcial")
abline(h = 0, lty = 2, col = "gray")

plot(modelogam, shade = TRUE, shade.col = "pink",
     seWithMean = TRUE, scale = 0, residuals = TRUE,
     pch = 16, cex = 0.4, col = "lightblue",
     select = 2, main = "s(imc) — efecto parcial")
abline(h = 0, lty = 2, col = "gray")

plot(modelogam, shade = TRUE, shade.col = "aquamarine",
     seWithMean = TRUE, scale = 0, residuals = TRUE,
     pch = 16, cex = 0.4, col = "purple",
     select = 3, main = "s(presion_sistolica) — efecto parcial")
abline(h = 0, lty = 2, col = "gray")

plot(modelogam, shade = TRUE, shade.col = "red",
     seWithMean = TRUE, scale = 0, residuals = TRUE,
     pch = 16, cex = 0.4, col = "black",
     select = 4, main = "s(glucosa_ayuno) — efecto parcial")
abline(h = 0, lty = 2, col = "gray")

par(mfrow = c(1, 1))
#En base a los reyltados se pueden notar distintas cosas en el modelo gam en primera el modelo explica el 81.8% de la devianza lo cual significa que reduce bastante el error lo que implica que las variables si logran explicar la varible riesgo alto y por otro lado tiene un buen ajuste de r cuadrada de 83.2%. Por otro lado las variables parametricas como fumador gracias al coeficiente se puede concluir que fumar si tiene un efecto significativo en el riesgo alto ademas de que es la varible con menor p-value indicando significancia, por otro lado se mustra como entre menor actividad existe mayor riesgo ya que tanto sendentario como moderadoson significativas y tienen coeficientes positivos, sin embargo el sendentario tiene un coeficiente mayor.Esto contrasta con el análisis individual, donde las variables categóricas no salieron significativas por sí solas. Sin embargo, esto puede pasar porque el GAM no analiza cada variable de forma aislada, sino dentro del conjunto completo de variables.

#En los terminos que son numericos, existe significancia en todos menos en el de glucosa, ademastodas menos glucosa muestran grados de libertan que les dan una flexibilidad a las relaciones que se parece a las relaciones establecidas anteriormente donde edad mustra la relacion mas simple en contraste con imc y presion

#En general en base al modelo se puede decir que las variables significativas en este caso son edad, IMC, presión sistólica,sexo, fumar y actividad física

#Arbol de decisiones

datose$riesgo_alto <- as.factor(datose$riesgo_alto)
datose$sexo <- as.factor(datose$sexo)
datose$fumador <- as.factor(datose$fumador)
datose$actividad_fisica <- as.factor(datose$actividad_fisica)
#entrenamiento 
set.seed(123)

train <- sample(1:nrow(datose), nrow(datose) / 2)
datose.test <- datose[-train, ]
tree.riesgo <- tree(riesgo_alto ~ edad + sexo + imc + presion_sistolica +glucosa_ayuno + actividad_fisica + fumador,data = datose,subset = train)

summary(tree.riesgo)
## 
## Classification tree:
## tree(formula = riesgo_alto ~ edad + sexo + imc + presion_sistolica + 
##     glucosa_ayuno + actividad_fisica + fumador, data = datose, 
##     subset = train)
## Variables actually used in tree construction:
## [1] "edad"              "imc"               "presion_sistolica"
## [4] "glucosa_ayuno"    
## Number of terminal nodes:  12 
## Residual mean deviance:  0.2553 = 60.77 / 238 
## Misclassification error rate: 0.06 = 15 / 250
#visualizacion
plot(tree.riesgo)
text(tree.riesgo, pretty = 0)

#proceso de podado
cv.riesgo <- cv.tree(tree.riesgo, FUN = prune.misclass)

plot(cv.riesgo$size, cv.riesgo$dev, type = "b",
     main = "Validación cruzada del árbol")

 #aqui el que tiene el nivel mas bajo de error es el nivel 7 ya que aqui se puede observar el "codo"

prune.riesgo <- prune.misclass(tree.riesgo, best = 7)

plot(prune.riesgo)
text(prune.riesgo, pretty = 0)

#comparacion con matrices
pred.riesgo <- predict(tree.riesgo,newdata = datose.test,type = "class")

pred.prune <- predict(prune.riesgo,newdata = datose.test,type = "class")

real.riesgo <- datose.test$riesgo_alto

tabla.riesgo <- table(Real = real.riesgo,Predicho = pred.riesgo)
tabla.prune <- table(Real = real.riesgo,Predicho = pred.prune)

print(tabla.prune)
##       Predicho
## Real   bajo alto
##   bajo  161   16
##   alto   11   62
print(tabla.riesgo)
##       Predicho
## Real   bajo alto
##   bajo  160   17
##   alto    8   65
accuracy.riesgo <- mean(pred.riesgo == real.riesgo)
accuracy.prune <- mean(pred.prune == real.riesgo)

error.riesgo <- 1 - accuracy.riesgo
error.prune <- 1 - accuracy.prune

comparacion.arboles <- data.frame(
  Modelo = c("Árbol sin podar", "Árbol podado"),
  Accuracy = c(accuracy.riesgo, accuracy.prune),
  Error = c(error.riesgo, error.prune))

print(comparacion.arboles)
##            Modelo Accuracy Error
## 1 Árbol sin podar    0.900 0.100
## 2    Árbol podado    0.892 0.108
#En base a este arbol de decisiones las variables que explican mejor o afectan mas el riesgo alto son la edad,la IMC,la presion sistotelica y glucosa.Por otro lado al evaluar el error y accuracy de ambos arboles se demostro como el arbol sin podar tiene un desepeño solo un poco mayor al arbol podado ademas de que el arbol podado reduce los falsos positivos o mejor dicho los bajos altos. Es entonces que el arbol podado es preferible ya que no se ajusta y de similar manera no deja que personas con resgo alto dean desapercividas