/scripts/cluster_papers.R

http://github.com/hpiwowar/alt-metrics_stats · R · 71 lines · 33 code · 17 blank · 21 comment · 0 complexity · d263d00af01b4c1b67687315a908d0e1 MD5 · raw file 

    

  1. (flip = t(dat.indep.stats.norm)
    2. 

  2. #flip.df = as.data.frame(flip[,sample(1:dim(flip)[2], 100, T)])
    3. 

  3. flip.df = as.data.frame(flip[,20000:20050])
    4. 

  4. 

  5. myhetcorr = hetcor.modified(flip.df, 
    6. 

  6. 							use="pairwise.complete.obs", 
    7. 

  7. 							std.err=FALSE, 
    8. 

  8. 							pd=FALSE, 
    9. 

  9. 							type="spearman")
    10. 

  10. mycor.unadjusted = myhetcorr$correlations
    11. 

  11. #write.table(mycor.unadjusted,"/Users/hpiwowar/statslink/mycor.unadjusted.txt",append=F,quote=F,sep="\t",row.names=T)
    12. 

  12. 

  13. # Are some correlations NA?
    14. 

  14. which(is.na(mycor.unadjusted))
    15. 

  15. #mycor.unadjusted[which(is.na(mycor.unadjusted))] = 1
    16. 

  16. 

  17. # Now fix the correlation matrix if it is not positive-definite
    18. 

  18. mycor = adjust.to.positive.definite(mycor.unadjusted)
    19. 

  19. 

  20. #display
    21. 

  21. library(gplots)
    22. 

  22. colorRange = round(range(mycor) * 15) + 16
    23. 

  23. colorChoices = bluered(32)[colorRange[1]:colorRange[2]]
    24. 

  24. heatmap.2(mycor, col=colorChoices, cexRow=0.9, cexCol = .9, symm = TRUE, dend = "row", trace = "none", main = "Thesis Data", margins=c(10,10), key=FALSE, keysize=0.1)
    25. 

  25. 

  26. showpanel <- function(col) {
    27. 

  27.   image(z=matrix(1:100, ncol=1), col=col, xaxt="n", yaxt="n" )
    28. 

  28. }
    29. 

  29. quartz()
    30. 

  30. showpanel(colorChoices)
    31. 

  31. showpanel(bluered(32))
    32. 

  32. 

  33. #pdf("heatmap.pdf", height=10, width=10)
    34. 

  34. #heatmap.2(mycor, col=colorChoices, cexRow=0.9, cexCol = .9, symm = TRUE, dend = "row", trace = "none", main = "Thesis Data", margins=c(10,10), key=FALSE, keysize=0.1)
    35. 

  35. #dev.off
    36. 

  36. 

  37. ############## FIRST ORDER ANALYSIS
    38. 

  38. 

  39. ##############  Determine number of First-Order factors
    40. 

  40. 

  41. # Determine Number of Factors to Extract
    42. 

  42. library(nFactors)
    43. 

  43. eigenvectors.1st <- eigen(mycor) # get eigenvalues
    44. 

  44. # this line takes a long time
    45. 

  45. aparallel.1st <- parallel(subject=nrow(dat.indep.stats.norm), var=ncol(dat.indep.stats.norm), rep=100, cent=.05)
    46. 

  46. scree.results.1st <- nScree(eigenvectors.1st$values, aparallel.1st$eigen$qevpea)
    47. 

  47. summary(scree.results.1st)
    48. 

  48. plotnScree(scree.results.1st) 
    49. 

  49. 

  50. # Pull out the "Optimal Coordinate"
    51. 

  51. # defined in nScree help page as 
    52. 

  52. # The optimal coordinates (OC) corresponds to an extrapolation of the preceding eigenvalue by a regression line between the eignvalue coordinates and the last eigenvalue coordinate
    53. 

  53. #number.factors.1st = scree.results.1st$Components$noc
    54. 

  54. number.factors.1st = 5
    55. 

  55. 

  56. 

  57. ##############  Do First-Order Factor Analysis
    58. 

  58. 

  59. # Maximum liklihood doesn't converge because too 
    60. 

  60. fit.ml = factanal(dat.indep.stats.norm, number.factors.1st, rotation="promax", covmat=mycor)
    61. 

  61. print(fit.ml, sort=TRUE)
    62. 

  62. 

  63. 

  64. # Use princip axis when maximum liklihood fails to converge:
    65. 

  65. library(psych)
    66. 

  66. fit.fa.1st = fa(mycor, number.factors.1st, fm="minres", rotate="promax", 
    67. 

  67.                 scores=FALSE, residuals=TRUE, n.obs=max(dim(dat.indep.stats.norm)))
    68. 

  68. 

  69. #to show the loadings sorted by absolute value
    70. 

  70. print(fit.fa.1st, sort=TRUE)
    71.