GitHub Documents

El presente cuaderno muestra una introducción al análisis estadístico base de información mediante el uso de R. Se muestra una caracterización inicial de los datos, luego se presenta un análisis estadístico base y se propondrá realizar un análisis exploratorio para la identificación de atípicos mediante el uso de diagramas de cajas y bigotes (boxplot) y una agrupación de los datos a escalas mensuales, anuales y multianuales, así como diferentes gráficas a lo largo de todo el proceso. Para el ejercicio de clase se utilizaron datos diarios de precipitación en la estación pluviométrica de Nataga [21050090] operada por el IDEAM y localizada en el municipio del mismo nombre en el dpto. del Huila, esta información fue descargada directamente de la plataforma DHIME.

Antes de iniciar con el tema en si, es recomendable empezar el código limpiando el espacio de trabajo y llamando de una vez todos las librerías que se vayan a necesitar. Si no se tiene instalada alguna de las librerías es necesario hacerlo con el comando “install.packages”. Así mismo, se selecciona el directorio de trabajo donde se tiene el archivo a usar.

rm(list = ls()) # Limpia el espacio de trabajo

# install.packages("xts")
# install.packages("lubridate")
# install.packages("moments")
# install.packages("ggplot2")
# install.packages("ggpubr")
# install.packages("RColorBrewer")

library(xts) # Series Temporales
library(lubridate) # Utilidades para manejo de fechas
library(moments) # Momentos estadísiticos
library(RColorBrewer) # colores para gráficos
library(ggplot2) # Plot
library(ggpubr) # Plot

setwd("D:\\Javeriana\\Analisis de Datos") # Carpeta de trabajo

Reconozcamos la Herramienta

En R podemos realizar operaciones básicas. Iniciemos con una suma, aplicar logaritmo y exponencial.

3 + 4
## [1] 7
log(10)
## [1] 2.302585
exp(-0.5)
## [1] 0.6065307

pero yo puedo almacenar esta información en variables, por ejemplo

x <- 3 +  4  
x     # x es un vector cuya primera componente es 7. Enseguida vamos con los vectores!
## [1] 7

Así como creamos x, también podemos crear otras variables. Intentemos ahora crear y y z.

y = 2 + 6
y
## [1] 8
z <- c(x,y) # la función c() se utiliza para crear vectores
z
## [1] 7 8

Ahora, iniciemos con medidas de tendencia central (media, max y min). El comando base de Rstudio es mean(). Asignemos nuevas variables a estos resultados.

mean(z)
## [1] 7.5
a <- mean(z)
a
## [1] 7.5
b <- max(z)
b
## [1] 8

Por favor creen una variable c que contenga el valor min del vector z.

Muy bien, ahora creemos una distribución normal de 100 datos y redondeemos las cifras a 2

d <- rnorm(100)
round(d, 2)
##   [1]  0.26  0.38 -0.07 -1.50  0.69 -1.30 -0.08  0.77  1.19  0.32 -0.23  0.00
##  [13]  0.06  0.75  0.26  0.30  1.39  0.84  0.77  0.98  0.55 -0.45  0.39  0.46
##  [25]  0.59 -0.98 -0.39  2.22  1.04 -0.74 -0.16  1.28 -0.85  0.52  1.20 -1.28
##  [37] -0.22 -0.72 -0.57 -0.51 -0.36  0.51  0.03 -0.02  0.19 -0.96 -1.86  0.86
##  [49] -1.39  0.05 -0.94  0.42 -1.49 -0.08  1.39  0.51  0.51  0.66  0.93 -0.10
##  [61] -0.57 -0.73  0.39  0.05  0.34  0.50  0.91  0.47 -1.87  0.63  0.76  1.06
##  [73] -1.22 -0.94 -0.31 -1.92 -1.50 -1.07 -0.75 -1.28 -0.83 -0.02  0.56 -0.63
##  [85] -1.34  0.10 -0.59 -2.01 -0.59  0.33 -0.12  0.31 -0.07  0.08 -1.01  2.04
##  [97]  0.63  1.34  0.90  0.67
set.seed(10)
e <- round(rnorm(100),2)
args(round)
## function (x, digits = 0) 
## NULL

y si hacemos nuestra primera gráfica

plot(d)

Histograma

En estadística, un histograma es una representación gráfica de una variable en forma de barras, donde la superficie de cada barra es proporcional a la frecuencia de los valores representados. Sirven para obtener una “primera vista” general, o panorama, de la distribución de la población, o de la muestra, respecto a una característica, cuantitativa y continua (como la longitud o el peso).

De esta manera ofrece una visión de grupo permitiendo observar una preferencia, o tendencia, por parte de la muestra o población por ubicarse hacia una determinada región de valores dentro del espectro de valores posibles (sean infinitos o no) que pueda adquirir la característica. Así pues, podemos evidenciar comportamientos, observar el grado de homogeneidad, acuerdo o concisión entre los valores de todas las partes que componen la población o la muestra, o, en contraposición, poder observar el grado de variabilidad, y por ende, la dispersión de todos los valores que toman las partes, también es posible no evidenciar ninguna tendencia y obtener que cada miembro de la población toma por su lado y adquiere un valor de la característica aleatoriamente sin mostrar ninguna preferencia o tendencia.

hist(d, col='orange', breaks=5, 
     ylab = "Frecuencia", xlab="Datos Distr. Normal", main = "Histograma ejemplo")

Pueden decirme dónde esta la media? podrían hacer un mejor análisis si amplian el número de clases?

Hablemos de los Tipos de Datos

El primero de ellos es el Vector. el vector es elemento más básico en R.Contiene elementos de la misma clase (son atómicos). Se crea con la función c(), que significa ‘concatenar’ o ‘combinar’.

x <- c(1,2,3,4)    # con la función c() creamos el vector x
class(x)           # la función class() devuelve el tipo de objeto
## [1] "numeric"
## [1] "numeric"

También tenemos otros tipos de datos, como los que vemos a continuación

creemos variables con los diferentes tipos de datos que nos permitan ver el ejercicio un poco más claro.

y <- c("a","b")
class(y)
## [1] "character"
## [1] "character"

z <- c(1L,2L,3L)   # escribimos L detrás del número para obligar a que sea entero
class(z)
## [1] "integer"
## [1] "integer"

w <- c(TRUE, F)    # en general, puede escribirse TRUE/FALSE o T/F
class(w)
## [1] "logical"
## [1] "logical"

t <- c(1+2i, 1+3i)
class(t)
## [1] "complex"
## [1] "complex"

k <- factor(x)    # convierte el vector x en un factor
k
## [1] 1 2 3 4
## Levels: 1 2 3 4
## [1] 1 2 3 4
## Levels: 1 2 3 4
class(k)
## [1] "factor"
## [1] "factor"

Caracterización Inicial

Nota Importante: este tipo de herramientas (Rsutdio, Python) están pensados para escritura en inglés. En este idioma no hay un uso frecuente de tildes y caracteres especiales, así que para preservar las tildes y caracteres especiales se debe especificar el tipo de codificación como “UTF-8”.

Ahora si, entremos en materia. Para el desarrollo de sus actividades académicas seguramente será muy útil relacionar información ambiental con fechas. Dado esto, el interés de este ejercicio de clase se centrará en las fechas y los valores. Con estas dos columnas se crea una variable tipo “xts”, las cuales se utilizan para el manejo de series temporales.

La mayoría de fuentes oficiales de información ambiental maneja las fechas en el formato estandar AAAA-MM-DD, si se tuviera otro formato es necesario especificarlo; si se deseara utilizar las fechas junto con la hora de adquisición se debería utilizar variables tipo POSIXct.

datos <- read.csv("excel.csv.csv", encoding = "UTF-8") # Lee el archivo del DHIME
datos <- as.xts(datos$Valor, order.by = as.Date(datos$Fecha)) # Organiza los valores como serie temporal

Cuando trabajamos información hidroclimatológica (precipitación, caudal, temperatura, velocidad del viento, humedad realtiva, brillo solar, entre otras), se debe tener especial cuidado en la fuente, dado que el IDEAM, por ejemplo, no incluye celdas cuando en ese día no se registró información, por lo que el número de días total de la serie y el número de registros varían.

cat("Número de días con registro: ", length(datos))
## Número de días con registro:  10944
cat("Número de días totales: ", max(index(datos))-min(index(datos))+1)
## Número de días totales:  10958

Para esto solucionar esto se utiliza un vector xts vacío pero cargado con la totalidad de los días que hay en la serie. La función cbind cuando se utiliza en variables xts hace la combinación de columnas emparejando las fechas.

fechas <- seq(min(index(datos)), max(index(datos)), by = "day") # Fechas completas
empty.vector <- xts(order.by = fechas) # Serie temporal vacia
datos <- cbind(datos, empty.vector) # Serie con fechas completas

cat("\fNúmero de días con registro: ", length(datos))
## Número de días con registro:  10958
cat("\nNúmero de días totales: ", max(index(datos))-min(index(datos))+1)
## 
## Número de días totales:  10958

Posteriormente, se hace la caracterización estadística inicial. Para el presente ejercicio se seleccionaron las estadísticas nombradas en la variable “nombres”. Al concluir la eliminación de valores atípicos y si se realizara llenado de datos faltantes estas estadísticas deben ser calculadas de nuevo.

# Estadísticas a calcular
nombres <- c("Media", "Desv. Estandar", "Asimetria", "Curtosis", "Máximo", "Mínimo", "Mediana","Coef. Variacíón", "Número Total de Datos", "Número de Datos Faltantes", "Porcentaje de Faltantes")

# Calculo de las estadísticas
carc <- vector(length = length(nombres)) # Vector vacio para caracterización estadística
carc[1] <- mean(datos, na.rm = T)
carc[2] <- sd(datos, na.rm = T)
carc[3] <- skewness(datos, na.rm = T)
carc[4] <- kurtosis(datos, na.rm = T)
carc[5] <- max(datos, na.rm = T)
carc[6] <- min(datos, na.rm = T)
carc[7] <- median(datos, na.rm = T)
carc[8] <- carc[2] / carc[1]
carc[9] <- sum(!is.na(datos))
carc[10] <- sum(is.na(datos))
carc[11] <- carc[10] / length(datos)

carc <- data.frame("Estadística" = nombres, Valor = carc)
carc
##                  Estadística        Valor
## 1                      Media 5.454130e+00
## 2             Desv. Estandar 1.097241e+01
## 3                  Asimetria 3.510595e+00
## 4                   Curtosis 2.024767e+01
## 5                     Máximo 1.320000e+02
## 6                     Mínimo 0.000000e+00
## 7                    Mediana 0.000000e+00
## 8            Coef. Variacíón 2.011762e+00
## 9      Número Total de Datos 1.094400e+04
## 10 Número de Datos Faltantes 1.400000e+01
## 11   Porcentaje de Faltantes 1.277605e-03

Dada la gran cantidad de datos es recomendable ayudarse continuamente de gráficas.

En este punto se hará una gráfica, donde se muestre la serie completa de datos de Precipitación en el tiempo. Para las siguientes gráficas se va a emplear la librería ggplot, para esto primero se debe organizar la información en un data frame, el cual se introduce como primera variable y posteriormente se especifican los nombres de las columnas que serán graficadas. Mediante el signo “+” se indica el tipo de geometría que se utilizará (“geom_…”) y se le da formato al gráfico.

# Plot de serie temporal mensual multianual
datos.df <- data.frame(t = index(datos), datos)
p <- ggplot(datos.df, aes(t, datos)) + geom_line() +
  labs(title = "Precipitación en la Estacion Nataga", x = "Fechas", y = "mm/h") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), text = element_text(size = 14))
p

A estos gráficos se les puede cambiar el color y el tipo de línea. Probemos con líneas azules!.

# Plot de serie temporal mensual multianual
datos.df <- data.frame(t = index(datos), datos)
p <- ggplot(datos.df, aes(t, datos)) + geom_line(colour = "blue") +
  labs(title = "Precipitación en la Estacion Nataga", x = "Fechas", y = "mm/h") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), text = element_text(size = 14))
p

y si ahora probamos con el grosor y tipo de la línea… por favor busquen en internet el comando que les pueda ayudar en esta labor.

# Plot de serie temporal mensual multianual
datos.df <- data.frame(t = index(datos), datos)
p <- ggplot(datos.df, aes(t, datos)) + geom_line(colour = "blue",lwd = 2, linetype = 2) +
  labs(title = "Precipitación en la Estacion Nataga", x = "Fechas", y = "mm/h") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), text = element_text(size = 14))
p

Bien ahora hablemos un poco más del análisis estadístico mediante el uso de boxplot. Recordemos…los gráficos de caja (boxplot) fueron originalmente desarrollados por M.E. Spear, permiten conocer cómo se distribuyen los datos dentro de una variable. A diferencia de los histogramas que requieren un tamaño de muestra de al menos 30 casos para ser útiles, los gráficos de caja pueden ser construidos con tan solo 5 casos y aportan más detalles acerca de las colas de la distribución. Representan la información que se observa en la figura:

Mediana. Valor que deja a la mitad de los casos por encima y a la otra mitad por debajo.

Primer Cuartil (Q1). El 25% de los casos se encuentran por debajo de este valor.

Tercer Cuartil (Q3). El 75% de los casos se encuentran por encima de este valor.

Rango Intercuartílico (RIC). Es la diferencia entre el tercer y el primer cuartil.

Límites Superior o Inferior (Ls o Li). Ls contiene los casos por encima de Q3 más 1,5 veces el rango intercuartílico o Li por debajo de Q1 – 1,5xRIC (Estilo de Tukey). Cuando los valores no son posibles en lugar de emplear la aproximación anterior se escogen los valores máximo o mínimo de la muestra (Estilo de Spears).

Los valores atípicos son aquellos que están más a allá de los límites inferior y superior. Cuando los valores atípicos están más allá de 3 veces el RIC en lugar del 1.5 son denominados valores extremos.

Para generar los gráficos de boxplot en una escala mensual se debe generar una columna que especifique el mes, para lo que se usa la función month y se especifica que se agrupara con esta columna.

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot() + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

Juguemos con los gráficos de boxplot…

Ejemplo 1. Resaltar la media (mean) dentro del grafico usando la función stat_summary:

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot() + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))+
  stat_summary(fun=mean, geom="point", color="red")
bp

Ejemplo 2. cambiar el color del borde y resaltar los datos atípicos:

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(color=4, outlier.colour="red") + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

y si queremos omitir los datos atípicos… cómo se les ocurre?

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(color=4, outlier.colour=NA) + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp

Ejemplo 3. Boxplot notches

# Plot boxplot mensual multianual
datos.df <- data.frame(t = index(datos), datos)
datos.df[,1] <- month(index(datos)) # Mes de adquisición del dato
bp <- ggplot(datos.df, aes(t, datos, group = t)) +
  geom_boxplot(notch=TRUE, notchwidth = 0.5,
               color="blue", alpha = 0.7,
               outlier.colour="red", outlier.size=0.7) + labs(x = "Meses", y = "mm/h") +
  scale_x_continuous(breaks = c(1:12)) +
  theme(text = element_text(size = 14))
bp
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.
## notch went outside hinges. Try setting notch=FALSE.

Qué otras formas de boxplot encuentran?

Para agrupar mensual o anualmente se utiliza la función aggregate, en donde se introduce un indice bajo el cual se desea que se agrupen los datos, se especifica la manera de agruparlos y cualquier elemento extra que se considere. Para obtener este indice a partir de las fechas diarias se “corta” la fecha con la función cut y se especifica a que resolución temporal se desear cortar. Este proceso cambia el tipo de variable, así que se debe volver a definirla como una serie temporal (xts).

Cuando la serie tiene numerosos datos faltantes se debe especificar una función que descarte un mes o año cuando a este le falta cierto porcentaje de datos.

## Agrupación
mensual <- aggregate(datos, by = cut(index(datos), "month"), sum, na.rm = T)
mensual <- as.xts(mensual, order.by = as.Date(index(mensual)))
head(mensual)
##            [,1]
## 1992-01-01  130
## 1992-02-01  129
## 1992-03-01   56
## 1992-04-01  108
## 1992-05-01  252
## 1992-06-01   34
anual <- aggregate(datos, by = cut(index(datos), "year"), sum, na.rm = T)
anual <- as.xts(anual, order.by = as.Date(index(anual)))
head(anual)
##            [,1]
## 1992-01-01 1517
## 1993-01-01 2236
## 1994-01-01 2636
## 1995-01-01 1787
## 1996-01-01 2508
## 1997-01-01 1971

En precipitación se debe tener en cuenta que para pasar de diarios a mensuales o anuales se suman pero para pasar a valores multianuales se promedia. Este tratamiento difiere de variable a variable.

mensual.ma <- aggregate(mensual, by = month(index(mensual)), mean, na.rm = T)
mensual.ma
##             
## 1  174.08333
## 2  186.98000
## 3  234.84000
## 4  217.41667
## 5  178.68333
## 6  107.44333
## 7   82.52000
## 8   56.96000
## 9   79.04333
## 10 185.04000
## 11 266.19000
## 12 220.46667
anual.ma <- mean(anual, na.rm = T)
anual.ma
## [1] 1989.667

Finalmente, se puede graficar el comportamiento anual promedio de la zona. Incluso se pueden mostrar varias variables en una sola gráfica para complementar la información.

temp <- c(30.0, 28.6, 27.9, 28.1, 29.7, 30.2, 31.0, 31.5, 30.7, 28.8, 27.9, 28.6)
mensual.ma <- data.frame(m = 1:12, m.abb = month.abb, ppt = mensual.ma, temp = temp)
p0 <- ggplot(mensual.ma, aes(m, ppt)) + geom_col(fill = "darkslategray3") +
  scale_x_continuous(breaks = 1:12, labels = mensual.ma$m.abb) +
  labs(title = "Precipitación Mensual Multianual en Nátaga", subtitle = sprintf("Precipitación total = %.f mm", anual.ma), x = "Meses", y = "mm") +
  theme(plot.title = element_text(hjust = 0.5))
p0

coef <- 0.2
p1 <- ggplot(mensual.ma, aes(m)) +
  geom_col(aes(y = ppt), fill = "darkslategray3") +
  geom_line(aes(y = temp/coef), color = "firebrick3") +
  scale_x_continuous(breaks = 1:12, labels = mensual.ma$m.abb) +
  scale_y_continuous(name = "Precipitación (mm)", sec.axis = sec_axis(~.*coef, name = "Temperatura (ºC)")) +
  labs(title = "Caracterización Mensual Multianual en Nátaga", subtitle = sprintf("Precipitación total = %.f mm \nTemperatura Media = %.1f ºC", anual.ma, mean(temp)), x = "Meses") +
  theme(plot.title = element_text(hjust = 0.5))
p1