Grafico de las Medias (Escala eje Y)

by Nacho

Si queremos graficar las medias respecto a un factor plotMeans no deja ajustar los ejes. Para poder ajustarlos este script permite crear la función pm.stg (que funciona como plotMeans) que si lo pertmite.

pm.stg<-function (response, factor1, factor2, error.bars = c("se",
"sd", "conf.int", "none"), level = 0.95, xlab =
deparse(substitute(factor1)),
    ylab = paste("mean of", deparse(substitute(response))), legend.lab
= deparse(substitute(factor2)),
    main = "Plot of Means", pch = 1:n.levs.2, lty = 1:n.levs.2,
    col = palette(), ylimits=range(response))
{
    if (!is.numeric(response))
        stop(gettextRcmdr("Argument response must be numeric."))
    xlab
    ylab
    legend.lab
    error.bars <- match.arg(error.bars)
    if (missing(factor2)) {
        if (!is.factor(factor1))
            stop(gettextRcmdr("Argument factor1 must be a factor."))
        valid <- complete.cases(factor1, response)
        factor1 <- factor1[valid]
        response <- response[valid]
        means <- tapply(response, factor1, mean)
        sds <- tapply(response, factor1, sd)
        ns <- tapply(response, factor1, length)
        if (error.bars == "se")
            sds <- sds/sqrt(ns)
        if (error.bars == "conf.int")
            sds <- qt((1 - level)/2, df = ns - 1, lower.tail = FALSE) *
                sds/sqrt(ns)
        sds[is.na(sds)] <- 0
        yrange <- if (error.bars != "none")
            c(min(means - sds, na.rm = TRUE), max(means + sds,
                na.rm = TRUE))
        else range(means, na.rm = TRUE)
        levs <- levels(factor1)
        n.levs <- length(levs)
        plot(c(1, n.levs), yrange, type = "n", xlab = xlab, ylab = ylab,
            axes = FALSE, main = main, ylim=ylimits)
        points(1:n.levs, means, type = "b", pch = 16, cex = 2)
        box()
        axis(2)
        axis(1, at = 1:n.levs, labels = levs)
        if (error.bars != "none")
            arrows(1:n.levs, means - sds, 1:n.levs, means + sds,
                angle = 90, lty = 2, code = 3, length = 0.125)
    }
    else {
        if (!(is.factor(factor1) | is.factor(factor2)))
            stop(gettextRcmdr("Arguments factor1 and factor2 must be factors."))
        valid <- complete.cases(factor1, factor2, response)
        factor1 <- factor1[valid]
        factor2 <- factor2[valid]
        response <- response[valid]
        means <- tapply(response, list(factor1, factor2), mean)
        sds <- tapply(response, list(factor1, factor2), sd)
        ns <- tapply(response, list(factor1, factor2), length)
        if (error.bars == "se")
            sds <- sds/sqrt(ns)
        if (error.bars == "conf.int")
            sds <- qt((1 - level)/2, df = ns - 1, lower.tail = FALSE) *
                sds/sqrt(ns)
        sds[is.na(sds)] <- 0
        yrange <- if (error.bars != "none")
            c(min(means - sds, na.rm = TRUE), max(means + sds,
                na.rm = TRUE))
        else range(means, na.rm = TRUE)
        levs.1 <- levels(factor1)
        levs.2 <- levels(factor2)
        n.levs.1 <- length(levs.1)
        n.levs.2 <- length(levs.2)
        if (length(pch) == 1)
            pch <- rep(pch, n.levs.2)
        if (length(col) == 1)
            col <- rep(col, n.levs.2)
        if (length(lty) == 1)
            lty <- rep(lty, n.levs.2)
        if (n.levs.2 > length(col))
            stop(sprintf(gettextRcmdr("Number of groups for factor2,
%d, exceeds number of distinct colours, %d."),
                n.levs.2, length(col)))
        plot(c(1, n.levs.1 * 1.4), yrange, type = "n", xlab = xlab,
            ylab = ylab, axes = FALSE, main = main, ylim=ylimits)
        box()
        axis(2)
        axis(1, at = 1:n.levs.1, labels = levs.1)
        for (i in 1:n.levs.2) {
            points(1:n.levs.1, means[, i], type = "b", pch = pch[i],
                cex = 2, col = col[i], lty = lty[i])
            if (error.bars != "none")
                arrows(1:n.levs.1, means[, i] - sds[, i], 1:n.levs.1,
                  means[, i] + sds[, i], angle = 90, code = 3,
                  col = col[i], lty = lty[i], length = 0.125)
        }
        x.posn <- n.levs.1 * 1.1
        y.posn <- sum(c(0.1, 0.9) * par("usr")[c(3, 4)])
        text(x.posn, y.posn, legend.lab, adj = c(0, -0.5))
        legend(x.posn, y.posn, levs.2, pch = pch, col = col,
            lty = lty)
    }
    invisible(NULL)
}
 
####################################
 
Esta función, pm.stg, permite ya definir los limites de y:
set.seed(1)
dv <- rnorm(100)
iv1 <- factor(sample(letters[1:2], 100, replace=T))
iv2 <- factor(sample(letters[14:17], 100, replace=T))
 
plotMeans(dv, iv1, iv2, error.bars="se") # los limites de Y los impone R
pm.stg(dv, iv1, iv2, error.bars="se", ylimits=c(-1,1)) # los limites de Y los imponemos nosotros

Gráficos múltiples de barras de error

by InesNaya

Aquí os dejo el código para hacer un gráfico múltiple de barras de error.
Son las estimaciones con sus desviaciones típicas de dos parámetros de la ecuación de crecimiento de von Bertalanffy para erizos de mar del submareal.

Tuve que hacerlo por partes para que los plots de L8 y K me salieran con el mismo ancho.
Aquí teneis los datos originales por si quereis hacer la prueba:
l8_k_estimates.txt
Y aquí va el código:

figure4 <- read.csv("l8_k_estimates.txt",header=T,  sep=" ")
 
pdf("figure4b.pdf", width=8, family="Times") 
op <- par(mfrow=c(1,3), mai=c(0.5,0,0.5,0.2))
 
# Se añaden los labels 
plot.new()
axis(2,at=seq(0.01,0.99, length.out=19), labels=figure4$name, tick= F, las=2, pos=c(1.1,1))
# Se añaden los estimadores de l8
dotchart(figure4$l8.value, pch=20, xlab="L8", color="grey20", xlim=c(min(figure4$l8.value)-max(figure4$l8.std),max(figure4$l8.value)+max(figure4$l8.std)))
segments(figure4$l8.value,1:19,figure4$l8.value+figure4$l8.std,1:19,col="grey20")
segments(figure4$l8.value,1:19,figure4$l8.value-figure4$l8.std,1:19,col="grey20")
# Se añaden los estimadores de k 
dotchart(figure4$k.value, pch=20, xlab="K", color="grey20", xlim=c(min(figure4$k.value)-max(figure4$k.std),max(figure4$k.value)+max(figure4$k.std)))
segments(figure4$k.value,1:19,figure4$k.value+figure4$k.std,1:19,col="grey20")
segments(figure4$k.value,1:19,figure4$k.value-figure4$k.std,1:19,col="grey20")
 
par(op)
dev.off()

Lorenz curves and Gini index

by mariafboan

La Lorenz curve y el Gini Index se usan mucho como indicadores de la concentración de la riqueza (miden el grado de concentración de la renta en un país). En la literatura pesquera ha sido utilizada unas pocas veces para ilustrar la concentración de la abundancia. Nosotros queremos utilizarlos para analizar la concentración del esfuerzo de pesca, lo que en principio no se ha hecho antes.

#PARA HACER CURVAS DE LORENZ E INDICES GINI:
#1. Hace falta tener cargado el package "ineq"
library(ineq)
Lorenz <- read.table(file = "C:\\R.gps\\Lorenz.txt",header = TRUE)
names(Lorenz)
summary(Lorenz)
 
#Loop que hizo Inés para hacer las curvas y estimar los indices de las 31 patches de erizo
#y guardar los primeros en un pdf (sin el loop tendría que hacer un código para cada uno de los 31 patches):
giniIndex <- NULL
pdf("GraficosLorentz.pdf", width=11.69, height=8.26)
op <- par(mfrow=c(3,7))
for(i in 1:31){
Mtemp <-  subset(Lorenz,Mancha==as.character(i), select=4)#Columna4=n.celda
plot(Lc(Mtemp[,1]), main=paste("Mancha", i, sep=" "))
giniIndex[i] <- ineq(Mtemp[,1])
}
par(op)
dev.off()
 
#Para guardar los índices:
write(giniIndex,file="clipboard",ncolumns=1)
#o
write.csv(giniIndex,file="ginis.csv")
#Resultó mejor la 2º pq da los datos en una unica columna con el nº 
#de patch asociado
 
##FIGURA DE LOS 3 PATCHES MÁS VISITADOS (7,15Y22):
 
M7<-subset(Lorenz,Mancha=="7")
M15<-subset(Lorenz,Mancha=="15")
M22<-subset(Lorenz,Mancha=="22")
Lc.7<-Lc(M7$n.celda)
Lc.15<-Lc(M15$n.celda)
Lc.22<-Lc(M22$n.celda)
plot(Lc.7,ylab="",xlab="",main="")
lines(Lc.15,lty="dashed")
lines(Lc.22,lty="dotted")

Predictions and/or confidence (or prediction) intervals on predictions

by MJuanJorda

Los intervalos de confianza de la funcion predict nos han dado muchos quebraderos de cabeza!
El codigo calcula los intervalos de confianza de los fitted values de los predict !!

library(nlme) 
fm1 <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject, data = Orthodont) 
 
plot(Orthodont)
newdat <- expand.grid(age=c(8,10,12,14), Sex=c("Male","Female")) 
 
newdat$pred <- predict(fm1, newdat, level = 0)
 
Designmat <- model.matrix(eval(eval(fm1$call$fixed)[-2]), newdat[-3]) 
predvar <- diag(Designmat %*% fm1$varFix %*% t(Designmat)) 
newdat$SE <- sqrt(predvar) 
 
newdat$SE2 <- sqrt(predvar+fm1$sigma^2)
 
library(ggplot2) 
pd <- position_dodge(width=0.4) 
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+ 
   geom_point(position=pd)+ 
  geom_linerange(aes(ymin=pred-2*SE,ymax=pred+2*SE), position=pd)
 
 
## prediction intervals 
ggplot(newdat,aes(x=age,y=pred,colour=Sex))+ 
   geom_point(position=pd)+ 
   geom_linerange(aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)

Fuente: Wikidot

Plots de Linea de Costa

by Nacho

Hay webs como http://rimmer.ngdc.noaa.gov/ que nos permiten extraer los datos de lat y long de la linea de costa de todo el mundo a escala 1:250000. Pudiendo seleccionar con zoom aquella que queremos con exactitud (para evitar archivos con muchos datos y dificiles de manejar). Se selecciona en formato Splus y nos da un archivo .dat que abrimos y guardamos en formato .txt.

Galicia<- read.table("Galicia.txt")
names(Galicia) <- c("Longitud","Latitud")
# Usando plot() la costa sale desproporcionada
plot(Galicia$Long, Galicia$Lat, type="l", xlab="Longitud",ylab="Latitud")
title(main="Galicia")
# Usar mejor "eqscplot()"
eqscplot(Galicia$Long,Galicia$Lat,type="l",xlab="Longitud",ylab="Latitud",tol=-0.10)
#tol premite ajustar los margenes para encajar el mapa
title(main="Galicia")

*Es necesario instalar y cargar el paquete MASS

DatosCoordenadas Galicia

Gráficos en R: Chuletas para ponerlos bonitos

by InesNaya

A la hora de poner bonitos los gráficos en R es buena idea tener a mano qué colores y símbolos podemos utilizar para hacerlo.

Os dejo aquí la lista de colores y la de símbolos con el script para que las podais crear vosotros mismos:

COLORES EN R
ColoursInR
Descargar versión en pdf

# Colores en R
pdf("ColoursInR.pdf", width=11.7, height=8.5)
plot.new()
legend(x="topleft",colors(), pch=15, col=colors(), bty="n", ncol=10, cex=0.5, pt.cex=1, inset=0)
dev.off()

SÍMBOLOS EN R (pch)
SymbolsInR
Descargar versión en pdf

# Símbolos en R
pdf("SymbolsInR.pdf", width=5.5, height=4.5)
plot.new()
legend(x="topleft", as.character(0:25), pch=c(0:25), bty="p", ncol=4, cex=1.5, pt.cex=1, inset=0)
dev.off()

Fuente original: aquí.

Leyenda en multiples plots

by Nacho

Si haceis multiples plots y quereis colocar la leyenda debajo de todos ellos o a un lado. Primero, podeis determinar los margenes o aumentar las columnas o filas con par(mfrow=c(3,2)). El comando que hay que añadir es xpd=TRUE

opar <-par(mar = c(20, 4, 4, 4)) #Con par guardamos las opciones gráficas que queremos
#en este caso  aumentamos los  margenes
plot(1:10)
lines(1:10)
par(xpd=TRUE)  #xpd para poder usar los margenes
legend(locator (1),lty=1, col="black", legend="straigh line")
#locator nos  permite colocar una o varias leyendas donde le  señalemos con el ratón
par(opar) #restauramos las opciones de graficos iniciales