En el último artículo hemos visto AWK como herramienta para el tratamiento masivo de datos, donde hemos demostrado que podía ser una alternativa parcial a otras herramientas comerciales como SAS, ACL o IDEA (muy utilizadas en el mundo de la Auditoría). AWK es una solución parcial dado que únicamente nos proporciona funcionalidades para hacer un tratamiento básico sobre los ficheros (p.ej. sumarizaciones, cruces de datos, operaciones matemáticas simples, etc), sin embargo no tenemos las características estadísticas que si podemos encontrar en las herramientas comerciales.
Afortunadamente para todas esas funcionalidades podemos apoyarnos en R: R es un lenguaje de programación que permite llevar a cabo análisis estadísticos avanzados, entendiendo por por estadística como la ciencia de recoger y analizar datos con el propósito de sacar conclusiones y tomar decisiones.
Con R también podemos realizar cálculos numéricos, aunque para esas tareas también vale la pena echar un vistazo a Octave.
Aparte de las potentes características de R, existen multitud de paquetes que amplían las funcionalidades: desde dedicados al análisis de datos psicológicos hasta financieros.
Al igual que AWK, R no ofrece una interfaz gráfica tan potente como las herramientas comerciales que comentavamos, aunque podemos apoyarnos en determinados mecanismos que nos harán la vida más fácil a la hora de tratar la información (p.ej. utilizando MySQL Query Browser) como veremos en las correspondientes secciones del artículo.
Para aprender a utilizar R mediante esta guía es muy recomendable replicar todos los ejemplos y visualizar los resultados directamente. Únicamente leyendo el artículo es más difícil entender el funcionamiento completo de la herramienta.
De antemano pido disculpas si soy inexacto o cometo algún tipo de error en las explicaciones estadísticas, no soy un experto en esa materia y estaré muy agradecido si detectáis incorrecciones y las hacéis llegar vía comentario o correo.
Finalmente, destacar que para la elaboración del artículo me he basado en diversos tutoriales que he encontrado por Internet, especialmente simpleR, con los cuales he ido aprendido a utilizar R de forma paralela a la redacción de esta guía.
Índice de contenidos
Estructuras complejas de datos
Distribuciones de probabilidad
Generación de números aleatorios
Generación de números secuenciales
Visualización de datos de una variable
Visualización de datos de dos variables
Identificar puntos en un gráfico
Identificar la distribución de los datos
Diferentes formas de entender la probabilidad
Probabilidades bayesienas o basadas en evidencias
T-Test: Datos con distribución normal
Una variable
Dos variables 26
Simulaciones sobre el T-Test
Wilcoxon Test: Datos con distribución diferente a la normal
ANOVA / Kruskal
Coeficiente de Pearson: Datos con distribución normal
Métodos no paramétricos: Datos con distribución diferente a la normal
Regresiones lineales múltiples
Técnicas para identificar fraude
Instalación
La instalación de R en Ubuntu GNU/Linux no puede ser más sencilla:
sudo apt-get install r-base-core
Para plataformas Windows, debemos descargar el instalador de la página del projecto R: http://www.r-project.org/
Introducción
Para interactuar con R se dispone de una potente linea de comandos:
$ R R version 2.6.2 (2008-02-08) Copyright (C) 2008 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R es un software libre y viene sin GARANTIA ALGUNA. Usted puede redistribuirlo bajo ciertas circunstancias. Escriba 'license()' o 'licence()' para detalles de distribucion. R es un proyecto colaborativo con muchos contribuyentes. Escriba 'contributors()' para obtener más información y 'citation()' para saber cómo citar R o paquetes de R en publicaciones. Escriba 'demo()' para demostraciones, 'help()' para el sistema on-line de ayuda, o 'help.start()' para abrir el sistema de ayuda HTML con su navegador. Escriba 'q()' para salir de R. >
Para salir de la consola basta con ejecutar "q()", entonces R nos permite grabar la sesión actual (historial de comandos y datos cargados).
Con todos los comandos que veremos podemos utilizar la función ‘help’ para obtener los detalles de su funcionamiento, por ejemplo: help(read.csv) o ?help
En caso de que no recordemos un comando exacto podemos utilizar "apropos(parte_del_comando)" para que el sistema liste posibles candidatos:
> apropos("rank")
[1] "dsignrank" "psignrank" "qsignrank" "rank" "rsignrank"
Estructuras simples de datos
R permite utilizar variables con tipado dinámico, es decir, no tenemos que declarar las variables y explicitar que tipo de valores van a contener. Los nombres de las variables puede ser cualquier cadena que empiece por un carácter. Veremos habitualmente el uso de puntos en lugar de subrayado/underscore para separar palabras y que las variables sean más autoexplicativas, por ejemplo transferencias.nacionales. Por tanto, si vemos esa nomenclatura no tiene porque tratarse de un método/atributo de un objeto (como pasa en otros lenguajes), sino que puede ser perfectamente una simple variable. Para declarar y asignar un valor a una variable se puede utilizar tanto ‘=’ como ‘<-‘:
contador = 0 contador < - 1
Una variable puede contener una estructura de datos simple o compleja.
Números
R permite el uso de numéricos enteros o decimales, estos últimos separados por puntos:
x = 1 y = 2.5
Sobre los valores numéricos se pueden realizar las habituales operaciones aritméticas, entre las que se encuentran: + * / – ^
z = (x + y ^ 2) - .5
Los números decimales entre 0 y 1 podemos especificarlos sin indicar el 0 (p.ej. 0.5 => .5).
Booleanos
Una variable booleana únicamente puede contener 2 valores: TRUE o FALSE (abreviadas también como T o F)
x = TRUE y = FALSE
Si utilizamos expresiones condicionales, su resultado será un valor booleano (p.ej. 17 %in% 1:100). Los operadores disponibles son los siguientes: < <= > >= == !=
Por otra parte, podemos combinar varios valores booleanos o expresiones condicionales mediante el uso de negaciones, “ands” y “or”: ! & |
Cadenas
Para el tratamiento de variable que contienen cadenas o de las cadenas directamente nos pueden ser de utilidad las siguientes funciones:
# Impresión
cat("Cadena!\n")
# Concatenación separada por nada
paste("Cadena", "!", sep="")
> s = c("Hello", " ", "World", "!")
> paste(s, collapse="")
[1] "Hello World!"
# Longitud de una cadena
> nchar("Hello World!")
[1] 12
# Substring de la posición 4 a la 6
> s = "Hello World"
> substring(s, 4, 6)
[1] "lo "
# Separa una cadena en campos, devuelve una lista
> s = "foo, bar, baz"
> x = strsplit(s, ", ")
> x[[1]] # Vector
[1] "foo" "bar" "baz"
> x[[1]][1] # Elemento 1
[1] "foo"
# Busca elementos que cumplan una expresión regular
> s = c("hola", "adios")
> grep("^a", s)
[1] 2
> grep("^a", s, value=T)
[1] "adios"
# O elementos que cumpla aproximadamente una expresión regular
> agrep ("abc", c("abbc", "jdfja", "cba"))
[1] 1
# Buscamos una expresión regular en una cadena
> x = regexpr("o", "Hello")
> x[1] # Posición donde empieza la cadena encontrada
> attr(x, "match.length") # Longitud de la cadena encontrada
# Reemplaza todos los ; por ,
> gsub(";", ",", "dato1;dato2;dato3")
[1] "dato1,dato2,dato3"
# Reemplaza solo el primer ; por ,
> sub(";", ",", "dato1;dato2;dato3")
[1] "dato1,dato2;dato3"
Fechas
R proporciona un tipo de variable especial para las fechas, para las cuales nos puede resultar de utilidad las siguientes funciones:
# Convertir cadena a fecha
> a = as.Date("01/02/03", format="%y/%m/%d")
# Diferencia de fechas en dias
> b = as.Date("01/02/03", format="%y/%d/%m")
> a - b
Time difference of -27 days
# Fecha actual + 15 dias
> Sys.Date() + 15
[1] "2009-03-08"
# Convertir fecha a cadena con un formato determinado
> format(Sys.Date(), format="%A, %d %B %Y")
[1] "sábado, 21 febrero 2009"
# Extraer el mes
> format( Sys.Date(), format="%m" )
[1] "02"
# Secuencia de fechas desde Enero 2005 a Julio 2005 mes a mes
> seq(as.Date("2005-01-01"), as.Date("2005-07-01"), by="month")
[1] "2005-01-01" "2005-02-01" "2005-03-01" "2005-04-01" "2005-05-01"
[6] "2005-06-01" "2005-07-01"
# Recorrer la secuencia
> dates = seq(as.Date("2005-01-01"), as.Date("2005-07-01"), by="month")
> for (i in 1:length(dates)) {
d = as.Date(as.character(dates[i])) # Hay que convertir o solo visualizaremos un numero
cat(paste(d, "\n", sep=""))
}
Estructuras complejas de datos
En R disponemos de las siguientes estructuras complejas de datos fundamentales:
- Vectores: 1
dimensión - Matrices/Arrays: 2
o más dimensiones - Data Frames: 2
dimensiones donde cada columna puede ser de un tipo diferente
(cadenas, números, fechas y booleanos) - Listas: pueden
contener cualquiera de las estructuras anteriores .
Para que R pueda trabajar con los datos, estos deben ser cargados en memoria. Para listar lo que actualmente existe en memoria utilizamos "ls()":
> datos = "test" > ls() [1] "datos"
Cuando queramos dejar de trabajar con un conjunto de datos y no los necesitemos en memoria, podemos eliminarlos mediante "rm(datos)".
Por otra parte, en todas las estructuras podemos utilizar las funciones:
> x = 5 > str(x) # Explica la estructura de x > summary(x) # Sumariza la información contenida en x
Vectores
Veamos las funciones relacionadas con los vectores:
- Creación
x = seq(-1, 1, by=.1) # Creamos un vector con números desde -1 al 1 sumando 0.1 names(x) = letters[1:length(x)] # Ponemos nombre a las coordenadas x["a"], x[1] x = c(1, 2, 3) # Creamos un vector con los datos indicados x = c(a=1, b=2, c=3) # Datos y nombres simultáneamente > rep(1,10) [1] 1 1 1 1 1 1 1 1 1 1 > rep(1:5,3) [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 > rep(1:5,each=3) # Repite del 1 al 5 cada elemento 3 veces [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5
- Filtros
> x["a"] # Elemento "a" > x[1] # Elemento 1 > x[5:10] # Elementos 5 al 10 [1] -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 > x[c(5,7:10)] # Elemento 5 y del 7 al 10 [1] -0.6 -0.4 -0.3 -0.2 -0.1 > x[-(5:10)] # Todos menos los elementos del 5 al 10 [1] -1.0 -0.9 -0.8 -0.7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 > x[ x>0 ] # Los valores superiores a 0 [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
- Operaciones
> x>0 # Devuelve vector con el resultado de aplicar la condición [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > x = rnorm(10) > sort(x) # Orden de menor a mayor [1] -1.4159893 -1.1159279 -1.0598020 -0.2314716 0.3117607 0.5376470 [7] 0.6922798 0.9316789 0.9761509 1.1022298 > rev(sort(x)) # Orden inverso [1] 1.1022298 0.9761509 0.9316789 0.6922798 0.5376470 0.3117607 [7] -0.2314716 -1.0598020 -1.1159279 -1.4159893 > o = order(x) # Orden de menor a mayor, pero devuelve el indice de los elementos > o [1] 3 1 9 6 4 7 8 10 2 5 > x[ o[1:3] ] # Asi podemos acceder a los valores [1] -1.415989 -1.115928 -1.059802 > x = sample(1:5, 10, replace=T) > unique(x) # Elimina duplicados (no necesita que el vector sea ordenado previamente) [1] 1 4 5 3 # Busquemos que elementos son igual a 0 y en que posición se encuentran: > x.test = c(1, 2, 0, 3, 2, 0, 0, 1) > x.test == 0 [1] FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE > which(x.test == 0) [1] 3 6 7 # Para concatenar dos vectores: > c(c(1,2,3), c(4,5,6)) [1] 1 2 3 4 5 6 # O sumar sus elementos: > c(c(1,2,3) + c(4,5,6)) [1] 5 7 9
- Cálculos
> min(x) # Mínimo > max(x) # Máximo > cummax(c(5,5,4,7,7,3)) # Máximo visto al llegar a cada elemento [1] 5 5 5 7 7 7 > cummin(c(5,5,4,7,7,3)) # Mínimo visto al llegar a cada elemento [1] 5 5 4 4 4 3 > sum(x) # Suma > length(x) # Longitud > round(x) # Redondeo de cada elemento > round(x, digits=2) # Redondeo de cada elemento a 2 decimales > rank(x) # Devuelve la posición de cada elemento si estuviese ordenado el vector de menor a mayor (p.ej. 9 15 7 => segundo tercero primero => 2 3 1)Rango > median(x) # Mediana: valor de la variable que deja el mismo número de datos antes y después que él, una vez ordenados estos > mean(x) # Media aritmetica: suma de todos dividida entre el número de elementos # Esperanza: # - Suma del producto de la probabilidad de cada suceso por el valor de dicho suceso. # - Por ejemplo, en un juego de azar el valor esperado es el beneficio medio. # - Si todos los sucesos son de igual probabilidad la esperanza es la *media aritmética*. > var(x) # Varianza: dispersión de una variable aleatoria respecto a su *esperanza* > sd(x) # Desviación estándar: media de las distancias que tienen los datos respecto de su media aritmética (raiz cuadrada de la varianza = sqrt(var(x))) > mad(x) # Desviación media normalizada respecto a la mediana (más resistente a outliers/valores extremos que la desviación estándar) > cor(x, x) # Correlación por coeficiente de Pearson: fuerza y dirección de una relación lineal entre dos variables aleatorias (obviamente la correlación del ejemplo es total y en la misma dirección dado que son los mismos números: 1)
- Cuartiles: en diversas ocasiones haremos mención a los cuartiles de un conjunto de datos, estos corresponden a los siguientes elementos:
- Q1: Los valores inferiores a Q1 representan el 25% de la población
- Q2: Igual a la media
- Q3: Los valores inferiores a Q3 representan el 75% de la población
> quantile(x) > summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. -2.70600 -0.56750 0.10400 0.08912 0.72590 2.65400 > IQR(x) # Diferencia entre el primer y tercer cuantil (Q1 y Q3). Más resistente a outliers/valores extremos que la media.
- Factores: vector con datos cualitativos
> x = factor( sample(c("Yes", "No", "Perhaps"), 20, replace=T) )
> table(x) # Sumarizamos (numero de ocurrrencias de cada factor)
x
No Perhaps Yes
5 8 7
> levels(x) # Listamos los valores posibles
[1] "No" "Perhaps" "Yes"
- Valores no disponibles (No Available)
> x = c(1,5,9,NA,2) > na.omit(x) # Limpia los valores NA > is.na(x) # Comprueba si es NA (TRUE/FALSE)
- Infinito
> log(0) [1] -Inf > is.infinite(log(0)) [1] TRUE
Matrices
En una matriz todos los datos deben ser del mismo tipo, es decir: cadenas, números o booleanos.
- Creación
> m = matrix( c(1,2,3,4), nrow=2 ) # Creamos una matriz numérica > m [,1] [,2] [1,] 1 3 [2,] 2 4 > m = cbind( c(1,2), c(3,4) ) # Identico a la anterior > m = rbind( c(1,2), c(3,4) ) # Idem
- Operaciones
# Producto (Fila1: 1*2 + 3*1, Fila2: 2*2 + 4*1 )
> m %*% matrix( c(2,1), nrow=2 )
[,1]
[1,] 5
[2,] 8
> det(m) # Determinante
> t(m) # Transposición
[,1] [,2]
[1,] 1 2
[2,] 3 4
> solve(m) # Inversa
[,1] [,2]
[1,] -2 1.5
[2,] 1 -0.5
> diag(c(1,2)) # Diagonal
[,1] [,2]
[1,] 1 0
[2,] 0 2
> diag(rep(1,2)) # Matriz identidad
[,1] [,2]
[1,] 1 0
[2,] 0 1
Data Frames
Hasta el momento hemos visto como trabajar con un vector o un par de vectores relacionados, pero también es habitual tener que hacer tratamiento sobre tablas que contienen más de una columna (más de un vector en definitiva). Para facilitar esta tarea, R dispone de los denominados data frames.
En un Data Frame cada columna puede ser de un tipo diferente: cadenas, números, fechas o booleanos. Por otra parte, esta estructura de datos permite realizar toda una serie de operaciones de gran utilidad como cruces entre 2 dataframes o sumarizaciones.
- Creación
> peso = c(75, 98, 67, 80)
> altura = c(1.75, 1.89, 1.61, 1.77)
> sexo = c("F","F","M","F")
> personas = data.frame(peso,altura,sexo) # Construye un dataframe con los 3 vectores y establece como nombre de columna el de las variables
> personas
peso altura sexo
1 75 1.75 F
2 98 1.89 F
3 67 1.61 M
4 80 1.77 F
> rm(peso) ; rm(altura); rm(sexo)
> df = data.frame(c(1, 2, 3), c(7, 3, 4)) # Sin nombres de columna
> df = data.frame(col1 = c(1, 2, 3), col2 = c(7, 3, 4)) # Con nombres de columna
- Acceso a elementos
# Muestra la columna 1
> df$col1
> df[,1]
> df[["col1"]]
> personas[2] # Columna 2
altura
Maria 1.75
Alicia 1.89
Pedro 1.61
Judith 1.77
> personas[,2] # Fila 2
[1] 1.75 1.89 1.61 1.77
> personas["Maria", "altura"]
[1] 1.75
> personas$altura
[1] 1.75 1.89 1.61 1.77
- Filas y columnas
> dim(df) # Numero de columnas y filas
> names(df) # Nombres de columnas
> row.names(df) # Nombres de filas
# Por defecto, las filas se enumeran del 1..n pero también podríamos asignar nombres
> row.names(personas)=c("Maria","Alicia","Pedro","Judith")
> personas
peso
altura sexo
Maria
75 1.75 F
Alicia
98 1.89 F
Pedro
67 1.61 M
Judith
80 1.77 F
# Acceso a las columnas como variables
> attach(df)
> col1
> dettach(df)
- Filtros
> personas[personas$sexo == 'F',]
peso
altura sexo
Maria
75 1.75 F
Alicia
98 1.89 F
Judith
80 1.77 F
> personas[personas$sexo == 'F' & personas$altura > 1.75,]
peso
altura sexo
Alicia
98 1.89 F
Judith
80 1.77 F
> attach(personas)
> personas[sexo == 'F' & altura > 1.75,]
> detach(personas)
- Agrupación y número de ocurrencias:
> attach(personas)
> table(sexo) # Contabiliza ocurrencias por sexo
sexo
F
M
3
1
> table(sexo, peso) # Contabiliza ocurrencias por sexo y peso
peso
sexo
67 75 80 98
F
0 1 1 1
M
1 0 0 0
> detach(personas)
- Cruces
> x = data.frame(id = c(1, 2, 3, 4), datos = c(70, 30, 40, 100))
> y = data.frame(id = c(1, 2, 3, 5), datos = c(.5, .3, .2, .1))
# Si no se especifica "by", el cruce se realiza por el
# resultado de intersect(names(x), names(y))
# Si queremos cruzar por campos con nombres diferentes en cada lado,
# podemos usar "by.x" y "by.y"
> merge(x, y, by=c("id")) # Inner Join
id datos.x datos.y
1
1 70 0.5
2
2 30 0.3
3
3 40 0.2
> merge(x, y, by=c("id"), all.x = TRUE) # Left Join
id datos.x datos.y
1
1 70 0.5
2
2 30 0.3
3
3 40 0.2
4
4 100 NA
> merge(x, y, by=c("id"), all.y = TRUE) # Right Join
id datos.x datos.y
1
1 70 0.5
2
2 30 0.3
3
3 40 0.2
4
5 NA 0.1
> merge(x, y, by=c("id"), all = TRUE) # Outer Join
id datos.x datos.y
1
1 70 0.5
2
2 30 0.3
3
3 40 0.2
4
4 100 NA
5
5 NA 0.1
- Sumarizaciones
> z = data.frame(id = c(1, 2, 2, 4, 5, 4), datos1 = c(70, 30, 40, 100,
22, 33), datos2 = runif(6))
> aggregate(z[,2:length(z)], by=list(z$id), FUN=sum) # FUN puede ser
# cualquier funcion de R (sum, mean...) o definida por nosotros
Group.1 datos1 datos2
1
1 70 0.6584459
2
2 70 0.4219876
3
4 133 0.8075231
4
5 22 0.5956466
- Otros tratamientos en conjunto
# Aplicar "mean" cogiendo las filas (1) como argumento > apply(z, 1, FUN=mean) [1] 23.870088 10.987341 14.264821 34.780961 9.259135 12.448966 # Aplicar "mean" cogiendo las filas (2) como argumento > apply(z, 2, FUN=mean) id datos1 datos2 3.00000 49.16667 0.63899 # Sumar 1 a todos los elementos > z + 1 id datos1 datos2 1 2 71 1.610265 2 3 31 1.962024 3 3 41 1.794464 4 5 101 1.342882 5 6 23 1.777406 6 5 34 1.346899
Listas
Una lista puede contener vectores, matrices y data frames:
- Creación
> h = list()
- Acceso a elementos
> h[["a"]] = 1 # Acceso a elemento mediante doble corchete [[]] > h[["b"]] = c(1, 2, 3) > h[["a"]] = NULL # Borra un elemento
Introducción manual de datos
Para introducir rápidamente datos en pequeños vectores podemos utilizar la funcion "c" tal y como hemos visto en la sección anterior:
> x = c(2,3,0,3,1,0,0,1) > x [1] 2 3 0 3 1 0 0 1
Esto nos generaría un vector denominado x con 8 valores. También podríamos utilizar la función scan() que nos facilita la vida a la hora de hacer copy/pastes:
> x = scan() 1: 3 4 1 1 3 4 3 3 1 3 2 1 2 1 2 3 2 3 1 1 1 1 4 3 1 26: Read 25 items
Esta función lee todos los elementos que le indiquemos hasta que se encuentra una linea en blanco.
Por otra parte, podemos crear y editar gráficamente un vector mediante:
> data.entry(x2=c(NA)) # Creamos el vector numérico x2 > data.entry(x3=c(“”)) # Creamos el vector de cadenas x2 > x = edit(as.data.frame(NULL)) # Creamos un dataframe
O editamos un/os vector/es ya existente/s:
> x = c(1, 2) > y = c(“a”, “b”) > data.entry(x) # Editamos x (variable ya existente) formateado en una tabla > data.entry(x, y) # Editamos x e y (variables existentes)
También podemos editar visualmente un dataframe mediante el comando “edit”:
> x = data.frame(c(1,2), c("a","b"))
> x = edit(x)
Librerías y datos de ejemplo
Como hemos comentado en la introducción, R dispone de una inmensa cantidad de librerías y paquetes con funciones y datos preparados para todo tipo de análisis. Ver listado completo en:http://cran.r-project.org/web/views/index.html
En Ubuntu GNU/Linux se almacenan los paquetes instalados en el directorio "/usr/lib/R/library/", y para instalar nuevos debemos o bien utilizar el sistema de paquetes (buscar paquetes con "aptitude search r-cran") o bien descargandolo y descomprimiendolo en el directorio anterior.
Para la plataforma Windows, los paquetes se almacenan habitualmente en el directorio donde se haya realizado la instalación de R (p.ej. "c:\Program files\R\"), en la carpeta library. En este caso, la interfaz gráfica Rgui nos permite instalar paquetes nuevos de CRAN automáticamente, o en su defecto tendremos que descargarlos y descomprimirlos en el directorio indicado.
Para listar las librerías disponibles podemos utilizar "library" sin argumento, y cuando necesitemos hacer uso de una de ellas tendremos que especificarlo para que sea cargada:
> library(MASS)
Por otra parte, tanto las librerías de R acostumbran a incorporar datos de ejemplo que podemos listar con el comando "data()", para después cargar un conjunto de datos:
> data(EuStockMarkets)
A partir de ese momento podremos acceder a los datos por su nombre, por ejemplo "EuStockMarkets".
Importar/exportar datos
Ficheros
En R es muy habitual trabajar con ficheros CSV (campos delimitados por ,) para así poder ser utilizado fácilmente en hojas de cálculo: Veamos como podríamos exportar a dicho formato:
# Exportar
> df = data.frame(runif(10), runif(10), runif(10))
> names(df) = c("dato1", "dato2", "dato3")
> write.table(df, file = "dataframe1.csv", sep = ",",
col.names = NA, qmethod = "double") # Con id/nombres de filas
> write.table(df, file = "dataframe2.csv", sep = ",",
row.names = FALSE, qmethod = "double") # Sin id/nombre de filas
La importación seria igual de sencilla:
# Importar de nuevo a R
> read.table("dataframe1.csv", header = TRUE, sep = ",",
row.names = 1) # Con id/nombres de filas
> read.table("dataframe2.csv", header = TRUE, sep = ",")
# Sin id/nombres de filas
La función read.table es muy versátil y permite multitud de opciones que nos permitirían cargar prácticamente cualquier fichero con campos delimitados por un signo (ver ayuda ?read.table).
En caso de que queramos cargar un fichero con campos de longitud fija (sin separación por signo), podemos utilizar la siguiente función:
# Importa una tabla de 3 columnas de 5, 10 y 3 caracteres respectivamente > read.fwf(ff, widths=c(5,10,3))
En ambos casos podemos utilizar el argumento colClasses para indicar el tipo de elemento que hay en cada columna. Si no se especifica R los intenta inferir, pero hay casos en los que nos interesará especificarlos. Por ejemplo:
> read.table("test.csv", header = TRUE, sep = ",",
colClasses="character")
> read.table("test.csv", header = TRUE, sep = ",",
row.names = 1, colClasses=c("Date", "character",
rep(10, "numeric")))
Veamos un ejemplo de fichero:
Mes Dia Edad Dato1 Dato2 Jan 13 25 15 115 Feb 15 32 24 226 Mar 15 24 34 228 Apr 31 52 63 420 May 16 34 29 208 Jun 31 42 75 492 Jul 24 34 67 436 Aug 15 34 47 316 Sep 13 55 37 277 Oct 29 54 68 525 Nov 20 87 82 577 Dec 17 35 61 401 Jan 21 36 64 620 Feb 26 58 80 652 Mar 24 75 70 495 Apr 21 70 74 514
Si el anterior fichero tiene por nombre "sample", ejecutariamos:
datos = read.csv("sample", header=TRUE, sep=" ",
quote="")
En lugar del nombre del fichero, podríamos indicar file.choose() para que R nos pregunte el fichero que queremos cargar.
Bases de Datos: MySQL
R almacena todos los datasets que cargamos o generamos en memoria, esto puede suponer una limitación si queremos trabajar con una cantidad de información que supera la memoria RAM de nuestro sistema (actualmente, esto significa que los datos deben ser realmente muy muy grandes como para que sea una limitación).
Podríamos trabajar importando y exportando continuamente ficheros CSV, pero no resulta eficiente a la larga. Como puntos de mejora, podemos apoyarnos en una BBDD para almacenar la información en disco y únicamente cargar a memoria aquello que necesitemos. Lo recomendable es utilizar R con MySQL, no obstante también existen librerías para acceder a otras BBDD (p.ej. SQLite pero se desaconseja por problemas de corrupción de datos y las propias limitaciones de la BBDD).
Para instalar la libreria rmysql en Ubuntu GNU/Linux:
apt-get install r-cran-rmysql
En Windows podemos hacerlo fácilmente también a través de la interfaz gráfica Rgui.
Para poder trabajar, necesitaremos también la BBDD Mysql. Para Ubuntu GNU/Linux podemos utilizar:
apt-get install mysql-client mysql-server mysql-query-browser mysql-admin
En Windows lo mejor es bajarnos la versión no instalable de MySQL y >MySQL GUI Tools.
Una vez bajado, lo descomprimimos y podemos iniciar la BBDD ejecutando bin/mysqld.exe, así como acceder a la misma mediante "bin/mysql.exe".
Para poder empezar a trabajar, debemos crear una BBDD llamada en MySQL que, por ejemplo, podemos llamar "R":
# mysql Welcome to the MySQL monitor. Commands end with ; or \g. Your MySQL connection id is 14 Server version: 5.0.51a-3ubuntu5.4 (Ubuntu) Type 'help;' or '\h' for help. Type '\c' to clear the buffer. mysql> CREATE DATABASE R; Query OK, 1 row affected (0.00 sec) mysql> quit Bye
Esta tarea también podemos hacerlo desde la interfaz gráfica MySQL Admin o el Query Browser, conectando contra localhost y usuario "root" (sin contraseña).
Par acceder a la BBDD desde R:
> Sys.setenv(MYSQL_HOME="C:/Users/xxx/Documents/NoBackup/Applications/MySQL/mysql-5.1.31-win32")
# Solo en Windows
> library(RMySQL)
Loading
required package: DBI
# Conectamos y seleccionamos la BBDD "R"
> con = dbConnect(dbDriver("MySQL"), dbname = "R",
user="root", host="localhost")
Únicamente en el caso de Windows tenemos que establecer la variable de entorno MYSQL_HOME apuntando al directorio donde hayamos descomprimido/instalado la MySQL.
Vamos a crear una tabla que contenga un dataframe de R y a continuación lo descargamos de la memoria:
# Creamos un data frame
> x = data.frame(c(1,2,3), c(5200, 1300, 3400), c("Juan",
"Judith", "Ana"))
> names(x) = c("id", "datos", "nombre")
# Creamos una tabla con el data frame anterior
> dbWriteTable(con, "Datos", x)
# Eliminamos el data frame de la memoria
> rm(x)
Ahora vamos a recuperar el dataframe leyendo la información de la BBDD:
# Recuperamos de la tabla el dataframe > x = dbReadTable(con, "Datos") # Listamos los campos de la tabla "Datos > dbListFields(con, "Datos")
También podemos construir dataframes a partir de consultas SQL, esto nos permite aprovechar la potencia de MySQL para hacer filtros, JOINs, sumarizaciones, etc (aunque todas esas operaciones también las podemos hacer con R sin depender de la BBDD):
# Construimos un dataframe a partir de una consulta > x = dbGetQuery(con, "SELECT id, nombre FROM Datos WHERE datos > 2000")
Otras operaciones:
# Listado de tablas
> dbListTables(con)
# Enviamos una sentencia DROP/CREATE/... que no devuelve información
# en este caso borramos la tabla Datos
> dbSendQuery("DROP TABLE Datos")
# Desconectamos
> dbDisconnect(con)
Finalmente, otro punto interesante de utilizar MySQL es que podemos usar el Query Browser gráfico y acceder al contenido de las tablas y mediante SQL realizar filtros rápidos, ordenación de datos, exportar a Excel/CSV, etc. Esto puede ser de gran utilidad para entender los datos de forma más rápida sin depender de la consola de R.
Distribuciones de probabilidad
Para cada distribución de probabilidad disponible R pone a nuestro servicio 4 tipos de funciones:
- r: generación de números aleatorios
- p: probabilidades
- d: función de densidad
- q: cuantiles
Por ejemplo, para la distribución uniforme tenemos las funciones runif, punif, dunif y qunif.
Generación de números aleatorios
Será habitual que utilicemos las funciones de generación de números aleatorios de R para construir vectores/matrices/data frames con números que siguen una distribución determinada.
En esta sección vamos a tratar con las funciones ‘r’ para la generación de números y las ‘d’ (densidad) para el dibujo de gráficas (que nos facilitarán entender las distribuciones).
Como gráficas generaremos histogramas y curvas para las distribuciones de números continuos y puntos para las discretas:
- Distribuciones de números continuos
# Uniforme (misma probabilidad para todos los números) > x=runif(100, 0, 1) > hist(x,probability=TRUE,col=gray(.9),main="uniform on [0,1]") > curve(dunif(x,0,1),add=T) # Densidad # Normal = Gaussiana > x=rnorm(100, mean=50, sd=25) # 100 números de media situada en 100 (en lugar de 0) y desviación estandar de 25 > hist(x,probability=TRUE,col=gray(.9),main="normal mu=0,sigma=1") > curve(dnorm(x, mean=50, sd=25),add=T) # Densidad # Exponencial > x=rexp(100,50) # 100 números con media en 50 > hist(x,probability=TRUE,col=gray(.9),main="exponential mean=2500") > curve(dexp(x,50),add=T) # Densidad ## Chi cuadrado con 1 grado de libertad > rchisq(100, 1) # Dibujamos Chi cuadrado con varios grados de libertad > curve(dchisq(x,1), xlim=c(0,10), ylim=c(0,.6), col='red', lwd=3) > curve(dchisq(x,2), add=T, col='green', lwd=3) > curve(dchisq(x,3), add=T, col='blue', lwd=3) > curve(dchisq(x,5), add=T, col='orange', lwd=3) ## Student's T > rt(100, 1) # Dibujamos varios grados de libertad > curve( dt(x,1), xlim=c(-3,3), ylim=c(0,.4), col='red', lwd=2 ) > curve( dt(x,2), add=T, col='blue', lwd=2 ) > curve( dt(x,5), add=T, col='green', lwd=2 ) > curve( dt(x,10), add=T, col='orange', lwd=2 ) > curve( dnorm(x), add=T, lwd=3, lty=3 ) ## Fisher's F > rf(100, 1, 1) > curve(df(x,1,1), xlim=c(0,2), ylim=c(0,.8), lty=2) > curve(df(x,3,1), add=T) > curve(df(x,6,1), add=T, lwd=3) > curve(df(x,3,3), add=T, col='red') > curve(df(x,6,3), add=T, lwd=3, col='red') > curve(df(x,3,6), add=T, col='blue') > curve(df(x,6,6), add=T, lwd=3, col='blue') ## Lognormal law > rlnorm(100) > curve(dlnorm(x), xlim=c(-.2,5), lwd=3, main="Log-normal distribution")
- Distribuciones de números discretos
# Muestreo con repetición y probabilidad uniforme x = sample(1:100, 10, replace=T) # 10 número del 1 al 100 # Binomial > n = 20 > p = 0.5 > x = rbinom(100, n, p) # 100 números: 20 intentos con una probabilidad de exito de cada intento del 50% (se suman el número de exitos) > hist(x,probability=TRUE,) # Dibujamos los puntos donde se encontraría la densidad de la distribución > xvals=0:n > points(xvals,dbinom(xvals,n,p),type="h",lwd=3) # Linea > points(xvals,dbinom(xvals,n,p),type="p",lwd=3) # Punto # Hypergeometric > x = rhyper(100, 300, 100, 100) > hist(x,probability=TRUE,) # Poisson > x = rpois(100, 1) > hist(x,probability=TRUE,) # Geometrica > x = rgeom(100, .5) > hist(x,probability=TRUE,) # Binomial negativa > x = rnbinom(100, 10, .25) > hist(x,probability=TRUE,)
Muestreo aleatorio
Las funciones de la sección anterior generan números aleatorios que pueden repetirse, sin embargo si queremos hacer un muestreo necesitamos una función que nos permita generar una serie de números aleatorios no repetidos (p.ej. Como cuando tenemos 10 bolas en un bombo y vamos extrayéndolas, el primer número que sale no se volverá a repetir). Para esto utilizaremos la función sample, la cual ya aparecía en la sección anterior pero con el argumento "replace=T" que hace variar su comportamiento. Por ejemplo cogemos 10 números entre 1 y 100 donde todos tienen las mismas probabilidades de ser elegidos:
> sample(1:100,10) [1] 53 99 63 92 60 87 54 25 51 73
Con esta función podemos extraer filas de una matriz, utilizando los números generados como indices:
> personas[sample(1:length(personas), 2),] peso altura sexo 1 75 1.75 F 3 67 1.61 M
Generación de números secuenciales
Si en lugar de números aleatorios queremos generar secuencias, podemos utilizar las siguientes notaciones:
> x = 1:10 [1] 1 2 3 4 5 6 7 8 9 10 > x = seq(1, 10, 2) # Números del 1 al 10 con incrementos de 2 [1] 1 3 5 7 9 > x = letters[1:11] [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k"
Probabilidades
Con las funciones ‘p’ y ‘q’ podemos responder preguntas del siguiente tipo:
- ¿Cual es la probabilidad que en una distribución normal con media en 0 y desviación estándar 1 un número sea inferior a 0.7?
> pnorm(.7, mean=0, sd=1) [1] 0.7580363
- ¿Y que sea superior a 0.7?
> pnorm(.7,lower.tail=F) [1] 0.2419637
- Que valor limita el área que representa un 75% de probabilidades (área de la derecha):
> qnorm(.75) [1] 0.6744898
Visualización de datos de una variable
- Gráficos tipo tarta: Genera tantas divisiones de la tarta como elementos en el vector.
> gastos = c(Gimnasio=60, Dietas=345, Alquiler=800, Transporte=200) > pie(gastos)
- Gráficos de barras: Visualiza tantas barras como elementos existen en el vector.
> gastos = c(Gimnasio=60, Dietas=345, Alquiler=800, Transporte=200) > barplot(gastos)
- Histogramas: Muestra un número determinado de barras, donde cada una representa un rango y la “altura” depende del número de elementos del vector que se encuentran en dicho rango. Nos permiten ver donde se concentran los números.
En el caso de pruebas de datos, un histograma nos puede servir para mostrar claramente donde se concentran las diferencias detectadas.
> transferencias=abs(rnorm(100, mean=100, sd=50)) > hist(transferencias) # En el eje Y se visualiza el número de # ocurrencias > hist(transferencias, probability=TRUE) # En el eje Y se visualizan # las probabilidades > hist(transferencias, breaks=4) # Los datos se agrupan en 5 (hay 4 divisiones)
- Histogramas con marcas de ocurrencia
> hist(transferencias) # En el eje Y se visualiza el número de ocurrencias > rug(transferencias) # Marcas con las ocurrencias
- Histogramas con linea de densidad: facilita la identificación de la distribución de probabilidad.
> hist(transferencias, probability=TRUE) # Para visualizar la densidad # se requiere que el histograma refleje la probabilidad > lines(density(transferencias), col="red", lwd=3) # Linea # con la densidad de probabilidad
- Boxplot: Visualiza los siguientes elementos de abajo a arriba (o izquierda a derecha)
- Mínimo
- Q1
- La mediana
- Q3
- Q3+1,5*IQR
- Puntos extremadamente separados
> transferencias.nacionales=abs(rnorm(100, mean=130, sd=50))
> transferencias.internacionales=abs(rnorm(100, mean=150, sd=60))
> transferencias.nacionales[1] = 310
> transferencias.nacionales[2] = 280
> transferencias.internacionales[1] = 350
> boxplot(transferencias.nacionales, transferencias.internacionales,
names=c("Nacionales", "Internacionales")), horizontal=TRUE)
- Histograma y boxplot simultáneos: Nos permite observar la distribución de probabilidad y los quantiles, media, etc.
"hist.and.boxplot" <-
function (x,...) {
## Plot both the histogram and the boxplot
## to show relationship easily.
op<-par(no.readonly=TRUE)
on.exit(par(op))
layout(matrix(c(1,2),2,1),heights=c(3,1))
par(mai=c(1,1,1,1)/2)
hist(x,xlab=FALSE,col=gray(0.95),yaxt='n',...)
rug(x)
boxplot(x,horizontal=TRUE)
}
> transferencias=abs(rnorm(100, mean=100, sd=50)) > hist.and.boxplot(transferencias)
- Stemplot: stemplot es un gráfico en modo texto similar al histograma pero invertido lateralmente, útil para conjuntos de datos pequeños:
> x = rnorm(100) > stem(x) The decimal point is at the | -3 | 0 -2 | -1 | 975443333222110 -0 | 998888666666555544333333322222211000 0 | 11222222223333344455667779 1 | 001223334444666677 2 | 016 3 | 4
-
Como se puede observar, en el ejemplo debemos considerar que la coma empieza en la | y por tanto tenemos un punto en -3 (-3.0), ninguno en -2, varios en -1 (-1.9, -1.7, -1.5, 2 puntos -1.4, etc.) y así sucesivamente. El gráfico es más representativo si inclinamos la cabeza hacia la derecha.
- Histograma y boxplot simultáneos:
"plot.hist.and.box" <-
function (x,...) {
## Plot both the histogram and the boxplot
## to show relationship easily.
op<-par(no.readonly=TRUE)
on.exit(par(op))
layout(matrix(c(1,2),2,1),heights=c(3,1))
par(mai=c(1,1,1,1)/2)
hist(x,xlab=FALSE,col=gray(0.95),yaxt='n',...)
rug(x)
boxplot(x,horizontal=TRUE)
}
> x = rnorm(100) > plot.hist.and.box(x)
Visualización de datos de dos variables
- Gráficos de barras: Imaginemos que tenemos dos variables clasificadas, la primera de ella corresponde a si la persona es fumadora o no y la segunda al número de horas que dedica al estudio (1 = menos de 5 horas, 2 = entre 5 y 10 horas, 3 = más de 3 horas):
> fumador = c("S","N","N","S","N","S","S","S","N","S")
> estudio = c(1,2,2,3,3,1,2,1,3,2)
> t = table(fumador,estudio) # Genera tabla de contingencias
> t
estudio
fumador 1 2 3
N 0 2 2
S 3 2 1
-
Como observamos, podemos utilizar de nuevo la función table() para contabilizar el número de apariciones de cada categoría. Si queremos conocer la probabilidad de caer en cada una de las parejas de categorías:
> options(digits=3) # Limitamos a 3 decimales como máximo > prop.table(t) # Todos los importes de la tabla anterior sumaran 1 estudio fumador 1 2 3 N 0.0 0.2 0.2 S 0.3 0.2 0.1
-
Para visualizar la información mediante un gráfico de barras:
> barplot(table(estudio,fumador),main="Horas de estudio segun fumador/no fumador",
+ beside=TRUE, # Separar las categorias en varias barras
+ legend.text=c("Menos 5 horas","5-10","Más de 10 horas"))
- Puntos en un plano: Dadas las coordenadas X e Y de un conjunto de puntos, dibujamos sobre un plano.
> punto.espacio = data.frame(x=c(1, 5, 10), y=c(200, 50, 100)) > puntos.espacio x y 1 1 200 2 5 50 3 10 100 > plot(puntos.espacio)
- Puntos en un plano con linea de regresión: Si entendemos que X es un valor que controlamos (p.ej. tabletas de chocolate que comemos) e Y un valor que medimos (p.ej. grado de felicidad), podemos utilizar una regresión lineal para predecir las Y para un valor desconocido de X (p.ej. cual sera mi grado de felicidad si como 20 tabletas de chocolate).
Veamos como podemos dibujar una linea en nuestro gráfico que corresponda con la regresión lineal:
> x.chocolate = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9) > y.felicidad = c(4, 7, 10, 6, NA, NA, NA, 1, 2, 1) > plot(x.chocolate, y.felicidad) > abline(lm(y.felicidad ~ x.chocolate), col="red") # lm = linear model (variable y como función lineal de x)
- Puntos en un plano con linea de regresión resistente: linea con capacidad de "resistir" a la influencia de un porcentaje determinado de datos muy dispersos (p.ej. reducir el efecto que tiene 2 datos puntuales extremos):
> library(MASS) # Cargamos un paquete que contiene la función rlm > x.chocolate = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9) > y.felicidad = c(4, 7, 10, 6, NA, NA, NA, 1, 2, 1) > plot(x.chocolate, y.felicidad) > abline(lm(y.felicidad ~ x.chocolate), col="red") # Linea de regresión > abline(rlm(y.felicidad ~ x.chocolate), col="red", lty=2) # Linea de regresión resistente
- Puntos en un plano con linea de regresión e histogramas de X e Y
"scatterplot" <-
function(x,y,...) {
def.par = par(no.readonly = TRUE)# save default, for resetting...
n<-length(x)
xhist = hist(x,sqrt(n), plot=FALSE)
yhist = hist(y, sqrt(n), plot=FALSE)
top = max(c(xhist$counts, yhist$counts))
xrange = c(min(x),max(x))
yrange = c(min(y),max(y))
nf = layout(matrix(c(2,0,1,3),2,2,TRUE), c(3,1), c(1,3), TRUE)
layout.show(nf)
par(mar=c(3,3,1,1))
plot(x, y, xlim=xrange, ylim=yrange, xlab="x", ylab="y",...)
abline(lm(y~x))
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0,col=gray(.95))
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes=FALSE, xlim=c(0, top),
space=0,col=gray(.95), horiz=TRUE)
par(def.par)#- reset to default
}
> x = abs(rnorm(100, mean=100, sd=50) ) > y = rexp(100, rate=.1) > scatterplot(x, y)
- Agrupa en categorías y visualiza boxplot (Mínimo, Q1, mediana, Q3, Q3+1,5*IQR, outliers):
> x = runif(10) # Valores > y = c(1, 1, 1, 2, 2, 1, 2, 3, 3, 2) # Categorias > boxplot(x ~ y)
Identificar puntos en un gráfico
Mostramos un gráfico a partir de unos datos:
> x = runif(100) > y = runif(100) > plot(x,y)
A continuación utilizamos ‘identify’ para indicar que queremos localizar 1 punto:
> identify(x, y, n=1) # Hacemos click en el gráfico y nos devuelve el indice del punto [1] 56
Para obtener los valores del punto, utilizamos el indice en los vectores usados para el gráfico:
> c(x[56],y[56]) [1] 0.2514029 0.9420260
Identificar la distribución de los datos
Determinadas pruebas estadísticas que podemos hacer implican que los datos se encuentran repartidos probabilisticamente según una distribución determinada. Por ejemplo, los T-Tests son útiles cuando tenemos datos que siguen una distribución normal (gaussiana), pero no cuando son uniformes. Por eso es importante conocer que herramientas tenemos para identificar la potencial distribución que describe mejor a nuestros datos.
Una primera validación podria ser la ejecución de la prueba Shapiro-Wilk que nos permite comprobar si una distribución es normal (ver sección "Pruebas estadísticas"). Pero también podemos apoyar nuestro análisis en gráficos que nos ayudaran a identificar o confirmar la distribución. La siguiente función que nos permite visualizar el histograma y densidad de unos datos conjuntamente con la densidad teorica de una distribución determinada:
hist.density = function(x, FUN=dnorm, ...) {
hist(x, probability=TRUE, breaks=20, col="light blue")
rug(jitter(x, 5))
# Densidad de x
points(density(x), type='l', lwd=3, col='red')
# Densidad teorica de una distribución concreta
f = function(t) {
# t => Graph's X
FUN(t, ... )
}
curve(f, add=T, col="red", lwd=3, lty=2)
}
Probamos a comparar un par de conjuntos de datos contra la distribución normal:
x = rnorm(100) hist.density(x, FUN=dnorm, mean=mean(x), sd=sd(x)) x = rexp(100) hist.density(x, FUN=dnorm, mean=mean(x), sd=sd(x))
Y ahora, contra la exponencial:
x = rnorm(100) hist.density(rnorm(100), FUN=dexp) x = rexp(100) hist.density(rexp(100), FUN=dexp)
Otra técnica que podemos utilizar es la visualización del Quartil-Quartil plot:
qq = function (y, ylim, quantiles=qnorm,
main = "Q-Q Plot", xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles", plot.it = TRUE, ...)
{
y = y[!is.na(y)]
if (0 == (n = length(y)))
stop("y is empty")
if (missing(ylim))
ylim = range(y)
x = quantiles(ppoints(n))[order(order(y))]
if (plot.it)
plot(x, y, main = main, xlab = xlab, ylab = ylab, ylim = ylim,
...)
# From qqline
y = quantile(y, c(0.25, 0.75))
x = quantiles(c(0.25, 0.75))
slope = diff(y)/diff(x)
int = y[1] - slope * x[1]
abline(int, slope, ...)
invisible(list(x = x, y = y))
}
La anterior función dibujara el QQ plot de los datos actuales y la linea que une el primer y tercer cuartil de los mismos. Las reglas para interpretar el gráfico son las siguientes:
1.- Si la parte izquierda se encuentra por ENCIMA de la linea => la distribución de los datos es MAYOR en la izquierda que la de la distribución comparada
2.- Si la parte izquierda se encuentra por DEBAJO de la linea => la distribución de los datos es MENOR en la izquierda que la de la distribución comparada
3.- Si la parte derecha se encuentra por ENCIMA de la linea => la distribución de los datos es MENOR en la derecha que la de la distribución comparada
4.- Si la parte derecha se encuentra por DEBAJO de la linea => la distribución de los datos es MAYOR en la derecha que la de la distribución comparada
Probamos a comparar unos datos contra la distribución uniforme, exponencial y normal:
qq(runif(100), quantiles=qunif) # La linea se solapa con los puntos qq(runif(100), quantiles=qexp) # Concentraciones diferentes qq(runif(100), quantiles=qnorm)
Diferentes formas de entender la probabilidad
Probabilidades frecuentista
Asociados a sistemas físicos aleatorios como la ruleta rusa, los dados o los átomos reactivos. En estos sistemas un tipo de evento tiende a ocurrir con una frecuencia persistente. Las probabilidades de frecuencia intentan describir justamente estas frecuencias estables.
Las 2 teorías principales:
- Teorías de la frecuencia probabilística: la probabilidad es la frecuencia relativa con la que se da un suceso, en el límite de un número muy grande de sucesos.
- Teorías de la propensión probabilística: la probabilidad es una propensión o tendencia física de un tipo de evento, por ejemplo que al lanzar muchas veces (teoría de los grandes números) una moneda al aire, el 50% de veces saldrá cara y el 50% cruz. Dado que según las teorías de la frecuencia probabilística no existe una frecuencia en el lanzamiento de una moneda, solo en un grupo de lanzamientos, las teorías de la propensión hablan de tendencia/propensión en un único lanzamiento.
La mayor parte de funciones de R que veremos en este documento parten de este enfoque.
Probabilidades bayesienas o basadas en evidencias
Los estadísticos de la escuela bayesiana aceptan la existencia e importancia de las probabilidades físicas, pero consideran que el cálculo de probabilidades basadas en evidencias también es válido y necesario. Estos entienden la probabilidad en base a 2 interpretaciones diferentes:
- Interpretación subjetiva: Se interpreta la probabilidad como el nivel de creencia de un individuo sobre una proposición.
- Interpretación objetiva: Las reglas de la estadística bayesiana pueden ser justificadas, racionalizadas e interpretadas como una extensión de la lógica aristotélica. Esta supone que la mente reproduce sólo la realidad, la existencia de las cosas tal y como son.
La inferencia bayesiana usa un estimador numérico del grado de creencia en una hipótesis aún antes de observar la evidencia y calcula un estimador numérico del grado de creencia en la hipótesis después de haber observado la evidencia. Veamos un ejemplo:
- Probabilidad subjetiva: El 80% de los correos que contienen la palabra "free" y "drugs" son SPAM.
- Nueva evidencia: Recibimos un email que contiene "free" y "drugs".
- Probabilidad: ¿El nuevo email es SPAM? Hay una probabilidad del 80% de que el mensaje sea SPAM.
- Después de observar la nueva evidencia: Comprobamos si el mensaje realmente es SPAM y se ajustan las probabilidades anteriormente calculadas.
La probabilidad subjetiva inicial puede estar basada en una creencia soportada por estudios sobre los correos basura que circulan por internet, pero no es una probabilidad empírica. No sabemos que tipo de correo recibimos en nuestro buzón, quizás trabajamos en el sector farmacéutico y por tanto es habitual recibir correos con la palabra "drugs". Sin embargo, la probabilidad bayesiana tiene la capacidad de "aprender" y ajustarse a las evidencias que se van observando. De hecho podríamos incluir fácilmente la evaluación de más variables (no solo la existencia de "drugs" y "free").
Por otra parte, los métodos bayesianos permiten incorporar de forma natural cualquier tipo de información previa. Si tuviésemos mensajes de SPAM que hemos recibido, podríamos calcular las probabilidades iniciales.
En esencia, los seguidores de la estadística tradicional sólo admiten probabilidades basadas en experimentos repetibles y que tengan una confirmación empírica mientras que los llamados estadísticos bayesianos permiten probabilidades subjetivas.
Es importante tener en cuenta que las probabilidades bayesianas abren un debate sobre si contribuyen con justificaciones válidas o simplemente creencias, dado que con la misma información las personas pueden tener diferentes interpretaciones en función de su grado de creencia previa.
Pruebas estadísticas
R nos pemite efectuar pruebas estadísticas, las cuales parten de 2 hipotesis:
- Hipotesis nula (H0): por ejemplo, el tabaco no incrementa el riesgo de cancer.
- Hipotesis alternativa (H1): por ejemplo, puede ser simétrica como "el tabaco incrementa/decrementa el riesgo de cancer" o asimétrica "el tabaco incrementa el riesgo de cancer". En la segunda ya estamos eliminando una posibilidad y por tanto debemos hacerlo con cuidado.
Una prueba nunca nos dirá si una hipotesis es cierta, sino que refutaran la hipotesis o fallaran en la refutación ("K. Popper: we never prove that something is true, we merely continuously try to prove it wrong and fail to do so").
Otra clasificación posible de las hipotesis en función de los parámetros que conocemos:
- Hipotesis simple: por ejemplo, queremos testear si una muestra procede de una población de la cual conocemos su distribución y su media y varianza.
- Hipotesis compuesta: si H1 es cierto, conocemos la distribución pero no la varianza.
En las pruebas estadísticas podemos cometer dos tipos de errores:
- Error del tipo I: Refutar incorrectamente la hipotesis nula
- Error del tipo II: No refutar la hipotesis nula, cuando se deberia refutar.
Cuando realicemos las pruebas, también obtendremos el intervalo de confianza en el cual podremos asegurar que se encuentra un parámetro (p.ej. media, varianza, etc.) de la población a la que pertenece la muestra que estamos analizando. El intervalo tiene además asociado un porcentaje de probabilidades de que sea correcto, habitualmente se trabaja con el 95%.
Finalmente, las pruebas siempre generarán un p-Value, que corresponde con la probabilidad de cometer el error de refutar incorrectamente la hipótesis nula (tipo I). Lo que significa que si p-value es muy pequeño, refutaremos la hipótesis nula.
Antes de realizar ninguna prueba se establece el limite (alpha) a partir del cual p-value refutará la hipótesis nula. Es habitual utilizar el limite 0.05, es decir, si es inferior a 0.05 entonces refutaremos la hipótesis.
Existen multitud de pruebas, entre las que podemos encontrar:
| Objetivo | Test paramétrico | Test no paramétrico |
| Comparar 2 medias | Student’s T test | Wilcoxon’s U test |
| Comparar más de 2 medias | ANOVA | Kruskal-Wallis test |
| Comparar 2 varianzas | Fisher’s F Test | Ansari-Bradley o Mood test |
| Comparar más de 2 varianzas | Barlett test | Fligner test |
Otros conceptos importantes:
- Robustez: un test es robusto si es valido cuando las suposiciones dejan de cumplirse
- Resistencia: un valor (media, mediana, varianza…) es resistente si no depende de valores extremos (outliers). Por ejemplo, la media no es resistente y la mediana si. Por ejemplo, como estimadores robustos tenemos:
- Ubicación (media): trimmed mean, median
- Dispersión (desviación estandar): IQR, MAD, Sn, L-moments, MCD
Prueba de distribución normal
El test Shapiro-Wilk nos permite probar si un conjunto de datos tiene una distribución normal (gausiana). Para ello definimos que la hipótesis nula es que la distribución es la normal y solo la refutaremos si el p-value es superior a 0.05:
shapiro.test(rnorm(100))$p.value [1] 0.4749011 # > 0.05 no podemos refutar shapiro.test(runif(100))$p.value [1] 0.005700384 # < 0.05 refutamos, no es normal
Prueba de proporción
Supongamos que hacemos un muestreo y seleccionamos 100 expedientes hipotecarios de una entidad y de estos, 42 no tienen adjunta documentación sobre el análisis del riesgo efectuado. Esto lo podríamos traducir en 3 afirmaciones diferentes (de menor a mayor exactitud):
- El 42% de los expedientes están incompletos
- El 42% de los expedientes están incompletos con un margen de error del 9%
- El 42% de los expedientes están incompletos con un margen de error del 9% y una confianza del 95%
La última afirmación que es la más correcta estadísticamente, para llegar a ella podemos utilizar el test de proporción de R:
> prop.test(42, 100, conf.level=0.95) 1-sample proportions test with continuity correction data: 42 out of 100, null probability 0.5 X-squared = 2.25, df = 1, p-value = 0.1336 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.3233236 0.5228954 sample estimates: p 0.42
Veamos otro ejemplo: En el mes de Diciembre la entidad dispone de 20 préstamos, y mediante una auditoría se detecta que 5 préstamos presentan diferencias entre lo que debería haberse liquidado por interés y lo liquidado realmente. Suponemos que las diferencias se encuentran distribuidas bajo la *normal* con la media en el 0, es decir la mayoría de liquidaciones de intereses no presentan diferencias y en caso contrario, éstas están más próximas al 0 que a diferencias más distantes (p.ej. es más probable que el sistema se equivoque de unos pocos euros que de miles de euros).
Información de los 5 préstamos:
> # Capital pendiente de los préstamos > capital = c(2000,2700,2100,2600,2900,3000,3300,2100,2700,2100,2100,2400,2500,2600,2400,2200,1900,2700,1800,2600) > # TAE (se debe dividir entre 12 para obtener el nominal mensual) > interes = c(0.042,0.037,0.038,0.040,0.049,0.044,0.028,0.042,0.045,0.036,0.031,0.042,0.038,0.041,0.039,0.047,0.045,0.036,0.032,0.040) > liquidacion.sistema = c(7.00,8.10,6.65,8.67,11.24,11.00,7.39,7.35,10.12,6.10,5.09,8.40,7.92,8.88,7.80,8.62,7.12,8.10,4.80,8.67) > liquidacion.teorica = round(capital * (interes/12), 2) > diferencias = liquidacion.teorica - liquidacion.sistema > diferencias.porcentuales = round(diferencias / liquidacion.teorica, 4) > prestamos = data.frame(capital, interes, liquidacion.sistema, liquidacion.teorica, diferencias, diferencias.porcentuales) > prestamos capital interes liquidacion.sistema liquidacion.teorica diferencias 1 2000 0.042 7.00 7.00 0.00 2 2700 0.037 8.10 8.32 0.22 3 2100 0.038 6.65 6.65 0.00 4 2600 0.040 8.67 8.67 0.00 5 2900 0.049 11.24 11.84 0.60 6 3000 0.044 11.00 11.00 0.00 7 3300 0.028 7.39 7.70 0.31 8 2100 0.042 7.35 7.35 0.00 9 2700 0.045 10.12 10.12 0.00 10 2100 0.036 6.10 6.30 0.20 11 2100 0.031 5.09 5.42 0.33 12 2400 0.042 8.40 8.40 0.00 13 2500 0.038 7.92 7.92 0.00 14 2600 0.041 8.88 8.88 0.00 15 2400 0.039 7.80 7.80 0.00 16 2200 0.047 8.62 8.62 0.00 17 1900 0.045 7.12 7.12 0.00 18 2700 0.036 8.10 8.10 0.00 19 1800 0.032 4.80 4.80 0.00 20 2600 0.040 8.67 8.67 0.00
Con esta información podemos afirmar que para el mes de Diciembre:
- El 25% de los préstamos (5 de 20) presentan diferencias en la liquidación de los préstamos.
- La liquidación de intereses presenta una diferencia media de 0.08 € (1,05%) para la totalidad de los préstamos.
> mean(diferencias) [1] 0.083 > mean(diferencias.porcentuales) [1] 0.0105
Sin embargo, las auditorias son anuales y abarcan 12 meses, por tanto lo que hemos hecho nosotros ha sido coger una muestra de la totalidad de liquidaciones que se efectúan durante todo un año. Así que si queremos extrapolar estas conclusiones podríamos decir:
- El 25% de las liquidaciones de préstamos (supongamos que durante todo el año hemos tenido los mismos 20) presentan diferencias.
- Estas diferencias presentan una media de 0.08 € (1,05%), por tanto la diferencia anual asciende a 19,20€ (0.08€ x 20 préstamos x 12 meses).
No obstante, estas afirmaciones no son exactas en términos estadísticos. Lo correcto seria en primer lugar definir un porcentaje de confianza de nuestra extrapolación (p.ej. 95% de probabilidades de que la extrapolación sea cierta) y a partir de ahí identificar el margen de error resultante. En el caso del número de liquidaciones, tenemos 5 liquidaciones que liquidan incorrectamente sobre 20 que hemos validado:
> prop.test(5, 20, conf.level=0.95) 1-sample proportions test with continuity correction data: 5 out of 20, null probability 0.5 X-squared = 4.05, df = 1, p-value = 0.04417 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.0959326 0.4941155 sample estimates: p 0.25
En el test de proporción hemos indicado que 5 de una muestra de 20 cumplen una hipótesis (liquidaciones incorrectas) y queremos saber el intervalo/margen de error para una confianza del 95%, es decir, un 95% de probabilidades de que nuestra extrapolación al resto de la población sea cierta. La función ha dado como resultado que entre el 9,59% y el 49,41% de los expedientes (margen de error del 39,81%) serán incompletos con una confianza del 95%.
Por tanto, podemos afirmar que el 25% de las liquidaciones son incorrectas con un margen de error del 39,81% y una confianza del 95%. Como se puede observar, el margen de error es considerable y podríamos disminuirlo bajando el grado de confianza (p.ej. 90%) o ampliando la muestra. Por ejemplo, si hubiésemos validado 2000 liquidaciones y hubiésemos identificado 500 erróneas:
> prop.test(500, 2000, conf.level=0.95) 1-sample proportions test with continuity correction data: 500 out of 2000, null probability 0.5 X-squared = 499.0005, df = 1, p-value < 2.2e-16 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.2312709 0.2697002 sample estimates: p 0.25
El porcentaje sigue siendo del 25% pero el margen de error ha disminuido considerablemente hasta llegar a un 3,84% (entre 23,12% y 26,97%). Como es obvio, entre más grande sea la muestra, más fiable serán los resultados y las posibles extrapolaciones.
En cuanto a la segunda extrapolación sobre las diferencias medias, tendremos que realizar la prueba que se describe en la siguiente sección
Pruebas sobre la media
Las pruebas sobre la media nos pueden resultar de utilidad para, entre otros, los siguientes objetivos:
- Tenemos unos resultados sobre una muestra referentes a la media (p.ej. diferencias/errores en el cálculo de las nóminas del mes de Junio) y queremos realizar una extrapolación a toda la población (p.ej. las nóminas de todo el año).
- Si tenemos 2 muestras diferentes y queremos saber si ambas pertenecen a la misma población, podemos utilizar las pruebas para validar si ambas extrapolan a la misma media.
T-Test: Datos con distribución normal
Una variable
En el ejemplo de la sección de "Pruebas de proporción" afirmábamos que "las diferencias presentan una media de 0.08 € (1,05%), por tanto la diferencia anual asciende a 19,20€ (0.08€ x 20 préstamos x 12 meses)". Esta extrapolación no es exacta y mediante el uso del T-Test y a partir de las diferencias encontradas vamos a extraer el margen de error de la media para una confianza del 95%:
> t.test(diferencias, conf.level = 0.95) One Sample t-test data: diferencias t = 2.2532, df = 19, p-value = 0.03625 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 0.005901257 0.160098743 sample estimates: mean of x 0.083
Podemos afirmar que la diferencia media sera de 0.08 € con un margen de error del 15,41% (0,01€ al 0,16€) y una confianza del 95%. Es decir, en el peor de los casos tendremos una diferencia anual de 38,4€ (0.16€ x 20 liquidaciones mensuales x 12 meses) y en el peor 2,4 €. De nuevo, para reducir el margen de error podemos decrementar el margen de confianza o incrementar la muestra.
Veamos otro ejemplo, generamos 10 valores de distribución normal pero con una media aleatoria desconocida:
n = 10 m = rnorm(1) x = m + rnorm(n) # Desplazamos la media
Partiendo de que conocemos que la población tiene una distribución normal y desconocemos su media, vamos a testear la muestra contra las siguientes hipótesis:
- Hipótesis nula (H0): La media es 0
- Hipótesis altenrativa (H1): La media no es 0
> t.test(x) One Sample t-test data: x t = -1.4372, df = 9, p-value = 0.1845 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: -1.2043513 0.2685511 sample estimates: mean of x -0.4679001
La probabilidad de falsar H0 incorrectamente es del 18,45%, como no es inferior a 0.05 no podemos refutar H0.
Probemos con una muestra más amplia:
n = 1000 # m = rnorm(1) # Mantenemos la misma desviación x = m + rnorm(n) # Desplazamos la media > t.test(x) One Sample t-test data: x t = -19.7442, df = 999, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: -0.6541373 -0.5358651 sample estimates: mean of x -0.5950012
Ahora tenemos un p-value muy inferior a 0.05, por tanto podemos refutar H0. Estos ejemplos nos sirven para demostrar la importancia del tamaño de las muestras dado que nos pueden hacer extraer conclusiones incorrectas (en nuestra ejecución m = -0.5964069 y por tanto H0 es realmente falso).
Veamos otro ejemplo: Un fabricante de coches dice que un determinado modelo tiene un consumo de 25 litros por cada 100 KM, se pregunta a un grupo de 10 usuarios y se obtiene una media de 22 litros con una desviación estándar de 1.5. ¿Soportan estos datos la afirmación del fabricante?
En este caso, la hipótesis nula seria que la media es igual a 25 y la alternativa que no es igual:
> consumo = rnorm(100, mean=22, sd=1.5) > t.test(consumo, conf.level=.95, mu=25) One Sample t-test data: consumo t = -17.6998, df = 99, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 25 95 percent confidence interval: 21.97223 22.58265 sample estimates: mean of x 22.27744
El p-value es inferior a 0.05 y por tanto tendriamos que refutar la hipotesis nula. Es decir, los datos no soportan la afirmación del fabricante de que el coche consume 25 litros por cada 100 KM.
Dos variables
Hasta ahora hemos evaluado una única muestra, veamos un ejemplo con un par de conjuntos de datos. Imaginemos que hemos realizado un experimento sobre el tiempo de recuperación de un paciente a partir de un nuevo medicamento y para ello hemos constituido 2 grupos:
> x = c(15, 10, 13, 7, 9, 8, 21, 9, 14, 8) > y = c(15, 14, 12, 8, 14, 7, 16, 10, 15, 12)
Al grupo X se le ha suministrado el fármaco, al grupo Y únicamente un placebo. Los números corresponden a los días de recuperación después de suministrar el fármaco/placebo.
Nuestra hipótesis nula será que ambos grupos tienen la misma media y varianza, mientras que la hipótesis alternativa es que el grupo de X (el del fármaco) tiene una media inferior (descartamos que el fármaco perjudique y provoque una media superior):
> t.test(x,y,alt="less",var.equal=TRUE)
Two Sample t-test
data: x and y
t = -0.5331, df = 18, p-value = 0.3002
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf 2.027436
sample estimates:
mean of x mean of y
11.4 12.3
P-value es superior a 0.05 y por tanto, no podemos refutar que tienen la misma media (el fármaco no presenta ninguna mejora).
En caso de que tengamos 2 conjuntos de datos relacionados, por ejemplo 2 profesores que han evaluado los mismos 10 exámenes, utilizaremos el argumento “paired”:
> x = c(3, 0, 5, 2, 5, 5, 5, 4, 4, 5)
> y = c(2, 1, 4, 1, 4, 3, 3, 2, 3, 5)
> t.test(x,y,paired=TRUE)
Paired t-test
data: x and y
t = 3.3541, df = 9, p-value = 0.008468
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.3255550 1.6744450
sample estimates:
mean of the differences
En este caso el p-value es inferior a 0.05 y por tanto podemos refutar la hipotesis nula: las medias son iguales y los profesores han evaluado igual.
Simulaciones sobre el T-Test
Hagamos una simulación de N casos para observar cuantos casos no refutaríamos H0 (media en 0) para muestras que tienen la media diferente de 0 y son de distribución normal:
N = 1000
n = 3
v = c()
for (i in 1:N) {
x = rnorm(n, mean=-25) # La media de la distribución esta en -25
v = append(v, t.test(x)$p.value)
}
sum(v>.05)/N # Porcentaje de casos donde no refutariamos H0 (media en 0)
[1] 0 # Nunca
No refutariamos nunca la hipótesis. En cambio, si estamos tratando con una población de distribución diferente a la normal y con media diferente de 0:
N = 1000
n = 3
v = vector()
for (i in 1:N) {
x = runif(n, min=-49, max=1) # La media de la distribución esta en -25
v = append(v, t.test(x)$p.value)
}
sum(v>.05)/N # Porcentaje de casos donde no refutariamos H0 (media en 0)
[1] 0.721 # Alto
El 72,10% de los casos no refutaría la hipotesis y por tanto seria errónea la deducción (la media realmente es diferente de 0). Como consecuencia, queda demostrado que el T-Test no es tan efectivo cuando trabajamos con datos de distribución diferente a la normal y en consecuencia es mejor utilizar otras alternativas.
Wilcoxon Test: Datos con distribución diferente a la normal
Si pensamos que la distribución no es normal (ver secciones de identificación de distribuciones y pruebas de normales) es mejor utilizar pruebas no paramétricas como Wilcoxon’s U test. Por ejemplo, un estudio sobre la duración de las llamadas telefónica arroja los siguientes minutos:
> x = c(12.8,3.5,2.9,9.4,8.7,.7,.2,2.8,1.9,2.8,3.1,15.8)
Si visualizamos los datos (comprobar también utilizando QQ plot o con el test Shapiro-Wilk):
> plot.hist.and.box(x) > stem(x) The decimal point is 1 digit(s) to the right of the | 0 | 01233334 0 | 99 1 | 3 1 | 6
Se observa una distribución con una cola larga lejos de ser una normal, por tanto el T-Test queda descartado y en su lugar utilizaremos el Wilcox test.
Como hipotesis nula a testear diremos que la mediana es 5, mientras que la alternativa es que la mediana es superior (descartamos que sea inferior => hipótesis compleja):
> wilcox.test(x,mu=5,alt="greater") Wilcoxon signed rank test with continuity correction data: x V = 39, p-value = 0.5156 alternative hypothesis: true location is greater than 5
El p-value es 0.5156 > 0.05 y por tanto no podemos refutar la hipótesis nula.
ANOVA / Kruskal
Los tests T y Wilcoxon nos permite analizar una o dos muestras de datos, pero si queremos evaluar más de 2 conjuntos necesitamos realizar un análisis de varianza (ANOVA). Veamos un ejemplo: Hemos realizado un examen a 24 personas y las correcciones se han repartido a 3 profesores. Las notas van de 1 a 5, siendo esta última la mejor:
x = c(4,3,4,5,2,3,4,5) y = c(4,4,5,5,4,5,4,4) z = c(3,4,2,4,5,5,4,4) scores = data.frame(x,y,z) boxplot(scores)
Queremos comprobar que todos los profesores han evaluado bajo los mismos criterios y que no se han producido grandes diferencias. Para ello realizaremos un análisis de las varianzas. Para esto, necesitamos que los datos estén en formato pila:
> scores.s = stack(scores) > scores .s values ind 1 4 x 2 3 x 3 4 x ... 9 4 y 10 4 y ... 17 3 z 18 4 z ...
Ejecutamos el test Oneway (suponemos que la distribución es normal) y comprobamos si todos los profesores tienen la misma media (hipótesis nula), asumiendo que la varianza será equivalente para todos:
> attach(scores.s) > oneway.test(values ~ ind, var.equal=T) One-way analysis of means data: values and ind F = 1.1308, num df = 2, denom df = 21, p-value = 0.3417
El p-value es de 0.3417 > 0.05 y por tanto no podemos refutar la hipotesis.
En caso de que la distribución no sea normal, utilizaríamos el Test Kruskal:
> kruskal.test(values ~ ind) Kruskal-Wallis rank sum test data: values by ind Kruskal-Wallis chi-squared = 1.9387, df = 2, p-value = 0.3793
Prueba sobre la varianza
Comparar varianzas de 2 conjuntos de datos:
- Paramétrico (distribución normal)
?fisher.test
- No paramétrico
?mood.test ?ansari.test
Comparar varianzas de más de 2 conjuntos de datos:
- Paramétrico (distribución normal)
?bartlett.test
- No paramétrico
?fligner.test
Prueba de independencia
Nos puede interesar validar si 2 muestras/variables son independientes entre si o si sus fluctuaciones se encuentran relacionadas, con este objetivo podemos hacer uso de las pruebas de independencia.
Coeficiente de Pearson: Datos con distribución normal
El coeficiente de correlación Pearson asume que la distribución de probabilidad de los datos es la normal (estadística paramétrica), si esta condición se cumple entonces las estimaciones son más precisas que otros métodos (ver siguiente sección) y viceversa.
La correlación siempre devuelve un valor entre -1 y 1, donde 0 significa que no hay ningún tipo de correlación, -1 que la correlación es absoluta pero en dirección opuesta y 1 la correlación es completa.
El cálculo del coeficiente de correlación de Pearson se lleva a cabo de la siguiente forma:
> x = rnorm(1000, 50, 10) # Partimos de datos aleatorios normales > y = rexp(1000, 50) # Y datos aleatorios exponenciales > cor(x, y) [1] -0.007434471 > cor.test(x, y) Pearson's product-moment correlation data: x and y t = -0.2349, df = 998, p-value = 0.8144 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.06939557 0.05458377 sample estimates: cor -0.007434471
En el ejemplo anterior, el coeficiente es de -0.007 y las probabilidades de que las variables sean independiente (hipotesis nula) son del 81,44% (p-value).
Métodos no paramétricos: Datos con distribución diferente a la normal
Cuando no estamos seguros de la distribución de probabilidad de los datos o se trata de una diferente a la normal, la mejor opción es utilizar métodos no paramétricos como el coeficiente de correlación de Spearman (rho), Kendall tau o chi cuadrado. Por ejemplo, coeficiente Spearman:
> cor(x, y, method="spearman") [1] 0.003597064 > cor.test(x, y, method="spearman") Spearman's rank correlation rho data: x and y S = 166066990, p-value = 0.9095 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 0.003597064
Y coeficiente Kendall:
> cor(x, y, method="kendall") [1] 0.002506507 > cor.test(x, y, method="kendall") Kendall's rank correlation tau data: x and y z = 0.1187, p-value = 0.9055 alternative hypothesis: true tau is not equal to 0 sample estimates: tau 0.002506507
Al igual que con el coeficiente de Pearson, cabe recordar que la correlación siempre devuelve un valor entre -1 y 1, donde 0 significa que no hay ningún tipo de correlación, -1 que la correlación es absoluta pero en dirección opuesta y 1 la correlación es completa.
En los ejemplos anteriores los coeficientes son de 0.004 (rho) y 0.003 (tau) respectivamente, y las probabilidades de que las variables sean independiente (coeficiente = 0) son del 90,95% y 90,55% (p-value).
Veamos un ejemplo rápido donde las variables no son independientes:
> cor.test(c(1, 2, 3, 4, 5), c(2, 3, 4, 5, 6), method="spearman") Spearman's rank correlation rho data: c(1, 2, 3, 4, 5) and c(2, 3, 4, 5, 6) S = 0, p-value = 0.01667 alternative hypothesis: true rho is not equal to 0 sample estimates: rho 1
Dado que los elementos del segundo vector son iguales que los del primero pero sumando 1, no existe independencia, por tanto rho es 1 y la probabilidad de que sean independientes ínfima: 1,67%.
Hagamos ahora el mismo test de independencia pero utilizando el método de chi cuadrado:
> chisq.test(data.frame(x, y)) Pearson's Chi-squared test data: data.frame(x, y) X-squared = 20.5816, df = 999, p-value = 1
El test nos muestra una probabilidad del 100% (p-value) de que la hipotesis nula sea cierta: las variables son independientes.
Análisis de regresiones
Regresiones lineales simples
Si entendemos que X es un valor que controlamos (p.ej. tabletas de chocolate que comemos) e Y un valor que medimos (p.ej. grado de felicidad), podemos utilizar una regresión lineal para predecir las Y para un valor desconocido de X (p.ej. cual sera mi grado de felicidad si como 20 tabletas de chocolate).
En primer lugar, a partir de un conjunto de X e Y conocidos inferimos los coeficientes de nuestro modelo: Y = B1*x + B0
> x = 1:100 > y = rexp(100) > lm.result = lm(y ~ x) # response ~ predictor > lm.result Call: lm(formula = y ~ x) Coefficients: (Intercept) x 1.130804 -0.003217
El modelo generado corresponde con la siguiente fórmula: y = -0.003217 * x + 1.130803
> coef(lm.result) (Intercept) x.chocolate 0.4395754 0.1466509
Para cada Y conocido (los que existen en el conjunto de datos inicial), podríamos coger sus correspondientes X y utilizar la formula deducida y obtener la Y’ que nos daría. Si restamos las Y reales con la Y’ calculadas tenemos los residuales:
> lm.res = resid(lm.result)
> summary(lm.res)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-5.449e-01 -2.725e-01 4.904e-02 3.342e-18 2.616e-01 4.947e-01
De los datos anteriores se desprende que la diferencia media es muy pequeña.
Para analizar si el modelo encaja con los datos, podemos visualizar ambos:
> plot(x,y) > abline(lm(y ~ x))
Y analizar los residuales mediante 4 gráficas:
> par(mfrow=c(2,2)) ; plot(lm.result)
- Residuals vs. Fitted: Muestra los valores usados de Y junto a los residuales (diferencia con Y’). Observando la distribución alrededor de la linea y=0 podríamos identificar tendencias que nos hagan pensar que existe una correlación entre los datos.
- Normal qqplot: Los residuales son normales si la gráfica se aproxima a una linea recta.
- Scale-Location: Muestra la raiz cuadrada de los residuales estandarizados. Los puntos más altos corresponden con los residuales más largos.
- Cook’s distance: Identifica los puntos que tienen mucha influencia sobre la linea de regresión.
Veamos un ejemplo donde el modelo encaja significativamente con los datos:
> xx = 1:100 > yy = xx*2 + 100 > plot(xx,yy) > abline(lm(yy ~ xx)) > par(mfrow=c(2,2));plot(lm(yy ~ xx))
Como hemos visto en el ejemplo inicial, el modelo generado correspondia con esta formula: y = -0.003217 * x + 1.130803
Los parámetros calculados son inferencias sobre los cuales se realizan un test de verificación y por tanto tienen asociado un p-value :
> coefficients(summary(lm.result))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.130804321 0.191432344 5.9070703 5.041288e-08
x -0.003217217 0.003291033 -0.9775707 3.306937e-01
Segundo ejemplo donde el modelo se ajusta más a la realidad:
> coefficients(summary(lm(yy~xx)))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100 4.002787e-15 2.498259e+16 0
xx 2 6.881441e-17 2.906368e+16 0
Con la función "predict" podemos obtener las Y para X determinadas:
> predict(lm.result,data.frame(x=c(1,2)))
1 2
1.127587 1.124370
Y también podemos generar para cada predicción el intervalo de confianza:
> predict(lm.result,data.frame(x=c(1,2)), interval="confidence")
fit lwr upr
1 1.127587 0.7533518 1.501822
2 1.124370 0.7557617 1.492978
Esto puede ser de utilidad observarlo gráficamente:
plot(x,y) abline(lm.result, col="red") curve(predict(lm.result,data.frame(x=x), level=.95, interval="confidence")[,3],add=T) curve(predict(lm.result,data.frame(x=x), level=.95, interval="confidence")[,2],add=T)
Regresiones lineales múltiples
Para predecir el precio de una vivienda, existen múltiples factores que influyen: número de habitaciones, barrio, etc. Por tanto, el modelo será diferente al visto en el anterior apartado:
Y = B0 + B1*X1 + B2*X2 + …
Si partimos de los siguientes datos (extraído del tutorial simpleR):
"homeprice" <-
structure(list(list = c(80, 151.4, 310, 295, 339, 337.5, 228.7,
245, 339, 43, 279, 599, 119, 289, 249, 178, 159, 289, 488, 376,
249, 275, 275, 459, 219, 359, 379, 189, 173), sale = c(117.7,
151, 300, 275, 340, 337.5, 215, 239, 345, 48, 262.5, 613, 119,
305, 249, 170, 153, 291, 450, 370, 245, 275, 272.5, 459, 230,
360, 370, 185, 185), full = c(1, 1, 2, 2, 2, 1, 3, 1, 1, 1, 1,
3, 2, 2, 1, 2, 1, 2, 3, 3, 1, 2, 2, 3, 1, 2, 2, 1, 1), half = c(0,
0, 1, 1, 0, 1, 0, 1, 2, 0, 2, 2, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0,
1, 1, 1, 0, 1, 1, 0), bedrooms = c(3, 4, 4, 4, 3, 4, 3, 3, 3,
1, 3, 4, 3, 3, 2, 3, 2, 3, 3, 3, 3, 4, 2, 5, 3, 3, 4, 4, 3),
rooms = c(6, 7, 9, 8, 7, 8, 7, 7, 7, 3, 8, 10, 7, 6, 4, 6,
7, 7, 7, 7, 8, 8, 6, 11, 7, 8, 8, 8, 7), neighborhood = c(1,
1, 2, 2, 3, 2, 1, 1, 2, 1, 1, 3, 1, 2, 2, 1, 1, 2, 3, 3,
2, 2, 2, 3, 2, 3, 3, 1, 2)), .Names = c("list", "sale", "full",
"half", "bedrooms", "rooms", "neighborhood"), row.names = c(NA,
-29L), class = "data.frame")
> homeprice
list sale full half bedrooms rooms neighborhood
1 80.0 117.7 1 0 3 6 1
2 151.4 151.0 1 0 4 7 1
3 310.0 300.0 2 1 4 9 2
4 295.0 275.0 2 1 4 8 2
5 339.0 340.0 2 0 3 7 3
...
Podemos observar el precio en función de las habitaciones de la vivienda y agrupados según barriadas:
> library(lattice);
> panel.lm = function(x,y) {
panel.xyplot(x,y)
panel.abline(lm(y~x))}
> attach(homeprice)
> xyplot(sale ~ rooms | neighborhood, panel= panel.lm,layout=c(3,1))
Si tenemos en cuenta los dormitorios y los baños completamente equipados:
> summary(lm(sale ~ bedrooms + full + neighborhood))
Call:
lm(formula = sale ~ bedrooms + full + neighborhood)
Residuals:
Min 1Q Median 3Q Max
-72.85 -43.11 -15.36 22.89 165.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -67.89 47.58 -1.427 0.1660
bedrooms 31.74 14.77 2.149 0.0415 *
full 28.51 18.19 1.567 0.1297
neighborhood 101.00 17.69 5.709 6.04e-06 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Residual standard error: 58.71 on 25 degrees of freedom
Multiple R-squared: 0.7764, Adjusted R-squared: 0.7496
F-statistic: 28.94 on 3 and 25 DF, p-value: 2.694e-08
Vemos que un dormitorio adicional supone 31 mil euros más, un baño completo 28 mil euros y estar situado en un barrio mejor 101 mil euros adicionales.
El precio de una vivienda sin habitaciones correspondería a (según la barriada):
> -67.89 + 101*(1:3) [1] 33.11 134.11 235.11
Las variables pueden ser combinadas de la siguientes forma:
|
Fórmula |
Interpretación |
| Y ∼ X1 |
Y modelado por X1 (visto en el apartado anterior) |
| Y ∼ X1 + X2 | Y modelado por X1 y X2 |
| Y ∼ X1 * X2 | Y modelado por X1, X2 y X1*X2 |
| (Y ∼ (X1 + X2)ˆ2) | Interacción en 2 sentidos |
| Y ∼ X1+ I((X2^2) | Y modelado por X1 y X2 al cuadrado |
| Y ∼ X1 | X2 | Y modelado por X1 condicionado en X2 (p.ej. precio piso (Y) según habitaciones (X1) pero la relación depende del barrio (X2 = Eixample, Gracia…), un modelo diferente por cada barrio) |
Estandarizar y escalar
Si queremos comparar dos variables aleatorias con escalas diferentes (media y desviación estándar diferentes) podemos estandarizarlas utilizando la función ‘scale’ que nos devolverá los datos con media 0 y desviación estándar 1:
> x = rnorm(5, mean=100, sd=16) > z = scale(x) > mean(z) # Casi 0 [1] 1.887162e-16 > sd(z) # 1 [1] 1
Técnicas para identificar fraude
Si necesitamos estudiar unos datos para identificar posibles fraudes (p.ej. blanqueo de capitales), podemos utilizar diversas técnicas que nos ayudaran a centrar nuestro trabajo:
- Ordenación de los datos y Agings: dado un conjunto de transacciones de un periodo, si se identifican fechas anteriores o futuras pueden ser indicios de potenciales fraudes.
- Clasificación: Si por ejemplo el 95% de las transacciones se encuentran destinadas a la misma persona, esto puede ser un indicio de fraude.
- Estadística:
A partir de un conjunto de valores, como por ejemplo:
> transacciones = round(rexp(100)*1000, 2)
Podemos utilizar como regla que un valor extremo (outlier) será todo aquel valor superior/inferior a la media mas/menos la desviación estándar multiplicada por 3. Ejemplo:
> transacciones[transacciones < mean(transacciones)-sd(transacciones)*3] numeric(0) > transacciones[transacciones > mean(transacciones)+sd(transacciones)*3] [1] 4621.22 5500.91 5049.35
Estas transacciones son potenciales anomalías que deben ser analizadas.
Otra forma de realizar la misma identificación de forma gráfica seria utilizando boxplot:
> boxplot(transacciones)
Existen otros métodos para la identificación de valores extremos como por ejemplo: Rosner test, Dixon test y Grubbs test, cuyas implementaciones para R pueden encontrarse en librerías disponibles en la red.
- Ley de Benford:
Es más probable que el dígito 1 sea el primero de una cantidad que el 9. Si existe algún tipo de fraude, es probable que desvirtúe esta distribución.
Definimos en R las siguientes funciones que definen una distribución discreta que sigue la ley de Benford:
# Función de densidad
dbenford = function(x){
log10(1 + 1/x)
}
# Función de distribución
pbenford = function(q){
cumprobs = cumsum(dbenford(1:9))
return(cumprobs[q])
}
# Función de los cuantiles
qbenford = function(p){
cumprobs = cumsum(dbenford(1:9))
cumprobs[9] = 1 # To fix a rounding error
quantiles = sapply(p, function(x) {10 - sum(cumprobs >= x)})
return(quantiles)
}
# Generación aleatoria de número según la distribución
rbenford = function(n){
sample(1:9, size = n, replace = TRUE, prob = dbenford(1:9))
}
También creamos una función que nos permita visualizar gráficamente un conjunto de datos y compararlos contra la distribución basada en la ley de Benford:
BenfordObsExp = function(x){
data = substitute(x)
n = length(x)
# Peel off the first digit
x = as.numeric(substring(formatC(x, format = 'e'), 1, 1))
obsFreq = tabulate(x, nbins = 9)
benFreq = n * dbenford(1:9)
plot(1:9, benFreq, xlim = c(1, 9), ylim = c(0, max(obsFreq, benFreq)), type = 'l',
main = paste(data, ": observed (red) and expected (Benford's Law)"),
xlab = "Digit", ylab = "Frequency")
axis(1, at = 1:9)
points(1:9, obsFreq, col = "red", pch = 16)
}
Visualizamos las transacciones para identificar anomalías:
> BenfordObsExp(transacciones)
Programación
Hasta el momento hemos visto como interactuar con la shell de R, pero también podemos generar ficheros con comandos para que sean ejecutados por R. Por ejemplo, creamos el fichero “test.R”:
x = runif(100)
y = runif(100)
z = x + y
cat("Imprimimos el listado generado:\n")
z
cat("Fin")
Y lo ejecutamos mediante:
$ R --slave --vanilla -f test.R Imprimimos el listado generado: [1] 0.8127675 0.2205953 0.8711474 1.1461206 1.5373148 0.8068898 1.1030845 [8] 0.3391666 1.2360769 1.3011322 1.5047748 0.5618684 0.9276672 0.8581834 [15] 1.6176204 0.8798410 0.8575242 0.8529046 1.7232819 0.7442834 1.1624653 [22] 1.5103152 1.2968150 1.2777382 0.6855236 0.8132157 1.3333253 1.2821548 [29] 1.3576682 1.1655971 0.3267455 1.1248876 0.2601012 0.3260605 0.7609038 [36] 0.5666418 0.8664088 1.2207019 1.1709423 1.4731985 1.4740750 1.5157082 [43] 1.2244977 1.2633403 1.6129252 1.6111450 0.8875049 0.9616655 1.0212954 [50] 0.3086974 0.8500047 0.9417143 1.7620811 0.9862873 1.2310904 0.7335394 [57] 1.2586577 0.9282816 0.6394448 0.9480084 1.0489753 0.1546561 0.8817316 [64] 0.8384202 0.9665596 1.9363211 1.1402878 0.5812365 0.5223378 1.5062012 [71] 0.4690346 1.1312263 1.1484181 1.5491576 1.9231288 0.8261800 1.6588865 [78] 0.2999772 0.5317558 1.0639816 1.6124205 1.0747384 0.8697748 0.7791368 [85] 1.0116984 0.7644074 1.5599921 0.8562209 1.5363938 1.1215349 1.1354001 [92] 0.9089423 1.7165493 0.8804664 0.6908375 0.1000780 1.3157824 1.3394445 [99] 1.1517646 0.6332460 Fin
También podríamos cargar y ejecutar el código del fichero "test.R" desde la terminal de R:
$ R -q
> source("test.R")
Codigo de ejemplo
Imprimimos el listado generado:
Fin
En este caso el comportamiento es ligeramente diferente, la variable Z no se imprime por pantalla pero si queda en memoria para que podamos trabajar con ella.
Estructuras de control
Condiciones:
if(...) {
...
} else {
...
}
x = rnorm(100) y = ifelse(x>0, 1, -1) z = ifelse(x>0, 1, ifelse(x<0, -1, 0))
Bucles:
for (i in 1:10) {
...
if(...) { next }
...
if(...) { break }
...
}
while(...) {
...
}
Definición de funciones
Podemos crear funciones donde el valor retornado corresponderá a la última expresión (como en Ruby, aunque también podemos utilizar return):
f = function(x, z=1) {
x^2 + x + z
}
Adicionalmente, la anterior función define 2 argumentos, el segundo puede no ser especificado y toma por defecto 1. A la hora de utilizar la función, si especificamos los nombres de los argumentos, estos pueden ir en
cualquier orden:
f(2) # x=2, z=1 f(3, 2) # x=3, z=2 f(z=2, x=3) #x=3, z=2
Cuando llamamos a la función, el valor retornado se imprime por pantalla, si queremos que lo devuelva igualmente pero que no se imprima podemos utilizar invisble:
f = function(x, z=1) {
invisible(x^2 + x + z)
}
r = f(2)
r
Podemos tener un número indefinido de argumentos si utilizamos "…":
f = function(x, ...) {
plot(x, ...)
}
f = function (...) {
l = list(...) # Put the arguments in a (named) list
for (i in seq(along=l)) {
cat("Argument name:", names(l)[i], "Value:", l[[i]], "\n")
}
}
Es importante tener en cuenta que una función siempre guarda toda la información en su ambito local, no puede modificar valores de variables globales.
Si queremos ver el código de una función, simplemente la llamamos sin parentesis ni argumentos:
> f
function(x, z=1) {
invisible(x^2 + x + z)
}
Depurar / Debug
Paso a paso:
f = function () {
x = rnorm(10)
y = rnorm(11)
x + y
}
debug(f)
f()
# n o ENTER para continuar al siguiente paso
# c para continuar hasta el siguiente cambio de contexto
# where llamadas
# Q termina
# Cualquier otra expresión será evaluada desde el entorno donde se llamó a f(), es decir, no podemos ver las variables locales
undebug(f)
Breakpoints:
f = function () {
x = rnorm(10)
y = rnorm(11)
browser()
x + y
}
f()
# Mismas posibilidades que con debug, excepto que aquí si podemos acceder a las variables locales (p.ej. x e y
Optimización / Profiling
# Tiempo de ejecución de una función > system.time(runif(100)) user system elapsed 0 0 0
Consultar ayuda:
?RProf
tengo una tarea de clase y no la puedo encontrar, sólo requiero una explicación sobre índice de anomalías
MUY Buen aporte… gracias!!
Yo como interfaces graficas para R recomiendo siempre :
-rcomander
-rstudio (muy buena: http://www.rstudio.com/ide/)
Excelente trabajo, felicito la iniciativa de difundir el conocimiento en la red
como puedo extraer los dígitos centrales de un numero entero, estoy en un casi de simulación y me falta esa parte.
Increíble. Muy bien hecho.