martes, 21 de agosto de 2012

Clase 2: Análisis unidimensional en R-project

Análisis unidimensional de los datos

Luego de familiarizarnos con el manejo básico de datos en R (post anterior), podemos comenzar a realizar el análisis unidimensional de los datos (de 1 dimensión o 1 variable). 
El estudio de una variable se divide en su análisis descriptivo, el inferencial y el probabilístico. En esta entrada nos basaremos en el primero de ellos.
El análisis descriptivo consta de tablas de frecuencias (absolutas, relativas, acumuladas), gráficos y cálculo de estadísticos o medidas que nos permitan resumir la variable en cuestión.


Tipo de variable.

El tipo de variable a estudiar nos determinará el tipo de análisis a seguir. Una variable puede ser cuantitativa (o numérica) o cualitativa (o categórica). Además, una variable cuantitativa puede ser continua (si sus valores toma números reales; e.g. el peso) o discreta (números enteros; e.g. el número de hijos), mientras que una variable  cualitativa puede clasificarse como nominal (categorías sin orden; e.g. el género) u ordinal (categorías con orden; e.g. estado de una enfermedad: leve, moderado, grave).

Tablas de frecuencia, gráficos, estadísticos e inferencia.

La clasificación de la variable de estudio determinará si la utilización de una tabla de frecuencias será con datos agrupados o no agrupados, el tipo de gráfico a seleccionar y qué inferencia podemos extraer de ella. Por ejemplo, las variables discretas, y cualitativas se representan mediante una tabla de frecuencias con datos no agrupados, mientras que para los datos continuos elaboramos intervalos y realizamos la tabla de frecuencia con datos agrupados en dichos intervalos. Ademas, si se trata de una variable continua elegiremos el histograma, mientras que para una variable discreta utilizaremos el diagrama de barras. Por su parte, si la variable es continua tendrá sentido calcular su media, mientras que si la variable es cualitativa este estadístico no lo podemos calcular (e.g. no se puede calcular la media del género).

Estadísticos descriptivos.

Los estadísticos se dividen principalmente en aquellos de centralidad o posición central (media, moda, mediana, etc.), en los de posición no cental (cuantiles: cuartiles, deciles y percentiles), de dispersión o variabilidad (varianza, desviación típica, error típico, coeficiente de variación, rango intercuartílico, etc.), y de forma (asimetría y curtosis).




Aquí les dejo una breve introducción a las principales funciones que hay en R para desarrollar cada uno de estos análisis. 
Por último les dejo algunas funciones para la manipulación de bases de datos.

Saludos!



###########################################################################
# An?lisis de datos con R #
# 2012 #
# Msc. Rosana Ferrero #
# http://statisticalecology.blogspot.com/ #
###########################################################################
  

###########################################################################
# Clase 2. Estad?stica descriptiva #
###########################################################################
# 0. Funciones de interés para manejar los datos #
# 1. Estad?stica unidimensional #
# a. Tablas de Frecuencia #
# b. Gráficos #
# c. Estadísticos #
# 2. Ejemplo: datos Iris #
# 3. Otras funciones #
###########################################################################
 
###########################################################################
# 0. Funciones de inter?s para manejar los datos #
###########################################################################
 
########## Leer los datos
 
### La funci?n "scan"
x=scan()
1: 6
2: 7
3: 3
4: 4
5: 8
6: 5
7: 6
8: 2
9:
 
### La funci?n "read.table". Producen un dataframe
data=read.table("c:\\temp\\regression.txt",header=T)
map=read.table("c:\\temp\\bowens.csv",header=T,sep=",")
data=read.table(file.choose(),header=T)
## opci?n seq=""
#\n newline
#\r carriage return
#\t tab character
#\b backspace
#\a bell
#\f form feed
#\v vertical tab
## chequear un archivo
file.exists("c:\\temp\\Decay.txt")
## leer datos desde archivos con formato no est?ndar. Producen un objeto lista no un dataframe
murders=scan("c:\\temp\\murders.txt", skip=1, what=list("","","",""))
murder.frame=as.data.frame(murders)
murders=read.table("c:\\temp\\murders.txt",header=T)
## leer archivos con diferentes n?meros/valores por l?nea
line.number=length(scan("c:\\temp\\rt.txt",sep="\n"))
(my.list=sapply(0:(line.number-1),function(x) scan("c:\\temp\\rt.txt",skip=x,nlines=1,quiet=T)))
## leer l?neas. Produce un objeto de clase character.
readLines("c:\\temp\\murders.txt",n=-1)
 
######## Ver ?tem 4 para obtener m?s funciones de inter?s
 
###########################################################################
# 1.a. Tablas de frecuencias #
###########################################################################
 
table(rpois(100,5)) # realiza tablas de frecuencias absolutas
table(rpois(100,5))/sum(rpois(100,5)) #realiza tablas de frecuencias relativas
table(cut(rpois(100,5),3)) #para variables num?ricas si queremos hacer una tabla de frecuencias absolutas agrupando los datos por intervalos
margin.table(tabla, 1 o 2 o ...) #realiza tablas marginales del objeto tabla
prop.table(tabla, 1 o 2 o ...) #tabla de frecuencias condicionadas por la dimensi?n indicada con 1 o 2 o... del objeto tabla
ftable(tabla) # tabla de frecuencias para 3 o m?s variables categ?ricas
 
m = matrix(1:4,2); m
margin.table(m,1) # totales por filas
margin.table(m,2) # totales por columnas
prop.table(m,1) # porcentaje por filas
prop.table(m,2) # porcentaje por columnas
 
data(Titanic)
Titanic
ftable(Titanic) # crea tablas de contingencia llanas; x= factores o lista o tabla
ftable(Titanic, row.vars = 1:2, col.vars = "Survived")
 
tabulate(c(2,3,3,5), nbins = 10) # cuenta el n?mero de veces que ocurre cada valor en un vector dado
 
# Tabla de frecuencias de 3 v?as (3 variables)
xtabs(~A+B+C, data) # tabla de frecuencias de 3 v?as
 
# Tabla de contingencia (Cross Tabulation) de 2 v?as (2 variables)
library(gmodels)
CrossTable(mydata$myrowvar, mydata$mycolvar)
 
requires(stats)
with(warpbreaks, table(wool, tension)) #distribuci?n de frecuencias simples
table(state.division, state.region)
with(airquality, table(cut(Temp, quantile(Temp)), Month)) #tabla de contingencia de 2 v?as
 
#La tabla ignora los valores perdidos, para incluir los NA como una categor?a para contar debemos incluir la opci?n exclude=NULL si la variable es un vector. Si la variable es un factor debemos crear un nuevo factor usando: newfactor = factor(oldfactor, exclude=NULL). (Obtenido de "QuickR" -http://www.statmethods.net/stats/frequencies.html-)
 
hist(x, nclass=, breaks=seq(0,100,by=5))
#por ejemplo:
(hist(velocidades, plot=FALSE, col='14',nclass=10,ylim=c(0,10), labels=TRUE,main= 'Velocidad autom?oviles', sub='Via el Volador', ylab='Frecuencias',xlab='Velocidad en Km/h'))
#$breaks (l?mites de clase)
#$counts (frecuencias absolutas)
#$intensities (frecuencias relativas)
#$mids (marcas de clase o puntos medios)
 
###########################################################################
# 1.b. Gr?ficos #
###########################################################################
 
x11() # activa un dispositivo gr?fico
dev.off() #cierra un dispositivo gr?fico
par() # selecciona ciertos par?metros gr?ficos que est?n definidos por defecto y pueden modificarse
demo(graphics) # demostraci?n de gr?ficos y c?digo de R
plot(x,y) # gr?fico de dispersi?n: x e y son vectores
plot(x) #si x es un vector num?rico realiza un gr?fico de sus elementos sobre el ?ndice, y si x es una serie temporal realiza un gr?fico de x frente al tiempo
plot(f) # si f es factor realiza un gr?fico de barras
plot(f, x) # si f es factor y x vector, realiza un diagrama de cajas para cada nivel del factor f
plot(hoja.datos) # realiza un gr?fico de dispersi?n para cada pareja de variables
 
## argumentos de las funciones gr?ficas
#type= "n" (nada) "p" (puntos) "l" (l?neas) "b" (both: punto+l?nea) "o" "s" "S" "h"
#xlab= ylab= (a?ade etiquetas a los ejes)
#main= sub= (a?ade t?tulo y subt?tulo)
#xlim=c(,) ylim=c(,) (selecciona los l?mites m?nimo y m?ximo para los ejes)
#add= T o F (solapa un gr?fico con otro ya existente o no)
#col= (indicar color). ver la opci?n colors()
 
## otras funciones gr?ficas
points(x,y, pch=n?) #a?ade puntos definidos por las coordenadas contenidas en los vectores x e y el aspecto indicado en pch. Se puede utilizar type
abline(h=, v=) # a?ade l?neas horizontales y/o verticales
abline(lm) o abline(a,b, lty=n?) #a?ade una recta y=a+b*x con el trazo indicado en lty
legend(x,y) # a?ade la leyenda en el punto indicado
title(main="" , sub="") #a?ade t?tulo y subt?tulo
axis() # modifica elementos referentes a los ejes como color, tickmarks, fuentes, etc
text(x,y,etiquetas) #a?ade etiquetas en las posiciones marcadas por x e y
lines()
curve()
 
#gr?ficos simples
hist(vector, nclass, breaks, probability, plot) #histograma
mid.age=c(2.5, 7.5, 13, 16.5, 17.5, 19, 22.5, 44.5, 70.5) ;
acc.count=c(28, 46, 58, 20, 31, 64, 149, 316, 103) ;
age.acc=rep(mid.age, acc.count) ; age.acc ;
brk=c(0, 5, 10, 16, 17, 18, 20, 25, 60, 80) ; hist(age.acc, breaks=brk)
A=hist(x); lines(c(min(A$breaks), A$mids, max(A$breaks)), c(0, A$counts, 0), type="l", col="red") #histograma con pol?gono de fercuencias
 
x=seq(100, 145, by=5) ; n=length(x) #funci?nd e distribuci?n emp?rica
plot(sort(x), (1: n)/n , type=?s?, ylim=c(0,1))
 
qqnorm(x) #gr?fico de probabilidad normal, en abscisas est?n los valores esperados de los cuantiles de la normal y en las ordenadas los valores de x
qqplot(x) #gr?fico cuantil-cuantil, en abscisas est?n lso cuantiles del vector x y en ordenadas los cuantiles del vector y
 
boxplot(vector1, vector2, plot=F o T) #diagrama de cajas como vectores
boxplot(formula, data=NULL, subset, plot=F o T) # diagrama de cajas m?ltiple
plot(factor, vector) #diagrama de cajas para cada nivel del factor f
 
stripchart(formula o x, method="overplot" o "jitter" o "stack") #gr?fico de puntos
stem(vector) #diagrama de tallo-hojas
 
plot(table()) #diagrama de barras para una variable num?rica discreta
plot(factor) #diagrama de barras para una variable categ?rica
barplot(vector o tabla, names.args=NULL) #cada valor representa la altura de una barra que podemos etiquetar opcionalmente con names.args
barplot(matriz o tabla, beside=T (adyacentes) o F (apiladas)) #diagrama de barras m?ltiples apiladas o no.
pie(x,labels=names(x), col) #diagrama de sectores
 
#gr?ficos m?ltiples por ventana
par(mfrow=c(filas, columnas)) #divide la pantalla gr?fica en tantas filas y columnas como se indique
split.screen()
#datos multivariantes
plot(data.frame) #realiza un gr?fico de dispersi?n de todos los pares de variables
pairs(A) # matriz de gr?ficos planos, un gr?fico de dispersi?n para cada pareja de columnas de A
coplot(a~b|c) #gr?ficos condicionales
persp() #representa gr?ficamente funciones de 2 variables
contour() #representa curvas de nivel de una funci?n de 3 variables
image() #representa un mapa de 2 variables
matplot() #representa una matriz frente a otra por columnas
matpoints o matlines # a?ade puntos o l?neas a un gr?fico (matrices)
 
#otros
library(lattice)
xyplot(circumference ~ age, groups=Tree, data=Orange)
 
# gr?ficos con barra de error: "The R Book" (Crawley, 2009)
# en diagrama de barras
error.bars=function(yv,z,nn)
{
xv=
barplot(yv,ylim=c(0,(max(yv)+max(z))),names=nn,ylab=deparse(substitute(yv)))
g=(max(xv)-min(xv))/50
for (i in 1:length(xv))
{
lines(c(xv[i],xv[i]),c(yv[i]+z[i],yv[i]-z[i]))
lines(c(xv[i]-g,xv[i]+g),c(yv[i]+z[i], yv[i]+z[i]))
lines(c(xv[i]-g,xv[i]+g),c(yv[i]-z[i], yv[i]-z[i]))
}
}
comp=read.table("c:\\temp\\competition.txt",header=T)
attach(comp)
names(comp) #[1] "biomass" "clipping"
se=rep(28.75,5)
labels=as.character(levels(clipping))
ybar=as.vector(tapply(biomass,clipping,mean))
error.bars(ybar,se,labels)
 
# en diagrama de puntos
xy.error.bars=function (x,y,xbar,ybar)
{
plot(x, y, pch=16, ylim=c(min(y-ybar),max(y+ybar)),xlim=c(min(x-xbar),max(x+xbar)))
arrows(x, y-ybar, x, y+ybar, code=3, angle=90, length=0.1)
arrows(x-xbar, y, x+xbar, y, code=3, angle=90, length=0.1)
}
x = rnorm(10,25,5)
y = rnorm(10,100,20)
xb = runif(10)*5
yb = runif(10)*20
xy.error.bars(x,y,xb,yb)
 
#guardar gr?ficos
pdf(file="f1,pdf", width=8, heigth=10)
plot(runif(10))
dev.off()
#otra forma
plot(runif(50))
dev.copy2eps() #el nombre del archivo por defecto es Rplot.eps y el formato es postcript
 
###########################################################################
# 1.c. Estad?sticos simples #
###########################################################################
 
#mean Media aritm?tica
#mean(datos, trim=.05) Media truncada o recortada al 5%
#median El percentil 0.5: la mediana
#min El m?nimo de una serie de n?meros
#max El m?ximo de una serie de n?meros
#quantile Los percentiles de una distribuci?n
#range M?nimo y m?ximo de un vector
#sum Suma aritm?tica
#var Cuasi-Varianza y cuasi-covarianza
#var(datos)*(length(datos)-1)/length(datos) Varianza
#sd Cuasi-desviaci?n t?pica
#sqrt(var(datos))*(length(datos)-1)/length(datos) Desviaci?n t?pica
#sd(datos)/abs(mean(datos)) Coeficiente de Variaci?n
#mad Desviaci?n media de la mediana
#skewness Asimetr?a
#kurtosis Curtosis o apuntalamiento
#summary Resumen de estad?sticas de una serie de datos
#fivenum Retorna los 5 n?meros de Tukey (de cajas de Tukey): minimum, lower-hinge, median, upper-hinge, maximum.
#cor Correlaci?n (admite uno o dos argumentos)
#cumsum Suma cumulativa de un vector
#prod El producto de los elementos de un vector
#sample Muestreo aleatorio (y permutaciones)
 
### Veamos algunos ejemplos:
mean(1:10) #media aritm?tica
exp(mean(log(abs(1:10)))) #media geom?trica
1/mean(1/1:10) #media arm?nica
median(1:10) #mediana
quantile(1:10) #cuartiles
quantile(1:10, probs=seq(0,1, by=1/10)) #deciles
weigthed.mean()
IRQ()
range()
cor()
tapply(x,f) #para obtener res?menes por grupos
100*sd(x)/mean(x) # CV coeficiente de variaci?n
 
#Estad?sticas simples:
x = seq(1:10)
x
# [1] 1 2 3 4 5 6 7 8 9 10
cumsum(x)
# [1] 1 3 6 10 15 21 28 36 45 55
median(x)
#[1] 5.5
 
 
###########################################################################
# 2. Ejemplo: datos Iris
###########################################################################
data(iris)
head(iris)
attach(iris)
 
table(Species)
summary(Sepal.Length) # resumen estad?stico de 1 variable
summary(iris) #resumen estad?stico de varias variables
tapply(Sepal.Length, Species, mean)
var(Sepal.Length)
cor(iris[1:4]) #matriz de correlaciones
 
plot(iris) #gr?fico de un dataframe
plot(Sepal.Length) # gr?fico de 1 variables
plot(Sepal.Length, Petal.Length) # gr?fico de 2 variables; ?dem plot(Petal.Length~Sepal.Length)
plot(Species[Sepal.Length>5.1]) # gr?fico de 1 factor
plot(Species, Sepal.Length) # gr?fico de 1 factor (diagrama de cajas)
hist(Sepal.Length) #histograma
pie(table(Species[Sepal.Length>5]))
qqline(Sepal.Length) # comparaci?n de un vector con los valores esperados normales
qqplot(Sepal.Length, Petal.Length) # comparaci?n de las distribuciones de dos vectores o variables
cols = rep(c("red","blue","green"), each=50)
pich = rep(c(21,22,23),each=50)
plot(Sepal.Length,Petal.Length, pch=pich, bg=colores)
 
library(ctest)
 
shapiro.test(Sepal.Length) # contraste de normalidad
mu=mean(Sepal.Length)
de=sd(Sepal.Length)
ks.test(Sepal.Length, "pnorm", mean=mu, sd=se) # contraste de normalidad
 
boxplot(Sepal.Length, Petal.Length)
t.test(Sepal.Length, Petal.Length) # igualdad de medias para dos muestras independientes sin hip?tesis sobre la varianza
var.test(Sepal.Length, Petal.Length) # igualdad de varianzas
wilcox.test(Sepal.Length, Petal.Length) # contraste para muestras no normales de Wilcoxon (o Mann-Whitney)
ks.test(Sepal.Length, Petal.Length) # Kolmogorov-Smirnov para comparar dos distribuciones continuas
 
################################################################################
# 3. Otras funciones #
################################################################################
 
######### Primer visi?n de los datos
worms=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/worms.txt",header=T)
attach(worms)
names(worms)
summary(worms)
#las funciones by y aggregate permiten resumir un dataframe seg?n los niveles del factor
by(worms, Vegetation, mean)
 
## Sub?ndice e ?ndices
worms[3, 5]
worms[14:19, 7]
worms[1:5, 2:3]
worms[3, ]
worms[, 3]
worms[, c(1,5)]
# Seleccionar filas aleatoriamente de un Dataframe
worms[sample(1:20,8),]
# Ordenando los Dataframes
worms[order(Slope), ]
worms[rev(order(Slope)), ]
worms[order(Vegetation, Worm.density)]
# Usar condiciones l?gicas para seleccionar filas de un Dataframe
worms[Damp == T,]
worms[Worm.density > median(Worm.density) & Soil.pH < 5.2,]
worms[,sapply(worms,is.numeric)]
worms[,sapply(worms,is.factor)]
# Drop filas
worms[-(6:15),]
worms[!(Vegetation=="Grassland"),]
worms[-which(Damp==F),]
 
### Omitir filas que contengan valores perdidos
data=read.table("c:\\temp\\worms.missing.txt",header=T)
data
na.omit(data)
new.frame=na.exclude(data)
data[complete.cases(data),]
apply(apply(data,2,is.na),2,sum)
 
#### Usar order y unique para eliminar la pseudoreplicaci?n
worms[rev(order(Worm.density)),][unique(Vegetation),]
 
#### Ordenar de manera compleja con direcciones mixtas
worms[order(Vegetation,-Worm.density),]
# solo podemos utilizar el signo de menos cuando ordenamos variables num?ricas, si no debemos utilizar la funci?n rank para hacer num?ricos los niveles de los factores.
worms[order(-rank(Vegetation),-Worm.density),]
 
### Un dataframe con nombres en las filas en lugar de n?meros
detach(worms)
worms=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/worms.txt",header=T)
#worms=read.table("c:\\temp\\worms.txt",header=T,row.names=1)
worms
 
##### Crear un dataframe
x=runif(10)
y=letters[1:10]
z=sample(c(rep(T,5),rep(F,5)))
new=data.frame(y,z,x)
new
 
y=rpois(1500,1.5)
table(y)
as.data.frame(table(y))
 
short.frame=as.data.frame(table(y)) # expandir un dataframe
long=as.data.frame(lapply(short.frame, function(x) rep(x, short.frame$Freq)))
long[,1]
 
##### Eliminar filas duplicadas en un Dataframe
dups=read.table("c:\\temp\\dups.txt",header=T)
dups
unique(dups)
dups[duplicated(dups),]
 
##### Juntar dos Dataframes: Merge
lifeform=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/lifeforms.txt",header=T); lifeform
flowering=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/fltimes.txt",header=T); flowering
merge(flowering,lifeforms) #contiene solo las filas que tienen entradas completas para ambos dataframes
(both=merge(flowering,lifeforms,all=T)) # si queremos incluir todas las especies poniendo NA donde no hay coincidencias
# cuando tenemos diferentes nombres en 2 dataframes que queremos unir:
(seeds=read.table("c:\\temp\\seedwts.txt",header=T))
merge(both,seeds,by.x=c("Genus","species"),by.y=c("name1","name2"))
 
##### Juntar dos Dataframes: Match
herb=read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/herbicides.txt",header=T)
herb
worms$hb=herb$Herbicide[match(worms$Vegetation,herb$Type)] ## agregar un herbicida recomendado
match(worms$Vegetation,herb$Type)
worms
 
##### Resumir los contenidos de un dataframe:
#summary: resume los contenidos de todas las variables
#aggregate: crea una tabla del tipo tapply
#by: crea funciones para cada nivel del factor especificado
aggregate(worms[,c(2,3,5,7)],by=list(veg=Vegetation),mean)
aggregate(worms[,c(2,3,5,7)],by=list(veg=Vegetation,d=Damp),mean)
 
x1=1:8; x2=11:18; x3=rep(c("a","i"),each=4); x4=rep(c("f","m"),4)
df=data.frame(y=x1,z=x2,age=x3,sex=x4); df
summary(df)
table(df$age,df$sex) # returns frequencies
table(df$age)
table(df$age)
 
## apply (to a matrix)
x=matrix(1:6, 2,3); x # 2 rows, 3 cols
sum(x) # 21
apply(x,1,sum) # 9 12 sum rows(1)
apply(x,2,sum) # 3, 7 11 sum cols(2)
 
## tapply (table apply)
tapply(df$y,df$age,mean) # tapply(vector,factor,function) for table apply
t=tapply(df$y,list(df$age,df$sex),mean); t
t["a","m"]; t[1,2] # 3 you can reference individual elements
 
# Unir y cortar vectores de un Dataframe
require(stats)
formula(PlantGrowth) # check the default formula
pg = unstack(PlantGrowth) # unstack according to this formula
pg
stack(pg) # now put it back together
stack(pg, select = -ctrl) # omitting one vector
subset(x, cond) #devuelve una selecci?n de x que cumple unas condiciones
split(x, f) #divide el vector o la hora de datos x en grupos definidos por los valores de f
cut(x, breaks, labels=T o F) # divide el rango del vector x en intervalos y codifica los elementos de x de acuerdo con el intervalo en el que caigan
Created by Pretty R at inside-R.org