jueves, 29 de septiembre de 2011

Diagramas de cajas (boxplot), cuantiles, IQR, valores atípicos, y más. Aplicación en R

Diagrama de caja o Caja con bigotes
Es un gráfico, basado en cuartiles, mediante el cual se visualiza un conjunto de datos. Está compuesto por un rectángulo, la "caja", y dos brazos, los "bigotes".
Se dibujará una caja entre las bisagras (hinges) superior e inferior con una línea en la mediana. Bigotes (whisker, una sola línea, no un cuadro) se extienden desde las bisagras.

Si preguntamos a R por los detalles:
?boxplot.stats


The two ‘hinges’ are versions of the first and third quartile, i.e., close to quantile(x, c(1,3)/4). The hinges equal the quartiles for odd n (where n <- code="" length="" x="">) and differ for even n. Whereas the quartiles only equal observations for n %% 4 == 1 (n = 1 mod 4), the hinges do so additionally for n %% 4 == 2 (n = 2 mod 4), and are in the middle of two observations otherwise.
The notches (if requested) extend to +/-1.58 IQR/sqrt(n). This seems to be based on the same calculations as the formula with 1.57 in Chambers et al. (1983, p. 62), given in McGill et al. (1978, p. 16). They are based on asymptotic normality of the median and roughly equal sample sizes for the two medians being compared, and are said to be rather insensitive to the underlying distributions of the samples. The idea appears to be to give roughly a 95% confidence interval for the difference in two medians.


A la izquierda observamos el gráfico que lanza R al pedirle:
boxplot(alltime.movies$Gross)

Mientras que a la derecha he generado un gráfico que explica cada estadístico utilizado. El código para este gráfico se observa debajo.



##############################################################
####### boxplot2 grafica un diagrama de cajas ################
####### identificando los estadísticos que intervienen #######
##############################################################
 
    install.packages("UsingR") 
    library(UsingR)
#creamos una nueva función        
    boxplot2<-function(x=Gross)
    {
      #calculo estadísticos
      stats=boxplot.stats(x)$stats #extreme of the lower whisker, the lower ‘hinge’, the median, the upper ‘hinge’ and the extreme of the upper whisker.
      iqr=stats[4]-stats[1] #rango intercuartílico
      f=fivenum(x) #minimum, lower-hinge, median, upper-hinge, maximum
      stats2<-c(f[1],stats,f[5]) #agrego los min y max del fivenum
 
      #grafico
      boxplot(x) #hago el gráfico de cajas
      abline(h=stats[1],lty=2,col="red") #agrego una línea roja en el lower whisker
      abline(h=stats[5],lty=2,col="red") #agrego una línea roja en el upper whisker
      text(rep(1.35,7), stats2, labels=c('minimum','lower whisker', 'lower hinge', 'median', 'upper hinge', 'upper whisker','maximum'), cex=0.6) ##agrego labels
      text(rep(.7,7), stats2, labels=round(stats2), cex=0.6) ##agrego los valores
      #para estos datos el mínimo coincide con el lower whisker
 
      text(1.35,f[2]-1.5*iqr , labels="lower notche", cex=0.6) ##agrego labels
      text(.7, f[2]-1.5*iqr, labels=round(f[2]-1.5*iqr), cex=0.6) ##agrego los valores
      abline(h=f[2]-1.5*iqr, lty=3,col="blue")
 
      text(1.35,f[4]+1.5*iqr , labels="upper notche", cex=0.6) ##agrego labels
      text(.7, f[4]+1.5*iqr, labels=round(f[4]+1.5*iqr), cex=0.6) ##agrego los valores
      abline(h=f[4]+1.5*iqr,lty=3,col="blue")
      #lower notche no sale en este gráfico porque es menor al mínimo   
    }
 
#ahora para cualquier dato podemos graficar el boxplot con la descripción
data(alltime.movies)
boxplot2(x=alltime.movies$Gross)
 
 
#################################################
#### EXPLICACIÓN ####
#################################################
# preguntar cómo funciona boxplot.stats y fivenum
 
boxplot.stats
      # function (x, coef = 1.5, do.conf = TRUE, do.out = TRUE) 
      # {
      #   if (coef < 0){stop("'coef' must not be negative")}
      #   nna <- !is.na(x)
      #   n <- sum(nna)
      #   stats <- stats::fivenum(x, na.rm = TRUE)
      #   iqr <- diff(stats[c(2, 4)])
      #   if (coef == 0){do.out <- FALSE}
      #   else {
            #     out <- if (!is.na(iqr)) {
            #       x < (stats[2L] - coef * iqr) | x > (stats[4L] + coef * iqr) #por fuera del notch
            #     }
            #     else !is.finite(x)
            #     if (any(out[nna], na.rm = TRUE)) 
            #       stats[c(1, 5)] <- range(x[!out], na.rm = TRUE)
      #   }
      #   conf <- if (do.conf) 
      #     stats[3L] + c(-1.58, 1.58) * iqr/sqrt(n)
      #   list(stats = stats, n = n, conf = conf, out = if (do.out) x[out & 
      #                                                                 nna] else numeric())
      # }
 
 
fivenum
      # function (x, na.rm = TRUE) 
      # {
      #   xna <- is.na(x)
      #   if (na.rm) 
      #     x <- x[!xna]
      #   else if (any(xna)) 
      #     return(rep.int(NA, 5))
      #   x <- sort(x)
      #   n <- length(x)
      #   if (n == 0) 
      #     rep.int(NA, 5)
      #   else {
      #     n4 <- floor((n + 3)/2)/2
      #     d <- c(1, n4, (n + 1)/2, n + 1 - n4, n)
      #     0.5 * (x[floor(d)] + x[ceiling(d)])
      #   }
      # }
 
 
 
########################################  
### PARA UNOS DATOS CUALQUIERA Y=1:10
########## ############################## 
 
### ¿Qué son los hinges (bisagras)?:
    #lower y upper hinges son versiones de los cuartiles 1 y 3, respectivamente.
    #cuando la muestra (o número de datos n=length(x)) es impar, son los cuartiles
    #cuando la muestra es par, difiere de los cuartiles. 
    #ej
    y=1:9 #si la muestra es impar, no hay problema, los cuartiles son números enteros
    quantile(y,c(.25,.75))  # 3 y 7
    y=1:10 #pero si la muestra es par, los cuartiles son números reales (salen con decimales), cuando en realidad "y" no tiene decimales
    quantile(y,c(.25,.75))  #3.25 y 7.75
 
### ¿Cómo se calculan los hinges?:
    #para solucionarlo hay 3 formas:
    # 0) la forma directa
    fivenum(y)
    fivenum(y)[2] # 3; lower hinges 
    fivenum(y)[4] # 8; upper hinges 
    # 1) redondeamos los datos
    round(quantile(y,c(.25,.75))) #3 y 8
    # 2) aplicamos la función floor((n + 3)/2)/2 , donde floor redondea el valor
    ny=length(y)
    n4 <- floor((ny + 3)/2)/2 
    dy <- c(1, n4, (n + 1)/2, n + 1 - n4, n)
    fdy<-0.5 * (y[floor(dy)] + y[ceiling(dy)]) #min, lower hinges, median, upper hinges, max
    fdy[2] # 3; lower hinges o fivenum(y)[2]
    fdy[4] # 8; upper hinges o fivenum(y)[4]
 
 
### ¿Qué son los whiskers (bigotes)?
    #upper whisker: máximo de los valores de x cuando eliminamos aquellos que están por encima de upper hinges+1.5*IQR (hinges+notches)  
    #lower whisker: mínimo de los valores de x cuando eliminamos aquellos que están por encima de lower hinges-1.5*IQR (hinges-notches)
    #donde, notches (muescas): son los valores hinges+-1.5*iqr/sqrt(n)
 
### ¿Cómo se calculan los whiskers (bigotes)?
    #0) forma rápida:
    boxplot.stats(y)$stats[c(1,5)]
    #1) paso a paso. primero calculamos los notches 
    IQR=fivenum(y)[4]-fivenum(y)[2]  #rango intercuartílico, diferencia entre el tercer y primer cuartil
    fivenum(y)[4]+1.5*IQR #notche superior; el /sqrt(n) lo usa solo para la salida de do.conf=T
    fivenum(y)[2]-1.5*IQR #notche inferior; el /sqrt(n) lo usa solo para la salida de do.conf=T
    #luego calculamos los whiskers
    outy <- if (!is.na(IQR)) {y < (fivenum(y)[2] - coef * IQR) | y > (fivenum(y)[4] + coef * IQR)} #por fuera de los notches 
    range(y[!outy], na.rm = TRUE) #whiskers
    #en estos datos de prueba no hay datos extremos o outliers (y=1:10)
    #outliers:  the values of any data points which lie beyond the extremes of the whiskers.
Created by Pretty R at inside-R.org




Ventajas
* Proporcionan una visión general de la simetría de la distribución de los datos; si la mediana no está en el centro del rectángulo, la distribución no es simétrica.
* Son útiles para ver la presencia de valores atípicos también llamados outliers.

Veamos entonces las definiciones de cada elemento:
  • Cuantiles (quantil): son colectivamente los cuartiles (Qi, i=1,...,3 divide los datos en 4 regiones), deciles (Di, i=1,...,9 divide los datos en 10 regiones) y percentiles (Pi, i=1,...,99 divide los datos en 100 regiones).
  • Rango intercuartílico (IQR): es la diferencia entre el tercer y el primer cuartil. IQR=Q3-Q1.
  • Valores atípicos (outliers): son observaciones numéricamente distantes del resto de los datos. Tomando como referencia el IQR, en un diagrama de cajas se considera un valor atípico el que se encuentra 1,5 veces esa distancia de uno de los cuartiles (atípico leve) o a 3 veces esa distancia (atípico extermo).
  1. atípicos leves: son aquellos que se encuentran en x < Q1-1.5*IQR o x > Q3 + 1.5*IQR
  2. atípicos extremos: son aquellos que se encuentran en x < Q1-3*IQR o x > Q3+3*IQR
  • Bisagras (hinges): son los mismos que los cuartiles a menos que al dividir el tamaño de la muestra por 4 el resto sea 3 (como el 39 / 4 = 9 I 3).
  1. bisagra inferior: es la mediana de la mitad inferior de los datos hasta e incluyendo la mediana.
  2. bisagra superior: es la mediana de la mitad superior de los datos hasta e incluyendo la mediana.
  • Resumen de 5-números (five number summary): consiste en el valor mínimo, la bisagra inferior, la mediana, la bisagra superior y el valor máximo. (aunque algunos libros utilizan los cuartiles en lugar de las bisagras).

Para observar cómo realizar estos gráficos en R, basta con pedirle ejemplos de la función "boxplot".



No hay comentarios:

Publicar un comentario

Libros para descargar (gratis) sobre Estadística Computacional