Posts Tagged ‘barras de error’

Grafico de las Medias (Escala eje Y)

Thursday, June 17th, 2010

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

Tuesday, June 8th, 2010

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()