Analysing a Regression Analysis Model

I'm playing around with hospital waiting lists trying to find out what factors affect waiting times. 

Using the Pearson Correlation of my features (which are the ethnic make up of a hospital waiting list and the time spent waiting), the data looks like this:

Pearson Correlation Heatmap: Ethnicity and waiting time

What if I normalise the rows?

Pearson Correlation Heatmap: normalized rows

Well that was silly. Of course there will be a negative correlation between ethnicities as the total needs to sum to 1.0.

Anyway, it was at this point I found the data was dirty beyond salvage due to upstream cleaning processes gone wrong. Having got a new data set, I tried again (this time ignoring the Pearson correlation between ethnicities):

Pearson correlation: ethnicities vs waiting list time

Note that this time, the data was standardized. It looks like waiting list time goes up for white people.

Inferential Linear Regression 

Putting Pearson correlation to one side, let's see a linear regression model trained on this data. Note, in Spark's linear regression algorithm, one is able to use ElasticNet which allows a mix between L1 (lasso) and L2 (ridge regression) regularization.

Using dummy encoding, the results look like this:

coefficient             category   p-value
2.9759572500367764 White 0.011753367452440378
0.607824511239789 Black 0.6505335882518613
1.521480345096513 Other 0.2828384545133975
1.2063242152220122 Mixed 0.4440061099633672
14.909193776961727 intercept  0.0

Hmm, those p-values look far from conclusive. OK, let's make the most populous group the zero-vector:

coefficient             category  p-value

-3.133385162559836 Mixed   0.01070291868710882

-1.7512494284036988 Black   0.10414898264964889

-1.4364038661386487 Other   0.08250720507990783

-2.3504004542073984 Asian   0.0006682319884454557

17.88573117661592 intercept 0.0

Now this is more interesting. The p-values are still a little too high to draw conclusions about two groups but it's starting to look like the waiting list size is lower if you are Asian or Mixed ethnicity.

Adding other columns makes the ethnicities at least look even more certain although the p-value for these new categories - including socioeconomic and age  -were themselves not particularly conclusive.

coefficient             category  p-value
-4.309164982942572 Mixed 0.0004849926653480718
-2.3795206765105696 Black 0.027868866311572482
-2.2066035364090726 Other 0.008510096871574113
-2.9196528536463315 Asian 2.83433623644580e-05
20.342211763968347 intercept 0.0

Now those p-values look pretty good and they're in agreement with my Pearson correlation. 

"A beautiful aspect of regression analysis is that you hold the other independent variables constant by merely including them in your model! [Alternatively,] omitting an important variable causes it to be uncontrolled, and it can bias the results for the variables that you do include in the model."[Regression Analysis, Jim Frost]

Work continues.

Statistical Covariance Part 2

I quoted in a previous post somebody who said that covariance matrices always have a determinant of 0. This isn't strictly true. What can I say? I'm not a very good mathematician.

The argument was that the means are always substracted from each row and with simple algebra you could demonstrate that the determinant is 0. But this ignores the fact that covariance is the expected value of two rows having their expected values subtracted. 

Say we have two different distributions from which we draw combinations. However, for reasons peculiar to the use case, certain combinations are not allowed. When we tabulate the probabilities for this state space, we'll have zeros in some cells. Even though the inner product of the corresponding two probability vectors might not be zero, the expected probability of the two together is.

Note another caveat. Bearing in mind that correlation is just covariance divided by the root of the product of both variances, "the correlation between A and B is only a measure of the strength of the linear relationship between A and B. Two random variables can be perfectly related, to the point where one is a deterministic function of the other, but still have zero correlation if that function is non-linear." [Daniel Shiebler's blog]

Positive Definite

A matrix is positive definite if x^T M x > 0 (an equivalent definition is "a symmetric matrix whose eigenvalues are all positive real numbers" - Coding the Matrix, Klein, definition 12.6.1) 

Note that all matrices that can be expressed as A^TA are positive semi definite (see this SO answer). The proof is simple: substitute A^TA for M above. 

x^T M x = x^T A^TA x = (x A)^T(A x)  

and any non-zero vector or real number multiplied by itself is positive. Since our covariance matrix can be expressed as A^TA it too is at least positive semidefinite (x^T M x ≥ 0). But we know it must also be positive definite as you can't invert M x = 0 for x ≠ 0.

Why this last statement is true can be explained at this elegant SO answer. Basically, if M x = 0 for x ≠ 0 then each row of M must be linarly dependent on each other for the equation 

M _i,j x_j = 0 ∀ i 

to hold. That is, if you give me (n-1) values in a row, I can tell you the value of the last one. Here's a simple Python/Numpy example where the last value is a row is just the sum of the first two. It's easy to see that when multipled with the vector [1, 1, -1] this matrix would be in the null space:

>>> M = np.asmatrix([[1, 2, 3], [4, 5, 9], [6, 7, 13]])
>>> np.linalg.det(M)
0.0

But if a matrix is linearly dependent, it's determinant must be 0. And if a matrix's determinant is 0, it cannot be inverted. QED.

What if my Matrix is Singular?

We add a tiny amount to make it non-singular. This is called conditioning. It sounds like a hack but it does have a basis in maths. In a frequentist interpretation, this is ridge regression. In a Bayesian interpretation, it's the prior.

Briefly, the argument goes that for ridge regression, we penalize large model parameters. So, instead of minimizing our error (θ X - y_actual) we minimize our error plus the penalty:

y_estimate = (θ X - y_actual)² + λ θ^T θ

by differentiating with respect to θ. Solve this equation and you'll see a [X^T X + λ I]^-1.

The argument for the Bayesian prior briefly goes like this: if we take a frequentist view and assume that the error in our data is Gaussian and plug y_estimate into it, we'll see our familiar equation for a Gaussian multiplied by e^λθTθ. Since the Bayesian posterior,  p(θ|Data,Model) must equal the frequentist probability, e^λθTθ is the only term that maps to p(θ|Model) simple because it's the only one with θ in it. Therefore, our conditioning manifests itself in the prior.

Full derivations appear in Jake VanderPlas' wonderful blog.

Ridge regression

"Ridge regression adds an additional λI to the matrix X^TX so that it's non-singular, and we can take the inverse of the whole thing" [1]

It adds bias. Take a matrix that has rank 1 (that is, each row adds no more information).

import numpy as np
from numpy.linalg import matrix_rank, det, inv, norm

mat_rank_1 = np.matrix('1    2    3;'
                       '10   20   30;'
                       '100, 200, 300')

We can demonstrate it does indeed have rank 1:

print mat_rank_1, "has rank ", matrix_rank(mat_rank_1)  # it is indeed 1

Now if were to use it, say, in linear regression, we'd have problems when we invert it:

mTm = np.dot(mat_rank_1.T, mat_rank_1)

print "\ndeterminant of rank 1 matrix multiplied by its transpose:", det(mTm)  # and not too surprisingly, the determinant is 0

The following blows up with "numpy.linalg.linalg.LinAlgError: Singular matrix":

inv(mmt)

So, we can add some bias:

l = np.eye(3, 3) * 0.1

wbTwb = np.dot(mat_rank_1.T, mat_rank_1) + l

print wbTwb, "has determinant", det(wbTwb), "and rank", matrix_rank(wbTwb)

print inv(wbTwb)

and this doesn't blow up.

"Ridge regression was originally developed to deal with the problem of having more features than data points. But it can also be used to add bias into our estimations. We can use the λ value to impose a maximum value on the sum of all our [weights]" [1]

John D Cook warns: "A little noise makes the system go from theoretically impossible to solve to theoretically possible to solve, but that may not be very useful in practice. It will let you compute a solution, but that solution may be meaningless... You use knowledge of your domain beyond the specific data at hand to guide how you change your matrix.

Ill Conditioned

He goes on: "The condition number of a matrix is the norm of the matrix times the norm of its inverse... A small change to a matrix might not change its norm much, but it might change the norm of its inverse a great deal. If that is the case, the matrix is called ill-conditioned because it has a large condition number.

"You can think of condition number as an error multiplier... Note that condition number isn’t limited to loss of precision due to floating point arithmetic. If you have some error in your input b, say due to measurement error, then you will have some corresponding error in your solution to Ax = b, even if you solve the system Ax = b exactly. If the matrix A is ill-conditioned, any error in b (rounding error, measurement error, statistical uncertainty, etc.) will result in a correspondingly much larger error in x."

In our case, it we can calculate this condition number thus:

print "norm of original matrix", norm(wbTwb)
print "norm of inverse of original matrix", norm(inv(wbTwb))
print "Condition number", norm(wbTwb) * norm(inv(wbTwb))

Which results in:

norm of original matrix 141414.1
norm of inverse of original matrix 14.1421356238
Condition number 1999897.38132

Yoiks. That modest λ=0.1 made our matrix quite ill conditioned.

Interestingly, making our original matrix rank 3 by avoid linear dependencies between rows:

mat_rank_1 = np.matrix('1    2    3;'

                       '20   10   30;'

                       '100, 300, 200')

reduces the condition number by an order of magnitude even for the same value of λ.

What's more, making λ much bigger (say 10) makes the original matrix less ill-conditioned by two orders of magnitude. But of course, now we're significantly distorting our data - the "meaningless" that Cook talks of.

An excellent article on how ill-conditioning can effect neural networks can be found here.

[1] Machine Learning in Action