Lambda Island

With the knowledge of transducers under your belt, it’s time to start analyzing some data. This episode provides a first introduction to the Kixi.stats statistical toolkit, by analyzing a data set and its distribution, with the goal of creating a predictive model through linear regression.

Transcript

In today’s episode we take our first steps into data science using the kixi.stats library, which implements over two dozen statistical functions, as well as many statistical distributions to sample from.

The author of kixi.stats is Henry Garner, who also wrote the Redux library I covered in the previous episode. Henry wrote the “Clojure for Data Science” book, and he has been a tremendous help in preparing this episode, spending more time to answer my naive statistics questions than I could reasonably have asked for.

This episode will out of necessity only scratch the merest surface of the vast field that is Data Science, but hopefully it’s enough to give you an appetite for more. If you do I recommend checking ou Henry’s book, “Clojure for Data Science”, which explains a lot of the ideas and principles in more depth. You can find a link in the show notes. In this episode and the next I’ll work through an example of Linear Regression, which is covered in much more depth in chapter 3 of the book.

The kixi.stats API is clustered around three namespaces:

kixi.stats.core contains a few dozen standard statistical functions, things like the mean, median, variance, and standard deviation.

kixi.stats.math defines some mathematical operations like the square root, logarithm, and exponent. These are used by kixi internally but can be handy in your own code as well, especially since they smooth over certain differences between Clojure and ClojureScript.

kixi.stats.distribution defines well known statistical distributions, like a Normal, Bernoulli, or Poission distribution, allowing you to play with the parameters that define the distribution, and sample from the result.

How easy you’ll find this episode will largely depend on your knowledge of statistics. That said, as always I’ll do what I can to make things generally understandable.

I’ve prepared a GitHub repo with everything you need to follow along with today’s episode, as always check out the show notes to find the link.

This project contains the data we’ll be working with, and has the right dependencies set up. It also contains a few utility namespaces you’ll be using later on.

$ git clone https://github.com/lambdaisland/ep43-data-science-kixi-stats.git

This project uses deps.edn, so you do need to make sure the Clojure CLI tools are installed on your system.

Among the dependencies you find kixi.stats, which is the star of today’s episode, as well as xforms and redux which I covered last time, and the thi.ng/geom library which you’ll use to visualize the data.

$ cat deps.edn
{:deps  {org.clojure/clojure     {:mvn/version "1.10.0-alpha6"}
         kixi/stats              {:mvn/version "0.4.3"}
         net.cgrand/xforms       {:mvn/version "0.18.2"}
         redux                   {:mvn/version "0.1.4"}
         org.clojure/data.csv    {:mvn/version "0.1.4"}
         thi.ng/geom             {:mvn/version "1.0.0-RC3"}}
 :paths ["src" "resources"]}

Under resources you’ll find the “abalone” data set, which contains measurements of physical features of 4000-odd abalones, a type of shelled sea creature.

The abalone namespace contains some code to read in and parse the data. graphs has functions for drawing common graphs like histograms and box-and-whiskers plots. The core namespace is where you can add your code.

$ tree
.
├── deps.edn
├── README.md
├── resources
│   ├── abalone.data
│   ├── abalone.names
│   └── README.md
└── src
    └── ep43_data_science
        ├── abalone.clj
        ├── core.clj
        └── graphs.clj

This data set contains things like the length, weight, diameter, and the number of “rings” in the shell of the abalone. These are kind of like year rings on a tree, they are a good indicator of the creature’s age, but counting them is hard, it involves slicing the shell in half, applying a dye to stain the rings, and then counting them under a microscope. What we’d like to find out is whether some of these other attributes are a good predictor for the amoung of rings, and hence for the age.

(ns ep43-data-science.core
  (:require [ep43-data-science.abalone :as abalone]
            [ep43-data-science.graphs :as g]
            [kixi.stats.core :as kixi]
            [kixi.stats.distribution :as kixi.dist]
            [kixi.stats.math :as m]
            [net.cgrand.xforms :as x]
            [redux.core :as redux]))

(rand-nth abalone/data)
;; => {:rings 9,
;;     :sex "I",
;;     :shell-weight 14.7,
;;     :diameter 54.0,
;;     :whole-weight 54.2,
;;     :viscera-weight 10.5,
;;     :length 71.0,
;;     :shucked-weight 28.5,
;;     :height 21.0}

By now you should have the github repo cloned, and the core namespace open in your editor. This is a hands-on episode, so if you haven’t done so now is a good time. Launch a connected REPL so you can evaluate forms.

Let’s try to get a better sense of the data.

We have data about 4177 abalones, 4177 data points, each represented as a map. Lenght measurements, like the length, height, and diameter are in milimeter, and have a milimeter precision, whereas weights are in grams, and have a precision of 1/10th of a gram, the ring count is always a whole number, and finally the sex, which can be M, F, or I, standing for male, female, and infant.

(count abalone/data)
;; => 4177

What you want to do first is get a general sense of the magnitude and shape of your data. What are the minimum and maximum values for each variable? Do certain values occur more often than others?

Let’s run this through kixi.stats to figure out some of the data’s characteristics. kixi.stats implements “reducing functions”, meaning you use them as the second argument to transduce. Let’s calculate the mean value for the rings.

Use a map transducer to pull out the :rings values, and use kixi/mean to reduce this “stream” of integers to a single value.

If this sample of 4000 abalones is representative for the population as a whole, then you can say that an average abalone has about ten rings.

(transduce (map :rings) kixi/mean abalone/data)
;; => 9.933684462532918

If you watched episode 42, “Inside Transducers”, then you know there’s no magic here . kixi/mean is simply a reducing function, still it’s worth asking why the API looks like this. Why not provide a mean function which simply takes a collection?

The reason is that kixi.stats is built for “awkward sized data”, data that’s small enough to process on a single machine, but too large to conveniently fit in memory. The streaming nature of transducers makes this easy, and you get the composability for free.

Let’s have a quick look under the hood to see how this works. mean in this case is an alias for the arithmetic-mean, which is the actual reducing function.

Remember that a reducing function has three arities: one to produce an initial value, one to process each “step”, and finally a “completing” arity that yields the final result. This reducing function uses a two-element vector to keep track of the total sum, and the count. As a final step it divides one by the other to yield the mean.

(def arithmetic-mean
  "Calculates the arithmetic mean of numeric inputs."
  (fn
    ([] [0.0 0.0])
    ([[^double s ^double c :as acc] e]
     (if (nil? e)
       acc
       (let [e (double e)]
         [(+ s e) (inc c)])))
    ([[s c]]
     (when-not (zero? c)
       (/ s c)))))

You can try that out by calling kixi/mean directly without any arguments, what you get is the initial accumulator.

You can thread this vector through calls to the two-argument step function to see the accumulator build up the total and the count.

Finally this gets passed to the one-argument completing step to yield the final result.

(kixi/mean)
;; => [0.0 0.0]

(-> (kixi/mean)
    (kixi/mean 5)
    (kixi/mean 7)
    (kixi/mean 12))
;; => [24.0 3.0]

(kixi/mean [24 3])
;; => 8

If you replace the completing step with the identity function then you can see that it computed the total for all 4177 abalones. This pattern of tracking multiple parameters and only computing the final result all the way at the end is used throughout kixi.stats. In fact some functions share the same underlying reducing function, only with a different completing step.

(transduce (map :rings) (completing kixi/mean identity) abalone/data)
;; => [41493.0 4177.0]

The arithmetic mean or “average” is broadly known and understood, but it can easily be distorted by outliers. Instead you can look at the median, which doesn’t have this problem.

The median is the value that splits the population in half, meaning that 50% of all abalones will have 9 rings or fewer, whereas 50% have 10 or more.

There’s something strange about this result though, the :rings values are all integers, so the median should be an integer too, but instead we got this decimal number. On top of that you might notice that when you run this code you get a different result. Even when I re-run this here I’ll get slightly different results each time.

To understand why that is first think about how you would normally compute the median. You would take all elements in your collection, sort them, and then take the middle one. This means keeping all data in memory, which is exactly what we’re trying to avoid.

Instead kixi.stats uses a special data structure under the hood called t-digest, which provides approximate answers about the shape of your data without having to keep everything in memory, so you trade off some accuracy for memory use.

(transduce (map :rings) kixi/median abalone/data)
;; => 9.125

For a better overview you can look at what statisticians call the “five number summary”: these are the minimum, first quartile, median, third quartile, and maximum.

:min and :max are the smallest and largest values in the data, the median is the value that lies smack in the middle. Right in the middle of the minimum and the median is the first quartile, it splits the lowest half of the data again in two, so that one quarter of all values are smaller than :q1, and three quarters are greater than :q1. Similarly :q3 is the midway point between median and maximum. You could say it’s the median of the top half of the data.

Kixi.stats also provides a sixth number: the IQR or “interquartile range”. This is simply the difference between Q1 and Q3. It follows by definition that half of all data falls in this range. We’re starting to learn some things about our data now. While abalones can have as few as 1 rings, or as many as 29, at least half of them have somewhere between 8 and 11, a fairly narrow range.

(transduce (map :rings) kixi/summary abalone/data)
;; => {:min 1.0, :q1 8.0, :median 9.282608695652174, :q3 11.0, :max 29.0, :iqr 3.0}

What you could do when you first get a data set like this is to look at the summary for each of the variables. Concatenate all the maps so you just get a big sequence of key-value pairs, then use xform’s by-key to get the summary for each variable. You want to split the data based on the map keys, and then get the value from the key-value pairs, and reduce these values to get the summary.

Notice how x/reduce turns the reducing function kixi/summary into a transducer, which gives you a lot of possibilities.

So now you have this big pile of numbers. You could take a moment to go through these to see if anything looks strange, but you could also visualize this information and see if anything jumps out.

(into {}
      (comp cat (x/by-key key val (x/reduce kixi/summary)))
      abalone/data)
;; =>
;; {:rings          {:min 1.0, :q1 8.0, :median 9.266666666666667, :q3 11.0, :max 29.0, :iqr 3.0},
;;  :sex            {:min nil, :q1 nil, :median nil, :q3 nil, :max nil, :iqr nil},
;;  :shell-weight   {:min 0.3,:q1 26.01370731707317,:median 46.70705294705295,:q3 65.77664576802508,:max 201.0,:iqr 39.76293845095191},
;;  :diameter       {:min 11.0, :q1 70.0, :median 85.0, :q3 96.0, :max 130.0, :iqr 26.0},
;;  :whole-weight   {:min 0.4,:q1 88.39417977855477,:median 159.9038757763975,:q3 230.66032011466794,:max 565.1,:iqr 142.26614033611315},
;;  :viscera-weight {:min 0.1,:q1 18.644969969969964,:median 34.12173913043479,:q3 50.57916666666666,:max 152.0,:iqr 31.934196696696695},
;;  :length         {:min 15.0, :q1 90.0, :median 109.0, :q3 123.0, :max 163.0, :iqr 33.0},
;;  :shucked-weight {:min 0.2,:q1 37.21798245614035,:median 67.18149867374008,:q3 100.43694513371933,:max 297.6,:iqr 63.218962677578986},
;;  :height         {:min 0.0, :q1 23.0, :median 28.0, :q3 33.0, :max 226.0, :iqr 10.0}}

These numbers, useful as they are, don’t really provide a sense of shape or scale. Just like a picture is worth a thousand words, a graph is worth a thousand numbers.

In the git repository you’ll find a graphs namespace with various plotting functions, based on the thi.ng/geom library. All graphs you’ll see in this episode are created with these functions. There’s quite a bit there already, in the example-graphs namespace you can see it in action.

I decided to use geom for this episode because it’s a really nice and powerful API, and because contrary to most alternatives it’s a pure Clojure library. It doesn’t wrap some Java thing, which also means it’s equally usable from Clojure and ClojureScript.

But I found that for quick data exploration the thi.ng/geom API is too low-level, I needed something that I could just throw some data at and see a result, and so I started writing my own graphing functions on top of thi.ng/geom. This turned out to be quite a rabbit hole, and that’s one of the reasons this episode took a bit longer than usual.

The result you can see in this graphs namespace. This might evolve into a proper library at some point, for now it’s simply included here so you can follow along with the episode.

You use it in three steps: you start with a map containing the data and various attributes and options, you run this through a graph function like histogram, box-plot, scatter-plot, or barcode-plot, and then render and possibly open the result.

g/open will use whatever program your operating system thinks is most appropriate for displaying SVGs.

(g/open (g/histogram {:data abalone/data :x :rings}))

To use a specific program like your browser use open-with. Open and open-with both render to a temporary file. You can use render-to-file to save out the result. If you have an image viewer that automatically refreshes when the file updates, then this can be a great way to work from the REPL and get quick feedback. You could even use Emacs for that, although I’ve had Emacs crash a few times when viewing SVGs, so your mileage may vary.

(g/open-with
 (g/histogram
  {:data abalone/data :x :rings})
 "firefox")

(g/render-to-file
 (g/histogram
  {:data abalone/data :x :rings})
 "graph.svg")

So what are we looking at? A histogram splits the domain, the x-axis, in equally sized parts called “bins”, and counts how many data points fall within each bin. So on the x-axis you have the amount of rings, whereas the y-axis are counts: how many abalones have about this many rings?

The amount of rings is discrete, it can be 1 or 5 or 7, but not 2.95 or 6.3, so it would make sense to have one “bin” for each number of rings in the dataset.

The minimum and maximum amount of rings are 1 and 29, so that means 29 distinct values, so 29 would make a good bin number. Mind you, not 28! If you’re familiar with the fence post analogy, then here we’re counting the fence posts, not the fences in between them.

(transduce (map :rings) (redux/juxt kixi/min kixi/max) abalone/data)
;;=> [1.0 29.0]

So let’s draw this histogram again but with 29 bins, now you can see clearly how prominent 9 and 10 ring abalones are, whereas for lower or higher number, the counts very quickly taper off.

Now a first question to ask is what kind of a data distribution this is. You do this partly by looking at the shape of this graph, and partly by understanding the underlying process that generated this data. You recognize a kind of bell shape in this graph, which could indicate a normal distribution, but things aren’t spread out symetrically around the median, instead there’s a much longer tail on the high end.

(g/open (g/histogram {:data abalone/data
                      :x :rings
                      :bin-count 29}))

That long tail on the high end should be a hint though that maybe this is a log-normal distribution. This means that when you plot the logarithm of the amount of rings, the result should be roughly normal.

Log-normal distributions are quite common in nature, they occur when multiple factors compound and multiply over time.

There’s an log function in the kixi.math namespace to calculate the “natural” logarithm. Use that to transform the value of :rings before building up the histogram.

Now that long tail is gone, but there seems to be holes in this curve. Remember that the ring count is discrete, it’s not continuous, only a limited number of distinct values occur over and over. It’s possible for some bins to fall “in between” two of these values, so they end up empty.

(g/open (g/histogram {:data abalone/data
                      :x (comp m/log :rings)
                      :bin-count 29}))

A trick you can use is add in a little bit of randomness, a bit of noise. Move each point a tiny bit to the left or the right, so that instead of being concentrated around a small number of discrete values, it looks more like a continuous variable. This is called jitter.

Now it starts to look more like that Gaussian bell-shaped curve you’d expect from a normal distribution. This is an important insight, it means that we might have more success predicting the logarithm of the rings, then the amount of rings directly.

(g/open (g/histogram {:data abalone/data
                      :x (comp m/log (g/jitter 0.5) :rings)
                      :bin-count 29}))

So far we’ve focused on the :rings, which is what statisticians call the “dependent” or “outcome” variable, this is the thing we’re trying to predict. In the next episode we’ll have a closer look at the “independent” input variables.

How are their values distributed? Do we need to apply some cleanup or data transformation? Which of them is most suitable for performing the linear regression we are after, and how good is the final model really able to predict the amount of rings?

Thanks for watching!

Data Science with Kixi.stats, part 1