library(qpcR)

### set all variables to NULL
TSSs <- NULL
RSSs <- NULL
REGSSs <- NULL
RSQ1s <- NULL
RSQ2s <- NULL
RSQ3s <- NULL
AICs <- NULL

mTSSs <- NULL
mRSSs <- NULL
mREGSSs <- NULL
mRSQ1s <- NULL
mRSQ2s <- NULL
mRSQ3s <- NULL
mAICs <- NULL

sTSSs <- NULL
sRSSs <- NULL
sREGSSs <- NULL
sRSQ1s <- NULL
sRSQ2s <- NULL
sRSQ3s <- NULL
sAICs <- NULL

### increase noise
outerLoop <- seq(0, 10, by = 0.5)
### number of simulations per outerLoop
innerLoop <- 1:500
### define sample size
SEQ <- seq(1, 30, length.out = 30)

                          
for (i in outerLoop) {
      X <- reps[1:30, 1]
      Y <- reps[1:30, 2]
      Y2 <- Y * 100
      ### generate true model
      mod <- pcrfit(cbind(X, Y2), 1, 2, l5())
      ### generate different sample size
      Y3 <- pcrpred(mod, SEQ, which = "y")       
      mod2 <- pcrfit(cbind(SEQ, Y3), 1, 2, l5())
      ### take fitted values      
      FITTED <- fitted(mod2)         
               
      for (j in innerLoop) { 
            ### add gaussian noise            
            NOISE <- sapply(FITTED, function(x) rnorm(1, x, x * (i/100)))
            ### generate noisy model
            noisemod <- pcrfit(cbind(SEQ, NOISE), 1, 2, l5())
            pcrplot(noisemod, main = paste(i, j))
            cat(i, j, "\n")
            ### calculate sum-of-squares
            RSS <- sum(residuals(noisemod)^2, na.rm = TRUE)
            REGSS <- sum((fitted(noisemod) - mean(fitted(noisemod)))^2, na.rm = TRUE)
            TSS <- sum((Y3 - mean(Y3, na.rm = TRUE))^2, na.rm = TRUE)
            RSSs <- c(RSSs, RSS)
            REGSSs <- c(REGSSs, REGSS)
            TSSs <- c(TSSs, TSS)
            AICs <- c(AICs, AIC(noisemod))
            
            ### calculate R-square variants
            RSQ1s <- c(RSQ1s, REGSS/TSS)
            RSQ2s <- c(RSQ2s, REGSS/(REGSS + RSS))
            RSQ3s <- c(RSQ3s, 1 - (RSS/TSS))
      }
      
      mRSSs <- c(mRSSs, mean(RSSs, na.rm = TRUE))
      mREGSSs <- c(mREGSSs, mean(REGSSs, na.rm = TRUE))
      mTSSs <- c(mTSSs, mean(TSSs, na.rm = TRUE))
      mRSQ1s <- c(mRSQ1s, mean(RSQ1s, na.rm = TRUE))
      mRSQ2s <- c(mRSQ2s, mean(RSQ2s, na.rm = TRUE))
      mRSQ3s <- c(mRSQ3s, mean(RSQ3s, na.rm = TRUE))
      mAICs <- c(mAICs, mean(AICs, na.rm = TRUE))
      
      sRSSs <- c(sRSSs, sd(RSSs, na.rm = TRUE))
      sREGSSs <- c(sREGSSs, sd(REGSSs, na.rm = TRUE))
      sTSSs <- c(sTSSs, sd(TSSs, na.rm = TRUE))
      sRSQ1s <- c(sRSQ1s, sd(RSQ1s, na.rm = TRUE))
      sRSQ2s <- c(sRSQ2s, sd(RSQ2s, na.rm = TRUE))
      sRSQ3s <- c(sRSQ3s, sd(RSQ3s, na.rm = TRUE))
      sAICs <- c(sAICs, sd(AICs, na.rm = TRUE))
      
      ### reset all variables
      RSSs <- NULL
      REGSSs <- NULL
      TSSs <- NULL
      RSQ1s <- NULL
      RSQ2s <- NULL
      RSQ3s <- NULL
      AICs <- NULL
}


print(mTSSs)
print(mREGSSs + mRSSs)
print((mREGSSs + mRSSs)/ mTSSs)
print(mRSQ1s)
print(mRSQ2s)
print(mRSQ3s)



write.table(rbind(mTSSs, mREGSSs + mRSSs, abs((mREGSSs + mRSSs) / mTSSs), mAICs, mRSQ1s, mRSQ2s, mRSQ3s),
            file = "clipboard", sep = "\t", col.names = FALSE, row.names = FALSE)

plot(1:length(outerLoop), (mREGSSs + mRSSs)/mTSSs, pch = 16, cex = 2, cex.lab = 1.4, cex.axis = 1.4)
plot(1:length(outerLoop), mAICs, pch = 16, cex = 2, cex.lab = 1.4, cex.axis = 1.4)
plot(1:length(outerLoop), mRSQ1s, pch = 16, cex = 2, cex.lab = 1.4, cex.axis = 1.4)
plot(1:length(outerLoop), mRSQ2s, pch = 16, cex = 2, cex.lab = 1.4, cex.axis = 1.4)
plot(1:length(outerLoop), mRSQ3s, pch = 16, cex = 2, cex.lab = 1.4, cex.axis = 1.4)