effectsplot.r

effectsplot = function(formula, data=NULL, main=NULL, xlab=NULL, ylab=NULL) {
  # author: Matthias C. M. Troffaes
  # date: 26 Nov 2009
  # license: GPLv3
  # usage:
  #
  #   effectsplot(y ~ x, data=..., main=..., xlab=..., ylab=...)
  #
  # where the last three arguments are optional, and you should mind that xlab
  # is a *vector* of two strings (first for group effect label, second for
  # residual label).


  # extract data frame from formula and data
  ep.frame = model.frame(formula, data=data)
  # set up labels
  if (is.null(main)) {
    main = "Effects Plot"
  }
  if (is.null(xlab)) {
    xlab = c(sprintf("Group Effect (%s)", names(ep.frame)[2]), "Residuals")
  }
  if (is.null(ylab)) {
    ylab = names(ep.frame)[1]
  }
  # get groups
  ep.groups = levels(as.factor(ep.frame[,2]))
  ep.num.groups = length(ep.groups)
  # get global effect
  ep.global.effect = mean(ep.frame[,1])
  # get group effects
  ep.group.effects = rep(0, length(ep.groups))
  for (i in 1:ep.num.groups) {
    ep.group.effects[i] = mean(ep.frame[,1][ep.frame[,2] == ep.groups[i]]) - ep.global.effect
  }
  # get residuals: we do this by copying data, substracting global
  # effect, and substracting group effects
  ep.residuals = ep.frame[,1] - ep.global.effect
  for (i in 1:ep.num.groups) {
    ep.residuals = ep.residuals - ep.group.effects[i] * (ep.frame[,2] == ep.groups[i])
  }

  # get plot y limits
  ymin = min(ep.residuals, ep.group.effects)
  ymax = max(ep.residuals, ep.group.effects)
  
  # plot group effects at x=1
  plot(rep(1, ep.num.groups), ep.group.effects, xlim=c(0.5,2.5), ylim=c(ymin, ymax), xaxp=c(1, 2, 1), main=main, xlab="", ylab=ylab, xaxt="n")

  # plot residuals at x=2
  boxplot(ep.residuals, add=TRUE, at=2)

  # annotate the x axis
  axis(1, at=c(1, 2), labels=xlab)

  # annotate group effects
  for (i in 1:ep.num.groups) {
    text(1, ep.group.effects[i], ep.groups[i], pos=2)
  }
}