Scores are still correlated after PCA I have 22 variables with more than 6000 observations. They are highly correlated. I know these data would work as great explanatory variables to a dichotomous event (present absent). Therefore, I intend to combine them via binary logistic regression and, so to avoid multicolinearity, I thought of "orthogonalizing" the original data using PCA. Then I could choose the main modes of variability (or even all the PC's) and use them as the explanatory variables in my regression, once PCs are orthogonal and independent, by definition. 
I am running this in MATLAB with function pca. The variables are nearly normally distributed and are first normalized between 0 +/- 2 standard deviations (around zero because it will be interesting to keep their signal for future analysis, I believe). I also chose not to center the data within the MATLAB function, once they're already normalized. 
Now here comes the pitfall. My first 2 PCs (scores) are correlated with r=0.7891! Any hints on this? Any suggestions on mistakes I may be making?

EDIT: In case it lightens th scenario, here is a biplot. Letting PCA center the data, I believe the cloud would be just displaced to around the origin, but still keeping its format / correlation, correct?

 A: Remember how PCA works: It starts by finding the single dimension (direction) of greatest variation in your dataset.  That becomes PC1.  Then it finds the direction of greatest variation that is at right angles to that.  It becomes PC2.  Etc.  
The important thing to recognize is that if your data aren't centered, the direction of greatest variation might well be from the origin of the space (e.g., $(0, 0)$ in a Cartesian plane) to the centroid (mean vector) of your data.  Every subsequent principal component is constrained by that first one:  They all have to be at right angles to it.  That means that the resulting PCs may not uncorrelate your data.  
Here is a quick illustration (written in R):  
library(MASS)   # we'll use this package
set.seed(7668)  # this makes the example exactly reproducible
X = mvrnorm(100, mu=c(-5, 0), Sigma=rbind(c(1,  1.6),   # here I generate data
                                          c(1.6,  4) ))

windows(height=4, width=7)
  layout(matrix(1:2, nrow=1))

  plot(X[,1], X[,2], xlim=c(-8, 1))
  abline(h=0, col="gray");  abline(v=0, col="gray")
  points(mean(X[,1]), 0, pch="*", cex=2, col="red")

  biplot(prcomp(X, center=FALSE))


