proyectR.Rmd

---
title: "PROYECTO_R"
output:
  html_document:
    df_print: paged
---
```{r}
library(ggplot2)

```

```{r}
#setwd("/Users/carlosmichelmourradiaz/documents/R_Project")
setwd("C:/Users/pauro/OneDrive/Documentos/R/Proyecto_R")
```
Por si quieren actualizar archivos:
```{r}
# Datos abiertos de Mexico
#download.file('http://datosabiertos.salud.gob.mx/gobmx/salud/datos_abiertos/datos_abiertos_covid19.zip', '../source/dataMx.zip')
#unzip('../source/dataMx.zip', exdir = '../source/')
```
Cargar datos:
```{r}
# RUTA PAU ../source/dataMx/210114COVID19MEXICO.csv
rawMxData <- read.csv('../source/dataMx/210114COVID19MEXICO.csv')
```
Volver columna de datos binaria:
```{r}
# Vivos
rawMxData$FECHA_DEF[rawMxData$FECHA_DEF != '9999-99-99'] <- 0
# Muertos
rawMxData$FECHA_DEF[rawMxData$FECHA_DEF == '9999-99-99'] <- 1
```

Filtrar por rangos de edad:
```{r}
# De 0-17
rawMxData$EDAD[rawMxData$EDAD >= 0 & rawMxData$EDAD < 18] <- 1
# De 18-44
rawMxData$EDAD[rawMxData$EDAD >= 18 & rawMxData$EDAD < 45] <- 2
# De 45-59
rawMxData$EDAD[rawMxData$EDAD >= 45 & rawMxData$EDAD < 60] <- 3
# De 60-99
rawMxData$EDAD[rawMxData$EDAD >= 60 & rawMxData$EDAD < 100] <- 4
# Ignorar >99
rawMxData$EDAD[rawMxData$EDAD >= 100] <- 5

# Ignorar centenarios para todos los datos
avoidCentenarians <- rawMxData$EDAD != 5

```

0 -> paciente vivo
1 -> paciente muerto
# Filtrar datos
Función por enfermedad:

```{r}
# FUNCION QUE REEMPLAZA 2 (FALSE) CON 0 (FALSE)
coMor <- function(comorbiditiesList, avoidCentenarians){
  comorbidities <-matrix(1:length(rawMxData$FECHA_ACTUALIZACION[avoidCentenarians]), ncol = 1)
  comorbidities <- cbind(comorbidities,as.numeric(rawMxData$FECHA_DEF[avoidCentenarians]), rawMxData$EDAD[avoidCentenarians])
  colnames(comorbidities) <- c('INDICE', 'DEFUNCION', 'EDAD')
  comorbiditiesAux <- sapply(comorbiditiesList, function(comorbidity){
    comorbidity <- comorbidity[avoidCentenarians]
    return(replace(rawMxData[[comorbidity]],rawMxData[[comorbidity]] == 2, 0))
    })
  comorbidities <- cbind(comorbidities, comorbiditiesAux)
  allDates <- apply(comorbidities[,4:length(comorbidities[1,])], MARGIN = 1, FUN = function(paciente){
    bools <- paciente <= 1
    if(all(bools)){
      return(TRUE)
    }else{
      return(FALSE)
    }
  })
  return(comorbidities[allDates,])
}
```
# MAPAS

# Instalacion paquetes en MacOs
```{r}
if (!require("devtools")) {
    install.packages("devtools")
}
# devtools::install_github("diegovalle/mxmaps")
```
```{bash}
brew install pkg-config
brew install gdal
```
```{r}
install.packages("sf")
install.packages("sf", configure.args = "--with-proj-lib=/usr/local/lib/")
install.packages("rgdal", configure.args = c("--with-proj-lib=/usr/local/lib/", "--with-proj-include=/usr/local/include/"))
library("mxmaps")
data("df_mxmunicipio")
```
# Parseo de datos
```{r}
df_mxmunicipio_2020$value <- NA
df_mxmunicipio_2020$cases <- NA
df_mxmunicipio_2020$cases_live <- NA
df_mxmunicipio_2020$cases_dead <- NA
df_mxmunicipio_2020$state_code<-as.numeric(df_mxmunicipio_2020$state_code)
df_mxmunicipio_2020$municipio_code <-as.numeric(df_mxmunicipio_2020$municipio_code)

estados <- sort(unique(rawMxData$ENTIDAD_RES))
for(entidad in estados){
  municipios <- sort(unique(rawMxData$MUNICIPIO_RES[rawMxData$ENTIDAD_RES == entidad]))
  for(municipio in municipios){
           if(municipio < 999){
           cases <- rawMxData$FECHA_DEF[rawMxData$ENTIDAD_RES == entidad & rawMxData$MUNICIPIO_RES == municipio]
           df_mxmunicipio_2020$cases[df_mxmunicipio_2020$state_code == entidad & df_mxmunicipio_2020$municipio_code == municipio] <- length(cases)
           df_mxmunicipio_2020$cases_live[df_mxmunicipio_2020$state_code == entidad & df_mxmunicipio_2020$municipio_code == municipio] <- length(cases[cases == 1])
           df_mxmunicipio_2020$cases_dead[df_mxmunicipio_2020$state_code == entidad & df_mxmunicipio_2020$municipio_code == municipio] <- length(cases[cases == 0])
           percentage <- (length(cases)/df_mxmunicipio_2020$pop[df_mxmunicipio_2020$state_code == entidad & df_mxmunicipio_2020$municipio_code == municipio])*100
           df_mxmunicipio_2020$value[df_mxmunicipio_2020$state_code == entidad & df_mxmunicipio_2020$municipio_code == municipio] <- percentage }}
           
df_mxmunicipio_2020$value[is.na(df_mxmunicipio_2020$value)] <- 0
```

# Mapas
```{r}
mxmunicipio_choropleth(df_mxmunicipio_2020, 
                       num_colors = 5,
                       title = "Porcentaje de población infectada con COVID-19",
                       legend = "%")
```
```{r}
mxmunicipio_choropleth(df_mxmunicipio_2020, num_colors = 1,
                       zoom = subset(df_mxmunicipio_2020, metro_area %in% 
                                       "Valle de México")$region,
                       title = "Porcentaje de poblacion por municipio en Valle de México",
                       show_states = FALSE,
                       legend = "%")
```
```{r}
mxmunicipio_choropleth(df_mxmunicipio_2020, num_colors = 1,
                       zoom = subset(df_mxmunicipio_2020, state_name %in% 
                                       c("Jalisco", "Morelos", "Hidalgo"))$region,
                       title = "Porcentaje de poblacion infectada por municipio en Jalisco, Hidalgo y Morelos",
                       show_states = FALSE,
                       legend = "%")
```

# COMORBILIDADES
```{r}
#generar la matriz con los datos de las comorbilidades indicadas
prueba <- coMor(c('OBESIDAD', 'DIABETES', 'NEUMONIA', 'EPOC', 'RENAL_CRONICA', 'ASMA', 'INMUSUPR', 'TABAQUISMO', 'HIPERTENSION'))
table(prueba[,4:6])
```


```{r}
#convertir la matriz en un dataframe para acceder a las columnas
#ignorar los datos que son desconocidos como aquellos cuyo valor es 98

coFrame <- data.frame('Muertos' = as.numeric(rawMxData$FECHA_DEF[rawMxData$OBESIDAD != 98]),
           'Obesos' = rawMxData$OBESIDAD[rawMxData$OBESIDAD != 98],
           'Edad' = rawMxData$EDAD[rawMxData$OBESIDAD != 98],
           'Diabetes' = rawMxData$DIABETES[rawMxData$OBESIDAD != 98],
           'Neumonia' = rawMxData$NEUMONIA[rawMxData$OBESIDAD != 98], 
           'Hipertension' = rawMxData$HIPERTENSION[rawMxData$OBESIDAD != 98],
           'Asma' = rawMxData$ASMA[rawMxData$OBESIDAD != 98],
           'EPOC' = rawMxData$EPOC[rawMxData$OBESIDAD != 98],
           'Inmunosupresion' = rawMxData$INMUSUPR[rawMxData$OBESIDAD != 98],
           'Tabaquismo' = rawMxData$TABAQUISMO[rawMxData$OBESIDAD != 98],
           'Renal_Cronica' = rawMxData$RENAL_CRONICA[rawMxData$OBESIDAD != 98])
```
0 -> sin comorbilidad
1 -> con comorbilidad
# Minimos cuadrados
```{r}
#glm -> General Linear Model
#utilizar la función glm para llevar a cabo una regresión lineal logistica
reLin <- glm(formula = Muertos~Obesos+Edad+Diabetes+Neumonia+Hipertension+Asma+EPOC+Inmunosupresion+Tabaquismo+Renal_Cronica, coFrame, family = binomial(link = 'logit'))
summary(reLin)
```

A continuación podemos ver que todas las variable a excepción de la diabetes son significativas por lo que procederemos a interpretar dichos coeficientes

# Prueba de Significancia

H0: El modelo no es significativo (B =B2=B3=...=Bk)
Ha: El modelo si es significativo

Estadistico de prueba

```{r}
#realizamos una comparación de las desviaciones
with(reLin,null.deviance-deviance)
```

Valor P del estadistico de Prueba

Valor $\alpha$ = 0.05

```{r}
print("Valor P")
with(reLin,pchisq(null.deviance-deviance,df.null-df.residual,lower.tail = FALSE))
```

### Conclusión 
Como alpha es mayor al valor P,se rechaza H0 por Ha ya que  podemos deducir que al menos una de las variables de  nuestro modelo pueden explicar los datos.


# Interpretación de Coeficientes

Podemos observar que hay una alta relación entre padecer obesidad y morir si se enferma de covid teniendo 67% mas de probabilidad, ademas se observo que los paciente que padecen de Neumonia poseen 2.8% mas de probabilidad.

Anteriormente descartamos Diabetes debido debido a que no se encontro una relación significante que pudieran explicar las muertes.

```{r}
# Preguntar significado de valor de EDAD
intercepts <- exp(coefficients(reLin))
intercepts
```

# Clasificación de observaciones

```{r}
#install.packages("InformationValue")
```


```{r}
library("InformationValue")
```

```{r}
#añadir columna que indique la probabilidad de muerte del paciente
coFrame$PREDICT <- reLin$fitted.values
coFrame
```

```{r}
#punto de corte maximiza la clasificación de las observaciones
#el punto de corte sera en base a muertos con respecto a probabilidad para que sea mas acertado
Punto_Corte <- optimalCutoff(coFrame$Muertos,coFrame$PREDICT,optimiseFor = "Both")
Punto_Corte
```

```{r}
#indica que porcentaje de las observaciones se clasifico correctamente
Obs_correctas <- 1-misClassError(coFrame$Muertos,coFrame$PREDICT,threshold = Punto_Corte)

#esta matriz indicara como estan clasificadas las observaciones en los cuadrantes
sensitivity <- confusionMatrix(coFrame$Muertos,coFrame$PREDICT,threshold = Punto_Corte)
sensitivity
# 1-1 -> Pacientes que el modelo predijo que no moririan y no murieron
# 1-0 -> Pacientes que el modelo predijo que no moririan y murieron
# 0-1 -> Pacientes que el modelo predijo que  moririan y no murieron
# 0-0 -> Pacientes  que el modelo predijo que no moriria y no murieron
```

```{r}
#un ~78% de las observaciones es muy buen resultado
Obs_correctas
```
## Gráficas
```{r}
# CODIGO EJEMPLO
#Obtener las clasificaciones de todos los genes con valor t > 3
diffCatGenes <- genes[abs(define_tTest(expValues[,,1], expValues[,,2])) > 3,2]

#Calcular porcentajes por clasificacion
percentage <- round(table(diffCatGenes)/sum(table(diffCatGenes))*100, 2)
percentage
#Componer los identificadores con porcentaje
labs <- paste(c('Activator', 'Ambivalent', 'Represor'), percentage)
labs <- paste(labs,"%",sep="")

#Generar la grafica de pie
colorGradient <- colorRampPalette(c('lightskyblue','yellow3'))
pie(table(diffCatGenes), labels = labs, main='Proportion of statistical values by category', col = colorGradient(3))

```

```{r}
# Pie porcentajes de comorbilidades
comobiditiesList <- c('Obesos', 'Diabetes', 'Neumonia', 'Tabaquismo')
numIndividuals <- length(coFrame$Muertos)
labs <- c('Muertos', 'Vivos')
color <- c('red', 'blue')

deadsPerComorbidity <- sapply(comobiditiesList, function(co){
  deadsPerComorb <- coFrame[[co]][coFrame$Muertos == 1 & coFrame[[co]] == 1]
  return(length(deadsPerComorb))
  })

alivePerComorbidity <- sapply(deadsPerComorbidity, function(x){
  total <- numIndividuals - x
  return(total)
})

percentageCalculator <- function(x, aliveOrDeadVector){
  percentage <- round((aliveOrDeadVector[x]/(deadsPerComorbidity[x]+alivePerComorbidity[x]))*100, 2)
  return(percentage)
}

# Calculo de porcentajes para las etiquetas
deadsPercentages <- percentageCalculator(1:4, deadsPerComorbidity)
alivePercentages <- percentageCalculator(1:4, alivePerComorbidity)

#Componer los identificadores con porcentaje
labs <-c( paste('Muertos',deadsPercentages),paste('Vivos',alivePercentages))
labs <- paste(labs,"%",sep="")
labs <- sapply(1:4, function(x){
  labsMatrix <- labs[seq(x,length(labs),4)]
  return(labsMatrix)
  })

layout(matrix(c(1,2,3,4), nrow=2, byrow=TRUE))
pie(c(deadsPerComorbidity[1], alivePerComorbidity[1]), labels =labs[,1], main = 'Mortalidad de Obesos', col = color, radius =1)
pie(c(deadsPerComorbidity[2], alivePerComorbidity[2]), labels = labs[,2], main = 'Mortalidad de Diabetes', col = color, radius = 1)
pie(c(deadsPerComorbidity[3], alivePerComorbidity[3]), labels = labs[,3], main = 'Mortalidad de Neumonia', col = color, radius = 1)
pie(c(deadsPerComorbidity[4], alivePerComorbidity[4]), labels = labs[,4], main = 'Mortalidad de Tabaquismo', col = color, radius = 1)

# Barplots de interceptos
barplot(intercepts[2:5], ylim = c(0,2)); abline(h=1, lty=2, col="red")

# Boxplot de acuerdo a la edad
boxplot(coFrame$Edad)

# Pie de falsos positivos
sensitivityData <- list(Vivos = coFrame$PREDICT[coFrame$Muertos== 1], Muertos = coFrame$PREDICT[coFrame$Muertos == 0])
densMuertos <- density(sensitivityData$Muertos)
densVivos <- density(sensitivityData$Vivos)
mi.rojo <- rgb(1, 0, 0, 0.3)
mi.verde <- rgb(0, 1, 0, 0.3)
plot(densMuertos, col='red', xlim = c(0.6, 1),ylim=c(0,140), main = 'Sensibilidad y específicidad', xlab = 'Valor predictivo', ylab='Frecuencia')
lines(densVivos, col='green')
polygon(densMuertos, col=mi.rojo, border='red')
polygon(densVivos, col=mi.verde, border='green')
abline(v = 0.93, col= 'blue')
text(0.8, 135, 'Predicción\nmuertos', cex = 0.9)
text(0.975, 135, 'Predicción\nvivos', cex = 0.9)
legend('topleft',c('Muertos','Vivos','Punto de corte'),fill=c('red','green', 'blue'))


# Diagrama de Venn 

```



### Pendientes
Redactar reporte estilo artículo
Gráficas: de acuerdo a la pregunta (ej. comparar prob entre categorías, combinación de condiciones, etc), comparar con estádisticas de acuerdo a las predicciones y categorías


IDEAS:
- Calcular probabilidad del usuario (input condiciones): shiny paquete para interfaz