#http://rtutorialseries.blogspot.com/2011/10/r-tutorial-series-exploratory-factor.html#http://www.statmethods.net/advstats/factor.html#http://www.statmethods.net/advstats/cluster.html

glut.factor.analyses.R

setwd("C:/Users/Sasha/Google Drive/LNGstuff/glut analyses")
data = read.csv("all.genesLabelled2.csv")

library(psych)
data$structure = NULL

#sample data format (after NULL):
#     GRIA1     GRIA2    GRIA3    GRIA4    GRIK1    GRIK2    GRIK3    GRIK4
#1 4.200910 13.198395 4.531563 4.283639 2.323050 1.680551 0.591786 1.653584
#2 1.193841  4.902818 1.869044 1.779108 1.620001 0.839528 0.138910 1.132033
#3 1.018170  3.949605 1.290182 0.966998 0.775067 0.623651 0.168454 1.008362
#4 3.091375  9.771480 3.656806 3.343708 1.854348 1.269660 0.636690 1.648468
#5 0.580231  3.027282 1.485677 0.929639 0.562741 0.335170 0.166034 1.145590
#6 0.409548  3.244448 1.216516 0.774867 0.827955 0.724840 0.115512 0.614748

corMat = cor(data)
solution <- fa(r = corMat, nfactors = 5, rotate = "varimax", fm = "pa") 

#customizing correlation matrix function to display diagonal values
my.cor.plot = 
function (r, numbers = FALSE, colors = TRUE, n = 51, main = NULL, 
                       zlim = c(-1, 1), show.legend = TRUE, labels = NULL, n.legend = 10, 
                       keep.par = TRUE, ...) 
{
    if (keep.par) 
        op <- par(no.readonly = TRUE)
    if (is.null(main)) {
        main <- "Correlation plot"
    }
    if (!is.matrix(r) & (!is.data.frame(r))) {
        if ((length(class(r)) > 1) & (class(r)[1] == "psych")) {
            if (class(r)[2] == "omega") {
                r <- r$schmid$sl
                nff <- ncol(r)
                r <- r[, 1:(nff - 2)]
            }
            else {
                r <- r$loadings
            }
        }
    }
    r <- as.matrix(r)
    if (min(dim(r)) < 2) {
        stop("You need at least two dimensions to make a meaningful plot")
    }
    r <- t(r)
    if (is.null(n)) {
        n <- dim(r)[2]
    }
    nf <- dim(r)[2]
    nvar <- dim(r)[1]
    if (is.null(labels)) {
        if (is.null(rownames(r))) 
            rownames(r) <- paste("V", 1:nvar)
        if (is.null(colnames(r))) 
            colnames(r) <- paste("V", 1:nf)
    }
    else {
        rownames(r) <- colnames(r) <- labels
    }
    max.len <- max(nchar(rownames(r)))/6
    if (is.null(zlim)) {
        zlim <- range(r)
    }
    if (colors) {
        gr <- colorRampPalette(c("red", "white", "blue"))
        colramp <- gr(n)
    }
    else {
        colramp <- grey((n:0)/n)
    }
    ord1 <- seq(nvar, 1, -1)
    if (nvar != nf) {
        r <- t(r)
    }
    r <- r[, ord1]
    MAR <- 5
    par(mar = c(MAR + max.len, MAR + max.len, 4, 0.5))
    if (show.legend) {
        layout(matrix(c(1, 2), nrow = 1), widths = c(0.9, 0.1), 
               heights = c(1, 1))
    }
    image(r, col = colramp, axes = FALSE, main = main, zlim = zlim)
    box()
    at1 <- (0:(nf - 1))/(nf - 1)
    at2 <- (0:(nvar - 1))/(nvar - 1)
    if (max.len > 0.5) {
        axis(2, at = at2, labels = colnames(r), las = 1, ...)
        axis(1, at = at1, labels = rownames(r), las = 2, ...)
    }
    else {
        axis(2, at = at2, labels = colnames(r), ...)
        axis(1, at = at1, labels = rownames(r), las = 1, ...)
    }
    at1 <- (0:(nf - 1))/(nf - 1)
    if (numbers) {
        rx <- rep(at1, ncol(r))
        ry <- rep(at2, each = ncol(r))
        text(rx, ry, round(r * 100), srt=-45)
    }
    if (show.legend) {
        leg <- matrix(seq(from = zlim[1], to = zlim[2], by = (zlim[2] - 
                                                                  zlim[1])/n), nrow = 1)
        par(mar = c(MAR, 0, 4, 3))
        image(leg, col = colramp, axes = FALSE, zlim = zlim)
        at2 <- seq(0, 1, 1/n.legend)
        labels = seq(zlim[1], zlim[2], (zlim[2] - zlim[1])/(length(at2) - 
                                                                1))
        axis(4, at = at2, labels = labels, las = 2, ...)
    }
    if (keep.par) 
        par(op)
}
my.cor.plot(corMat,numbers=TRUE,colors=TRUE,n=51,main=NULL,labels=NULL)