R es un lenguaje y ambiente para estadistica computacional y gráficas. Es un proyecto GNUsimilar al ambiente y lenguaje S, el cual fue desarrollado en los Laboratorios Bell (anteriormente AT&T, ahora “Lucent Technologies”) por John Chambers y sus colegas. R puede ser considerado como una implementación diferente de S. Hay algunas diferencias importantes, pero la mayoría del código escrito para S corre en R.
R provée una amplia variedad de estadísticas (modalción lineal y nonlineal, pruebas clásicas de estadística, análisis de series de tiempo, clasificación, agrupamiento, ...) y técnicas gráficas, y es altamente extensible. El lenguaje S es amenudo un vehículo a escoger para investigación en metodologías estadísticas, y R es “Open Source”.
Una de la fortalezas de R es la facilidad de realizar gráficas bien diseñadas con calidad de publicación, incluyendo símbolos matemáticos y fórmulas donde se necesiten.
Ejercicios:
x <- rnorm(10)
y <- rnorm(10)
plot(x,y)
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="red", bg="yellow", bty="l", tcl=0.4, main="Como personalizar un grafico en R", las=1, cex=1.5)
opar <- par() ##abre una ventana en R64
par(bg="lightyellow", col.axis="blue", mar=c(4,4,2.5,0.25))
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="red", bg="yellow", bty="l", tcl=-0.25, las=1, cex=1.5) ## pone las figuras, pone ejes, escalas y dibuja fondo
title("Como personalizar un grafico en R (bis)", font.main=3, adj=1) ## pone titulo.
##se puede sobreescribir las 2 ultimas lineas y obtenemos fig5b
par(bg="lightgrey", col.axis="blue", mar=c(4,4,2.5,0.25))
plot(x,y, xlab="Diez numeros al azar", ylab="Otros diez numeros", xlim=c(-2, 2), ylim=c(-2,2), pch=22, col="cyan", bg="pink", bty="l", tcl=-0.25, las=1, cex=1.5)
par(bg="lightgrey", mar=c(2.5, 1.5, 2.5, 0.25))
plot(x,y, type="n", xlab="", ylab="", xlim=c(-2, 2), ylim=c(-2,2), xaxt="n", yaxt="n")
rect(-3, -3, 3, 3, col="cornsilk")
points(x, y, pch=10, col="red", cex=2)
axis(side=1, c(-2, 0, 2), tcl=-0.2, labels=FALSE)
axis(side=2, -1:1, tcl=-0.2, labels=FALSE)
title("Como personalizar un grafico en R (ter)", font.main=4, adj=1, cex.main=1)
mtext("Diez números al azar", side=1, line=1, at=1, cex=0.9, font=3)
mtext("Otros diez números", line=0.5, at=-1.8, cex=0.9, font=3)
mtext(c(-2, 0, 2), side=1, las =1, at=c(-2,0,2), line=0.3, col = "blue", cex=0.9)
mtext(-1:1, side=2, las=1, at=-1:1, line=0.2, col = "blue", cex=0.9)
LATTICE y GRID
LATTICE Y GRID…
library(lattice)
n <- seq (5, 45, 5)
> n
[1] 5 10 15 20 25 30 35 40 45
y <- factor(rep (n,n), labels=paste("n =", n))
> y
[1] n = 5 n = 5 n = 5 n = 5 n = 5 n = 10 n = 10 n = 10 n = 10 n = 10 n = 10
[12] n = 10 n = 10 n = 10 n = 10 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15
[23] n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 15 n = 20 n = 20 n = 20
[34] n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 20
[45] n = 20 n = 20 n = 20 n = 20 n = 20 n = 20 n = 25 n = 25 n = 25 n = 25 n = 25
[56] n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25
[67] n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 25 n = 30 n = 30
[78] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30
[89] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 30
[100] n = 30 n = 30 n = 30 n = 30 n = 30 n = 30 n = 35 n = 35 n = 35 n = 35 n = 35
[111] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35
[122] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35
[133] n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 35 n = 40 n = 40 n = 40
[144] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[155] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[166] n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40 n = 40
[177] n = 40 n = 40 n = 40 n = 40 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[188] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[199] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[210] n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45 n = 45
[221] n = 45 n = 45 n = 45 n = 45 n = 45
Levels: n = 5 n = 10 n = 15 n = 20 n = 25 n = 30 n = 35 n = 40 n = 45
densityplot(~x|y)
data(quakes)
mini <- min(quakes$depth)
maxi <- max(quakes$depth)
int <- ceiling((maxi-mini)/9)
inf <- seq(mini, maxi, int)
quakes$depth.cat <- factor(floor(((quakes$depth - mini) / int)), labels=paste(inf, inf + int, sep="-"))
xyplot(lat ~ long | depth.cat, data = quakes)
data(InsectSprays)
aov.spray <- aov(sqrt(count) ~ spray, data = InsectSprays)
El argumento principal (y obligado) de aov() es una formula que especifica la respuesta en el lado izquierdo del símbolo ~ hy la variable explicativa en el lado derecho. La opción data = InsectSprays indica que las variables deben ser buscadas en el marco de datos InsectSprays. Est sintaxis es equivalente a:
aov.spray <- aov(sqrt(InsectSprays$count) ~ InsectSprays$spray)
o también, si conocemos el numero de las columnas de las variables:
aov.spray <- aov(sqrt(InsectSprays[, 1]) ~ InsectSprays[,2])
Los resultados del análisis no se muestran ya que son asignada a un objeto llamado aov.spray. Usaremos otras funciones para extraer los resultados, como print() para ver un resumen breve del análisis (sobre todo los parámetros) y sumary() para ver mas detalles (incluyendo pruebas estadísticas):
aov.spray
Call:
aov(formula = sqrt(count) ~ spray, data = InsectSprays)
Terms:
spray Residuals
Sum of Squares 88.43787 26.05798
Deg. of Freedom 5 66
Residual standard error: 0.6283453
Estimated effects may be unbalanced
summary(aov.spray)
Df Sum Sq Mean Sq F value Pr(>F)
spray 5 88.44 17.688 44.8 <2e-16 ***
Residuals 66 26.06 0.395
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Recordemos que escribir el nombre del objeto como un comando, es similar al comando print(aov.spray). Para ver una representación gráfica de los resultados podemos usar plot() o termplot(). Antes de escribir plot(aov.spray) dividiremos la ventana gráfica en cuatro partes de tal manera que las cuatro gráficas diagnosticas se dibujan en la misma ventana. Los comandos son:
> opar <- par()
> par(mfcol = c(2,2))
> plot(aov.spray)
y se obtiene la grafica 13.
y con la instrucción
termplot(aov.spray, se=TRUE, partial.resid=TRUE, rug=TRUE)se obtiene el grafico 14
Funciones como aov o lm devuelven una lista donde los diferentes elementos corresponden a los resultados del análisis. Podemos ver la estructura del objeto devuelto por aov() con la instrucción:
str(aov.spray, max.level = -1)
List of 13
- attr(*, "class")= chr [1:2] "aov" "lm"
Otra manera de ver esta estructura es visualizando los nombres del objeto, con la instrucción:
names(aov.spray)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "contrasts" "xlevels" "call" "terms"
[13] "model"
Los elementos se pueden exraer como en cualquier otra lista, con la instrucción:
aov.spray$coefficients
(Intercept) sprayB sprayC sprayD sprayE sprayF
3.7606784 0.1159530 -2.5158217 -1.5963245 -1.9512174 0.2579388
summry() también crea una lista la cual, en el caso de aov(), es simplemente una tabla de pruebas:
str(summary(aov.spray))
List of 1
$ :Classes ‘anova’ and 'data.frame': 2 obs. of 5 variables:
..$ Df : num [1:2] 5 66
..$ Sum Sq : num [1:2] 88.4 26.1
..$ Mean Sq: num [1:2] 17.688 0.395
..$ F value: num [1:2] 44.8 NA
..$ Pr(>F) : num [1:2] 6.33e-20 NA
- attr(*, "class")= chr [1:2] "summary.aov" "listof"
> names(summary(aov.spray))
NULL
Las funciones genéricas también se denominan métodos. Esquematicamente, se construyen como methodfoo, donde foo es la función de análisis. En el caso de summary, podemos ver las funciones que se usan en este método.
apropos("^summary")
[1] "summary" "Summary" "summary.aov"
[4] "summary.aovlist" "summary.connection" "summary.data.frame"
[7] "Summary.data.frame" "summary.Date" "Summary.Date"
[10] "summary.default" "Summary.difftime" "summary.factor"
[13] "Summary.factor" "summary.glm" "summary.infl"
[16] "summary.lm" "summary.manova" "summary.matrix"
[19] "summary.mlm" "Summary.numeric_version" "Summary.ordered"
[22] "summary.POSIXct" "Summary.POSIXct" "summary.POSIXlt"
[25] "Summary.POSIXlt" "summary.srcfile" "summary.srcref"
[28] "summary.stepfun" "summary.table" "summaryRprof"
Podemos ver la diferencia para este método comparando como actúa en una regresión lineal contra un análisis de varianza:
x <- y <- rnorm(5);
mod <- lm(y ~ x)
names(mod)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model"
names(summary(mod))
[1] "call" "terms" "residuals" "coefficients"
[5] "aliased" "sigma" "df" "r.squared"
[9] "adj.r.squared" "fstatistic" "cov.unscaled"
Los objetos devueltos por aov(), lm(), summary(), … son listas, pero no se visualizan como las listas "comunes y corrientes" que hemos visto anteriormente. De hecho, son métodos sprint de estos objetos (recuerde que escribir el nombre de un objeto equivale a usar print()):
apropos("^print")
[1] "print" "print.anova"
[3] "print.AsIs" "print.by"
[5] "print.coefmat" "print.condition"
[7] "print.connection" "print.data.frame"
[9] "print.Date" "print.default"
[11] "print.density" "print.difftime"
[13] "print.DLLInfo" "print.DLLInfoList"
[15] "print.DLLRegisteredRoutines" "print.factor"
[17] "print.family" "print.formula"
[19] "print.ftable" "print.function"
[21] "print.glm" "print.hexmode"
[23] "print.hsearch" "print.infl"
[25] "print.integrate" "print.libraryIQR"
[27] "print.listof" "print.lm"
[29] "print.logLik" "print.NativeRoutineList"
[31] "print.noquote" "print.numeric_version"
[33] "print.octmode" "print.packageInfo"
[35] "print.POSIXct" "print.POSIXlt"
[37] "print.proc_time" "print.restart"
[39] "print.rle" "print.simple.list"
[41] "print.srcfile" "print.srcref"
[43] "print.summary.table" "print.summaryDefault"
[45] "print.table" "print.terms"
[47] "print.ts" "print.warnings"
[49] "printCoefmat"
Todos estos métodos de print permiten cierta visualización dependiendo del análisis.
Funcion que realiza el análisis numérico del modelo de Ricker. Este modelo se usa ampliamente en ecología, particularmente en el estudio demográfico de peces. El objetivo es simular el modelo con respecto a la tasa de crecimiento r y el numero inicial de individuos en la población N_0 (la capacidad de carga K es fijada en 1), los resultados se mostraran como una gráfica del numero de individuos en funcion del tiempo. Agregaremos una opción ver los resultados en los últimos pasos de la simulación (por default los resultados son graficados).
> ricker <- function(nzero, r, K=1, tiempo=100, desde=0, hasta=tiempo)
+ { N<- numeric(tiempo+1)
+ N[1] <- nzero
+ for(i in 1:tiempo) N[i+1] <- N[i]*exp(r*(1-N[i]/K))
+ Tiempo <- 0:tiempo
+ plot(Tiempo, N, type="l", xlim = c(desde, hasta))
+ }
> layout(matrix(1:3, 3, 1))
> ricker(0.1, 1); title("r=1")
> ricker(0.1, 2); title("r=2")
> ricker(0.1, 3); title("r=2")
Bibliografia:
No hay comentarios:
Publicar un comentario