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
## Single term deletions
## 
## Model:
## riesgo_alto ~ poly(edad, 3)
##               Df Deviance    AIC    LRT  Pr(>Chi)    
## <none>             239.37 247.37                     
## poly(edad, 3)  3   609.16 611.16 369.79 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##   Modelo      RMSE
## 1 Lineal 0.2775445
## 2  Poli2 0.2773615
## 3  Poli3 0.2754523
## 4  Poli4 0.2754523

SEXO

#Basándose en el hecho de que la variable es categórica, solo se pueden usar modelos paramétricos. En este caso, se usa el modelo glm ya que la variable es categórica.
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 muestra que el cambio de sexo no es estadísticamente significativo para explicar el riesgo alto se utiliza entonces un chisq test para observar si realmente no hay relación estadística
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 realmente no es significativa al mostrar un p.value mayor a 0.05

#Por tanto se utilizara un modelo con controles 
modelosexocontrol <- glm(
  riesgo_alto ~ sexo + edad + imc + presion_sistolica + glucosa_ayuno,
  data = datose,
  family = binomial
)
summary(modelosexocontrol)
## 
## Call:
## glm(formula = riesgo_alto ~ sexo + edad + imc + presion_sistolica + 
##     glucosa_ayuno, family = binomial, data = datose)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -34.649822   4.225409  -8.200 2.40e-16 ***
## sexoM               0.782722   0.408770   1.915   0.0555 .  
## edad                0.285364   0.030969   9.215  < 2e-16 ***
## imc                 0.319718   0.054676   5.847 4.99e-09 ***
## presion_sistolica   0.060106   0.014359   4.186 2.84e-05 ***
## glucosa_ayuno       0.010369   0.009391   1.104   0.2695    
## ---
## 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: 164.52  on 494  degrees of freedom
## AIC: 176.52
## 
## Number of Fisher Scoring iterations: 8
#Aun con el modelo de controles sexo demuestra no ser estadisticamente significativa 

#representación visual 
boxplot(riesgo_cardiovascular ~ sexo,
        data = datose,
        col = "aquamarine",
        main = "Riesgo cardiovascular ~ Sexo")

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
#significancia 
 
drop1(modeloimc, test = "Chisq")
## Single term deletions
## 
## Model:
## riesgo_alto ~ poly(imc, 3)
##              Df Deviance    AIC  LRT  Pr(>Chi)    
## <none>            584.96 592.96                   
## poly(imc, 3)  3   609.16 611.16 24.2 2.269e-05 ***
## ---
## 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, type = "response"))^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 0.4467027
## 2       LOESS 0.4409212
## 3      Spline 0.4475874
#Aquí el modelo que reduce el error y por ende modela mejor la relación es el modelo LOESS, sin embargo es de relevancia mencionar el span chico que se está utilizando por lo cual el modelo es más flexible y sensible a variaciones locales. Por ello, si el objetivo fuera construir un modelo más predictivo y menos dependiente de la forma específica de la muestra el polinomial será el mejor modelo a utilizar.

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, test = "Chisq" )
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
#significancia 
 
drop1(modelopr, test = "Chisq")
## Single term deletions
## 
## Model:
## riesgo_alto ~ poly(presion_sistolica, 3)
##                            Df Deviance    AIC    LRT  Pr(>Chi)    
## <none>                          590.39 598.39                     
## poly(presion_sistolica, 3)  3   609.16 611.16 18.766 0.0003055 ***
## ---
## 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), type = "response")

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, type = "response"))^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 0.4481671
## 2       LOESS 0.4440639
## 3      Spline 0.4484260
#Aquí el modelo que reduce el error y por ende modela mejor la relación es el modelo LOESS, sin embargo es de relevancia mencionar el span chico que se está utilizando por lo cual el modelo es más flexible y sensible a variaciones locales. Por ello, si el objetivo fuera construir un modelo más predictivo y menos dependiente de la forma específica de la muestra el polinomial será el mejor modelo a utilizar

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
#significancia 
 
drop1(modeloga, test = "Chisq")
## Single term deletions
## 
## Model:
## riesgo_alto ~ poly(glucosa_ayuno, 4)
##                        Df Deviance    AIC    LRT Pr(>Chi)
## <none>                      607.12 617.12                
## poly(glucosa_ayuno, 4)  4   609.16 611.16 2.0392   0.7285
#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), type = "response")

predsplinega <- predict(modelospga,x = xrangega)

#comparativa
plot(datose$glucosa_ayuno, datose$riesgo_num,
     pch = 16,
     col = "lightblue",
     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, type = "response"))^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 0.4564547
## 2       LOESS 0.4498901
## 3      Spline 0.4573785
#Aquí el modelo que reduce el error y por ende modela mejor la relación es el modelo LOESS, sin embargo es de relevancia mencionar el span chico que se está utilizando por lo cual el modelo es más flexible y sensible a variaciones locales. Por ello, si el objetivo fuera construir un modelo más predictivo y menos dependiente de la forma específica de la muestra el polinomial será el mejor modelo a utilizar.

Fumador

#Basándose en el hecho de que la variable es categórica solo se pueden usar modelos paramétricos en este caso se usa el modelo glm ya que la variable es categórica

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 muestra que el cambio de si es fumador o no no es estadísticamente significativo para explicar el riesgo alto se utiliza un chisq test para observar si realmente no hay relación estadística
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
#Se usa un modelo con controles 
modelofumcontrol <- glm(
  riesgo_alto ~ fumador + edad + imc + presion_sistolica + glucosa_ayuno,
  data = datose,
  family = binomial
)
summary(modelofumcontrol)
## 
## Call:
## glm(formula = riesgo_alto ~ fumador + edad + imc + presion_sistolica + 
##     glucosa_ayuno, family = binomial, data = datose)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -44.39908    6.08703  -7.294 3.01e-13 ***
## fumadorsi           2.67219    0.61302   4.359 1.31e-05 ***
## edad                0.36281    0.04540   7.992 1.33e-15 ***
## imc                 0.44432    0.07414   5.993 2.06e-09 ***
## presion_sistolica   0.06666    0.01598   4.172 3.02e-05 ***
## glucosa_ayuno       0.01613    0.01040   1.552    0.121    
## ---
## 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: 143.24  on 494  degrees of freedom
## AIC: 155.24
## 
## Number of Fisher Scoring iterations: 8
#bajo este modelo la variable fumador sí es estadísticamente significativa al tener un p.value menor a 0.05, lo que indica que ser fumador sí afecta el riesgo cardiovascular manteniendo las demás variables constantes

#representacion visual 
boxplot(riesgo_cardiovascular ~ fumador,
        data = datose,
        col = "lightblue",
        main = "Riesgo cardiovascular ~ Fumador")

Actividad fisica

#Basándose en el hecho de que la variable es categórica solo se pueden usar modelos paramétricos en este caso se usa el modelo glm ya que la variable es categórica 

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 muestra que el cambio de si hace o no actividad física no es estadísticamente significativo para explicar el riesgo alto se utiliza un chisq test para observar si realmente no hay relación estadística
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 realmente significativo al mostrar un p.value mayor a 0.05

#Se usa un modelo con controles 
modeloactcontrol <- glm(
  riesgo_alto ~ actividad_fisica + edad + imc + presion_sistolica + glucosa_ayuno,
  data = datose,
  family = binomial
)
summary(modeloactcontrol)
## 
## Call:
## glm(formula = riesgo_alto ~ actividad_fisica + edad + imc + presion_sistolica + 
##     glucosa_ayuno, family = binomial, data = datose)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -37.414766   4.657622  -8.033 9.51e-16 ***
## actividad_fisicamoderado     1.136279   0.567196   2.003 0.045142 *  
## actividad_fisicasedentario   2.025430   0.608720   3.327 0.000877 ***
## edad                         0.317457   0.036116   8.790  < 2e-16 ***
## imc                          0.338119   0.056209   6.015 1.80e-09 ***
## presion_sistolica            0.055662   0.014599   3.813 0.000137 ***
## glucosa_ayuno                0.010658   0.009961   1.070 0.284622    
## ---
## 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: 155.73  on 493  degrees of freedom
## AIC: 169.73
## 
## Number of Fisher Scoring iterations: 8
#bajo este modelo la variable de actividad física sí es estadísticamente significativa y muestra cómo ser sedentario tiene un mayor efecto en tener riesgo cardiovascular, al esta variable tener el mayor coeficiente de todos los estados de actividad. En general este modelo muestra cómo tener actividad física reduce el riesgo cardiovascular en tanto todas las demás variables sean constantes. 

#representacion visual 
boxplot(riesgo_cardiovascular ~ actividad_fisica,
        data = datose,
        col = "blue",
        main = "Riesgo cardiovascular ~ Actividad física")

#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 visual de 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.52
## s(imc)               9.00 3.23    1.01    0.61
## s(presion_sistolica) 9.00 3.68    0.95    0.13
## s(glucosa_ayuno)     9.00 1.00    0.97    0.22
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 resultados 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 sí logran explicar la variable riesgo alto y por otro lado tiene un buen ajuste de r cuadrada de 83.2%. Por otro lado las variables paramétricas como fumador gracias al coeficiente se puede concluir que fumar sí tiene un efecto significativo en el riesgo alto además de que es la variable con menor p-value indicando significancia, por otro lado se muestra como entre menor actividad existe mayor riesgo ya que tanto sedentario como moderado son significativas y tienen coeficientes positivos, sin embargo el sedentario 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 términos que son numéricos, existe significancia en todos menos en el de glucosa, además todas menos glucosa 

#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 árbol de decisiones las variables que explican mejor o afectan más el riesgo alto son la edad, el imc, la presión sistólica y glucosa. Por otro lado al evaluar el error y accuracy de ambos árboles se demostró como el árbol sin podar tiene un desempeño solo un poco mayor al árbol podado además de que el árbol podado reduce los falsos positivos o mejor dicho los bajos altos. Es entonces que el árbol podado es preferible ya que no se ajusta y de similar manera no deja que personas con riesgo alto sean desapercibidas.