Evaluating the distribution of a continuous variable in a two dimensional space

I have performed a Principal Component Analysis on a set of hydrological indices. Those hydrological indices are derived from the discharge of some rivers (e.g. how long the river needs to get back to its mean discharge after a flood). This gave me two principal components that account for 80 % of the variance in those indices. When I plot them against each other I get the following plot.

Every point of this plot refers to a certain river catchment and the two axes are the first and second principal component. The next step is to use this plot and overlay it with some attributes of those same river catchment like the area of the catchment. (Yellowish colors indicates high values, reddish colors indicate lower values).

Or the aridity of the catchment.

When looking at those two images it becomes obvious that the aridity of the catchment follows a pattern, while the area of the catchment does not.

So, my question is: How can I quantify differences in the pattern of the aridity in comparison to the area?

For my first try I seperated the continuous variable of area and aridity and split them into categories. Then I evaluated the "clusterlikeness" of the categories with the Calinski-Harabaz Score. This got me some good results. It has the drawback though that the amount of bins I split the area and the aridity into is arbitrary, which I want to avoid.

So, is there a way I can evaluate the distribution of the pattern of the continuous variables (aridity, area or other) directly?

Basically I just want a number that tells me reliably that the distribution of the aridity is less random then the distribution of the area.

• PCs should be uncorrelated, but your scatter plot shows clear correlation. Something is wrong. – amoeba Nov 20 '18 at 14:11
• What could be the cause of this? I simply used my six hydrological indices for the PCA and told sklearn I wanted as many principal components that 80 % of the variance is explained. – F. Jehn Nov 20 '18 at 15:05
• I also used the StandartScaler from sklearn. – F. Jehn Nov 20 '18 at 15:11
• Seems like the visual impression is missleading. I just checked it and PC 1 and PC 2 have a correlation of basically 0. In the pictures a few datapoints are missing. I cropped it to make the picture smaller. – F. Jehn Nov 20 '18 at 15:14
• Well, I would not expect the correlation between PC1 and PC2 to be anything else than zero ;-) – January Nov 21 '18 at 15:15

One thing that comes to my mind is PCA regression.

Another thing that I myself like to use for strictly exploratory purposes (and which is frowned upon by serious people, because it amounts to data dredging) is do the reverse: test the regression of first $$N$$ interesting components on the variables of interest (aridity, area) and see which components have the highest $$r^2$$. You can then plot $$r^2$$ versus component number and see "how far back" your variables of interest are.

However in most of my applications, the variables of interest are truly independent, and the variables which go into PCA are dependent. For example, a variable of interest can be "sex", and the PCA variables are measurements of transcripts in blood which are dependent on sex, but not the other way round, therefore I regression of PCA components over the variables of interest is somewhat justified.

If I understand your problem correctly, your situation is reverse: the variables you use for coloring are the dependent variables, and the geological indices which you insare the independent variables.

However, the tools recommended in this response may be of use.

Code example in R:

## simulate a PCA
set.seed(321)
x <- rnorm(200)
y <- rnorm(200)
z <- rnorm(200)
xx <- replicate(10, rnorm(200, x))
yy <- replicate(10, rnorm(200, y))
zz <- replicate(10, rnorm(200, z))
data <- cbind(xx, yy, zz)
pca <- prcomp(data, scale.=T)

## calculate the r²
getR2 <- function(pca, x) {
m <- lm(pca$$x ~ x) ms <- summary(m) sapply(ms, function(i) i$$r.squared)
}

r2 <- getR2(pca, x)
plot(r2, type="b", xlim=c(1,5), xlab="Component", ylab="r²",
ylim=c(0, 1), bty="n")

## randomize the predictor n times
n <- 100
x.rand <- replicate(n, sample(x))
r2.rand <- apply(x.rand, 2, function(x) getR2(pca, x))
apply(r2.rand, 1, function(r2) lines(1:n, r2, col="grey"))


This last bit gives us an idea about whether the observed r² values are significant. But note that this is not proper significance testing, this is just an exploratory heuristic. Below is the result; on the left, the code run as above; on the right, the same code, but data is now cbind(yy, yy, zz) (no x involved):

The right plot shows what happens if you don't have any PC which can be explained by the independent variables.

As you can see on the left plot, we have "interesting" components 2 and 3. If we plot these components coloring by the value of x, we get:

library(tagcloud) # for smoothPalette
library(pca3d)    # for pca2d
pal <- c("black", "purple", "red", "orange", "yellow")
cols <- smoothPalette(x, palfunc=colorRampPalette(pal))
pca2d(pca, components=c(2,3), col=cols)


So, how does that all help you? Well, you were asking for a number that tells you that one independent variable, area, is "random" (which I take to mean: not influencing the dependent variables that go into the PCA), while the other, aridity, is not. However, you can't be sure that, for example, aridity can explain the variance in components 1 and 2, while area can explain the variance in components 3 and 4. The calculations above allow you to:

• visually identify the components which can be explained by the independent variables;
• calculate the fraction of the variance of the given component (r²) which can be explained by the independent variable;
• get a rough idea (again, exploratory mode only) whether this fraction is unlikely to have arisen by chance.
• I would rather say the variables I use for coloring a independent. The behaviour of a river (which is captured by the hydrological indices) is dependent on the climate (like aridity) and geology of the river catchment. I mainly want some number that indicates that the distribution of the aridity is less random than the distribution of the area. – F. Jehn Nov 20 '18 at 14:08
• In case someone comes across this and wonders about PCA regression. I found a Python package, which is mentioned here: stats.stackexchange.com/a/375781/227551 – F. Jehn Nov 21 '18 at 9:59
• Tried the PCR. Does not deliver any reasonable results. – F. Jehn Nov 21 '18 at 13:44
• Did you try the other approach (trying to find the PC components which can be explained by your independents)? I can give you the R code for this. – January Nov 21 '18 at 15:16
• Take a look at fig 1 in this paper, the left side shows the kind of plots I am talking about. – January Nov 21 '18 at 15:18