My question: When doing SVD analysis, you can extract the maximum contributors to each component from the V matrix - does the sign of the the component matter?
Background
I'm currently enjoying the Coursera class on Data Analysis by Jeff Leek. It's a great course for filling in my knowledge gaps. In his lectures on clustering he proposed that one way to use the SVD is to find the most important principal components, and then look at which variable contributes the most to that component by choosing the maximum row-value from the corresponding column in the V matrix. My problem is that he focuses on the maximum value instead of the absolute maximum.
Update: Simulating a pattern
I guess the following simulation answers my question but I would appreciate some feedback
To create a possible scenario I've tried to simulate a pattern where there are a few columns completely random.
Here is the simulation code:
set.seed(12345);
dataMatrix <- matrix(rnorm(400),nrow=40)
colnames(dataMatrix) <-
c(paste("Pos.3:", 1:3, sep=" #"),
paste("Neg.15:", 4:5, sep=" #"),
paste("No pattern:", 6:8, sep=" #"),
paste("Pos.15:", 9:10, sep=" #"))
for(i in 1:40){
# flip a coin
coinFlip <- rbinom(1,size=1,prob=0.5)
# if coin is heads add a common pattern to that row
if(coinFlip){
cols <- grep("Pos.3", colnames(dataMatrix))
dataMatrix[i, cols] <- dataMatrix[i, cols] + 3
}
}
for(i in 1:40){
# flip a coin
coinFlip1 <- rbinom(1,size=1,prob=0.5)
coinFlip2 <- rbinom(1,size=1,prob=0.5)
# if coin is heads add a common pattern to that row
if(coinFlip1){
cols <- grep("Neg.15", colnames(dataMatrix))
dataMatrix[i, cols] <- dataMatrix[i, cols] - 15
}
if(coinFlip2){
cols <- grep("Pos.15", colnames(dataMatrix))
dataMatrix[i,cols] <- dataMatrix[i,cols] + 15
}
}
This generates a simple heatmap with an obvious pattern (the column names indicate the pattern)
After I run the matrix through the svd() function I do a barplot of the V column to examine the values:
svd_out <- svd(scale(dataMatrix))
library(lattice)
key <- simpleKey(rectangles = TRUE, space = "top", points=FALSE,
text=c("Positive", "Negative"))
key$rectangles$col <- c("steelblue", "darkred")
barchart(as.table(svd_out$v[,1]),
horizontal=FALSE, col=ifelse(svd_out$v[,1] > 0,
"steelblue", "darkred"),
ylab="Impact value",
xlab="SVD - percentage explained by V column",
scales=list(x=list(rot=55, labels=colnames(dataMatrix), cex=1.1)),
key = key)
In the plot above the first V column indicates a strong impact from the patterned variables in both directions. The plot below shows the second V column and here the maximum value is a column without a pattern - if we used the absolute value we would select a patterned column.
To conclude: In the lecture slides this line:
maxContrib <- which.max(svd_out$v[,2])
should probably be:
maxContrib <- which.max(abs(svd_out$v[,2]))
Old example
An example based on R code that was used in the lectures
I've used a dataset from the course first assignment together with the Hmisc, lattice and mice package for exploring the issue. You can load the dataset here (although the data needs some data munging):
http <- "https://spark-public.s3.amazonaws.com/dataanalysis/loansData.rda"
con <- url(http)
load(con)
When looking at the first column vector of svd$v very few values are negative:
numvars <- names(loansData)[sapply(loansData, is.numeric)]
# Don't use the outcome variable in any clustering/svd
numvars <- numvars[numvars %nin% c("interest_rate")]
library(mice)
imp <- mice(loansData[, numvars])
c_imp <- complete(imp)
svd_out <- svd(scale(c_imp))
perc_explained <- svd_out$d^2/sum(svd_out$d^2)
barchart(as.table(svd_out$v[,1]),
horizontal=FALSE, col=ifelse(svd_out$v[,1] > 0,
"steelblue", "darkred"),
ylab="Percentage explained",
xlab="SVD - percentage explained by V column",
scales=list(x=list(rot=55, labels=label(loansData[, numvars]))))
The same for the second column:
When I selected the maximum contributor and any variable with at least 90 % of the maxcontributor function using the which.max(abs()) the result seems about right:
and when I do with just the which.max() it looks rather suspicious:
Without using the absolute value http://s4.postimage.org/k8vu343uj/Svd_non_abs_example.png
As we see many maximum contributors are in multiple columns, while this may happen the amount of repetitiveness is not something that I would expect.
Here is the function that I've created to get the plot and the variables of interest:
getSvdMostInfluential <- function(mtrx, quantile,
show_selection=TRUE,
varnames=NULL,
similarityThreshold = 1){
svd_out <- svd(scale(mtrx))
perc_explained <- svd_out$d^2/sum(svd_out$d^2)
cols_expl <- which(cumsum(perc_explained) < quantile)
# Select the variables of interest
vars <- list()
for (i in 1:length(perc_explained)){
v_abs <- svd_out$v[,i]
maxContributor <- which.max(v_abs)
similarSizedContributors <- which(v_abs >= v_abs[maxContributor]*.9)
if (any(similarSizedContributors %nin% maxContributor)){
maxContributor <- similarSizedContributors[order(v_abs[similarSizedContributors], decreasing=TRUE)]
}
vars[[length(vars) + 1]] <- maxContributor
}
if (show_selection){
require(lattice)
# Create transition colors
selected_colors <- colorRampPalette(c("darkgreen", "#FFFFFF"))(length(perc_explained)+2)[1:length(cols_expl)]
nonselected_colors <- colorRampPalette(c("darkgrey", "#FFFFFF"))(length(perc_explained)+2)[length(cols_expl)+1:length(perc_explained)]
names <- unlist(lapply(vars, FUN=function(x){
if (is.null(varnames)){
varnames <- colnames(mtrx)
}
paste(varnames[x], collapse="\n")
}))
las <- 2
m <- par(mar=c(8.1, 4.1, 4.1, 2.1))
on.exit(par(mar=m))
rotation <- 45 + (max(unlist(lapply(vars, length)))-1)*10
if (rotation > 90)
rotation <- 90
p1 <- barchart(perc_explained ~ 1:length(perc_explained),
horiz=FALSE,
ylab="Percentage explained",
xlab="SVD - percentage explained by V column",
col=c(selected_colors, nonselected_colors),
key=list(text=list(c("Selected", "Not selected")),
rectangles=list(col=c("darkgreen", "#777777"))),
scales=list(x=list(rot=rotation, labels=names)))
print(p1)
}
return(unique(unlist(vars)))
}
getSvdMostInfluential(c_imp, 0.8, varnames=label(loansData[, numvars]))
I've posted this question on the course forums but didn't get an answer. Note, this is not homework.
