a b/PCA_regression.R
1
#  Principal component regression analysis plots
2
lmmatrix<- function (y, cov){
3
  p <- matrix(NA, ncol = ncol(cov), nrow = ncol(y))
4
  for (j in 1:ncol(cov)) {
5
    x <- cov[, j]
6
    for (i in 1:ncol(y)) {
7
      fit <- summary(lm(y[, i] ~ x, na.action = na.omit))
8
      f <- fit$fstatistic
9
      p[i, j] <- pf(f["value"], f["numdf"], f["dendf"], lower.tail = FALSE)
10
    }
11
  }
12
  colnames(p) <- names(cov)
13
  return(p)
14
}
15
16
lmmatrix2<- function (y, cov){
17
  p <- matrix(NA, ncol = ncol(cov), nrow = ncol(y))
18
  #for (j in 1:ncol(cov)) {
19
  #x <- cov[, j]
20
  for (i in 1:ncol(y)) {
21
    fit <- anova(lm(y[, i] ~ ., cov, na.action = na.omit))
22
    p[i, ] <- fit$`Pr(>F)`[1:ncol(cov)]
23
  }
24
  #}
25
  colnames(p) <- names(cov)
26
  return(p)
27
}
28
29
plotp2 <- function (p, yaxis, xmax, title){
30
  plot(1, xlim = c(0, xmax), ylim = c(0, length(yaxis) + 1), type = "n", bty = "n", axes = FALSE, xlab = "Principal Component", ylab = "", main = title)
31
  axis(1, at = c(1:xmax), pos = 0.5, las = 1, lwd = 3)
32
  for (i in 1:length(yaxis)) {
33
    text(0.3, i, yaxis[i], xpd = TRUE, adj = 1)
34
  }
35
  for (i in 1:ncol(p)) {
36
    for (j in 1:nrow(p)) {
37
      pp <- p[j, i]
38
      colcode <- "white"
39
      if (pp <= 1e-09) {
40
        colcode = "darkred"
41
      }
42
      else if (pp <= 1e-04) {
43
        colcode = "red"
44
      }
45
      else if (pp <= 0.01) {
46
        colcode = "orange"
47
      }
48
      else if (pp <= 0.05) {
49
        colcode = "pink"
50
      }
51
      polygon(c(j - 0.5, j - 0.5, j + 0.5, j + 0.5), c(i - 0.5, i + 0.5, i + 0.5, i - 0.5), col = colcode, border = NA)
52
    }
53
  }
54
  legend("topright", c("<0.05", "<0.01", "<10E-5", "<10E-10"), col = c("pink", "orange", "red", "darkred"), pch = 15, pt.cex = 2, bty = "o", horiz = TRUE, xpd = TRUE)
55
}