PK %B—C²ßDCÊ/ Ê/ # aprenderhaciendo-latest/inicio.html
Contenidos:
1 2 3 | # free memory
rm(list = ls())
gc()
|
1 2 3 4 | a <- ts(1:30, frequency=12, start=c(2011,3))
print(a)
str(a)
attributes(a)
|
1 | plot(AirPassengers)
|
1 2 3 4 5 6 7 8 9 10 11 | # decompose time series
apts <- ts(AirPassengers, frequency=12)
f <- decompose(apts)
# seasonal figures
f$figure
plot(f$figure, type="b", xaxt="n", xlab="")
# get names of 12 months in English words
monthNames <- months(ISOdate(2011,1:12,1))
# label x-axis with month names
# las is set to 2 for vertical label orientation
axis(1, at=1:12, labels=monthNames, las=2)
|
1 | plot(f)
|
1 2 3 4 5 6 7 8 | fit <- arima(AirPassengers, order=c(1,0,0), list(order=c(2,1,0), period=12))
fore <- predict(fit, n.ahead=24)
# error bounds at 95% confidence level
U <- fore$pred + 2*fore$se
L <- fore$pred - 2*fore$se
ts.plot(AirPassengers, fore$pred, U, L, col=c(1,2,4,4), lty = c(1,1,2,2))
legend("topleft", c("Actual", "Forecast", "Error Bounds (95% Confidence)"),
col=c(1,2,4), lty=c(1,1,2))
|
1 2 3 4 5 6 | library(dtw)
idx <- seq(0, 2*pi, len=100)
a <- sin(idx) + runif(100)/10
b <- cos(idx)
align <- dtw(a, b, step=asymmetricP1, keep=T)
dtwPlotTwoWay(align)
|
1 2 3 4 5 | sc <- read.table("./data/synthetic_control.data", header=F, sep="")
# show one sample from each class
idx <- c(1,101,201,301,401,501)
sample1 <- t(sc[idx,])
plot.ts(sample1, main="")
|
1 | set.seed(6218)
|
1 2 3 4 5 6 7 8 9 10 11 12 | n <- 10
s <- sample(1:100, n)
idx <- c(s, 100+s, 200+s, 300+s, 400+s, 500+s)
sample2 <- sc[idx,]
observedLabels <- rep(1:6, each=n)
# hierarchical clustering with Euclidean distance
hc <- hclust(dist(sample2), method="average")
plot(hc, labels=observedLabels, main="")
# cut tree to get 6 clusters
rect.hclust(hc, k=6)
memb <- cutree(hc, k=6)
table(observedLabels, memb)
|
1 2 3 4 5 6 7 8 | library(dtw)
distMatrix <- dist(sample2, method="DTW")
hc <- hclust(distMatrix, method="average")
plot(hc, labels=observedLabels, main="")
# cut tree to get 6 clusters
rect.hclust(hc, k=6)
memb <- cutree(hc, k=6)
table(observedLabels, memb)
|
1 2 3 4 5 6 7 8 9 10 | classId <- rep(as.character(1:6), each=100)
newSc <- data.frame(cbind(classId, sc))
library(party)
ct <- ctree(classId ~ ., data=newSc,
controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5))
pClassId <- predict(ct)
table(classId, pClassId)
# accuracy
(sum(classId==pClassId)) / nrow(sc)
plot(ct, ip_args=list(pval=FALSE), ep_args=list(digits=0))
|
1 2 3 4 5 6 7 8 9 | library(wavelets)
wtData <- NULL
for (i in 1:nrow(sc)) {
a <- t(sc[i,])
wt <- dwt(a, filter="haar", boundary="periodic")
wtData <- rbind(wtData, unlist(c(wt@W, wt@V[[wt@level]])))
}
wtData <- as.data.frame(wtData)
wtSc <- data.frame(cbind(classId, wtData))
|
1 2 3 4 5 6 7 | # build a decision tree with DWT coefficients
ct <- ctree(classId ~ ., data=wtSc,
controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5))
pClassId <- predict(ct)
table(classId, pClassId)
(sum(classId==pClassId)) / nrow(wtSc)
plot(ct, ip_args=list(pval=FALSE), ep_args=list(digits=0))
|
1 2 | # fix seed to get a fixed result in the chunk below
set.seed(100)
|
1 2 3 4 5 6 7 | k <- 20
# create a new time series by adding noise to time series 501
newTS <- sc[501,] + runif(100)*15
distances <- dist(newTS, sc, method="DTW")
s <- sort(as.vector(distances), index.return=TRUE)
# class IDs of k nearest neighbors
table(classId[s$ix[1:k]])
|
1 2 3 | # free memory
rm(list = ls())
gc()
|
1 | set.seed(8953)
|
1 2 3 | iris2 <- iris
iris2$Species <- NULL
(kmeans.result <- kmeans(iris2, 3))
|
1 | table(iris$Species, kmeans.result$cluster)
|
1 2 3 4 | plot(iris2[c("Sepal.Length", "Sepal.Width")], col = kmeans.result$cluster)
# plot cluster centers
points(kmeans.result$centers[,c("Sepal.Length", "Sepal.Width")], col = 1:3,
pch = 8, cex=2)
|
1 2 3 4 5 6 | library(fpc)
pamk.result <- pamk(iris2)
# number of clusters
pamk.result$nc
# check clustering against actual species
table(pamk.result$pamobject$clustering, iris$Species)
|
1 2 3 | layout(matrix(c(1,2),1,2)) # 2 graphs per page
plot(pamk.result$pamobject)
layout(matrix(1)) # change back to one graph per page
|
1 2 | pam.result <- pam(iris2, 3)
table(pam.result$clustering, iris$Species)
|
1 2 3 | layout(matrix(c(1,2),1,2)) # 2 graphs per page
plot(pam.result)
layout(matrix(1)) # change back to one graph per page
|
1 | set.seed(2835)
|
1 2 3 4 | idx <- sample(1:dim(iris)[1], 40)
irisSample <- iris[idx,]
irisSample$Species <- NULL
hc <- hclust(dist(irisSample), method="ave")
|
1 2 3 4 | plot(hc, hang = -1, labels=iris$Species[idx])
# cut tree into 3 clusters
rect.hclust(hc, k=3)
groups <- cutree(hc, k=3)
|
1 2 3 4 5 | library(fpc)
iris2 <- iris[-5] # remove class tags
ds <- dbscan(iris2, eps=0.42, MinPts=5)
# compare clusters with original class labels
table(ds$cluster, iris$Species)
|
1 | plot(ds, iris2)
|
1 | plot(ds, iris2[c(1,4)])
|
1 | plotcluster(iris2, ds$cluster)
|
1 2 3 4 5 6 7 8 9 10 11 12 | # create a new dataset for labeling
set.seed(435)
idx <- sample(1:nrow(iris), 10)
newData <- iris[idx,-5]
newData <- newData + matrix(runif(10*4, min=0, max=0.2), nrow=10, ncol=4)
# label new data
myPred <- predict(ds, iris2, newData)
# plot result
plot(iris2[c(1,4)], col=1+ds$cluster)
points(newData[c(1,4)], pch="*", col=1+myPred, cex=3)
# check cluster labels
table(myPred, iris$Species[idx])
|
This post examines the features of R Markdown using knitr in Rstudio 0.96. This combination of tools provides an exciting improvement in usability for reproducible analysis. Specifically, this post
#. provides a basic example of producing console output and plots using R Markdown; #. highlights several code chunk options such as caching and controlling how input and output is displayed; #. demonstrates use of standard Markdown notation as well as the extended features of formulas and tables; and #. discusses the implications of R Markdown.
This post was produced with R Markdown. The source code is available here as a gist. The post may be most useful if the source code and displayed post are viewed side by side. In some instances, I include a copy of the R Markdown in the displayed HTML, but most of the time I assume you are reading the source and post side by side.
To work with R Markdown, if necessary:
Install R Install the lastest version of RStudio (at time of posting, this is 0.96) Install the latest version of the knitr package: install.packages(“knitr”)
To run the basic working example that produced this blog post:
Open R Studio, and go to File - New - R Markdown
If necessary install ggplot2 and lattice packages: install.packages(“ggplot2”); install.packages(“lattice”)
Paste in the contents of this gist (which contains the R Markdown file used to produce this post) and save the file with an .rmd extension
Click Knit HTML
`{r } set.seed(1234) library(ggplot2) library(lattice) `
## Basic console output To insert an R code chunk, you can type it manually or just press Chunks - Insert chunks or use the shortcut key. This will produce the following code chunk:
Pressing tab when inside the braces will bring up code chunk options.
The following R code chunk labelled basicconsole is as follows:
`{r } x <- 1:10 y <- round(rnorm(10, x, 1), 2) df <- data.frame(x, y) df `
The code chunk input and output is then displayed as follows:
`{r basicconsole} x <- 1:10 y <- round(rnorm(10, x, 1), 2) df <- data.frame(x, y) df `
## Plots Images generated by knitr are saved in a figures folder. However, they also appear to be represented in the HTML output using a [data URI scheme]( http://en.wikipedia.org/wiki/Data_URI_scheme). This means that you can paste the HTML into a blog post or discussion forum and you don’t have to worry about finding a place to store the images; they’re embedded in the HTML.
### Simple plot Here is a basic plot using base graphics:
`{r } plot(x) `
`{r simpleplot} plot(x) `
Note that unlike traditional Sweave, there is no need to write fig=TRUE.
### Multiple plots Also, unlike traditional Sweave, you can include multiple plots in one code chunk:
`{r } boxplot(1:10~rep(1:2,5)) plot(x, y) ``{r multipleplots} boxplot(1:10~rep(1:2,5)) plot(x, y) `
### ggplot2 plot Ggplot2 plots work well:
`{r ggplot2ex} qplot(x, y, data=df) `
### lattice plot As do lattice plots:
`{r latticeex} xyplot(y~x) `
Note that unlike traditional Sweave, there is no need to print lattice plots directly.
## R Code chunk features ### Create Markdown code from R The following code hides the command input (i.e., echo=FALSE), and outputs the content directly as code (i.e., results=asis, which is similar to results=tex in Sweave).
`{r , results='asis', echo=FALSE} cat("Here are some dot points\n\n") cat(paste("* The value of y[", 1:3, "] is ", y[1:3], sep="", collapse="\n")) ``{r dotpointprint, results='asis', echo=FALSE} cat("Here are some dot points\n\n") cat(paste("* The value of y[", 1:3, "] is ", y[1:3], sep="", collapse="\n")) `
`{r , results='asis', echo=FALSE} cat("x | y", "--- | ---", sep="\n") cat(apply(df, 1, function(X) paste(X, collapse=" | ")), sep = "\n") `
`{r createtable, results='asis', echo=FALSE} cat("x | y", "--- | ---", sep="\n") cat(apply(df, 1, function(X) paste(X, collapse=" | ")), sep = "\n") `
### Control output display The folllowing code supresses display of R input commands (i.e., echo=FALSE) and removes any preceding text from console output (comment=”“; the default is comment=”##”).
`{r echo=FALSE, comment="", echo=FALSE} head(df) `
`{r echo=FALSE, comment="", echo=FALSE} head(df) `
### Control figure size The following is an example of a smaller figure using fig.width and fig.height options.
`{r , fig.width=3, fig.height=3} plot(x) `
`{r smallplot, fig.width=3, fig.height=3} plot(x) `
### Cache analysis Caching analyses is straightforward. Here’s example code. On the first run on my computer, this took about 10 seconds. On subsequent runs, this code was not run.
If you want to rerun cached code chunks, just [delete the contents of the cache folder](http://stackoverflow.com/a/10629121/180892)
## Basic markdown functionality For those not familiar with standard [Markdown](http://daringfireball.net/projects/markdown/), the following may be useful. See the source code for how to produce such points. However, RStudio does include a Markdown quick reference button that adequatly covers this material.
### Dot Points Simple dot points:
and numeric dot points:
and nested dot points:
### Equations Equations are included by using LaTeX notation and including them either between single dollar signs (inline equations) or double dollar signs (displayed equations). If you hang around the Q&A site [CrossValidated](http://stats.stackexchange.com) you’ll be familiar with this idea.
There are inline equations such as $y_i = alpha + beta x_i + e_i$.
And displayed formulas:
$$frac{1}{1+exp(-x)}$$
knitr provides self-contained HTML code that calls a Mathjax script to display formulas. However, in order to include the script in my blog posts I [took the script](https://gist.github.com/2716053) and incorporated it into my blogger template. If you are viewing this post through syndication or an RSS reader, this may not work. You may need to view this post on my website.
### Tables Tables can be included using the following notation
A | B | C — | — | — 1 | Male | Blue 2 | Female | Pink
### Hyperlinks
### Images Here’s an example image:
![image from redmond barry building unimelb](http://i.imgur.com/RVNmr.jpg)
### Code Here is Markdown R code chunk displayed as code:
`{r} x <- 1:10 x `
And then there’s inline code such as x <- 1:10.
Hay un mito en la tecnociencia según la cual las herramientas “no son buenas ni malas, sino dependen de como se utilicen”. Sin embargo, las herramientas en efecto condicionan el quehacer, e incluso el sentido de la investigación cientÃfica. Este curso pretende enriquecer el abanico de las alternativas disponibles.
Se podrÃa intentar ser polÃticamente correcto y decir que queda a juicio de cada quién decir si R es mejor o peor que otras herramientas estadÃsticas. Sin embargo, desde el punto de vista de cualquier investigador riguroso el uso de R mejora con creces las funcionalidades que ofrece una hoja de cálculo, y en términos éticos y prácticos trabajar con software libre y de código abierto, basado en estándares abiertos, es sin duda una opción superior.
Muchos profesionales de la informática una vez graduados evaden las tareas de programación, de manera similar a como muchas personas evaden las matemáticas tanto como pueden. Esto se debe a que tanto una como la otra son, por lo general, muy mal enseñadas. En la actualidad, la programación es una destreza básica y fundamental al alcance de cualquier persona.
Estas notas buscan son una presentación de R basada en ejemplos, en lugar de hacer repaso pormenorizado de todos los conceptos y elementos del entorno R, se trabajará en base a demostraciones del flujo de trabajo usual para muchos usuarios de R con lo cuál se espera que cada persona gane comprensión de los retos y oportunidades que ofrece el uso de una herramienta de este tipo.
R es un lenguaje de programación, que se ejecuta sobre un entorno que ofrece capacidades de visualización y una gran cantidad de funciones y métodos de estadÃstica computacional, y mecanismos para extender sus funcionalidades. Por lo tanto, el correcto manejo de R implica desarrollar destrezas de programación.
R cuenta con el conjunto más extenso de métodos estadÃsticos disponibles en cualquier software estadÃstico, ya sea libre o privativo, además de potentes y flexibles capacidades de graficación. Sin embargo, el diseño de R en tanto que lenguaje de programación, presenta debilidades que han motivado crÃticas incluso por sus usuarios más fieles, un ejemplo de esto es el libro The R inferno.
La versión de R para Windows se puede conseguir desde: http://cran.r-project.org/bin/windows/base/
Descargar e instalar el archivo R-X.X.X-win.exe desde el enlace donde dice: Download R X.X.X for Windows, donde X.X.X es la versión de R más reciente.
Generalmente la versión disponible en los repositorios de las distintas distribuciones está algo desactualizada, esto no es conveniente porque se pueden tener incompatibilidad con los paquetes propios de R.
Para resolver esto se configuran los repositorios según la información que aparece en: http://cran.r-project.org/bin/linux/
Se generan paquetes de las últimas versiones de R para Debian, Redhat, SUSE y Ubuntu, si utiliza una distribución que no es compatible con ninguno de estos repositorios, puede compilar la última versión de R desde las fuentes http://cran.r-project.org/sources.html
R También está empaquetado para MacOSX Snow Leonard y superior
RStudio es un entorno de desarrollo integrado para R, disponible para Linux (Debian/Ubuntu/Fedora/OpenSUSE), Windows y MacOSX.
Es un entorno diseñado para programación y en este sentido es lo mejor disponible para R, aparte de ser de código abierto por lo que está siendo continuamente mejorado.
Existen otros editores orientados a usuarios que quieren aplicar técnicas estadÃsticas a datos, este enfoque termina siendo limitado porque con mucha frecuencia se requiere modificar la estructura de los datos, experimentar con las técnicas, aplicar las técnicas de manera iterativa, generar reportes automáticos o configurar en detalle los graficos, en cualquiera de estos casos es preferible un entorno de programación como RStudio.
Desde el sitio: http://rstudio.org/download/desktop Descargar e instalar la versión más reciente
Ofrece cuatro zonas personalizables: Editor, Consola, Datos y Ayuda.
Ilustración 1: Interfaz de RStudio
La primera vez que se ejecuta el programa la zona del Editor aparece desactivada, se activará en el momento que comience a editar un script de R.
Se comienza con el modo interactivo (en la consola) haciendo pruebas y utilizándolo como manual de ayuda.
La serie de comandos que devuelven los resultados esperados se guardan y modifican (en el editor) como scripts para crear rutinas reproducibles y automatizables.
En el espacio Datos aparecen los objetos del espacio de trabajo actual, y el histórico de comandos utilizados.
En el espacio Ayuda se encuentran por defecto el explorador de archivos, el gestor de gráficos, el gestor de paquetes, y la interfaz de ayuda.
Todos los espacios están conectados y los cambios realizados en un espacio con frecuencia afectan a los demás.
Se recomienda que cada proyecto tenga una estructura de archivos con las carpetas: (puede traducir los nombres si lo desea): data, documents, graphics , images, notes, code. También puede ser otra estructura adaptada a las preferencias del investigador, cualquiera es mejor que tener todos los archivos sueltos en una misma carpeta.
Ilustración 2: Estructura de un proyecto de análisis de datos
También es recomendable, si se utiliza RStudio, utilizar proyectos con el fin de gestionar los entornos de trabajo de manera que cada uno esté asociado a un espacio de trabajo, un historial, y un código bien diferenciado.
Empecemos por abrir en RStudio el archivo simpleRprogram.R.
A continuación se va a ejecutar el programa lÃnea por lÃnea usando el botón Run que está en la parte superior derecha del editor.
1 2 3 4 5 6 7 8 9 10 11 12 13 | setwd("[ruta a]/myRfolder/")
mydata <- read.csv("mydata.csv")
mydata
mydata$workshop <- factor(mydata$workshop)
summary(mydata)
plot(mydata$q1, mydata$q4)
myModel <- lm(q4 ~ q1 + q2 + q3, data=mydata)
summary(myModel)
anova(myModel)
plot(myModel)
|
La función setwd() (establecer directorio de trabajo) indica a R cual va a ser el directorio que va a tomar como referencia para acceder a los archivos.
La expresión mydata <- read.csv("mydata.csv") utiliza la función read.csv() para leer el archivo separado por comas mydata.csv. El resultado que es un objeto dataframe lo asigna mediante el operador <- a la variable mydata.
Un objeto dataframe es una estructura de datos de tipo tabla en la que cada columna es un campo y cada fila un registro de forma similar a como se almacenan los datos en una base de datos o en una hoja de cálculo.
Como estamos ejecutando el archivo lÃnea por lÃnea, y por lo tanto estamos en el modo interactivo al ejecutar la lÃnea 3 donde aparece mydata mostrará en la consola el contenido del objeto.
1 2 3 4 5 6 7 8 9 10 | > mydata
workshop gender q1 q2 q3 q4
1 1 f 1 1 5 1
2 2 f 2 1 4 1
3 1 f 2 2 4 3
4 2 3 1 NA 3
5 1 m 4 5 2 4
6 2 m 5 4 5 5
7 1 m 5 3 4 4
8 2 m 4 5 5 5
|
Al leer del archivo el campo workshop se considera de tipo numérico. Para que R lo interprete como un campo categórico se utiliza la función factor().
En la siguiente lÃnea, la función summary() devuelve las estadÃsticas descriptivas básicas de cada uno de los campos. Nótese como el campo Workshop muestra un conteo de la ocurrencia de cada categorÃa “1” y “2”.
La función plot() genera un gráfico de dispersión entre los valores en los campos q1 y q4, para acceder a estos valores se utiliza la notación mydata$q1 que hace referencia a los valores del campo q1 como un vector numérico.
A continuación, se quiere construir un modelo lineal, en este caso como la respuesta q4 se explica mediante las contribuciones de q1, q2 y q3. Para lo cual se escribe la fórmula q4 ~ q1 + q2 + q3. El argumento data=mydata indica a R de cual conjunto de datos toma las variables.
El resultado se asigna a la variable myModel en este caso es un objeto de la clase lm. Esta es una estructura que contiene los componentes generados durante el ajuste: coefficients, residuals, fitted.values, entre otros.
A este modelo lineal se le puede aplicar un análisis de varianza mediante la función anova(), y finalmente se genera un conjunto de gráficos predefinidos para el modelo lineal con plot(myModel). Hay muchas funciones como plot() y summary() que responden de forma distinta dependiendo del argumento de entrada.
1 | summary(iris)
|
1 2 | quantile(iris$Sepal.Length)
quantile(iris$Sepal.Length, c(.1, .3, .65))
|
1 2 | var(iris$Sepal.Length)
hist(iris$Sepal.Length)
|
1 | plot(density(iris$Sepal.Length))
|
1 2 | table(iris$Species)
pie(table(iris$Species))
|
1 | barplot(table(iris$Species))
|
1 2 3 4 | cov(iris$Sepal.Length, iris$Petal.Length)
cov(iris[,1:4])
cor(iris$Sepal.Length, iris$Petal.Length)
cor(iris[,1:4])
|
1 | aggregate(Sepal.Length ~ Species, summary, data=iris)
|
1 | boxplot(Sepal.Length~Species, data=iris)
|
1 2 | with(iris, plot(Sepal.Length, Sepal.Width,
col=Species, pch=as.numeric(Species)))
|
1 | plot(jitter(iris$Sepal.Length), jitter(iris$Sepal.Width))
|
1 | pairs(iris)
|
1 2 | library(scatterplot3d)
scatterplot3d(iris$Petal.Width, iris$Sepal.Length, iris$Sepal.Width)
|
1 2 | ## library(rgl)
## plot3d(iris$Petal.Width, iris$Sepal.Length, iris$Sepal.Width)
|
1 2 | distMatrix <- as.matrix(dist(iris[,1:4]))
heatmap(distMatrix)
|
1 2 3 | library(lattice)
levelplot(Petal.Width~Sepal.Length*Sepal.Width, iris, cuts=9,
col.regions=grey.colors(10)[10:1])
|
1 2 | filled.contour(volcano, color=terrain.colors, asp=1,
plot.axes=contour(volcano, add=T))
|
1 | persp(volcano, theta=25, phi=30, expand=0.5, col="lightblue")
|
1 2 | library(MASS)
parcoord(iris[1:4], col=iris$Species)
|
1 2 | library(lattice)
parallelplot(~iris[1:4] | Species, data=iris)
|
1 2 | library(ggplot2)
qplot(Sepal.Length, Sepal.Width, data=iris, facets=Species ~.)
|
1 2 3 4 5 6 7 8 9 10 11 | ## # save as a PDF file
## pdf("myPlot.pdf")
## x <- 1:50
## plot(x, log(x))
## graphics.off()
## #
## # Save as a postscript file
## postscript("myPlot2.ps")
## x <- -20:20
## plot(x, x^2)
## graphics.off()
|