library(readxl)
library(randomForest)
library(ggplot2)
randomForest <- read_excel("~/Desktop/randomForest.xlsx")

set.seed(123)
forest <- randomForest(NUE~., data = randomForest, importance = TRUE, ntree = 300, nPerm = 600)
forest

importance_.scale <- data.frame(importance(forest, scale = TRUE), check.names = FALSE)
importance_.scale

importance_.scale <- importance_.scale[order(importance_.scale$'%IncMSE', decreasing = TRUE), ]

importance_.scale$Z_name <- rownames(importance_.scale)
importance_.scale$Z_name <- factor(importance_.scale$Z_name, levels = importance_.scale$Z_name)
p <- ggplot(importance_.scale, aes(Z_name, `%IncMSE`)) +
  geom_col(width = 0.5, fill = '#FFC068', color = NA) +
  labs(title = NULL, x = NULL, y = 'Increase in MSE (%)', fill = NULL) +
  theme(panel.grid = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = 'black')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(expand = c(0, 0), limit = c(-3, 8))

p

p <- p +
  annotate('text', label = 'SOM', x = 10, y = 6, size = 4) +
  annotate('text', label = sprintf('italic(R^2) == %.2f', 44.93), x = 10, y = 5.5, size = 3, parse = TRUE)

p

library(rfPermute)
set.seed(123)
rfP <- rfPermute(NUE~., data = randomForest, importance = TRUE, ntree = 300, nrep = 600, num.cores = 1)
rfP

importance_.scale <- data.frame(importance(rfP, scale = TRUE), check.names = FALSE)
importance_.scale

importance_.scale.pval <- (rfP$pval)[ , , 2]
importance_.scale.pval
plotNull(rfP, scale = TRUE, preds = 'YN') 
plotNull(rfP, scale = TRUE, preds = 'YP') 

plot(rfp.importance(rfP, scale = TRUE))

importance_.scale <- importance_.scale[order(importance_.scale$'%IncMSE', decreasing = TRUE), ]

library(ggplot2)
importance_.scale$Z_name <- rownames(importance_.scale)
importance_.scale$Z_name <- factor(importance_.scale$Z_name, levels = importance_.scale$Z_name)
p <- ggplot() +
  geom_col(data = importance_.scale, aes(x = Z_name, y = `%IncMSE`), width = 0.8, fill = '#FFC068', color = NA) +
  labs(title = NULL, x = NULL, y = 'Increase in MSE (%)', fill = NULL) +
  theme(panel.grid = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = 'black')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(expand = c(0, 0), limit = c(-3, 7)) + coord_flip()

p

for (Z in rownames(importance_.scale)) {
  importance_.scale[Z,'%IncMSE.pval'] <- importance_.scale.pval[Z,'%IncMSE']
  if (importance_.scale[Z,'%IncMSE.pval'] >= 0.05) importance_.scale[Z,'%IncMSE.sig'] <- ''
  else if (importance_.scale[Z,'%IncMSE.pval'] >= 0.01 & importance_.scale[Z,'%IncMSE.pval'] < 0.05) importance_.scale[Z,'%IncMSE.sig'] <- '*'
  else if (importance_.scale[Z,'%IncMSE.pval'] >= 0.001 & importance_.scale[Z,'%IncMSE.pval'] < 0.01) importance_.scale[Z,'%IncMSE.sig'] <- '**'
  else if (importance_.scale[Z,'%IncMSE.pval'] < 0.001) importance_.scale[Z,'%IncMSE.sig'] <- '***'
}

p <- p +
  geom_text(data = importance_.scale, aes( x = Z_name , y = `%IncMSE`, label = `%IncMSE.sig`), nudge_y = 1)

p

p <- p +
  annotate('text', label = 'Grain yield', x = 10, y = 6, size = 4) +
  annotate('text', label = sprintf('italic(R^2) == %.2f', 48.71), x = 9.3, y = 6, size = 3, parse = TRUE)

p

library(A3)
set.seed(123)
otu_forest.pval <- a3(NUE~., data = randomForest, model.fn = randomForest, p.acc = 0.001, model.args = list(importance = TRUE, ntree = 500))
otu_forest.pval

p <- p +
  annotate('text', label = sprintf('italic(P) < %.3f', 0.001), x = 8.6, y = 6, size = 3, parse = TRUE)

p