AIC.vs.Rsq <- function()
{
      library(qpcR)       
      
      ### define different definitions of R-square
      RSS <- function(x) sum(residuals(x)^2, na.rm = TRUE)
      REGSS <- function(x) sum((fitted(x) - mean(fitted(x)))^2, na.rm = TRUE)
      TSS <- function(x) sum((x$data[, 2] - mean(x$data[, 2], na.rm = TRUE))^2, na.rm = TRUE)
      RSQ1 <- function(x) REGSS(x)/TSS(x)
      RSQ2 <- Rsq
      RSQ3 <- function(x) 1-(RSS(x)/TSS(x))
      
      m <- multdrc(F1.1 ~ Cycles, data = reps[10:25, ], fct = l5())

      ## define loop
      nsim <- 2000
      ## define noise
      noise <- 0.01

      ### define model list
      modList <- list(l5(), l4(), l3(), b5(), b4(), b3(), w4(), w3(), baro5())

      ## define result matrices and variables
      aic <- rsq1 <- rsq2 <- rsq3 <- loglik <- rss <- matrix(nrow = length(modList), ncol = nsim)
      m.aic <- m.rsq1 <- m.rsq2 <- m.rsq3 <- m.loglik <- m.rss <- matrix(nrow = length(modList), ncol = length(noise))
      s.aic <- s.rsq1 <- s.rsq2 <- s.rsq3 <- s.loglik <- s.rss <- matrix(nrow = length(modList), ncol = length(noise))
      counter <- 1

      ## take fitted values as optimal curve
      optDat <- fitted(m)
      lenDat <- length(optDat)

      ## loop through increasing noise
      for (j in noise) {
            ## iterate nsim simulations
            for (i in 1:nsim) {
                  ## generate data with increasing constant noise for each datapoint
                  #simDat <- sapply(optDat, function(x) rnorm(1, x, j))
                  ## if desired generate data with noise that is modelled from the value
                  simDat <- sapply(optDat, function(x) rnorm(1, x, x * j))
                  ## build models for each iteration and coerce to list
                  resList <- lapply(modList, function(x) try(multdrc(simDat ~ I(10:25), fct = x), silent = TRUE))
                  pcrplot(resList[[1]], main = paste(j, i))
                  cat(j, i, "\n")
                  ## generate measures  
                  aic[, i] <- sapply(resList, function(x) as.numeric(try(AIC(x), silent = TRUE)))
                  rsq1[, i] <- sapply(resList, function(x) as.numeric(try(RSQ1(x), silent = TRUE)))
                  rsq2[, i] <- sapply(resList, function(x) as.numeric(try(RSQ2(x), silent = TRUE)))
                  rsq3[, i] <- sapply(resList, function(x) as.numeric(try(RSQ3(x), silent = TRUE)))
                  loglik[, i] <- sapply(resList, function(x) as.numeric(try(logLik(x), silent = TRUE)))
                  rss[, i] <- sapply(resList, function(x) as.numeric(try(resVar(x), silent = TRUE)))
            }
            ## average all measures
            m.aic[, counter] <- apply(aic, 1, function(x) mean(x, na.rm = TRUE))
            m.rsq1[, counter] <- apply(rsq1, 1, function(x) mean(x, na.rm = TRUE))
            m.rsq2[, counter] <- apply(rsq2, 1, function(x) mean(x, na.rm = TRUE))
            m.rsq3[, counter] <- apply(rsq3, 1, function(x) mean(x, na.rm = TRUE))
            m.loglik[, counter] <- apply(loglik, 1, function(x) mean(x, na.rm = TRUE))
            m.rss[, counter] <- apply(rss, 1, function(x) mean(x, na.rm = TRUE))
            
            s.aic[, counter] <- apply(aic, 1, function(x) sd(x, na.rm = TRUE))
            s.rsq1[, counter] <- apply(rsq1, 1, function(x) sd(x, na.rm = TRUE))
            s.rsq2[, counter] <- apply(rsq2, 1, function(x) sd(x, na.rm = TRUE))
            s.rsq3[, counter] <- apply(rsq3, 1, function(x) sd(x, na.rm = TRUE))
            s.loglik[, counter] <- apply(loglik, 1, function(x) sd(x, na.rm = TRUE))
            s.rss[, counter] <- apply(rss, 1, function(x) sd(x, na.rm = TRUE))

            counter <- counter + 1
      }
      return(list(m.aic = m.aic, s.aic = s.aic, m.rsq1 = m.rsq1, s.rsq1 = s.rsq1, 
                  m.rsq2 = m.rsq2, s.rsq2 = s.rsq2, m.rsq3 = m.rsq3, s.rsq3 = s.rsq3,
                  m.loglik = m.loglik, s.loglik = s.loglik, m.rss = m.rss, s.rss = s.rss, 
                  resList = resList))
}

### run function
res <- AIC.vs.Rsq()

### plotting setup  (Figure 1 & Supplemental Figure 1)
bp <- barplot(res$m.aic[, 1], ylim = c(34, 40), col = rainbow(9))
par(mar = c(5,5,5,5)); plot(res$m.rsq3[, 1], type = "p",  pch = 16, col = rainbow(9), lwd = 4, cex = 4, xlab = "", ylab = "", ylim = c(0.975, 0.99), axes = FALSE); lines(res$m.rsq3[, 1], col = 1, lwd = 2); axis(side = 2, line = 1)
pcrplot(res$resList[[1]], lwd = 1.5, xlim = c(12, 17), xlab = "Cycles", ylab = "RF"); sapply(2:9, function(x) pcrplot(res$resList[[x]], add = T, lwd = 1.5, col = rainbow(9)[x]))
residMat <- sapply(res$resList, function(x) residuals(x)); matplot(10:25, residMat, type = "l", lwd = 2, col = rainbow(9), lty = 1, xlab = "Cycles", ylab = "Residual value")
plot(1:9, rep(1,9), pch = 16, col = rainbow(9), cex = 4, axes = F)
write.table(res$m.rsq3, file = "clipboard", sep = "\t", row.names = F, col.names = F)

### calculate p-values from chi-square of likelihood ratio
pval.LR <- function(x, y, DFx, DFy)
{
      stat <- abs(2*(x-y))
      deltaDF <- abs(DFx - DFy)
      pval <- 1 - pchisq(stat, deltaDF)
      return(pval)
}
pval.LR(res$m.loglik[1,], res$m.loglik[2,], 5, 4) # taking l5 and l4 as example

### calculate Akaike weights ###
aw <- akaike.weights(res$m.aic)$weights
write.table(aw, file = "clipboard", sep = "\t", row.names = F, col.names = F)