4

The area under the curve can be calculated using the trapz function of the pracma package in R. iAUC is much more accurate in many instances, particularly in biology. However there is no, to me known, R function to calculate this. The trap function follows:

Example:

a <- c(1, 4, 5, 6)
b <- c(2, 4, 6, 8)
plot(a, b)
trapz(a, b)

Answer: 21

Does anyone know how to use R to calculate the iAUC? A guide can be found here.

Thanks in advance.

Update In response to Cheng-Yow's answer below: Another example

# Recall that trapz-function is available in the pracma package
# Here is the function kindly provided by Cheng-Yow Timothy Lin
library(pracma) # should be loaded already
traps <- function(a, b) {
     for (i in 1:length(b)-1){
          iAUC[i] <-((b[i]+b[i+1])/2)*(a[i+1]-a[i])
     }
     return(sum(iAUC))
}

# Data example
a <- c(-15,  -5,   1,   2,   3,   5,   8,  10,  15,  20,  30,  45,  60,  90, 120)
b <- c(50.20604,  49.59338,  47.39944,  56.38831,  69.43493,  73.92512,  61.92072,  67.92632, 115.45669,
195.03242, 322.15894, 291.30094, 289.20284, 238.70562, 156.23798)

traps(a, b)

trapz(a, b)

They yield the same results- is this the incremental area under the curve? It is not how they have explained the trapz-function...

Thanks for any enlightenment!

Adam Robinsson
  • 1,651
  • 3
  • 17
  • 29

1 Answers1

4
    a <- c(1,2,3,4,5)
b <- c(1,2,-1,-2,2)
plot(a,b)
lines(a, b)
abline(b[1],0)
iAUC <- vector()

for (i in 1:(length(a)-1)){

  if((b[i+1]-b[1] >= 0) && (b[i]-b[1] >= 0))
  {iAUC[i] <-((b[i]-b[1]+b[i+1]-b[1])/2)*(a[i+1]-a[i])} 

  else if((b[i+1]-b[1] < 0) && (b[i]-b[1] >= 0))
  {iAUC[i] <-(b[i]-b[1])*((b[i]-b[1])/(b[i]-b[i+1])*(a[i+1]-a[i])/2)} 

  else if((b[i+1]-b[1] >= 0) && (b[i]-b[1] < 0))
  {iAUC[i] <-(b[i+1]-b[1])*((b[i+1]-b[1])/(b[i+1]-b[i])*(a[i+1]-a[i])/2)} 

  else if((b[i]-b[1] < 0) && (b[i+1]-b[1] < 0))
  {iAUC[i] <- 0}
}
sum(iAUC)
iAUC
Cheng-Yow Timothy Lin
  • 93
  • 8
  • Many many thanks Cheng-Yow for a nice function. Please see my update in the question above, as a response to your function. – Adam Robinsson Aug 12 '15 at 19:53
  • Hi Adam, one of the key differences is that my code uses a loop, which is slow when your data gets larger. The function in the pracma package is vectorized so it will perform much faster. Try running this after your snippet above! `a <- rep(a,1000) b <- rep(b,1000) plot(a,b) ptm <- proc.time() traps(a, b) proc.time()-ptm ptm <- proc.time() trapz(a, b) proc.time()-ptm`Also, don't forget to initalze iAUC vector with `iAUC<-vector()` in your function 'traps'" – Cheng-Yow Timothy Lin Aug 12 '15 at 20:48
  • Many thanks @Cheng-Yow, but does your traps function provide the incremental AUC (I'm sorry for repeating a possibly stupid question)? – Adam Robinsson Aug 12 '15 at 21:08
  • Because it provides the same answer as pracma, which does not state that it calculates incremental AUC... =) – Adam Robinsson Aug 12 '15 at 21:09
  • I just reread your word document that you provided. Sorry, according to Woelver, the pracma and the function i wrote does not calculate iAUC. it just calculates AUC. iAUC is baseline corrected, which means that instead of taking the areas about the origin, y=0, it takes y=(baseline measure, or first measure) I will edit my answer to reflect these changes. Sorry for the confusion. – Cheng-Yow Timothy Lin Aug 12 '15 at 21:55
  • many many thanks for the answers You have provided, I really appreciate Your time. However, copying and pasting the code directly into R, trying to execute it, failed. I tried editing it (although that's beyond me) but couldn't figure it out. Did anyone else have troubles executing it? Moreover, the reason for the: "for (i in 1:4){" is it because there are four change-points? It would be lovely if the code was general, so that it could be used for data sets with varying number of x's and y's. Thanks in advance. – Adam Robinsson Aug 16 '15 at 08:10
  • Applying the code to the vector I provided (a and b from the original post) returns an implausible value. I apologize if I've misunderstood the code. – Adam Robinsson Aug 16 '15 at 08:11
  • I did a quick check, appears to be correct now! A million thanks! – Adam Robinsson Aug 18 '15 at 08:14