Why set cannot be a functor

This is well known but keeps tripping me up when I think of an example so here it is.

Set demands that its elements implement Eq (or equals() or whatever is semantically equivalent). List has no such requirement. A bogus equals implementation will lead to Set violating the functor laws. 

Imagine an implementation that looked something like:

class OddEquals(val id: Int, val label: String) {

  override def equals(a: Any): Boolean =

    if (a.isInstanceOf[OddEquals]) {

      val other = a.asInstanceOf[OddEquals]

      id == other.id

    } else false

  override def hashCode(): Int = id

}

which would be a fairly typical piece of code for objects whose identity depends solely on their primary key.

Now, image two services. One finds an entity given any related string, we'll call it f. I'll fake it like this:

  type F = String => OddEquals

  val f: F = x => new OddEquals(x.length, x)

The other service simply extracts the entity's name:

  type G = OddEquals => String

  val g: G = _.label

The functor law says:

xs.map(x => g(f(x))) == xs.map(f).map(g))

But if xs is Set("hello", "world"), the left-hand side is the same as the input but the right-hand side is Set("hello").

Compare this to List("hello", "world") and you'll see that both left- and right-hand side are equal (that is, same as xs)

Python's Garbage

Python's GC appears to be hugely different to the GC in Java we know and love. I have some PySpark code that collects data from a data set back in the driver and then processes it in Pandas. The first few work nicely but ultimately PySpark barfs with memory issues and no amount of increasing the driver memory seems to solve the problem. This is odd, because each collect itself only retrieves a modest amount of data. After that, I'm not interested in it.

"Reducing memory usage in Python is difficult, because Python does not actually release memory back to the operating system. If you delete objects, then the memory is available to new Python objects, but not free()'d back to the system.

"As noted in the comments, there are some things to try: gc.collect (@EdChum) may clear stuff, for example. At least from my experience, these things sometimes work and often don't. There is one thing that always works, however, because it is done at the OS, not language, level.

Then if you do something like

import multiprocessing
result = multiprocessing.Pool(1).map(huge_intermediate_calc, [something_])[0]

Then the function is executed at a different process. When that process completes, the OS retakes all the resources it used. There's really nothing Python, pandas, the garbage collector, could do to stop that." [StackOverflow]

However, this does not work well for turning PySpark DataFrames to Pandas data frames. Spark is lazy so the DataFrame cannot be serialized and given to the new process (you cannot serialize, for instance, an open network port).

Ultimately, I had to cut my losses and return from StatsModels to Spark's logistic regression implementation.

In addition to this, we have the horrors of what happens with a Spark UDF when written in Python:

"Each row of the DataFrame is serialised, sent to the Python Runtime and returned to the JVM. As you can imagine, it is nothing optimal. There are some projects that try to optimise this problem. One of them is Apache Arrow, which is based on applying UDFs with Pandas." [Damavis]

Yuck. 

Logistic Regression in Action

Once again, I'm looking at patient data waiting times. This time, it's the A&E department. I'm not so much interested in how long they wait but rather the binary condition of whether they waited more or less than four hours.

Logistic Regression with Spark 

I used Spark's GeneralizedLinearRegression. Hyperparameter tuning did not produce much difference in output.

Splitting the data 90/10 for train and test gave me the following output on the test curves:

Prediction where true value is 1

The above shows the predictions for patients that we know really did wait more than 4 hours for treatment. The below shows predictions where we really know they waited less:

Predicition where true value is 0

Clearly, the model has spotted something but the overlap between the two states is large. What should the threshold be to best differentiate one from the other? Fortunately, I'm not trying to predict but to infer.

You can evaluate how good the model is with some simple PySpark code:

from pyspark.ml.evaluation import BinaryClassificationEvaluator
evaluator = BinaryClassificationEvaluator(labelCol=dependent_col,                                                    rawPredictionCol=prediction_col)
auroc = evaluator.evaluate(predictions, {evaluator.metricName: "areaUnderROC"})
auprc = evaluator.evaluate(predictions, {evaluator.metricName: "areaUnderPR"})

Evaluating with Spark's BinaryClassificationEvaluator gave me an area under the ROC curve of 0.79. Hardly stellar but not terrible.

"A no-skill classifier will be a horizontal line on the plot with a precision that is proportional to the number of positive examples in the dataset. For a balanced dataset this will be 0.5... Precision-recall curves (PR curves) are recommended for highly skewed domains where ROC curves may provide an excessively optimistic view of the performance." [MachineLearningMastery]

Probabilities, Odds, and Odds-Ratios

"Odds are the probability of an event occurring divided by the probability of the event not occurring. An odds ratio is the odds of the event in one group ... divided by the odds in another group." [Making sense of odds and odds ratios]

"The [linear regression] model is specified in terms of K−1 log-odds or logit transformations (reflecting the constraint that the probabilities sum to one)... The choice of denominator is arbitrary in that the estimates are equivariant under this choice." [Elements of Statistical Learning]

K here is the number of groups. So, for the case where we're looking at true/false outputs, K=2 and the equations in ESL become the simple equation in the following section.

Model output

An example: "given estimated parameters and values for x₁ and x₂ , the model might predict y = 0.5, but the only meaningful values of y are 0 and 1.  It is tempting to interpret a result like that as a probability; for example, we might say that a respondent with particular values of x₁ and x₂ has a 50% chance. But it is also possible for this model to predict y = 1.1 or y = −0.1, and those are not valid probabilities. Logistic regression avoids this problem by expressing predictions in terms of odds rather than probabilities." [Think Stats, 2nd Edition, Downey]

In a two-state output, Downey's equation would be:

y = log(odds) = log(p / (1-p)) =  β₀ + β₁ x₁ + ...

Odds(-ratios) are used when reporting effect size. "when X, the predictor is continuous, the odds ratio is constant across values of X.  But probabilities aren’t" because of that non-linear relationship above between odds and p.

"It works exactly the same way as interest rates.  I can tell you that an annual interest rate is 8%.  So at the end of the year, you’ll earn $8 if you invested $100, or $40 if you invested $500.  The rate stays constant, but the actual amount earned differs based on the amount invested. Odds ratios work the same." [TheAnalysisFactor] Since there is no bound on x₁ etc, you could easily have a value greater than 1.0 or less than 0.0 which just doesn't make sense for probabilities.

Interpreting the coefficients

I have 36 features and the p-values of approximately 30 of them are all less than 0.05. This seems suspiciously good even if sometimes their corresponding coefficients are not particularly significant. "If that difference between your coefficient and zero isn't interesting to you, don't feel obligated to interpret it substantively just because p < .05." [StackOverflow]

"How does one interpret a coefficient of 0.081 (Std. Error = 0.026)? ... Incorporating the standard error we get an approximate 95% confidence interval of exp(0.081 ± 2 × 0.026) = (1.03, 1.14).

"This summary includes Z scores for each of the coefficients in the model (coefficients divided by their standard errors)... A Z score greater than approximately 2 in absolute value is significant at the 5% level." [Elements of Statistical Learning]

Intepreting the intercept

If we take the example here [QuantifyingHealth] where we're looking at the relationship between tobacco consumption (x₁) and heart disease (y) over a 10 year period, then the interpretation of the intercept's coefficient depends on several conditions. 

For example, say that the intercept's coefficient is -1.93. This leads us to believe the probability of heart disease for a non-smoker (x₁=0) is 0.13 by plugging -1.93 in the above equation then rearranging.

If the smoking consumption is standardized, x₁=0 "corresponds to the mean lifetime usage of tobacco in kg, and the interpretation becomes: For an average consumer of tobacco, the probability of having a heart disease in the next 10 years is 0.13" [ibid].

If we're only interested in smokers not the general population, then x₁≠0. "Let’s pick the maximum [623kg] as a reference and calculate the limit of how much smoking can affect the risk of heart disease... We get P = 0.22. The interpretation becomes: the maximum lifetime consumption of 623 kg of tobacco is associated with a 22% risk of having heart disease".

Logistic Regression Notes

Just a few, miscellaneous notes on logistic regression models that I'm using Spark to build for me.

Interpreting Coefficients

The coefficients for one-hot encoded features can simply be compared to show the strength of the effect. "Things are marginally more complicated for the numeric predictor variables. A coefficient for a predictor variable shows the effect of a one unit change in the predictor variable." [DisplayR]

"How does one interpret a coefficient of 0.081 (Std. Error = 0.026) for tobacco, for example? Tobacco is measured in total lifetime usage in kilograms, with a median of 1.0kg for the controls and 4.1kg for the cases. Thus an increase of 1kg in lifetime tobacco usage accounts for an increase in the odds of coronary heart disease of exp(0.081) = 1.084 or 8.4%. Incorporating the standard error we get an approximate 95% confidence interval of exp(0.081 ± 2 × 0.026) = (1.03, 1.14)." [1]

Computational Problems

How do we avoid problems with singular matrices? "it's computationally cheaper (faster) to find the solution using the gradient descent in some cases... it works even when the design matrix has collinearity issues." [StackExchange]

"Switching to an optimizer that does not use the Hessian often succeeds in those cases. For example, scipy's 'bfgs' is a good optimizer that works in many cases"[StackOverflow]

Generalized Linear Models

"Where the response variable has an exponential family distribution [Gaussian, Bernoulli, gamma etc], whose mean is a linear function of the inputs ... this is known as a generalized linear model, and generalizes the idea of logistic regression to other kinds of response variables." [Machine Learning-A Probabilistic Perspective - Murphy]

"Software implementations can take advantage of these connections. For example, the generalized linear modeling software in R (which includes logistic regression as part of the binomial family of models) exploits them fully. GLM (generalized linear model) objects can be treated as linear model objects, and all the tools available for linear models can be applied automatically." [1]

Z-Scores

Z-Scores are "coefficients divided by their standard errors. A nonsignificant Z score suggests a coefficient can be dropped from the model.  Each of these correspond formally to a test of the null hypothesis that the coefficient in question is zero, while all the others are not (also known as the Wald test). A Z score greater than approximately 2 in absolute value is significant at the 5% level." [1]

Although this sounds like the job of a p-score, note that they are different (but related) to z-scores. A "p-value indicates how unlikely the statistic is. z-score indicates how far away from the mean it is. There may be a difference between them, depending on the sample size." [StackOverflow]

You can see quite clearly that jiggling sample size can profoundly change the numbers of the metrics. It's very useful to be able to quickly create fake but realistic data and play around and watch this effect. Given a large amount of data, you might want to reduce the p-value cut off since "in large samples is more appropriate to choose a size of 1% or less rather than the 'traditional' 5%." [SO]

[1] Elements of Statistical Learning

Miscellaneous Spark

Feature names in Spark ML

Why isn't this better known? This is a great way to extract the features names from a Dataset after the features have been turned into vectors (including even one-hot encoding):

    index_to_name = {}
    for i in df.schema[OUTPUT_COL].metadata["ml_attr"]["attrs"]:
        meta_map = df.schema[OUTPUT_COL].metadata["ml_attr"]["attrs"][i]
        for m in meta_map:
            index = m['idx']
            name = m['name']
            index_to_name[index] = name 

where OUTPUT_COL is the column containing a Vector of encoded features.

Non-deterministic Lists

I've seen this Spark-related bug twice in our code base: even though lists have an order, collect_list "is non-deterministic because the order of collected results depends on the order of the rows which may be non-deterministic after a shuffle." [StackOverflow]

The solution in the SO answer is to run the collect over a Window that performs an orderBy.

Journeys in Data Engineering

I'm currently helping a huge health provider with its data issues. Here are some things I've learned:

Don't bikeshed. There's no point wondering whether you've captured the correct ethnic description ("White? Irish White? Eastern European White?") when there are bigger fish to fry. For instance, hospitals were putting the patients' personal information in the primary key field. This sent the data officer apoplectic until it was cleansed from the files before being sent downstream. But as a result, the downstream keys were not unique. The data scientists consuming this data aggregated it and came to conclusions unwittingly based on one individual having visited the hospital over three million times.

Don't proactively go looking for data quality issues. They'll come looking for you. Continue to build your models but be extremely circumspect. This is a more efficient process than spending time wondering if the data looks good.

Just because the data looks good, it doesn't mean it's usable. How often is the data updated? Is it immutable or is it frequently adjusted? In short, there's an orthogonal axis in data space separate to data quality and it's a function of time. Perfect data that becomes old is (probably) no longer perfect. Perfect data that becomes incomplete is (probably) no longer perfect.

Give a thought to integration tests. Platforms like Palantir are pretty under-developed in this area (answer to my requests for an integration test platform: "we're working on it"). So, you may need to write bespoke code that just kicks the tires of your data every time you do a refresh. 

Remember that documentation is the code. I've had a nice experience using ScalaTest with its Given, When, Thens. It ensured that when running integration tests, the code generated the documentation so the two would never fall out of synch. This is, of course, much harder (impossible?) when running in a walled garden like Palantir.

Stick with a strongly typed language. This can be hard since pretty much all data scientists use Python (and therefore PySpark) but there is nothing worse than waiting several minutes for your code to run only to found out that a typo trips you up. In a language that is at least compiled, such problems would not occur. I'd go so far to say that Python is simply not the correct tool for distributed computing.

Python tools like Airflow are made much easier by using Docker since Python library version management is a pain in the arse. Far better for each application to have its own Docker container.

Never forget that your schema is your data contract as much as an API is your coding contract. Change it, and people downstream may squeal.

Persisting intermediate data sets hugely help debugging. 

Finally, don't use null to indicate something for which you know the meaning. If, for instance, when a record is in a certain state, it's better to say value X equals a flag to indicate that. If you use null, somebody reading that data doesn't know if the record is in this state or if the data is just missing.

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.