Acerca de la autora: Soy estudiante ecuatorina de la Carrera de Ingeniería en Botecnolgía de la Escuala Politécnia del Ejército. En esta página quiero mostrar los conocimientos adquiridos en el Curso de R, dictado por Karen Calva de AsoiMat, EPN . Y motivar a desarrollar habilidades en el tratamiento de datos con herramientas que actualmente están al alcance de muchos, como R:)

A. Base de datos

Figura 1.Factores de riesgo cardiovasular

Figura 1.Factores de riesgo cardiovasular

a1. Describir la base de datos y las varibles.

La base de datos “pacientes” consiste en un conjunto de 70 observaciones referentes a los factores de riesgo cardiovascular. Donde se toman en cuenta las siguientes variables:

a2. Leer la base de datos, nombres de columnas adecuados, definir tipo de varibles y etiquetas a las tipo factor.

colnames(datos)[colnames(datos) == "GENERO"] <-"SEXO"   
#levels(datos$SEXO) <- c(T,F)
#levels(datos$SEXO) <- c("Masculino","Femenimo")
levels(datos$SEXO) <- list("Masculino" = "Hombre", "Femenino" = "Mujer")
datos$IMC<-as.double(datos$IMC)
datos$TAD<-as.numeric(datos$TAD)

datos$IMC[is.na(datos$IMC)==TRUE] <-ceiling(mean(datos$IMC,na.rm = TRUE))
datos$TAD[is.na(datos$TAD)==TRUE] <-ceiling(mean(datos$TAD,na.rm = TRUE))


print(datos)  
##    PACIENTE EDAD COLESTEROL IMC TAD      SEXO
## 1         1   42        292  57  13 Masculino
## 2         2   64        235  50  11 Masculino
## 3         3   47        200  34   9 Masculino
## 4         4   56        200  38   8  Femenino
## 5         5   54        300  58   3 Masculino
## 6         6   48        215  22   6 Masculino
## 7         7   57        216  13   2  Femenino
## 8         8   52        254  42   7 Masculino
## 9         9   67        310  27   4 Masculino
## 10       10   46        237  16   7  Femenino
## 11       11   58        220  34   7 Masculino
## 12       12   62        233  46   8  Femenino
## 13       13   49        240  44  11 Masculino
## 14       14   56        295  18  12 Masculino
## 15       15   63        310   2  12 Masculino
## 16       16   64        268  49  11  Femenino
## 17       17   67        243  24  10  Femenino
## 18       18   49        239  17   8  Femenino
## 19       19   53        198  41   8  Femenino
## 20       20   59        218   2  10 Masculino
## 21       21   65        215  26   7 Masculino
## 22       22   67        254  47   4 Masculino
## 23       23   49        218  35  10 Masculino
## 24       24   53        221  31   9  Femenino
## 25       25   57        237  33  11  Femenino
## 26       26   47        244  25  10  Femenino
## 27       27   58        223  30   7  Femenino
## 28       28   48        198  36  10  Femenino
## 29       29   51        234  41   9 Masculino
## 30       30   49        175  45   9 Masculino
## 31       31   68        230  51   7  Femenino
## 32       32   58        248  15   8  Femenino
## 33       33   54        218  39  12  Femenino
## 34       34   59        285  55   3 Masculino
## 35       35   45        253  29   8 Masculino
## 36       36   53        187  23   9 Masculino
## 37       37   43        208  43   5 Masculino
## 38       38   57        246  12   9 Masculino
## 39       39   64        275  19  12 Masculino
## 40       40   43        218   7   8 Masculino
## 41       41   47        231  38   8 Masculino
## 42       42   58        200  37  11 Masculino
## 43       43   58        214  39   8  Femenino
## 44       44   48        230  28   7  Femenino
## 45       45   62        280  40   3  Femenino
## 46       46   54        198  12   5 Masculino
## 47       47   67        285  52  12 Masculino
## 48       48   68        201  14   9 Masculino
## 49       49   55        206   6   5  Femenino
## 50       50   50        223  21   8  Femenino
## 51       51   53        290  59  12  Femenino
## 52       52   63        315  53   3  Femenino
## 53       53   60        220  48   9  Femenino
## 54       54   46        230   9   8 Masculino
## 55       55   45        175  18   7 Masculino
## 56       56   53        213  19   7  Femenino
## 57       57   59        220  11   5 Masculino
## 58       58   62        287  59  12 Masculino
## 59       59   60        290  60  11 Masculino
## 60       60   62        209  10   8 Masculino
## 61       61   58        290  54   9 Masculino
## 62       62   57        260  53  12 Masculino
## 63       63   49        202  10   9  Femenino
## 64       64   61        214   5  11  Femenino
## 65       65   52        231   8   8  Femenino
## 66       66   59        280  56   3 Masculino
## 67       67   50        220  32   7 Masculino
## 68       68   46        233  20   8 Masculino
## 69       69   44        215   4   7 Masculino
## 70       70   60        202   3   5  Femenino

B. Análisis exploratorio

1b. Mostral la estructura (str) y resumen de la base de datos (summary) (mínimo, media, máximo, desviación estándar, primer cuartil de cada variable numérica y la frecuencia en el caso de variables categóricas)

ESTRUCTURA

#datos[datos==","]<- " "
str(datos)
## 'data.frame':    70 obs. of  6 variables:
##  $ PACIENTE  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ 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  57 50 34 38 58 22 13 42 27 16 ...
##  $ TAD       : num  13 11 9 8 3 6 2 7 4 7 ...
##  $ SEXO      : Factor w/ 2 levels "Masculino","Femenino": 1 1 1 2 1 1 2 1 1 2 ...

RESUMEN

summary(datos)
##     PACIENTE          EDAD         COLESTEROL         IMC       
##  Min.   : 1.00   Min.   :42.00   Min.   :175.0   Min.   : 2.00  
##  1st Qu.:18.25   1st Qu.:49.00   1st Qu.:214.2   1st Qu.:16.25  
##  Median :35.50   Median :56.00   Median :230.0   Median :31.50  
##  Mean   :35.50   Mean   :55.24   Mean   :236.8   Mean   :30.77  
##  3rd Qu.:52.75   3rd Qu.:60.00   3rd Qu.:254.0   3rd Qu.:44.75  
##  Max.   :70.00   Max.   :68.00   Max.   :315.0   Max.   :60.00  
##       TAD                SEXO   
##  Min.   : 2.000   Masculino:41  
##  1st Qu.: 7.000   Femenino :29  
##  Median : 8.000                 
##  Mean   : 8.157                 
##  3rd Qu.:10.000                 
##  Max.   :13.000

DESVIACIONES ESTÁNDAR

sd(datos$EDAD)
## [1] 7.092424
sd(datos$COLESTEROL)
## [1] 34.60467
sd(datos$IMC)
## [1] 17.25862
sd(datos$TAD)
## [1] 2.684554

2b. genere diagramas de caja para variables continuas y diagramas de barras para variables discretas, describir resultados.

CONTINUAS 
par(mfrow=c(1,3))
boxplot(datos$COLESTEROL,col=heat.colors(1,alpha = 0.2), ylab="mg/dL",main= "Colesterol")


boxplot(as.double(datos$IMC),col=topo.colors(3,alpha = 0.2), ylab="Unidades",main= "Índice de masa corporal (IMC)")
boxplot(as.integer(datos$TAD),col=terrain.colors(8,alpha = 0.2), ylab="mmHg",main= "Tensión Arterial Distólica (TAD)")

En el diagrama de caja de la variable “Colesterol” hay un dato atípica en el extremo superior; mientras que, en la variable “TAD” hay un dato atípico en el extremo inferior.

En cuanto a la dispersión se observa que en “Colesterol” los datos comprendidos entre el 25% y el 50% de la población está más dispersa que entre el 50% y el 75%. En cambio sucede lo contrario con los valores de la Variable “TAD”. Y en la la varible “IMC”, la dispersión es similar en ambos lados de la mediana.

DISCRETAS 
hist(datos$EDAD,breaks = "Sturges",col=topo.colors(6,alpha = 0.4), xlab="Años",ylab = "Número de pacientes",main= "Edad \n",freq=T,labels = T,ylim = c(0, 20) )

El rango de edad más frecuente es entre 55-60 años, en tanto que los rangos entre 40-45 y 65-70 son los menos frecuentes dentro de la población de pacientes.

3b. Calcule la correlación entre la variable dependiente y cada una de las variables explicativas (numéricas).

library(readr)
library(corrplot) 
library(mlbench) 
library(ggplot2)
library(Amelia) 
library(plotly) 
library(reshape2) 
library(caret)  
library(caTools) 
library(dplyr)
cor(select(datos,-PACIENTE,-SEXO))
##                   EDAD COLESTEROL        IMC         TAD
## EDAD        1.00000000 0.33144210 0.16302266 -0.02182395
## COLESTEROL  0.33144210 1.00000000 0.43064453  0.06076694
## IMC         0.16302266 0.43064453 1.00000000  0.08993568
## TAD        -0.02182395 0.06076694 0.08993568  1.00000000
corrplot(cor(select(datos,-PACIENTE,-SEXO)))

De acuerdo a los datos de correlación el 2Colesterol" presenta más realción con las variables “IMC”(0.43), “Edad”" (0.16), que con “TAD” (0.06). Consecuentemente el “IMC” tiene mayor relación con las variables “Colesterol” (0.43) y “Edad” (0.16), sin embargo con “TAD”(0.08) la correlación es menor.

4b. Muestre el efecto de las variables independientes con respecto a la variable dependiente.

par(mfrow=c(2,2))
datos%>%select(c(COLESTEROL, EDAD,IMC,TAD)) %>%
  melt(id.vars = "COLESTEROL") %>%
  ggplot(aes(x = value, y = COLESTEROL, colour = variable)) +
  geom_point(alpha = 0.7) +
  stat_smooth(aes(colour = "black")) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  labs(x = "VARIABLES", y = "COLESTEROL ") +
  theme_minimal()

C. Análisis ANOVA

Considere una variable categórica y realice un análisis ANOVA (como el revisado en clase).

Mujeres<- datos[,3][datos$SEXO=="Femenino"]
length(Mujeres)
## [1] 29
Hombres<- datos[,3][datos$SEXO=="Masculino"]
length(Hombres)
## [1] 41
Hombres<-sample(Hombres,29)
length(Hombres)
## [1] 29

Generamos dos variables: COLESTEROL es la variable respuesta y SEXO es la variable factor :

colesterol <- c(Mujeres,Hombres)
sexo <- as.factor(c(rep(c("Mujeres","Hombres"), each =29))) #factor
Exploramos los datos de la muestra: 
boxplot(colesterol ~ sexo, col = c("khaki1", "lightblue1"), ylab = "Índice de colesterol (mg/mL)") 

tapply(colesterol, sexo, mean)
##  Hombres  Mujeres 
## 245.1379 230.7931

Tabla de ANOVA

fm = aov( lm(colesterol ~ sexo) ) 
summary(fm)
##             Df Sum Sq Mean Sq F value Pr(>F)
## sexo         1   2984    2984    2.58  0.114
## Residuals   56  64768    1157

Elementos generados en el ANOVA:

names(fm)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

:. Las medias son estdísticamente iguales: 0.593>0.05

Prueba de post-hoc

qf(0.05, 1, 56, lower.tail = F)
## [1] 4.012973

Valores del estadístico > 4.012973estarán incluidos en la región de rechazo. En nuetro caso F=0.288 es menor que el valor crítico obtenido.Se acepta la H0, entonces las medias seran estadísticamente iguales.

Intervalo de confianza de la media del nivel de colesterol en Mujeres y Hombres, con un nivel de confianza del 95%:

Mujeres

media <- mean(colesterol[sexo =="Mujeres"]) 
valor_t <- pt(0.05/2, 56) 
sp <- sqrt(1169) 
ee  <- valor_t * (sp/ sqrt(29))

Límites de intervalo:

#Limite uperior
media + ee 
## [1] 234.0307
#Limite inferior
media - ee
## [1] 227.5555

Hombres

media <- mean(colesterol[sexo =="Hombres"]) 
valor_t <- pt(0.05/2, 56) #prueab T, depende del grado de confianza
sp <- sqrt(1169)  #desviación típica de la varianza muestral común,46 esti de la sd en anov  arri
ee  <- valor_t * (sp/ sqrt(29))  #error de estimación , 6 n de cada grupo
media#PREDICIION DEL NUMERO DE INSECTOS DE UN COLOR A APRTI DE LA MEUSTRA DADA
## [1] 245.1379

Límites de intervalo:

#Limite uperior
media + ee 
## [1] 248.3755
#Limite inferior
media - ee
## [1] 241.9004

Prueba de TUKEY HSD

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(colesterol ~ sexo))
## 
## $sexo
##                      diff       lwr      upr     p adj
## Mujeres-Hombres -14.34483 -32.23589 3.546232 0.1138607

El valor de p-adjunto=0.593424 crorrobora que las medias son estadísticamente iguales, ya que este valor se tiende más a 1 que a 0.

plot(intervals,mean.labels = TRUE, connect = FALSE) # si cercano a 0 son diferentse, si son

Si un intervalo no contiene cero, las medias correspondientes son significativamente diferentes.Este caso muestra lo contrario.

Validación del modelo ANOVA

A partir de los residuos del modelo comprobaremos si el modelo ANOVA es adecuado. Los supuestos que se deben cumplir son tres: independencia, homocedasticidad y normalidad

Independencia

plot(fm$residuals)

Si los rediduos son buenos, las predicciones son buenas; mientras mas dispersos mejor, porque no estarían correlacionados.

Normalidad

summary(fm$residuals)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -70.138 -25.052  -7.793   0.000  14.448  84.207
boxplot(fm$residuals)

hist(fm$residuals)

qqnorm(fm$residuals) 
qqline(fm$residuals)

La distribución tiende a ser Normal, aunque hay un dato atípico.

Prueba de Shapiro-Wilk para Normalidad

El test de Shapiro-Wilk indica si hay evidencia suficiente para rechazar la hipótesis nula (normalidad de los residuos)

shapiro.test(fm$residuals)  
## 
##  Shapiro-Wilk normality test
## 
## data:  fm$residuals
## W = 0.96303, p-value = 0.07451

Los residuos son normales pues el p-valor=0.03727 es menor a 0,05.

Homocedasticidad

Los gráficos y descriptivos nos informan si se verifica la igualdad de varianzas en los grupos descritos:

boxplot(fm$residuals~sexo, col = c("khaki1", "lightblue1")) 

desviaciones <- tapply(fm$residuals, sexo, sd)
desviaciones
##  Hombres  Mujeres 
## 39.01806 28.12012

Comparando la desviación máxima con la mínima obtenemos una orientación sobre la falta de homocedasticidad (>2 aproximadamente)

max(desviaciones) / min(desviaciones) #no hay hocedaticidad, las varinzas no osn iguales    
## [1] 1.387549

:.Las varianzas son estadísticamente iguales, pues se cumple que la relacion entre la desviación máxima y mínima es menor a 2.

Prueba de Barlett

El test de Bartlett indica que no tenemos evidencia suficiente para rechazar la hipótesis nula (las varianzas son iguales).

bartlett.test(fm$residuals ~ sexo)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  fm$residuals by sexo
## Bartlett's K-squared = 2.8999, df = 1, p-value = 0.08859

Dado que el p-valor=0.1593, menor que 0,05; no se encuentra evidecnai suficiente para rechazar la H0.

Kruskal-Wallis y pruebas post-hoc*

¿Qué hipótesis contrasta el test ANOVA?

Ho: las medias son iguales en todas las poblaciones

Ha: hay alguna media distinta

¿Qué hipótesis contrasta la prueba de Kruskal-Wallis?

Ho: la variable respuesta es la misma en todas las poblaciones valoradas

Ha: la variable respuesta es mayor en alguna de las poblaciones

Cuando no se cumplen las hipótesis exigidas por el modelo ANOVA, es posible utilizar la prueba no paramétrica Kruskal-Wallis: ¿hay diferencias significativas entre las poblaciones?

kruskal.test(colesterol, sexo)#ANOVA mulñtivariable
## 
##  Kruskal-Wallis rank sum test
## 
## data:  colesterol and sexo
## Kruskal-Wallis chi-squared = 2.1164, df = 1, p-value = 0.1457

Con un nivel de significación alfa = 0.05. Se acepta la Ha debido a que el p-valor=0.6916, mayor a 0,05 . Es decir que la variable respuesta es mayor en alguna de las poblaciones

Bajo la Ho el estadístico de contraste H del test de Kruskal-Wallis se distribuye como una Chi-cuadrado de grados de libertad (I-1) (donde I es el número de grupos que disponemos). Así obtenemos el cuantil buscado:

qchisq(0.05, 2-1, lower.tail = F)
## [1] 3.841459

Valores del estadístico > 3.841459 estarán incluidos en la región de rechazo. En nuetro caso 0.15738 es menor que el valor crítico obtenido.

D. Construir el modelo y predicción

1d. Considerando los cálculos anteriores genere el modelo de regresión lineal múltiple que mejor se ajuste a los datos

Train Y Test Data

set.seed(123) 
split <- sample.split(datos,SplitRatio =0.75) #Porcentaje  de trainimg data: 75%

train <- subset(datos,split==TRUE) 
test <- subset(datos,split==FALSE)

Entrenando el modelo

#Primer modelo
#model <- lm(COLESTEROL ~   EDAD +  IMC + TAD + SEXO, data = train)

#Modelo corregido

model <- lm(COLESTEROL ~   -1 + EDAD +  IMC , data = train)

summary(model)
## 
## Call:
## lm(formula = COLESTEROL ~ -1 + EDAD + IMC, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.543 -19.796   4.582  24.404  89.187 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## EDAD   3.8438     0.1943  19.780   <2e-16 ***
## IMC    0.7258     0.3265   2.223   0.0313 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.73 on 45 degrees of freedom
## Multiple R-squared:  0.9768, Adjusted R-squared:  0.9757 
## F-statistic: 946.6 on 2 and 45 DF,  p-value: < 2.2e-16

Ecuación del modelo 

model$coefficients
##     EDAD      IMC 
## 3.843842 0.725823
paste("COLESTEROL=",round(model$coefficients[1],2)," + ",
      round(model$coefficients[2],2),names(model$coefficients[2])
      )
## [1] "COLESTEROL= 3.84  +  0.73 IMC"

2d. Realice un análisis detallado de los residuos

Visualizacion del modelo

Regresión lineal de los residuos(diferencia entre el valor observado de la variable dependiente y el valor predicho de y))

res <- residuals(model)

# Convertir residuos en un DataFrame
# 
res <- as.data.frame(res)
ggplot(res,aes(res)) +  geom_histogram(fill='green3',alpha=0.5)

par(mfrow=c(2,2))
plot(model)

Predicciones

Se prueba el modelo prediciendo en un conjunto de datos de prueba.

test$predicted.COLESTEROL <- predict(model,test)

pl1 <-test %>% 
  ggplot(aes(COLESTEROL,predicted.COLESTEROL)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +  #cacl el intrv de confianza=(meda-desv s),tendencia datos
  xlab('Valor real de COLESTEOL') +
  ylab('Preciddión de COLESTEROL')+
  theme_bw()

ggplotly(pl1)

Evaluación del modelo

nuestros valores predichos.

error <- test$COLESTEROL-test$predicted.COLESTEROL
rmse <- sqrt(mean(error)^2) #RAIZ MEDIA CUADRATICA, si RMSE es cercano a 1 el odelo es bueno
rmse
## [1] 1.964042

3d. Analice la significancia de las variables y los parámetros individuales.