library(qpcR)
library(lmtest)
library(lawstat)
library(gplots)
library(lattice)
CheckPeriodApp <- function(x, cycle = 20, loess.span = 0.1, peak.window = 5, plot = TRUE, ...) {
SEL <- x[, 2]
POS <- x[, 1]
## make quadratic model
nSample <- 1L:length(SEL)
LM <- lm(SEL ~ I(nSample) + I(nSample^2))
FITTED <- fitted(LM)
RESID <- residuals(LM)
COEFS <- as.numeric(summary(LM)$coefficients[2:3, 4])
LOESS <- loess(RESID ~ nSample, span = loess.span)
RUNS <- runs.test(RESID)
BOX <- Box.test(RESID, lag = 1, type = "Ljung-Box")
ACF <- acf(RESID, lag = length(RESID), plot = FALSE)
FP <- findpeaks(as.numeric(ACF$acf), bw = peak.window)
PERIOD <- diff(FP$xmax)
## return statistics as list and for "textplot"ing
lmStat <- c(signif(COEFS[1], 3), signif(COEFS[2], 3))
names(lmStat) <- c("lin", "quad")
randStat <- c(signif(RUNS$p.value, 3), signif(BOX$p.value, 3))
names(randStat) <- c("Runs test", "Ljung-Box test")
periodStat <- c(round(mean(PERIOD, na.rm = TRUE), 1), round(sd(PERIOD, na.rm = TRUE), 1))
names(periodStat) <- c("Period (Mean)", "Period (SD)")
list(SEL = SEL,
POS = POS,
FITTED = FITTED,
COEFS = COEFS,
RESID = RESID,
LOESS = LOESS,
RUNS = RUNS,
BOX = BOX,
PERIOD = PERIOD,
FP = FP,
ACF = ACF,
#output for users
user = list(lmStat = lmStat, randStat = randStat, periodStat = periodStat))
}
## plot data with fitted line from quadratic model
## and show coefficients
plotFit <- function(PeriodApp) {
COL <- rainbow(length(PeriodApp[["FITTED"]]))
plot(PeriodApp[["SEL"]], col = COL, pch = 16, xlab = "Sample id", ylab = "Cq")
lines(PeriodApp[["FITTED"]], col = 1, lwd = 2)
#title(main = paste("lin:", signif(PeriodApp[["COEFS"]][1], 3), " quad:",
# signif(PeriodApp[["COEFS"]][2], 3)), line = -1.5, cex.main = 1.5)
}
## plot residuals with Loess fit
## and do Runs test & Ljung-Box Q-test
plotRes <- function(PeriodApp) {
COL <- rainbow(length(PeriodApp[["FITTED"]]))
plot(PeriodApp[["RESID"]], col = COL, pch = 16, xlab = "Sample id", ylab = "Residual")
lines(fitted(PeriodApp[["LOESS"]]), lwd = 2)
# title(main = paste("Runs test:", signif(PeriodApp[["RUNS"]][["p.value"]], 3),
# " Ljung-Box test:", signif(PeriodApp[["BOX"]][["p.value"]], 3)),
# line = -1.5, cex.main = 1.5)
}
## plot autocorrelation
plotAc <- function(PeriodApp) {
BP <- barplot(PeriodApp[["ACF"]][["acf"]][, , 1], col = "darkgrey", space = 2, ylab = "Autocorrelation of residuals",
ylim = c(min(PeriodApp[["ACF"]][["acf"]][, , 1]) * 1.05, max(PeriodApp[["ACF"]][["acf"]][, , 1])*1.2))
xPos <- BP[, 1]
points(xPos[PeriodApp[["FP"]][["xmax"]]], 0.1 * mean(PeriodApp[["FP"]][["ymax"]]) + PeriodApp[["FP"]][["ymax"]], pch = 25,
col = "darkred")
# title(main = paste("Estimated Periodicity:", round(mean(PeriodApp[["PERIOD"]], na.rm = TRUE), 1), "\u00B1",
# round(sd(PeriodApp[["PERIOD"]], na.rm = TRUE), 1)), line = -1.5, cex.main = 1.5)
}
plotHm <- function(PeriodApp) {
COL <- colorRampPalette(c("darkblue", "white", "darkred"))
pos <- as.character(PeriodApp[["POS"]])
pos_letter <- substr(pos, 0, 1)
pos_number <- vapply(pos, function(i)
substr(i, 2, nchar(i)), "a", USE.NAMES = FALSE)
MAT <- matrix(PeriodApp[["RESID"]], nrow = length(unique(pos_letter)), ncol = length(unique(pos_number)),
dimnames = list(unique(pos_letter), unique(pos_number)), byrow = TRUE)
levelplot(t(MAT[nrow(MAT):1, ]), col.regions = COL(100))
}
## Taken from https://rtricks.wordpress.com/2009/05/03/an-algorithm-to-find-local-extrema-in-a-vector/
findpeaks <- function(vec, bw = 1 , x.coo = c(1:length(vec)))
{
pos.x.max <- NULL
pos.y.max <- NULL
pos.x.min <- NULL
pos.y.min <- NULL
for(i in 1:(length(vec) - 1)) {
if((i + 1 + bw) > length(vec)) sup.stop <- length(vec) else sup.stop <- i + 1 + bw
if((i - bw) < 1) inf.stop <- 1 else inf.stop <- i - bw
subset.sup <- vec[(i + 1):sup.stop]
subset.inf <- vec[inf.stop:(i - 1)]
is.max <- sum(subset.inf > vec[i]) == 0
is.nomin <- sum(subset.sup > vec[i]) == 0
no.max <- sum(subset.inf > vec[i]) == length(subset.inf)
no.nomin <- sum(subset.sup > vec[i]) == length(subset.sup)
if(is.max & is.nomin){
pos.x.max <- c(pos.x.max,x.coo[i])
pos.y.max <- c(pos.y.max,vec[i])
}
if(no.max & no.nomin){
pos.x.min <- c(pos.x.min,x.coo[i])
pos.y.min <- c(pos.y.min,vec[i])
}
}
return(list(xmax = pos.x.max, ymax = pos.y.max, xmin = pos.x.min, ymin = pos.y.min))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.