Grupo R : RMyP

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

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)