#LAB ESTADISTICOS PRINCIPALES
#######################################
#Este laboratorio incluye los siguientes tópicos. 
##1.Media
##2.Moda
##3.Mínimo, máximo y rango
##4.Suma y conteo
##5.Varianza y desviación
##5.Cuantiles 
##6.Correlación y Causación
##7.Tablas de contingencia 
#######################################
#Cargando librerías 
rm(list=ls())
require(DescTools)
require(psych)
require(modeest)
require(dplyr)
require(ggplot2)
require(Epi)
require(MASS)
require(corrplot)
require(reshape2)
require(scales)
require(Epi)
require(pubh )
require(ggmosaic)

#########MEDIA DE UN SET DE DATOS ###########
#Existen varios tipos de media, la media aritmética es la mas usada, pero no la única.
#La media aritmética es el promedio simple de los valores de la serie, es utilizado en 
#una serie de cálculos como se verá más adelante. 
#No confundir con la mediana que es el valor del medio.
#La función en R es 
#mean(x, trim = 0, na.rm = FALSE, .)

#Vamos a usar el dataset midwest, contiene una serie de estadísticas interesantes 
#de la población de los USA . 
#Mas información en https://www.rdocumentation.org/packages/ggplot2/versions/3.3.2/topics/midwest
#Cargamos el dataset con : 
data(midwest)
#Para referirnos a los campos (columnas)  como variables podemos usar el comando:
attach(midwest)
#Para obtener la clase de objeto que es usaremos:
class(midwest)
#El dataset es de tipo tibble, una subclase de dataframe pero tiende a "quejarse" mucho 
#prefiero verlo como data.frame, para esto uso un cast
mid <- as.data.frame(midwest)
class(mid)
#ahora si, veamos que columnas tiene
names(mid)
#veamos cuantas filas tiene el dataset 
NROW(mid)
#Primero me interesa saber cual es la media de los estados respecto
#a personas con educación superior:
mean(mid$percollege)
#Note el signo de $ indica la columna del dataset

########MEDIA PONDERADA O MEDIA CON PESOS ############### 
#La media aritmética ponderada es similar a una media aritmética, excepto que
#en lugar de que cada uno de los puntos de datos contribuya igualmente al
#promedio final, algunos puntos de datos contribuyen más que otros. 
#La noción de media ponderada juega un papel importante  en la estadística descriptiva.
#La función es : 
#Mean(x, weights = NULL, trim = 0, na.rm = FALSE, ...)
#Veamos un ejemplo con una estadística de nacimiento de niños desde 1901
data(B.dk)
View(B.dk)
#Nos interesa saber cual es el promedio de nacimientos 
mean(B.dk$m)
#Y de niñas
mean(B.dk$f)
#Pero ahora nos interesa saber  cuantos niños varones nacerán en enero 2019 , podría
#tomar el promedio mensual del 2009
nac2018 <- B.dk[B.dk$year==2018,3 ] # el 3 representa la tercera columna (niños)
round(mean(nac2018))
#Pero el problema es que esta fórmula toma el promedio mensual de todo el año 2018,
#sin embargo , es más probable que se acerque a los valores de los últimos meses del 2018
#Una forma de hacerlo sería una media ponderada.
pesos <- c(1,1,1,1,2,2,2,2,3,3,4,4)

#La media ponderada reside en el paquete DescTools

round(Mean(nac2018, weights=pesos) ) #note Mean con mayúscula 

#############MEDIA GEOMÉTRICA ############################
#Otro tipo de media es la media geométrica, ésta toma en cuenta el promedio multi anual 
# y es especialmente usada en ambientes donde la incidencia multi anual es significativa. 
#La media geométrica se calcula como la raiz n de la multiplicación del conjunto de datos
#Esta función reside en el paquete psych

data(gmortDK)  # tasas y causas  de mortalidad en Dinamarca 
View(gmortDK)
#Por ejemplo calcular el porcentaje de mortalidad  por problemas cerebrovasculares (r7)
#para los períodos 88 al 92 (ver significado de variable per código 88)
data_per88 <- gmortDK[gmortDK$per==88 , 13]  #13 indica la columna 13 (r7) que es causas cerebrovasculares
#la sintaxis es : geometric.mean(x,na.rm=TRUE)
#Note que debido a la formula que usa no permite que la data tenga valores de 0
geometric.mean(data_per88)


########### MÍNIMO, MÁXIMO Y RANGO #############
max(mid$percollege)
min(mid$percollege)
#Nos llama la atención el mínimo, me interesa saber en que estado se 
#produce este mínimo. Utilizaremos subsetting.
mid[mid$percollege==min(mid$percollege),"state" ]
mid[mid$percollege>= 40,"state" ]
mid[mid$percollege <= 20,"state" ]
# Subsetting nos permite ver exactamente filas y columnas de interés.
#El formato básico es : 
#dataset[filas, columnas]
#Pero podemos poner también condicionales como en el ejemplo: 
#Si nos interesa una cantidad de filas y/o colunas podemos usar vectores en su lugar
vector_cols<- c("state", "poptotal", "popdensity","percblack","percwhite","percasian", "percother","percamerindan") 
mid[mid$percollege==min(mid$percollege),vector_cols]
#Podemos concluir que el mínimo se obtiene en un estado donde la población
#IndioAmericana es mayoritaria.


#La función range me da  los dos estadísticos anteriores 
range(mid$percollege)
#Note que el resultado es un vector numérico

##############LA MODA ####################
#Hay veces que la media no es suficiente para explicar el comportamiento de los datos.

moda_densidad<- mlv(mid$popdensity, method="mfv")
moda_densidad
#Nos da un valor, pero en realidad no nos dice mucho más

#Vamos a usar la data nickel, Los datos se refieren a la afectación  por exposición al níquel en 
#fundiciones  de níquel, y su análisis de un grupo de trabajadores en el sur de Gales 
#Por ejemplo se requiere saber cual es el valor más repetido en dicho grupo. 
#Este es el concepto de MODA. La moda es el valor más repetido de una serie de datos .
# La función para  la moda se halla en el paquete modest y su sintaxis es : 
#mlv(x, na.rm = FALSE, method="mfv")
#Analicemos la densidad poblacional de estos estados, nos interesa saber
#la moda

data(nickel)
View(nickel)
moda_niquel <- mlv((nickel$age1st), method="mfv")
moda_niquel  #la mayoría no tenían exposición al níquel 
#Note que nos esta diciendo que existen como 3 picos  donde se logran máximos

######### FUNCIONES DE DISPERSIÓN VARIANZA Y DESVIACIÓN #########
#Un tema que frecuentemente  queremos encontrar es que tan dispersos están
#los datos.Esto nos ayuda a entender que tan normal es una nueva observación
#cuando esta se compara con las observaciones anteriores. 
#La primera función para este propósito es la varianza 

#Varianza es una medida de dispersión que representa la variabilidad de una 
#serie de datos respecto a su media. 
#Una serie de datos podría tomar un número infinito de  valores, pero en la
#practica cuando se analiza un parámetro  de interés tienden a estar 
#más o menos dispersos dentro de un rango. 
#La varianza indica esta medida de dispersión. 

#La formula de la dispersión es la diferencias de las observaciones con la media
#elevada al cuadrado divido por la cantidad de observaciones. 


#Ejemplo : Nos preguntamos que tan dispersa está población (densidadd poblacional)
#habría que revisar de cuantos estados tenemos información
unique(mid$state)
var(mid$popdensity)

#calculamos manualmente  a ver si coincide, note la corrección de Bessel 
var1 <- (sum( (mid$popdensity - mean(mid$popdensity) )^2)  /(NROW(mid)-1))
var1
sqrt(var1)
sd(mid$popdensity)
#podemos ver gráficamente la distribución de la densidad de la población 
ggplot(mid , aes(x=popdensity)) +
  geom_density()

#Nos sorprende este número tan alto, veamos los extremos 
max(mid$popdensity)
min(mid$popdensity)
#Entonces concluimos que nuestra población está muy dispersa, 
#tenemos estados muy poblados y otros muy poco
sd(mid$popdensity)
#no se preocupe mucho si no entiende la lógica de la función, 
#veremos más en detalle cuando veamos  distribución

#Queremos analizar la distribución de población blanca por condado


############CUANTILES ###################
#Los desgloses anteriores están bien, pero en estadística generalmente se usa
#la revisión por cuantiles. 
#Un cuantil es aquel punto que divide la distribución de una variable en 
#intervalos, por tanto, no es más que una técnica estadística para separar
#los datos de una distribución.


vect_cuantiles <- c(0.25 , 0.5, 0.75 )
quantile(mid$percwhite, probs=vect_cuantiles)
#esto nos indica que el cuantil más bajo  ya es de 94% de población blanca 

#O la distribución de adultos 
#vect_cuantiles <- c(0.10 , 0.25, 0.50, .9 )
quantile(mid$percwhite, probs=vect_cuantiles)
quantile(mid$percblack)

#Podríamos ver estos cuantiles en forma  gráfica 

ggplot(mid, aes(x=percwhite)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=4) + coord_flip()


###########COUNT , LENGTH Y NROW ##########
#Comúnmente nos interesa saber cuantas observaciones cumplen con cierta condición
#Existen algunas funciones que son importante entenderlas 

#Nos interesa saber cuántos condados tiene población de adultos 
#mayores en situación pobreza
count(mid[mid$percelderlypoverty>=10,])
#Note que el subsetting nos devuelve un data frame y aquí count funciona bien
count(mid[mid$percelderlypoverty>=10, "county"])
#Aquí count() no funciona porque el subsetting nos ha devuelto un vector, en este caso usar length 
length(mid[mid$percelderlypoverty>=10, "county"])
#o NROW
NROW(mid[mid$percelderlypoverty>=10, "county"])
#Por último summary nos entrega todos los estadísticos anteriores
summary(mid$poptotal)




#CORRELACIÓN 
#=====================================
#Correlación es la medidas de relación que tiene dos variable 
#Existen diferente metodologías para el ´calculo de correlación 
#Métodos paramétricos  y métodos no paramétricos
#Dentro de los métodos paramétricos está la correlación de Pearson que veremos a continuación 

#Vamos a usar un dataset de  carros,
#primero cargamos  la data 
data(mtcars)
#Examinemos la data 
names(mtcars)
class(mtcars)
#Usamos subsetting para  que no salga una matriz de todas las columnas 
#contra todas las columnas. Seleccionamos solo las priemras 4 
#El estudiante podria seleccionar otras de su interés 

cor(mtcars[,1:4])
#Analice los resultados. Indique cuales son las variables que de alta correlación


#Veamos otro ejemplo con el dataset economics
#Cargamos la data 
data(economics)
#Analizamos la data 
class(economics)
NROW(economics)
names(economics)
head(economics)
# las columnas son : fecha , 
#date Month of data collection
#pce personal consumption expenditures, in billions of dollars, http://research.stlouisfed.org/fred2/series/PCE
#pop total population, in thousands, http://research.stlouisfed.org/fred2/series/POP
#psavert personal savings rate, http://research.stlouisfed.org/fred2/series/PSAVERT/
# uempmed median duration of unemployment, in weeks, http://research.stlouisfed.org/fred2/series/UEMPMED
#unemploy number of unemployed in thousands, http://research.stlouisfed.org/fred2/series/UNEMPLOY

#queremos ver si estas variables están relacionadas 

cor(economics$pce, economics$psavert)
#Que indica este valor? 
cor(economics[,c(2, 4:6)])
#Que indica esta matriz?
cor_economics <-cor(economics[,c(2, 4:6)])
cor_economics
#esto se ve mejor en un mapa de calor ,pero para esto habrá que  poner en un
#una sola columna , como está ggplot no lo puede visualizar adecuadamente 

econ_melted <- melt(cor_economics, varname=c("X","Y"), value.name="Correlacion")
econ_melted
class(econ_melted)
#ahora le ploteamos como heat map
#antes le ordenamos para que nos salga bonito el gráfico
econ_melted <- econ_melted[order(econ_melted$Correlacion),]
econ_melted <- econ_melted[order(econ_melted$Correlacion, decreasing = TRUE),]

ggplot(econ_melted, aes(x=X , y=Y)) + 
  geom_tile(aes(fill=Correlacion))+
  scale_fill_gradient2(low=muted("steelblue"), mid="white" , high="red")+
  theme_minimal() +
  labs(x=NULL , y=NULL)

#Un problema recurrente al momento de correlacionar son los 
#datos con NA, hay opciones para esto 

#cor(df ,use="_______")
# donde el parámetro puede ser :
#everything : usa las filas completas no las NAs
#all.abbs todas las filas deben estar completas o dará un error
#na.or.complete usa todo pero retornar varios  NAs
#pairwise.complete.obs  las 2 filas relacionadas deben estar presentes y sin nas

#VEAMOS UN EJEMPLO DE CORRELACION  ENTRE MATRICES 

mat1 <- c(9,9,NA,3,NA,5,8,1,10,4)
mat2 <- c(2,NA,1,6,6,4,1,1,6,7)
mat3 <- c(8,4,3,9,10,NA,3,NA,9,9)
mat4 <- c(10,10,7,8,4,2,8,5,5,2)
mat <- cbind(mat1,mat2,mat3,mat4)
cor(mat)
cor(mat , use="everything")
cor(mat , use="all.obs")
cor(mat , use="complete.obs")
cor(mat , use="na.or.complete")
cor(mat , use="pairwise.complete.obs")


##########################TABLAS DE CONTINGENCIA  ###################
#Una tabla de contingencia es una tabla en la que se resume la distribución de 
#frecuencias conjuntas de dos variables. En esta tabla, cada fila representa una 
# variable X (o una categoría de dicha variable) y cada columna representa una categoría
#de la variable Y ,cada celda de la tabla representa la frecuencia con la que aparecen
#conjuntamente las categorías de la fila y la columna en la que está ubicada.

#Las tablas de contingencia también son conocidas como tablas cruzadas, 
#tablas de frecuencias de doble entrada o tablas de frecuencias bidimensionales.

#Las tablas de frecuencia son similares a las obtenidas con un aggregate(), pero nos 
#provee informaciónn adicional 

#Existen varios paquetes que proveen esta funcionalidad, éstos se diferencian en 
#cuanto al formato de salida, opciones e indicadores que proveen 


######Usando base::table

#Creamos un primer vector 
X <- sample(c("Si", "No"), 10, replace = TRUE)

# Generamos datos aleatorios para la variable Y
Y <- sample(c("Europa", "America", "Africa"), 10, replace = TRUE)

#Creamos un data.frame para poder visualizar
df <- data.frame(sino= X, continente=Y)
#En base a la visualización de este df cuente las ocurrencias de la variable 
#"sino"  para cada continente 

#Ahora vemos la misma información como table de contingencia en base al paquete 
#base, no es muy explicativo pero hace su trabajo 

tabla <- table(X, Y)
View(tabla)

######Usando la función Epi::stat.table
#Veamos un ejemplo, usaremos un dataset que contiene data sobre la Oncocercosis en 
#Sierra León
#Cargamos la data 
data(Oncho)
#Veamos las primera columnas 
View(Oncho)
Oncho %>% head()
#veamos que clase de objeto es 
class(Oncho)  # es una tibble 
#Convertimos a data.frame 
oncho <- as.data.frame(Oncho)
#entendamos el dataset 
names(oncho)
View(oncho)

#Su sintaxis es : stat.table(index, contents = count(), data, margins = FALSE)
#Donde: 
#index es la columna(2) por las cuales se clasificará y contará  su frecuencia 
#y el conteo se hace con la función count()


#Veamos la distribución en base al area 
Epi::stat.table(area, contents = count(), Oncho, margins = TRUE)
#Veamos la distribución  en base al área y al grupo de edades
stat.table(list(area,agegrp), contents = count(), Oncho, margins = TRUE)

### Usando pubh cross_tbl 

#Pubh también tiene su propia función cross_tbl para  tablas de contingencias
#Nos interesa saber cual es la relación  de infecciones en la sabana vs. la selva
#como no nos interesa los otras columnas  vamos a seleccionar solo lo que nos interesa
#La función pubh::cross_tbl permiten crear una tabla de frecuencias más completa 
#Su sintaxis es : 
#cross_tbl( data,by,head_label = " ",bold = TRUE,show_total = TRUE,p_val = FALSE,pad = 3,method = 2,)
#Donde :
#Data es el data frame bajo estudio 
#by es la variable(s) categórica(2) (factores) a examinar 
#p_val es una variable lógica que indica si se debe presentar los valores p 
#method 1= calcular media y sd , 2= medias y IQR
#
oncho %>%
  dplyr::select(mf, area) %>% 
  cross_tbl(by = "area") %>%
  theme_pubh(2)

#Que nos indica esta tabla: Que en la proporción  de infectados y no infectados 
#es similar, pero en la selva la situación se complica y hay muchos más infectados.

#Veamos en forma gráfica 
onchosel <- oncho %>% dplyr::select(mf, area)
onchotab <- as.data.frame(table(onchosel$mf, onchosel$area))
names(onchotab) <- c("Infected", "Area" , "Freq")
ggplot(data=onchotab) +
  geom_mosaic(aes(weight=Freq, x=product(Infected, Area)),color="blue")   +
  xlab("Area") +
  ylab("Infected")




#Correlación no implica causación
#PROBLEMA 1:
#Dado el dataset mtcars encuentre relaciones entre variables que indiquen 
#una correlación aunque esto no indique causación 


#PROBLEMA 2
#Dado el siguiente dataset, cargue la data, entienda la data, 
housing <- read.csv(url("https://www.jaredlander.com/data/housing1.csv"))
#Encuentre las correlaciones entre las variables numérica, indique
#que relaciones son fuertes y cuales son débiles, indique si hace sentido 


#TAREA MODULO 2
data(mtcars)
data(economics)


#cONTESTE LAS SIGUIENTES PREGUNTAS 
# La media de millas por galón  mpg de los vehículos
mean(mtcars$mpg)

#-       Cuál es el  máximo hp   de los vehículos 
max(mtcars$hp)

#-       Cuál es el tiempo mínimo de cuarto de milla logrado 
min(mtcars$qsec)

#-       Indique los vehículos  que se hallan en el primer cuartil de mpg
valorq <- quantile(mtcars$mpg , probs=c(0.25) )
mtcars[mtcars$mpg < valorq , ]

#-       Cuál es la varianza y desviación estándar de disp.
var(mtcars$disp)
sd(mtcars$disp)

# -       Cual es la media y la media geométrica de psavert para los años 2010-1015 (deberá investigar sobre la fórmula de la media geométrica)
eco <- economics[economics$date> "2010-01-01" , 4]
mean(economics$psavert)
exp(mean(log(eco$psavert)))

#PARTE 2
cor(economics[,2:6])

#-       Cual es la covarianza entre el número de desempleados y el tiempo promedio que esta desempleado
cov(economics$uempmed, economics$unemploy)

#-       Cual fue el porcentaje más alto y más bajo de desempleados y en que fechas ocurrió. Para esto deberá crear una nueva variable en el dataset con el siguiente snippet completando la última sentencia  						


economics$pdes <- economics$unemploy/economics$pop


minimo <- min(economics$pdes)
maximo <- max(economics$pdes) 

economics[economics$pdes==minimo , ]$date
economics[economics$pdes==maximo , ]$date