1

I have plotted the following data with ggplot2 and then add 3rd order polynomial lines for each group consulting the post: Adding a 3rd order polynomial and its equation to a ggplot in r Data:

>lai.se1    
    DOS DAS N   LAI sd  se  ci
    D1  31  24  1.5879167   0.42763230  0.08729008  0.18057328
    D1  84  24  4.3241667   0.32478644  0.06629675  0.13714529
    D1  113 24  3.7037500   0.34151596  0.06971165  0.14420954
    D1  132 24  2.9704167   0.33386380  0.06814966  0.14097832
    D1  160 24  0.1879167   0.09868611  0.02014422  0.04167149
    D2  35  24  1.7679167   0.18551876  0.03786886  0.07833770
    D2  82  24  3.7670833   0.38212767  0.07800148  0.16135836
    D2  108 24  3.4104167   0.31431747  0.06415978  0.13272463
    D2  126 24  2.7879167   0.35024189  0.07149283  0.14789418
    D2  146 24  0.1950000   0.08836682  0.01803780  0.03731404
    D3  37  24  1.3179167   0.16378616  0.03343271  0.06916083
    D3  83  24  3.5233333   0.29256982  0.05972057  0.12354140
    D3  94  24  3.1604167   0.28257326  0.05768002  0.11932022
    D3  113 24  2.4587500   0.44131535  0.09008312  0.18635113
    D3  134 24  0.2758333   0.09536733  0.01946677  0.04027009

p<-ggplot(lai.se1, aes(x=DAS, y=LAI, colour=DOS)) + 
  geom_errorbar(aes(ymin=LAI-sd, ymax=LAI+sd), colour ="black", size =.3, width=1,
                position=position_dodge(.9)) +
  geom_point(size=1, shape=21, fill=FALSE)+ theme_bw()
p + stat_smooth(method="lm", se=TRUE, fill=NA,             ## to add polynomial lines
                formula=y ~ poly(x, 3, raw=TRUE))

Plot: enter image description here

## Add equation in the plot
lm_eqn = function(lai.se1){
  m=lm(y ~ poly(x, 3), lai.se1)#3rd degree polynomial
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
                   list(a = format(coef(m)[1], digits = 2),
                        b = format(coef(m)[2], digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq))
}


p + annotate("text", x=0.5, y=15000, label=lm_eqn(lai.se1), hjust=0, size=8, 
             family="Times", face="italic", parse=TRUE)


Error:
Error in model.frame.default(formula = y ~ poly(x, 3), data = lai.se1,  : 
  variable lengths differ (found for 'poly(x, 3)')

But as I tried to place the equations using similar function mentioned the above linked post, it gave an error of variable length. I am not really good at writing new functions in r and need your help to solve it. Please help.

Community
  • 1
  • 1
Barun
  • 185
  • 1
  • 2
  • 10
  • 1
    What command gave the error? What is the error `exactly`? It'd be a good place for people trying to help to start. So far your question is not clear (or rather there is no question except that you would like to get help). – Arun Jan 29 '13 at 13:27
  • 1
    Dear Arun, I want to put the respective equations and r^2 values for those lines in the same plot. that's my problem. I am adding the code and the generated error. – Barun Jan 29 '13 at 13:57
  • (+1) for the nice and complete edit! – Arun Jan 29 '13 at 15:04

1 Answers1

1

Your code is right, except for two things.

  • First, the function lm_eqn function does not undertand what is y and x. So, you'll have to pass the corresponding column names of lai.se1 as follows:

    lm_eqn = function(lai.se1){
        m=lm(LAI ~ poly(DAS, 3), lai.se1) #3rd degree polynomial
        eq <- substitute(italic(LAI) == a + b %.% italic(DAS)*","~~italic(r)^2~"="~r2, 
                 list(a = format(coef(m)[1], digits = 2), 
                 b = format(coef(m)[2], digits = 2), 
                 r2 = format(summary(m)$r.squared, digits = 3)))
                 as.character(as.expression(eq))
    }
    
  • 2) Just the order has to be reversed a bit (due to stat_smooth as annotate just adds geoms to a plot. Otherwise, the stat_smooth output will be replaced by annotate ):

    p <- ggplot(lai.se1, aes(x = DAS, y = LAI, colour = DOS))
    p <- p + geom_errorbar(aes(ymin = LAI-sd, ymax = LAI+sd), 
             colour = "black", size = .3, width = 1, 
             position = position_dodge(.9)) 
    p <- p + stat_smooth(method = "lm", se = TRUE, fill = NA, 
             formula = y ~ poly(x, 3, raw = TRUE))   
    p <- p + geom_point(size = 1, shape = 21, fill = FALSE) + theme_bw()
    p + annotate("text", x=0.5, y=7.5, label=lm_eqn(lai.se1), hjust=0, size=8, 
                 family="Times", face="italic", parse=TRUE)
    p
    

Edit: Proposed solution after OP's comment. How about this?

require(plyr)
lm_eqn = daply(lai.se1, .(DOS), function(w) {
    m = lm(LAI ~ poly(DAS, 3), w)
    eq <- substitute(italic(LAI) == a + b %.% italic(DAS)*","~~italic(r)^2~"="~r2, 
             list(a = format(coef(m)[1], digits = 2), 
             b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
             as.character(as.expression(eq))
    })


p <- ggplot(lai.se1, aes(x = DAS, y = LAI, colour = DOS))
p <- p + geom_errorbar(aes(ymin = LAI-sd, ymax = LAI+sd), 
         colour = "black", size = .3, width = 1, 
         position = position_dodge(.9)) 
p <- p + stat_smooth(method = "lm", se = TRUE, fill = NA, 
         formula = y ~ poly(x, 3, raw = TRUE))   
p <- p + geom_point(size = 1, shape = 21, fill = FALSE) + theme_bw()
p <- p + annotate("text", x=0.5, y=5.5, label=lm_eqn[1], hjust=0, size=8, 
             family="Times", face="italic", parse=TRUE)
p <- p + annotate("text", x=0.5, y=6.5, label=lm_eqn[2], hjust=0, size=8, 
             family="Times", face="italic", parse=TRUE)
p <- p + annotate("text", x=0.5, y=7.5, label=lm_eqn[3], hjust=0, size=8, 
             family="Times", face="italic", parse=TRUE)
p

Basically, you generate all equations and then annotate that many times (you can probably loop them if you have more...)

ggplot2_edit2

Arun
  • 116,683
  • 26
  • 284
  • 387
  • Thanks for the help. But what about 3 separate equations for D1, D2 and D3? – Barun Jan 29 '13 at 15:27
  • check out the edit. don't forget to up-vote and **mark** as answered (by checking the tick mark) if you have got it answered! – Arun Jan 29 '13 at 15:36
  • 1
    thanks a ton Arun. just want little clarification: is lm_eqn[1] is for D1 and so on? – Barun Jan 29 '13 at 15:52
  • For this case, Yes. As a general answer, you can check by doing `print(x$DOS)` inside the equation function. The function splits your data.frame by DOS and creates an equation for each unique value of DOS. So, depending on the levels of this factor (DOS), you'll get the index. If you require it in the order D1, D2, D3, then you can reorder the factor as: `lai.se1$DOS = factor(lai.se1$DOS, levels=c("D1", "D2", "D3"), ordered = TRUE)` – Arun Jan 29 '13 at 15:55