Primitivas en GNU R
Índice
- 1. Instalación
- 2. Básicas
- 3. Tipos de Datos
- 4. Vectores
- 4.1. Básico
- 4.2. Calculo vectorizado. Equivale a un bucle implicito
- 4.3. Generacion de secuencias regulares
- 4.4. Ordenacion de vectores
- 4.5. Operadores comparativos: <, <=, >, >=,
=, !
- 4.6. Operadores logicos: &,|, !, xor()
- 4.7. Identificacion y sustitucion de valores perdidos. Funciones is.na() y which()
- 4.8. Indeterminaciones e infinito
- 4.9. Manipulacion de vectores de caracteres
- 4.10. Indexacion de vectores mediante variables de caracteres. Funcion names()
- 5. Factores
- 6. Matrices
- 7. Listas
- 8. Hojas de datos (data frames)
- 9. Cargar datos nativos de R
- 10. Paquetes R
- 11. Escritura de datos en un archivo
- 12. Lectura de datos desde un archivo
- 13. Gráficos
- 13.1. Gestion de ventanas graficas
- 13.2. Funcion plot(). Es la función gráfica basica.
- 13.3. Plot a Png
- 13.4. Funcicion coplot(). Matriz de graficos de dispersion separados segun factores
- 13.5. Funcion pairs(). Ejemplo con datos simulados de los "Big Five" de personalidad
- 13.6. Funcion hist(). Histogramas
- 13.7. Funcion boxplot()
- 13.8. Funciones qqnorm() y qline
- 13.9. Catalogo de nombres de color
- 13.10. Parametros gráficos. Funcion par()
- 13.11. Representacion de funciones
- 14. Muestreo y probabilidad
- 15. Definicion de funciones
- 16. Descriptivos
- 17. Distribuciones
- 18. Test estadísticos
- 18.1. Distribución Binomial
- 18.2. Test Binomial
- 18.3. Distribución F
- 18.4. Bondad de ajuste multinomial
- 18.5. Test de Chi
- 18.6. Test Wilcox
- 18.7. Ecuación de la Regresión Logística
- 18.8. Test de significado para la regresión logística
- 18.9. Mann-Whitney test: comparar si hay diferencias entre medianas muestras independientes
- 19. Ecuación de Regresión Logística
- 20. Análisis de Varianza
- 21. Análisis de Clúster Jerárquico
- 22. Licencia
1 Instalación
1.1 Debian/Ubuntu
En Debian/Ubuntu encontraremos paquetes R empaquetados con el sufijo r-base o r-cran, buscando con apt-cache search. Es posible usar metapaquetes dependiendo de la versión como r-recommended
1.2 R Commander
Para las personas que vienen adictas a SPSS y similares, puede ser necesario instalarse R para empezar a cuestionar utilizarlo. Por favor, no olviden que la estadística es programación no manejo de un interfaz. Para instalarlo también es posible usar paquetería Debian como
$ sudo apt-get install r-cran-rcmdr
Y después ejecutar:
$ R > library("Rcmdr")
1.3 Paquetería R
La ventaja de R frente a otras opciones para hacer estadística es el enorme número de paquetes. R tiene facilidades propias si su sistema operativo no tiene soporte, o no le apetece utilizarlo.
> install.packages("snow")
2 Básicas
2.1 Configuración del Entorno
$ R # Para entrar en el entorno R una vez instalado # es necesario ejecutar la R desde la consola getwd() # Devuelve el # directorio de trabajo setwd("/home/davidam/public_html/software/R") # Cambia el directorio # de trabajo (debe # existir previamente) list.files() # Lista el contenido # del directorio # de trabajo. Tambien ls() source("/home/davidam/public_html/software/R/ejercicios-descriptivos.R") q() # Salir del entorno
2.2 Ayuda on-line en R
help (aov) ?mean help.search('mean') ??'standard deviation' help.start() # abre una ventana del navegador example (mean) # ejemplos del uso de mean
2.3 Operaciones Aritméticas
5+7 # se evalua y el resultado se imprime # por pantalla pero no se guarda 5**7 5-7 5/7 5^7 # elevado a 51 %% 7 # modulo (resto de una division entera) 51 %/% 7 # division entera
2.4 Asignacion de valores a variables
x <- 1 # asigna el valor 1 a la variable x x # lo escribe por pantalla y <- x^2 x = 1 # una alternativa 1 -> x # otra alternativa
2.5 Funciones
q # Muestra informacion de la funcion ?q # Algunas funciones elementales: log(), log10(), # exp(), sqrt(), sin(), cos(), tan() # Mas funciones basicas: c(), max(), min(), # pmax(), pmin(), range(), length(), sort() # sum(), rowSums(), prod(), mean(), colMeans(), # var(), cumsum(), etc. c(2,69,5,47,3647) sum(1,3,5,7) raiz64 <- sqrt(64) # el resultado de una funcion # puede guardarse en una variable print(x) print(pi,digits=16); print(pi,digits=5)
2.6 Listar y Eliminar Objetos del Espacio de Trabajo
ls() rm(x,raiz64) ls() rm(list=ls()) ls()
3 Tipos de Datos
3.1 Booleanos
verdad <- T mentira <- F edad <- c(27,30,29); is.numeric(edad) is.integer(edad) is.double(edad) # tambien existen is.complex(), # is.logical(), is.character()... typeof(edad) # devuelve el tipo del objeto # introducido como argumento
3.2 Arrays
calle <- c("c/ rue del percebe","28013");
3.3 Conversiones de Tipo
as.character(verdad) as.numeric(verdad) # doble precision por defecto y <- c(0,1,27,14.9,NA,-3.2);y # ojo, trunca en vez de redondear as.logical(y) as.character(y) as.integer(y) y<-c("hola", "adios","true","TRUE", "T", "F","129",NA ,"6");y as.logical(y) as.numeric(y)
4 Vectores
4.1 Básico
x <- c(1,4,6,3,7); x x[3] # devuelve el valor guardado en 3a posicion x[-4] # devuelve todo excepto el valor # guardado en 4a posicion x[1] <- 35 # guarda el valor 35 en la primer posicion length(x) # devuelve la longitud de x x[11] <- 99; x
4.2 Calculo vectorizado. Equivale a un bucle implicito
x <- c(1,2,3,4,5) y <- c(5,4,3,2,1) x+y y <- c(1,2,3);y suma <- x+y; suma # recicla el vector mas corto suma2 <- x-1 # con escalares parece mas natural x <- c(0.0001, 0.001, 0.01, 0.1, 1) log10(x) # muchas funciones tambien se # aplican vectorizadamente
4.3 Generacion de secuencias regulares
1:10 1:10-1 # el operador ':' tiene maxima prioridad 1:{10-1} #funcion seq(from, to, by, length) x <- seq (-3,3,0.5);x # valores desde -3 hasta 3 a intervalos de 0.5 seq (-3,3,length=10) # 10 valores equiespaciados desde -3 hasta 3 seq (-3, by=0.5,length=10) # 10 valores a intervalos de 0.5 desde -3 seq(along=x) # Genera la secuencia 1, 2, 3,..., length(x) #funcion rep(x, times, each) x <- 1:4 rep(x,times=3) # repite el contenido de x tres veces rep(x,each=3) # repite cada elemento de x tres veces rep(x,times=2,each=3) rep(c(0,1),times=c(4,3)) # 'times' puede ser un vector. Aqui ya no cabe 'each' x <- seq(19,6,-3);x rep(x,1:length(x)) # 'times' puede estar implicito times={1:length(x)}
4.4 Ordenacion de vectores
x <- c(20,80,30,50,0) order (x, decreasing=F) # devuelve las posiciones del vector ordenadas segun su contenido sort (x, decreasing=F) # devuelve el contenido del vector ordenado rank(x) # devuelve el orden de cada posicion segun su contenido min(x) which.min(x) # equivale a which(x == min(x)) x <- c(1,1,3:1,1:4,3); y <- c(0,9:1) x.ord <- order(x,y) # ordena los valores de x y, en caso de empate, utiliza los valores de y en las posiciones correspondientes. Ambos, x e y deben tener la misma longitud x.ord;x;y
4.5 Operadores comparativos: <, <=, >, >=, =, !
x <- c(1:10); x valor.verdad <- x>5; valor.verdad x <- 1:10 x[x >= 5] <- 20; x # condicional implicito x[x == 1] <- 25; x x <- 1:3; y <- 1:3 x == y # equivalente a identical(x, y) o all.equal(x, y, tol=0) y <- y+0.001; y identical(x, y) all.equal(x, y, tol=0.001)
4.6 Operadores logicos: &,|, !, xor()
a<-c(TRUE,TRUE,FALSE,FALSE) b<-c(F, T, F, T) a&b a|b !(a&b) (!a)|(!b)
4.7 Identificacion y sustitucion de valores perdidos. Funciones is.na() y which()
rm(list=ls()) x <- c(3,6,4,2,8,9) print (x); length(x) x[8:10] <- 3;x is.na(x) !is.na(x) which(is.na(x)) x[is.na(x)]<-999;x # codifica como 999 los valores perdidos x==NA # la expresion logica x == NA sa un resultado muy distinto de is.na(x)
Entonces si tenemos data$q005 como un frame con algunos datos perdidos y queremos solo los datos no perdidos podemos hacer
data5 <- data$q_005[!is.na(data$q_005)]
4.8 Indeterminaciones e infinito
x <- c(0,7,8); x/x 1/x -1/x is.nan(x) is.nan(x/x) is.nan(1/x) # Hay que tener cuidado con los NaN porque cualquier operacion con un NaN resulta en un NaN
4.9 Manipulacion de vectores de caracteres
# Concatena objetos en un vector de caracteres. Funcion paste(..., sep = " ", collapse = NULL) juntar <- paste("Una ", "frase ", "cualquiera", collapse ="");juntar v1<-c("A","B") v2<-2:3 codigos <- paste(v1,v2, sep = "");print(codigos) codigos <- paste(v1,v2, sep = ".");print(codigos) x <- paste(LETTERS[1:5]);x x <- paste(LETTERS[1:5], collapse="");x # Concatena e imprime. Funcion cat(... , file = "", sep = " ", fill = F, labels = NULL, append = F) verano <- month.abb[7:9]; verano cat(verano) # el resultado no puede guardarse en ua variable cat(verano,"\n") cat(verano, sep=',') # concatena separando con comas e imprime cat(verano, sep=';', fill=3) cat("Estaciones:","\t","Moncloa","\n", "\t", "\t","Aluche","\n") #Todo junto: cat() y paste() x<-2/3; cat(paste("resultado", signif(x,2), sep=" : " ),"\n")
4.10 Indexacion de vectores mediante variables de caracteres. Funcion names()
edad <- c(12,22,15,16,10) names(edad) # por defecto no se asignan nombres names(edad) <- paste("suj#", sep = "",length(edad):1); edad names(edad) edad["suj#2"] # devuelve la edad de Sujeto #2 (almacenado en la penultima posicion) edad[length(edad)-1] names(edad) <- NULL; edad # elimina los nombre asignados
5 Factores
Son vectores para datos categoricos. Permiten prescindir de la codificacion numerica y referirse a los niveles mediante nombres
estudios <- c(1,3,1,1,3,4,3) factor(estudios) nivel.estudios <- factor(estudios, levels=1:4,labels=c("primarios", "secundarios", "superiores", "doctorado"), ordered=T) nivel.estudios factor(c(5,2,2,4,5,4,3,3,1), 2:5, exclude=4) # los valores 1 y 4 se consideran valores perdidos # OJO. Internamente R siempre se asigna 1 al primer nivel del factor, 2 al segundo etc. levels(nivel.estudios) # levels() extrae los posibles niveles de un factor as.numeric(nivel.estudios) # recodifica el factor numericamente codigo.postal <- factor(c('28011', '28044', '28011','28013', '28013','28023')) equipo.futbol <- factor(c('VAL', 'VAL', 'FCB','VAL', 'FCB','ATM')) equipo.futbol # ordena los factores alfabeticamente as.numeric(equipo.futbol) as.character(equipo.futbol) # convierte en cadenas de caracters no en factor de nuevo
6 Matrices
Suponen la generalizacion de vectores a 2D. Todos los elementos de una matriz deben ser del mismo tipo. Una convencion relativamente extendida es cmenzar con mayuscula el nombre de una matriz
6.1 Función matrix()
X <- 1:12 dim(X) <- c(4,3); X # los elementos se organizan por columnas matrix(1:12,nrow=3,ncol=4,byrow=T) matrix(1:12,nrow=4,byrow=T) Mi.matriz <- matrix(1:12,3,4,F); Mi.matriz tamano <- dim(Mi.matriz) # otro uso de dim(); devuelve las dimensiones de la matriz rownames(Mi.matriz) <- LETTERS[1:tamano[1]]; Mi.matriz colnames(Mi.matriz) <- paste("Var",1:tamano[2], sep=""); Mi.matriz dimnames(Mi.matriz)
6.2 Concatenación de matrices
X1 <- c(3,7,5) X2 <- c(8,3,1) Xx <- cbind(X1,X2);Xx Yy <- rbind(X2,X1);Yy Zz <- cbind(X1,Xx);Zz
6.3 Indexacion de matrices
X <- matrix(c(1,4,12,15),2,2);X X[1,2] # elemento guardado en la 1a fila, 2da columna X[1, ] # todos los elementos de la primera fila X[ ,2] # todos los elementos de la segunda columna X[3] # para pensar un poco... Mi.matriz['B',] #tambien se pueden usar los nombres (se han asignado)
6.4 Operaciones con matrices
X <- matrix(c(1,4,12,15),2,2); X Y <- matrix(1:4,2,2); Y X+Y X-Y X%*%Y # producto matricial X*Y # producto elemento por elemento t(X) # traspuesta det(X) # determinante X.inv <- solve (X) # inversa de X (siempre que X sea cuadrada no singular, claro) X.inv X%*%X.inv # Comprueba el resultado. Ojo a los errores de redondeo # En general, solve(a,b) es una funcion que resuelve a %*% x = b para x, donde b puede ser un vector o una matriz. Si no se explicta se asume que es la matriz identidad y la funcion devuelve la inversa de a A <- matrix(c(1,4,12,15),2,2); A B <- matrix(c(5,2),2,1);B X <- solve(A,B); X A%*%X # comprueba que el resultado es efectivamente B
6.5 Algunas funciones que operan sobre filas o columnas completas
X <- matrix(c(1,4,12,15),2,2); X rowSums(X) colSums(X) rowMeans(X) colMeans(X)
6.6 Apply
apply(X,1,sum) # suma por filas; equivale a rowSums(X) apply(X,2,mean) # media por columnas; equivale a colMeans(X) apply(X,1,sd)
6.7 Matrices multidimensionales (arrays)
A <- array(1:24, c(3, 4, 2)); A dimnames(A) <- list(c("fila1", "fila2", "fila3"), c("col1", "col2", "col3", "col4"), c("capa 1", "capa 2"));A
7 Listas
Son una especie de contenedores generales donde pueden mezclarse todo tipo de componentes (objetos de cualquier tipo y cualquier longitud)
Son unos objetos poco estructurados y, por tanto, muy flexibles.
Muchas funciones nativas de R devuelven el resultado en forma de lista
rm(list=ls()) mis.num <- seq(1.0, 2.0, 0.1); mis.num2 <- 2:4 Mi.matriz <- matrix(1:12,3,4); Mi.matriz mis.caracteres <- paste(LETTERS[1:5]);mis.caracteres mis.logicos <- mis.num > 1.65; mis.logicos lista1 <- list(mis.num, mis.num2, Mi.matriz, mis.caracteres,mis.logicos) lista1
7.1 Indexacion, nombres y atributos
length(lista1) # devuelve el numero de componentes de la lista str(lista1) # devuelve informacion sobre la estructura de la lista a <- lista1[[1]] # devuelve el 1er objeto de la lista en forma de vector/matriz (y nombre excluido) a; typeof(a) b <- lista1[1] # devuelve una sublista compuesta por los elementos de la 1a entrada de la lista (nombre incluido) b; typeof(b) nomb <- c("reales","enteros", "matriz", "caracteres", "logicos") names(lista1) <- nomb; # asigna nombres a los componenes de la lista names(lista1) # devuelve los nombres de los componentes (si los hay) lista1[["reales"]] lista1$logicos # equivalente a lista1[["logicos"]] y a lista[[5]] # Otro modo de definir una lista que incluye nombres para los componentes lista2 <- list(A=mis.num, B=mis.num2); lista2 names(lista2)
7.2 Attach
Las funciones attach() y detach () "cargan" listas y simplifican las referencias.
attach(lista1, warn.conflicts = T)
8 Hojas de datos (data frames)
Son un tipo especial de lista:
- con estructura tabular
- donde las columnas pueden ser de distinto tipo
- suelen utilizarse como bases de datos donde cada fila representa una unidad de observacion y cada columna una variable
var1 <- seq(150,700,50); var2 <- 6:17; var3 <- c(10, 35, 17); var4 <- 2:8 data.frame(var1, var2) data.frame(var1, var3) # recicla var3 data.frame(var1, var4) # mensaje de error porque la longitud var1 no es multiplo de la de var4 var5 <- factor(rep(c(2,1,1,1), times=3), levels=1:2, labels=c("UsaTwitter", "NoUsaTwitter"));var5 mis.datos <- data.frame(var5, var1, var2, var3);mis.datos # por defecto los nombres de las columnas corresponden a los nombres de los objetos # para nombrar las filas (no siempre necesario) row.names(mis.datos) <- paste("suj",1:max(length(var1),length(var2),length(var3)), sep='');mis.datos # Equivale a rownames(mis.datos)
8.1 Algunos atributos
rownames(mis.datos); colnames(mis.datos) dim(mis.datos) # dimensiones str(mis.datos) # estructura attributes(mis.datos) # clase y nombres de fila y columna
8.2 Particiones de una hoja de datos segun un factor, funcion split()
Uso: split (x, factor). unsplit(x, factor)
trozos <- split (mis.datos, mis.datos$var5); trozos # devuelve una lista "trozos" con los valores de var1 segregados segun var5 typeof(trozos) junto.otravez <- unsplit(trozos, mis.datos$var5); junto.otravez # el 2do argumento funciona como factor, no es necesario que este definido como tal. # Ademas se pueden seleccionar solo algunas variables (ej, var2 y var3) trozos <- split (c(mis.datos$var1,mis.datos$var2), mis.datos$var3); trozos # lo mismo podria conseguirse con una sucesion de instrucciones como esta: twitter.si <- mis.datos[mis.datos$var5=="UsaTwitter",];twitter.si
8.3 Filtrado segun condiciones logicas. Función subset()
Uso: subset(x, subset, select)
subset(mis.datos, var3 >= 15) subset(mis.datos, var3 >= 15, select=c(var1,var2)) subset(mis.datos, var3 >= 15, select=-var3)
subset(data$q_005, (data$q_005 > 1900) & (data$q_005 < 2013))
8.4 Transformaciones de datos, funcion transform()
Uso: transform()
mis.datos2 <- transform (mis.datos, log.var1= log(var1)); mis.datos2
8.5 Edición manual de bases de datos
R dispone de una utilidad de aspecto sejemante a una hoja de calculo para:
- introducir datos manualmente
- modificar unos pocos elementos de bases de datos existentes
edit(data.frame()) # abre el editor con una 'hoja de datos' vacia # los huecos se rellenaran solos con NA # doble click en "var?" para cambiar el nombre de la variable # ¿que ha pasado?, ¿donde estan los datos? respuestas <- edit(data.frame()) # Es esencial asignarle un nombre a la hoja de datos respuestas respuestas2 <- edit(respuestas) # guarda los cambios en una nueva hoja de datos (y conserva la antigua) respuestas;respuestas2 fix(respuestas) # sobreescribe los cambios. Es equivalente a respuestas <- edit(respuestas)
9 Cargar datos nativos de R
La instalacion basica de R incluye el paquete 'datasets'.
data() # lista todas las bases de datos disponibles en datasets' help(women) # informacion sobre ese conjunto de datos concreto data(women, package="datasets") # carga esos datos # Si el paquete de procedencia esta cargado basta simplemente: women data(women) # Equivalente data('women')
Muchos paquetes adicionales tambien vienen con sus datos de ejemplo
10 Paquetes R
Cargar un paquete. Requiere instalarlo previamente bien:
- descargandolo en www.r-project.org y siguiendo las instrucciones
- utilizando el menu windows de RGui y eligiendo un repositorio
- con la función install.packages("package")
install.packages("adagio") library (foreign) # o mediante el menu windows de RGui library(help=foreign) # informacion sobre el paquete 'foreign'
10.1 Ruta de busqueda (search path).
Puede haber muchos objetos con el mismo nombre en distintos paquetes. La ruta de busqueda explicita el orden en que se van a recorrer los paquetes para buscar objetos partes deOrden en el que R
search() # muestra la ruta de busqueda # Cuidado con los nombres de las variables del espacio de trabajo: Tienen prioridad rm(list=ls()) women attach(women) mean(weight) weight <- 0; mean(weight) mean(women$weight)
10.2 ACERCA de attach()
attach() inserta el paquete en la ruta de busqueda de R de modo que todos los objetos del paquete estan accesibles mediante sus nombres Uso: attach(what, pos = 2, name = deparse(substitute(what)), warn.conflicts = TRUE) Precauciones:
- Cuidado con el enmascaramiento (nombres duplicados)
- attach no proporciona actualizacion dinamica de los objetos
10.2.1 Retirar paquetes
detach (package:ISwR) detach() # elimina el paquete que haya en 2a posicion plot(c(1,2,5,6), c(2,5,6,7)) # representa un diagrama de dispersion (vease la seccion 'Graficos') detach() # elimina el nuevo 2do paquete plot(c(1,2,5,6), c(2,5,6,7))
11 Escritura de datos en un archivo
11.1 Guardar en un archivo de texto plano una hoja de datos. Funcion write.table()
rm(list=ls()) data(women) write.table(women, "../datos/mujeres.dat", append=T) # Si la carpeta datos no existe previamente dara un mensaje de error ?write.table write.table(women, "../datos/mujeres.dat", append=T, row.names=F, col.names= F, sep='', dec=',') write.table(women, "../datos/mujeres.dat", append=F,quote=F)
11.2 Guadar en texto plano cualquier objeto. Funcion write(). Es mas general, pero no funciona bien con data frames
x<-seq(1,3,0.1) write(x,"../datos/numeros.dat", ncolumns=1,append=T) write(x,"../datos/numeros.dat", ncolumns=length(x),append=F)
12 Lectura de datos desde un archivo
12.1 Información completa en el manual "R Data Import/Export" accesible desde la ayuda
Es posible cargar datos desde un fichero de texto y hacer algún cálculo.
?read.table rm(list=ls()) mis.datos <- read.table("../datos/mujeres.dat", header=T, check.names=T, sep=" ",na.strings=999, comment.char = "#") edit(mis.datos) median(mis.datos$weight) colMeans(mis.datos) attach(mis.datos) sum(weight >=125)/length(weight) # calcula la proporcion de individuos que pesan al menos 125 lb
También es posible cargar datos desde un fichero de texto y manipular datos.
rm(list=ls(all=TRUE)) b=read.table("Oster.txt", header=T) dim(b) summary(b) names(b)# obtengo los nombres de las cabeceras de columna head(b)# Quiero ver los 6 primeros valores tail(b) # Los seis últimos attach(b)# anclamos la base de datos para trabajar con ella View(b)# Veo mis datos edit(b)
12.2 Existen variantes de read.table() para formatos algunos formatos habituales: read.csv, read.csv2, read.delim …
Tambien hay multiples paquetes como 'foreign' para leer archivos creados por programas como SPSS, SAS, Minitab, etc., aunque resulta mas comodo usar esos programas para guardar los datos en texto plano y usar read.table para importar a R
library(foreign) # En primer lugar hay que cargar el paquete mi.encuesta.spss <- read.spss("../datos/SurveyStatisticsI.sav", to.data.frame=T) edit(mi.encuesta.spss)
Mas informacion sobre el uso de foreign en el pdf descargable junto con el paquete.
Paquete para importar y exportar datos de excel: ‘xlsReadWrite’ (no incluido en el sistema base; hay que importarlo)
12.3 Funcion scan().
Mas flexible pero mas complicada de utilizar. Permite especificar el modo de las variables
?scan mis.datos <- scan("../datos/mujeres.dat", what = list(character(), numeric(),numeric()), skip=1);mis.datos mis.datos2 <- scan("../datos/mujeres.dat", what = list(0,0,0), skip=3)
13 Gráficos
R tiene un sistema de graficos muy completo que permite producir graficos con calidad de publicacion Hay múltiples funciones, parametros y paquetes graficos. Se van aprendiendo a medida que se van utilizando Para una buena introduccion: R for Beginners de E. Paradis (citado en la documentacion)
13.1 Gestion de ventanas graficas
x11() # abre una nueva ventana grafica (Unix y Wndows) windows(width=15, height=15, title='Figura 1') # abre una nueva ventana grafica (Windows) quartz() # abre una nueva ventana grafica (Mac OS) pdf("prueba.pdf") # abre un pdf de nombre "prueba.pdf" y redirige el output gráfico al archivo pdf postscript("prueba.eps") # abre un archivo de postscript encapsulado de nombre "prueba.eps" y redirige el output gráfico dev.list() # lista todos los dispositivos grafico abiertos. Los nums son los identificadores dev.cur() # devuelve el identificador del dispositivo grafico activo dev.set(2) # establece como dispositivo gráfico activo el dispositivo num 2 dev.off() # cierra el dispositivo grafico activo; dev.off(3) cierra el dispositivo grafico num 3
13.2 Funcion plot(). Es la función gráfica basica.
Segun los argumentos puede generar un diagrama de dispersion, de barras, caja y bigotes etc.
x <- 1:15; y <- 30:16 # con datos cuantitativos produce diagramas de dispersion plot (x,y, pch=2) plot(women) # puede usarse una hoja de datos o matriz. OJO solo toma las 2 primeras columnas plot(y) # representa los valores de 'y' en funcion de su posicion en el vector # plot (factor) produce un diagrama de barras f <- factor(rep(c(1,2,2,1,1,3), each=2), levels=1:3, labels=c("No fumador", "Fumador", "NS/NC")) plot(f) # tambien existe la funcion barplot para crear diagramas de barras # plot (var.cuant, factor) produce diagramas de caja organizados según los niveles del factor entorno <- factor(rep(c(1,2,2), each=5), levels=1:2, labels=c("Urbano", "Rural")) plot(entorno,y)
13.3 Plot a Png
Se usa la función png para almacenar la salida en un fichero de imagen.
x <- 1:15; y <- 30:16 # con datos cuantitativos produce diagramas de dispersion png(filename="/tmp/figure.png", height=295, width=300, bg="white") plot (x,y, pch=2) dev.off()
13.4 Funcicion coplot(). Matriz de graficos de dispersion separados segun factores
datos2 <- data.frame(entorno,women) coplot(women$weight ~ women$height | entorno) # representa peso en funcion de altura separando segun entorno socioeconomico<- factor(c(3,3,2,2,2,2,1,1,1,1,1,1,1,2,2), levels=1:3, labels=c("bajo","medio","alto")) plot(socioeconomico) coplot(women$weight ~ women$height | entorno+socioeconomico)
13.5 Funcion pairs(). Ejemplo con datos simulados de los "Big Five" de personalidad
pers <- matrix(rnorm(1000), ncol = 5) # genera 1000 datos con distribucion normal. Vease seccion 'Muestreo y probabilidad' colnames(pers) <- c("extraversion", "estabilidad", "apertura", "responsabilidad", "amabilidad") pairs(pers) # Nótese la diferencia con el resultado de plot() x11(); plot(pers) # usa solo las 2 primeras variables
13.6 Funcion hist(). Histogramas
pers <- as.data.frame(matrix(rnorm(1000), ncol = 5)) colnames(pers) <- c("extraversion", "estabilidad", "apertura", "responsabilidad", "amabilidad") cor(pers, method='pearson') # matriz de correlaciones. Otros metodos: 'spearman' y 'kendall' round(cor(pers),3) # la misma matriz de antes, pero con solo 3 decmales attach(pers) hist(extraversion) x11(); hist(extraversion, labels=T, col= "gray", xlim=c(-4,4), ylab="Frecuencia absoluta", sub= "Figura 1. Ejemplo de histograma", main= "Extraversión") hist(amabilidad, col= "red", border="white", add=T, xlim=c(-4,4), sub= "Figura 2. Otro ejemplo de histograma", main= c("Extraversión Amabilidad superpuestos")) x11(); hist(extraversion,freq=F, labels=T, col= "light blue",xlim=c(-4,4),ylab="Frecuencia relativa", las=1,sub= "Figura 1b. Ejemplo de histograma", main= "Extraversión") hist(amabilidad, col= "gray", border="white", xlim=c(-4,4), sub="", main= "Figura 3. Otro ejemplo de histograma cambiando el color") pdf("ejemplos_histograma.pdf") hist(extraversion, nclas=3, col= "light blue", sub= "Figura 3. Ejemplo de histograma", main= "Extraversión") hist(extraversion, breaks=seq(-4,4,0.5),col= "pink",main= "Figura 3. Ejemplo de histograma") dev.off()
13.7 Funcion boxplot()
boxplot(extraversion, main= "Diagrama de caja de extraversión") text(locator(5),"más outliers", adj=0) boxplot(pers)
13.8 Funciones qqnorm() y qline
qqnorm(extraversion, col="red") qqline(extraversion)
Tambien existe qqplot(), que permite utilizar otras distribuiones distintas de la normal
13.9 Catalogo de nombres de color
colors()
13.10 Parametros gráficos. Funcion par()
Produce cambios permanentes en el dispositivo gráfico activo Permite contral cada elemento del grafico final
?par par.actual <- par() # guarda los valores actuales (para poder reestablecerlos posteriormente) par(font.main=4, las=1, xaxs= "r", col.lab="blue", col="skyblue3",lwd=3) hist(extraversion,xlim=c(-4,4)) plot(extraversion, amabilidad) # par(par.actual) # reestablece los valores guardados par(mfrow=c(2,1)) # distribuye el espacio grafico e dos filas
13.11 Representacion de funciones
curve(x^2) curve(x^2, from=-3, to=3) # los cambios introducidos a traves de par no tienen efecto en la nueva ventana # tambien podria haberse hecho con plot x11() x <- seq(-3,3,0.25) y <- x^2 plot(x,y) plot (x,y, type="o") # 'p': puntos, "l": lines, 'b': puntos conectados mediante lineas, etc. # un ultimo ejemplo hist(extraversion, freq=F) curve(dnorm(x), add=T) detach(pers)
14 Muestreo y probabilidad
14.1 Muestreo
datos <- c(1,5,9,8,52,2,4) sample(datos,3) # toma 3 elementos de 'datos'. Por defecto, muestrea SIN reposición dado <- 1:6 # para simular el lanzamiento de 5 dados hace falta muestreo CON reposición sample(dado,5,replace=T) # sample() tambien permite incluir un vector de pesos si los sucesos no son equiprobables
Otro ejemplo: extracción simultánea de dos cartas de una baraja española.
baraja <- paste(rep(c("As",2:7,"Sota","Caballo","Rey"),4),rep(c("Oros","Copas","Espadas","Bastos"),each=10));baraja sample(baraja,2)
Generacion de numeros pseudo-aleatorios.
# set.seed() # inicializa el generador de aleatorios con una semilla concreta (reproductibilidad) set.seed(-98)
14.2 Distribuciones de probabilidad
Funciones de generacion de numeros pseudo-aleatorios conforme a una distrbucion dada Estructura: rfunc(n, par1, par2,…)
datos.gauss <- rnorm (100, mean=5, sd=2) # genera 100 numeros ~N(5,2) runif (4,-10,10) # genera 4 numeros ~U(-10, 10)
Tambien existen rbinom, rpois, rt, rf, rchisq, rbeta, rgamma, etc en ocasiones conviene mantener la semilla estable.
Funciones de densidad. Devuelven la ordenada de la funcion (de densidad) de probabilidad para un vector de valores de la V.A.
Estructura: dfunc(n, par1, par2,…, log=FALSE)
x1 <- seq(-5,5,0.01) gauss <-dnorm(x1) plot (x1,gauss,"l", lwd=2, col="blue") student <- dt(x1, df=4) x11(); plot (x1,student,"l", lwd=2,col="green") x2 <- seq(0,10,0.01) chi2 <-dchisq(x2,df=3) x11(); plot (x2,chi2,"l",lwd=2,col="red")
Funciones de distribucion probabilidad Estructura: pfunc(n, par1, par2,…, lower.tail = TRUE, log.p = FALSE)
pnorm(0.7) # devuelve la prob de que una V.A.con distribución N(0,1) tome un valor menor o igual que 0.7 pnorm(0.7, mean=10, sd=0.5) pchisq(3.84, df=3, lower.tail=F) # devuelve la probabilidad de que una V.A.~ Chi2 con 3gl tome un valor mayor que 3.84
Funciones que devuelven cuantiles Estructura: qfunc(n, par1, par2,…, lower.tail = TRUE, log.p = FALSE)
qnorm(0.025) # devuelve el cuantil 2.5 (valor de la V.A. que deja por debajo el 2.5% de la distribución) qf(0.975, df1=4, df2=53) # se puede usar para obtener los puntos de corte de la región de rechazo en un contraste de hipótesis
15 Definicion de funciones
Estructura: nombre.funcion <- function(argumento1, argumento2){ instruccion1; instruccion2; …; return(variable.resultado)} Las llaves no son necesarias si la funcion se define en una sola linea
Ejemplo: Definicion de una funcion que calcula la media
media <- function(datos){ # pueden darse valores por defecto a los argumentos numerador <-sum(datos) denominador <- length(datos) xm <- numerador/denominador # podriamos llamar "media" a la variable de retorno return(xm)} # devuelve al programa principal el VALOR de xm (pero la variable xm es local y sale)
Utilizacion de la funcion una vez definida (y cargada en memoria)
x <- rnorm(100) media(x) # la llamada se hace mediante el nombre de la función mi.media <- media(datos=x);mi.media
Recomendaciones:
- crear un directorio especifico para funciones propias
- guardar cada funcion en un script independiente
Cargar una funcion en el espacio de trabajo, funcion source()
source("../funciones/descriptivos.r") # carga la funcion 'descriptivos' que esta en el subdirectorio 'funciones' resultados <- descriptivos(x) # devuelve una LISTA de elementos que se almacenan en la variable 'resultados'
Para crear funciones minimamente complejas hace falta aquirir unos minimos conocimientos de programación.
16 Descriptivos
Cuando vamos a analizar unos datos lo primero que debemos mirar normalmente son los descriptivos, veamos cómo:
library(foreign) # descriptivos de paises. Ejercicio 2.1 paises.csv <- read.csv('paises.csv', sep=';', na.strings=999) media <- mean(paises.csv$Porcentaje) mediana <- median(paises.csv$Porcentaje) desviacionTipica <- sd(paises.csv$Porcentaje) varianza <- var(paises.csv$Porcentaje) cuantil <- quantile(paises.csv$Porcentaje) todos <- summary(paises.csv$Porcentaje) max <- max(paises.csv$Porcentaje) min <- min(paises.csv$Porcentaje) amplitud <- max - min tr=c(3,4,5,3,4,5,4,3,2,3,12,11,3,4,89) mean(tr,trim=5/100) # media recortada al 5% IQR(tr) # rango intercuartil stem(tr) # diagrama de tallos y hojas (ver http://www.estadisticaparatodos.es/taller/graficas/tallos_hojas.html)
En nuestro ejemplo utilizaremos los datos de paises.csv
17 Distribuciones
17.1 Binomial
> dbinom(4, size=12, prob=0.2) [1] 0.1329 > pbinom(4, size=12, prob=0.2) [1] 0.92744
17.2 Poisson
> ppois(16, lambda=12) # lower tail [1] 0.89871 > ppois(16, lambda=12, lower=FALSE) # upper tail [1] 0.10129
17.3 Distribución Contínua Uniforme
> runif(10, min=1, max=3) [1] 1.6121 1.2028 1.9306 2.4233 1.6874 1.1502 2.7068 [8] 1.4455 2.4122 2.2171
17.4 Exponencial
> pexp(2, rate=1/3) [1] 0.48658
17.5 Normal
> pnorm(84, mean=72, sd=15.2, lower.tail=FALSE) [1] 0.21492
17.6 Chi cuadrado
> qchisq(.95, df=7) # 7 degrees of freedom [1] 14.067
17.7 T de Student
> qt(c(.025, .975), df=5) # 5 degrees of freedom [1] -2.5706 2.5706
17.8 F distribución
> qf(.95, df1=5, df2=2) [1] 19.296
18 Test estadísticos
18.1 Distribución Binomial
Si hay doce coches cruzando un puente por minuto de media, encuentra la probabilidad de tener dieciseis o más coches cruzando el puente en un minuto particular.
> ppois(16, lambda=12) # lower tail [1] 0.89871 > ppois(16, lambda=12, lower=FALSE) # upper tail [1] 0.10129
La probabilidad es del 10.1%
18.2 Test Binomial
Supongamos que tenemos un juego de mesa que depende del lanzamiento de un dado (ver dado) y asigna una especial importancia al lanzamiento de un 6. En un juego en particular, el dado es lanzado 366 veces, y el 6 sale 51 veces. Si el dado es justo, esperamos que el 6 salga 366/6 = 61 veces. ¿Es la proporción del 6 significativamente más alta de lo que se esperaría por casualidad, sobre la hipótesis nula de un dado justo?
binom.test(51, 336, (1/6))
scipy.stats.binom_test(51, 336, 1.0/6)
18.3 Distribución F
El 95 percentil de la distribución F con (5, 2) grados de libertad es 19.296
> qf(.95, df1=5, df2=2) [1] 19.296
18.4 Bondad de ajuste multinomial
Dadas unas estadísticas de fumar en un campus, determine si los datos de ejemplo en la encuesta soportan un nivel de significación 0.5
> library(MASS) > levels(survey$Smoke) > smoke.freq = table(survey$Smoke) > smoke.freq > smoke.prob = c(.045, .795, .085, .075) > chisq.test(smoke.freq, p=smoke.prob) Chi-squared test for given probabilities data: smoke.freq X-squared = 0.1074, df = 3, p-value = 0.991
Como el p-valor 0.991 es mayor que el nivel de significación 0.05, no rechazamos la hipótesis nula que los datos de ejemplo en la encuesta soportan las estadísticas de fumar.
18.5 Test de Chi
Establezcamos la hipótesis del hábito de fumar de los estudiantes es independiente de su nivel de ejercicio en un nivel de significancia de 0.05.
> library(MASS) # load the MASS package > tbl = table(survey$Smoke, survey$Exer) > tbl # the contingency table Freq None Some Heavy 7 1 3 Never 87 18 84 Occas 12 3 4 Regul 9 1 7 > chisq.test(tbl) Pearson’s Chi-squared test data: table(survey$Smoke, survey$Exer) X-squared = 5.4885, df = 6, p-value = 0.4828 Warning message: In chisq.test(table(survey$Smoke, survey$Exer)) : Chi-squared approximation may be incorrect
Como el p-valor es 0.4828 es mayor que el nivel de significado 0.05, no rechazamos la hipótesis nula que el hábito de fumar es independiente del nivel de ejercicio de los estudiantes.
18.6 Test Wilcox
Sin asumir que los datos tienen distribución normal, chequea en el nivel de significación 0.05 si los rendimientos entre 1931 y 1932 en los datos tienen distribuciones de datos idénticas.
> library(MASS) # load the MASS package > head(immer) Loc Var Y1 Y2 1 UF M 81.0 80.7 2 UF S 105.4 82.3 > wilcox.test(immer$Y1, immer$Y2, paired=TRUE) Wilcoxon signed rank test with continuity correction data: immer$Y1 and immer$Y2 V = 368.5, p-value = 0.005318 alternative hypothesis: true location shift is not equal to 0 Warning message: In wilcox.test.default(immer$Y1, immer$Y2, paired = TRUE) : cannot compute exact p-value with ties
Las distribuciones son NO idénticas.
18.7 Ecuación de la Regresión Logística
Para un coche de 120 km/h y un peso de 2800, la probabilidad de estar ajustado es del 64%.
> am.glm = glm(formula=am ~ hp + wt, data=mtcars, family=binomial) > newdata = data.frame(hp=120, wt=2.8) > predict(am.glm, newdata, type="response")
18.8 Test de significado para la regresión logística
En el nivel de significado 0.05, decidimos si cualquiera de las variables independientes en el modelo de regresión logística de la transmisión del vehículo en los datos asigna mtcars es estadísticamente insignificante.
> am.glm = glm(formula=am ~ hp + wt, + data=mtcars, + family=binomial) > summary(am.glm) Call: glm(formula = am ~ hp + wt, family = binomial, data = mtcars) Deviance Residuals: Min 1Q Median 3Q Max -2.2537 -0.1568 -0.0168 0.1543 1.3449 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 18.8663 7.4436 2.53 0.0113 * hp 0.0363 0.0177 2.04 0.0409 * wt -8.0835 3.0687 -2.63 0.0084 ** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 43.230 on 31 degrees of freedom Residual deviance: 10.059 on 29 degrees of freedom AIC: 16.06 Number of Fisher Scoring iterations: 8
18.9 Mann-Whitney test: comparar si hay diferencias entre medianas muestras independientes
boxplot(Peso ~ Manejo) wilcox.test(Peso~Manejo) boxplot(Ostertagia ~ Manejo) plotmeans(Ostertagia~Manejo,col="coral", connect=T, ylab="Carga parasitaria (Ostertagia)", xlab="Manejo", data=b)# plot means and CI95% Grafica mitjana
19 Ecuación de Regresión Logística
Usando la ecuación de regresión logística, estime la probabilidad de un vehículo siendo ajustado con una transmisión manual si tiene un motor de 120 y pesa 2800.
> am.glm = glm(formula=am ~ hp + wt, data=mtcars, family=binomial) > newdata = data.frame(hp=120, wt=2.8) > predict(am.glm, newdata, type="response") 1 0.64181
Para un automóvil con un motor de 120 hp y 2800 lbs de peso, la probabilidad de estar ajustado con una transmisión manual es del 64%.
En el nivel de significado 0.05, decide si cualquiera de las variables independientes en el modelo de regresión logística de la transmisión del vehículo en los datos mtcars es estadísticamente insignificante.
> am.glm = glm(formula=am ~ hp + wt, data=mtcars, family=binomial) > summary(am.glm) Call: glm(formula = am ~ hp + wt, family = binomial, data = mtcars) Deviance Residuals: Min 1Q Median 3Q Max -2.2537 -0.1568 -0.0168 0.1543 1.3449 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 18.86630 7.44356 2.535 0.01126 * hp 0.03626 0.01773 2.044 0.04091 * wt -8.08348 3.06868 -2.634 0.00843 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 43.230 on 31 degrees of freedom Residual deviance: 10.059 on 29 degrees of freedom AIC: 16.059 Number of Fisher Scoring iterations: 8
20 Análisis de Varianza
20.1 Diseño Completamente Aleatorio
> df1 = read.table("fastfood-1.txt", header=TRUE); df1 Item1 Item2 Item3 1 22 52 16 2 42 33 24 3 44 8 19 4 52 47 18 5 45 43 34 6 37 32 36 > r = c(t(as.matrix(df1))) # response data > r [1] 22 52 16 42 33 ... > f = c("Item1", "Item2", "Item3") # treatment levels > k = 3 # number of treatment levels > n = 6 # observations per treatment > tm = gl(k, 1, n*k, factor(f)) # matching treatments > tm [1] Item1 Item2 Item3 Item1 Item2 ... > av = aov(r ~ tm) > summary(av) Df Sum Sq Mean Sq F value Pr(>F) tm 2 745 373 2.54 0.11 Residuals 15 2200 147
21 Análisis de Clúster Jerárquico
En los datos, se puede ejecutar la matriz distancia con hclust y plotear un dendograma que muestre una relación jerárquica entre vehículos.
> d <- dist(as.matrix(mtcars)) # find distance matrix > hc <- hclust(d) # apply hirarchical clustering > plot(hc) # plot the dendrogram
Para un conjunto de datos con 4000 elementos, coge hclust 2 minutos para finalizar el trabajo en un AMD.
> test.data <- function(dim, num, seed=17) { + set.seed(seed) + matrix(rnorm(dim * num), nrow=num) + } > m <- test.data(120, 4500)
22 Licencia
Este documento está bajo una Licencia Creative Commons Reconocimiento Unported 3.0