Re: Equation and R2 on graph

399 views
Skip to first unread message

Winston Chang

unread,
Jan 18, 2013, 1:28:58 PM1/18/13
to Oscar Branson, ggplot2
Hi Oscar - 

Here's another way of doing it, for linear models. Maybe you can find some ideas from it. In the code below, I show two methods: one with manually-created expressions, and another with expressions automatically generated from the lm object.


# Given a model, predict values of yvar from xvar
# This supports one predictor and one predicted variable
# xrange: If NULL, determine the x range from the model object. If a vector with
# two numbers, use those as the min and max of the prediction range.
# samples: Number of samples across the x range.
# ...: Further arguments to be passed to predict()
predictvals <- function(model, xvar, yvar, xrange=NULL, samples=100, ...) {
  # If xrange isn't passed in, determine xrange from the models.
  # Different ways of extracting the x range, depending on model type
  if (is.null(xrange)) {
    if (any(class(model) %in% c("lm", "glm")))
      xrange <- range(model$model[[xvar]])
    else if (any(class(model) %in% "loess"))
      xrange <- range(model$x)
  }
  newdata <- data.frame(x = seq(xrange[1], xrange[2], length.out = samples))
  names(newdata) <- xvar
  newdata[[yvar]] <- predict(model, newdata = newdata, ...)
  newdata
}



model <- lm(mpg ~ wt, mtcars)

# First generate prediction data
pred <- predictvals(model, "wt", "mpg")
sp <- ggplot(mtcars, aes(x=wt, y=mpg)) + geom_point() + geom_line(data=pred)

# === Plot with manual annotation ===
# Use x=Inf and y=Inf to place in corner, then use hjust and vjust to adjust
# position of text label inside plot area.
sp + annotate("text", label="r^2 == 0.75", parse = TRUE, x=Inf, y=Inf, hjust=1.2, vjust=1.5)


# =================================
# Auto-generate expression from model object
eqn <- as.character(as.expression(
  substitute(italic(y) == a + b * italic(x) * "," ~~ italic(r)^2 ~ "=" ~ r2,
    list(a = format(coef(model)[1], digits=3),
         b = format(coef(model)[2], digits=3),
         r2 = format(summary(model)$r.squared, digits=2)
))))
eqn
# "italic(y) == \"37.3\" + \"-5.34\" * italic(x) * \",\" ~ ~italic(r)^2 ~ \"=\" ~ \"0.75\""


# === Plot with generated expression ===
sp + annotate("text", label=eqn, parse=TRUE, x=Inf, y=Inf, hjust=1.1, vjust=1.5)


-Winston





On Fri, Jan 18, 2013 at 10:49 AM, Oscar Branson <ob...@cam.ac.uk> wrote:
I often find that I want to place the equation of a line drawn with stat_smooth or geom_line on a ggplot graph.  As far as I can tell, there is no standard way to do this.  Granted, there is the problem that the fit functions are very variable, so it's hard to format the output correctly in all cases, but there must be some way around this...

I am no expert at R coding, and have come up with a rather 'messy' solution, that combines fragments of strings together to produce a character string that can be parsed by ggplot's annotate command.  At the moment, it can deal with linera (y=mx+c) and polynomial (ns(x,order), from the splines package) functions.

I realise that my solution is rather crude, but I am not aware of a better way of going about it.  Any thoughts/comments/suggestions would be very welcome!

Thanks!

......Code follows......

set.seed(1000)

ts=data.frame(x=seq(0,10,0.2),y=seq(0,10,.2)^2-0.1*seq(0,10,.2)^3+rnorm(length(seq(0,10,.2)),0,0.4))

# have a look at the data - I want to put the equation of the line, and the R2 value at the top of the graph.
p=ggplot(ts) + geom_point(aes(x=x,y=y)) + stat_smooth(aes(x=x,y=y),method='lm',formula=y~ns(x,3))
p

# Equation for getting the coefficients and R2 values from a fit.  Parameters are
# df = data frame
# x.var = string denoting the x data in the df
# y.var = string denoting the y data in the df
# signif.figs = the number of significant figures you want the coefficients and R2 to.
# eq.plot = T/F, whether to plot the equation, or just the R2.  F = just R2
# model.type = 'lin' or 'poly', for the two types of model I have designed it for
# order = if using 'poly', the order of the polynomial.

lm_eqn = function(df,x.var,y.var,signif.figs,eq.plot=T,model.type,order){
  if(missing(x.var) | missing(y.var) | class(x.var)!='character' | class(y.var)!='character') stop('x.var and y.var must be the names of the columns you want to use as x and y as a character string.' )
  if(missing(model.type)) stop("model.type must be 'lin' (linear y~x model) or 'poly' (polynomial y~ns(x,order) model, generated by splines package).")
  if(model.type=='poly' & missing(order)) stop("order must be specified if poly method is used.")
  
  if(eq.plot==T) {
    # Linear y=mx+c equation
    if(model.type=='lin') {
      fit = lm(df[[y.var]] ~ df[[x.var]]);
      eq <- substitute(italic(y) == c + m %.% italic(x)*","~~italic(r)^2~"="~r2, 
                       list(c = signif(coef(fit)[1], signif.figs), 
                            m = signif(coef(fit)[2], signif.figs), 
                            r2 = signif(summary(fit)$r.squared, signif.figs)))
      as.character(as.expression(eq));
    }
    # polynomial expression generated with the ns(x,order) function [splines package]
    if(model.type=='poly') {
      fit = lm(df[[y.var]] ~ ns(df[[x.var]],order));
      
      base = gsub('!c!',signif(coef(fit)[1],signif.figs),"italic(y) == !c! + ")
      element.1 = "!m! %.% italic(x)~"
      element.2 = " + !m! %.% italic(x)^!o!~"
      element.r2 = gsub('!r2!',signif(summary(fit)$r.squared,signif.figs),"~~italic(r)^2~\"=\"~!r2!")
      eq=""
      
      for(o in 1:(order)) {
        if(o==1) {
          if(coef(fit)[(o+1)]<0) tmp=gsub("[+]","",base) else tmp=base
          eq=paste(tmp,gsub('!m!',signif(coef(fit)[(o+1)],signif.figs),element.1),sep="")
        }
        if(o>1) {
          if(coef(fit)[(o+1)]<0) tmp=gsub("[+]","",element.2) else tmp=element.2
          eq=paste(eq,gsub('!o!',o,gsub('!m!',signif(coef(fit)[(o+1)],signif.figs),tmp)),sep="")
        }
        if(o==(order)) eq=paste(eq,"\",\"",element.r2,sep="")
      }
    }
  }
  if(eq.plot==F) {
    # Linear y=mx+c equations
    if(model.type=='lin') {
      fit = lm(df[[y.var]] ~ df[[x.var]]);
      eq <- substitute(italic(r)^2~"="~r2, 
                       list(r2 = signif(summary(fit)$r.squared, signif.figs)))
      as.character(as.expression(eq));
    }
    # polynomial expression generated with the ns() function [splines package]
    if(model.type=='poly') {
      fit = lm(df[[y.var]] ~ ns(df[[x.var]],order));
      
      eq = gsub('!r2!',signif(summary(fit)$r.squared,signif.figs),"italic(r)^2~\"=\"~!r2!")
    }
  }
  return(eq)
}

# To implement:
eqn = annotate("text",x=5,y=16,label=lm_eqn(ts,'x','y',3,T,'poly',3),parse=T)
p + eqn


--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: https://github.com/hadley/devtools/wiki/Reproducibility
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

expression.png
auto-expression.png

Bryan Hanson

unread,
Jan 18, 2013, 1:40:47 PM1/18/13
to Oscar Branson, ggp...@googlegroups.com
If you want to stick with lm or rlm then you can use this function:


Bryan
Reply all
Reply to author
Forward
0 new messages