{"cells":[{"metadata":{"_uuid":"d181b3b8ee5952298529cf26aa652dddc7ced5ee"},"cell_type":"markdown","source":"I am attempting to do factor analysis here based on variables/characteristics from the given data."},{"metadata":{"_uuid":"7f91f4011d6d1e25c77f6430aca427462d39acce","_execution_state":"idle","trusted":true},"cell_type":"code","source":"# Reading the input files\nlibrary(tidyverse) # metapackage with lots of helpful functions...\nlist.files(path = \"../input\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8c3cf6b471ac57090b87366180912938381d474c"},"cell_type":"code","source":"# Create a dataframe by reading the input file\ndf <- read.csv('../input/training_set_metadata.csv')\ndf <- df[,c(2:11)] # Select only the required columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8741bdb67c8e5f914a66ccb89072b075f4bb25a9"},"cell_type":"code","source":"# Replacing NA values with 0\ndf[is.na(df)] <- 0\nhead(df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"95eaad84339d61d00ea13c2ba4b8731cac3a73e0"},"cell_type":"markdown","source":"**Find out which variables can be considered for Factor Analysis:**\n\nWhich of these variables are relevant? Kaiser-Meyer-Olkin (KMO) test can be performed to answer this question. Function KMO() from psych package checks if a variable is suitable for factor analysis"},{"metadata":{"trusted":true,"_uuid":"64b58f97ec6ba57ab2d94eb154e3f2556b610354"},"cell_type":"code","source":"library(psych)\ndf_corr <- cor(scale(df)) # Create a correlation matrix\nKMO(df_corr) # Kaiser-Meyer-Olkin factor adequacy","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d4a00a28f8936c77a18175eef8db341ca599dd50"},"cell_type":"markdown","source":"MSA (measure of sampling adequacy) is a measure for exclusion of variables. If MSA < 0.5 the variable should be dropped. Variables with MSA > 0.6 are suitable, variables with MSA > 0.8 very well suited for factor analysis. The result tells us to drop \"ra\",\"decl\",\"gal_l\", gal_b\", \"ddf\" variables"},{"metadata":{"trusted":true,"_uuid":"53162f0c26cf8bc8541cda644df33f105675556c"},"cell_type":"code","source":"# Create a new dataset with only qualified variables\ndf1 <- df[,c(\"hostgal_specz\",\"hostgal_photoz\",\"hostgal_photoz_err\",\"distmod\",\"mwebv\")]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"69cb4bb13ec86fddfffdf8c7b4535e17c81b129c"},"cell_type":"markdown","source":"Now, we will create a correlation matrix with rest of the variables"},{"metadata":{"trusted":true,"_uuid":"4dd6622f6b75027f90f2a49cf480ff24dc6368f7"},"cell_type":"code","source":"library(corrplot)\ndf1_corr <- cor(df1) # Create a correlation matrix\ncol <- colorRampPalette(c(\"#BB4444\", \"#EE9988\", \"#FFFFFF\", \"#77AADD\", \"#4477AA\"))\ncorrplot(round(df1_corr, 2), method=\"color\", col=col(200),  \n         type=\"upper\", order=\"hclust\", \n         addCoef.col = \"black\", # Add coefficient of correlation\n         tl.col=\"black\", tl.srt=45, #Text label color and rotation\n         diag=FALSE) # hide correlation coefficient on the principal diagonal\nround(df1_corr, 2) # Correlation matrix in table form","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4ad27318afe9a8044422fc0a627d900b45cff82a"},"cell_type":"markdown","source":"Factor analysis will be performed with fa() function from psych package. As parameters are passed: a correlation matrix, the number of factors, and a rotation method."},{"metadata":{"trusted":true,"_uuid":"c08e89ca5e66d518666b65c93803c13a03c5434e"},"cell_type":"code","source":"nfactors <- 2\nnvars <- dim(df1_corr)[1]\nfactors <- fa(r = df1_corr, nfactors = nfactors, rotate = \"Varimax\")\nfactors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b89ffeb00225245b1ab09a141a212b9cf837d89b"},"cell_type":"code","source":"library(ggplot2)\n# Plot Eigenvalues / Represented Variance\neigenvalues <- data.frame(factors$e.values)\ncolnames(eigenvalues) <- c(\"Values\")\neigenvalues$Number <- 1:nrow(df1_corr)\n\neigenvalues$RepresentedVariance <- NA\nfor (i in 1:nrow(df1_corr)) {\n    eigenvalues$RepresentedVariance[i] <- sum(eigenvalues$Values[1:i])/sum(eigenvalues$Values) * \n        100\n}\neigenvalues$RepresentedVariance_text <- paste(round(eigenvalues$RepresentedVariance, \n    0), \" %\")\n\ne1 <- ggplot(eigenvalues, aes(Number, y = Values), group = 1)\ne1 <- e1 + geom_bar(stat = \"identity\")\ne1 <- e1 + geom_line(aes(y = Values), group = 2)\ne1 <- e1 + xlab(\"Number [-]\")\ne1 <- e1 + ylab(\"Eigenvalue [-]\")\ne1 <- e1 + geom_hline(aes(yintercept = 1), col = \"red\")\ne1 <- e1 + geom_text(aes(label = RepresentedVariance_text), nudge_y = 0.2)\ne1 <- e1 + ggtitle(\"Eigenvalues and explained Variance\")\ne1 <- e1 + theme_bw()\ne1 <- e1 + scale_x_continuous(breaks = seq(1, 10, 1))\ne1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0323b3c67fb51082378bdde41a003d2aa7e91259"},"cell_type":"code","source":"library(dplyr)\nlibrary(tidyr)\nloadings_mat <- as.data.frame(matrix(nrow = nvars, ncol =nfactors))\nloadings_mat$Variable <- colnames(df1)\nfor (i in 1:nfactors) {\n  for (j in 1:nvars) {\n    loadings_mat[j, i] <- factors$loadings[j, i]  \n  }\n}\ncolnames(loadings_mat) <- c(\"Factor1\",\"Factor2\",\"Variable\")\nloadings_mat_gather <- loadings_mat %>% gather(\"Factor\", \"Value\", 1:nfactors)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e30b1e947c3efd7cbd4afdfeb5d35258b8452035"},"cell_type":"code","source":"loadings_mat","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7c15a29dce32323a91ec3bdefc892cf40c99483f","scrolled":false},"cell_type":"code","source":"g1 <- ggplot(loadings_mat_gather, aes(Variable, abs(Value), fill=Value))\ng1 <- g1 + facet_wrap(~ Factor, nrow=1)\ng1 <- g1 + geom_bar(stat=\"identity\")\ng1 <- g1 + coord_flip()\ng1 <- g1 + scale_fill_gradient2(name = \"Loading\", \n                       high = \"blue\", mid = \"white\", low = \"red\", \n                       midpoint=0, guide=F) \ng1 <- g1 + xlab(\"Variable\")  # improve x-axis label\ng1 <- g1 + ylab(\"Factor Loading\")  #improve y-axis label\ng1 <- g1 + ggtitle(\"Factors\")\ng1 <- g1 + theme(axis.text=element_text(size=10),\n        axis.title=element_text(size=12, face=\"bold\"))\ng1 <- g1 + theme(plot.title = element_text(size=12))\ng1 <- g1 + theme_bw(base_size=12)\ng1","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"226583815c06f3ea6c9dd79f7ad7c9bef00f569e"},"cell_type":"markdown","source":"Finally, we can name the factors by looking at the variables that are loaded into.\n\n**Factor 1:**  (All about the wavelength of the light)\n1. hostgal_specz: *the spectroscopic redshift of the source*\n2. hostgal_photoz:*The photometric redshift of the host galaxy of the astronomical source*\n3. distmod:*The distance to the source calculated from hostgal_photoz and using general relativity.*\n\n**Factor 2:** (The uncertainty)\n\nhostgal_photoz_err: *The uncertainty on the hostgal_photoz based on LSST survey projections*\n\nThe remaining variables might have been explained in other factors since we have chosen only 2 factors we could see only these factors. We can also see that **mwebv** is negatively loaded on two factors.\nWe can still try with different number of factors to see if it makes more sense for this data. Hope this gives an idea on treating the variables.\n"}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}