10  Introduction to Linear Regression

last update: 2024-12-29

In this tutorial, we introduce the fundamentals of linear regression, guided by the In Depth: Linear Regression chapter of the Python Data Science Handbook by Jake VanderPlas.

10.1 Setup

(ns noj-book.linear-regression-intro
  (:require
   [tech.v3.dataset :as ds]
   [tablecloth.api :as tc]
   [tablecloth.column.api :as tcc]
   [tech.v3.datatype.datetime :as datetime]
   [tech.v3.dataset.modelling :as ds-mod]
   [fastmath.ml.regression :as reg]
   [scicloj.kindly.v4.kind :as kind]
   [fastmath.random :as rand]
   [scicloj.tableplot.v1.plotly :as plotly]))

10.2 Simple Linear Regression

We begin with the classic straight-line model: for data points \((x, y)\), we assume there is a linear relationship allowing us to predict \(y\) as \[y = ax + b.\] In this formulation, \(a\) is the slope and \(b\) is the intercept, the point where our line would cross the \(y\) axis.

To illustrate, we’ll use Fastmath and Tablecloth to create synthetic data in which the relationship is known to hold with \(a=2\) and \(b=-5\).

For each row in the dataset below, we draw \(x\) uniformly from 0 to 10 and compute \(y = ax + b\) plus an extra random noise term (drawn from a standard Gaussian distribution). This noise is added independently for every row.

(def simple-linear-data
  (let [rng (rand/rng 1234)
        n 50
        a 2
        b -5]
    (-> {:x (repeatedly n #(rand/frandom rng 0 10))}
        tc/dataset
        (tc/map-columns :y
                        [:x]
                        (fn [x]
                          (+ (* a x)
                             b
                             (rand/grandom rng)))))))

simple-linear-data

_unnamed [50 2]:

:x	:y
8.43507195	11.67923879
0.81211567	-3.17546710
9.96872807	14.31359196
6.58413935	9.77055203
8.72384167	11.93271843
8.82741547	12.21489527
9.67620945	12.87133681
2.65301585	0.74672126
2.32064128	-1.47945017
9.86531067	14.31804683
…	…
5.96804714	6.72500358
4.19868422	3.37834290
9.68653297	15.13284509
9.99031448	15.62983675
2.79830503	0.82522248
5.90288877	7.60069629
4.31125832	4.07367578
1.04621530	-2.99750308
7.73808432	10.75323893
9.59535217	14.61437006
3.01903963	-0.12733051

Let’s plot these points using Tableplot’s Plotly API.

(-> simple-linear-data
    plotly/layer-point)

10.2.1 Regression using Fastmath

We can now fit a linear model to the data using the Fastmath library.

(def simple-linear-data-model
  (reg/lm
   ;; ys - a "column" sequence of `y` values:
   (simple-linear-data :y)
   ;; xss - a sequence of "rows", each containing `x` values:
   ;; (one `x` per row, in our case):
   (-> simple-linear-data
       (tc/select-columns [:x])
       tc/rows)
   ;; options
   {:names ["x"]}))

(type simple-linear-data-model)

fastmath.ml.regression.LMData

simple-linear-data-model

{:model :ols,
 :intercept? true,
 :offset? false,
 :transformer nil,
 :xtxinv
 #object[org.apache.commons.math3.linear.BlockRealMatrix 0x18151ca5 "BlockRealMatrix{{0.0638723276,-0.0095397869},{-0.0095397869,0.0020743721}}"],
 :intercept -4.964537061157466,
 :beta [2.000289714793488],
 :coefficients
 [{:estimate -4.964537061157466,
   :stderr 0.21955910709391505,
   :t-value -22.611392106973344,
   :p-value 3.051477665988989E-27,
   :confidence-interval [-5.405990233233432 -4.5230838890815]}
  {:estimate 2.000289714793488,
   :stderr 0.039567482629912616,
   :t-value 50.553878635718135,
   :p-value 0.0,
   :confidence-interval [1.9207339589460928 2.0798454706408833]}],
 :offset
 (0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0),
 :weights
 (1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0
  1.0),
 :residuals
 {:weighted
  (-0.2288118051025041
   0.1646033370538591
   -0.6622152012023861
   1.5649028724865257
   -0.552955272824196
   -0.47795603930751085
   -1.519348370934047
   0.4044580010172647
   -1.1568679865547737
   -0.4508955749140302
   1.4375408809486063
   -1.711874286904794
   -0.9979149288597808
   0.3171121984869545
   0.4601400770850823
   -0.4342369380532425
   -0.01169480717427529
   0.4178230074820686
   0.8018558976196826
   1.1925890310562846
   -0.4838746090352446
   -0.09677068144628986
   0.5463642777595328
   0.586452744773883
   0.3062822283496889
   0.2512066842279186
   1.2104310042922668
   0.45191869576873545
   -0.016456330861003954
   -0.36356807598719776
   -0.8735918995507279
   -1.1517151634895462
   1.3798285110698778
   -2.0877086514697996
   -1.0954086114728767
   0.12509653583696978
   1.3934709984611722
   0.8757770089612711
   -1.2042344490388253
   -0.24828267020289374
   -0.05570488775390148
   0.7215098717731898
   0.6108505053678357
   0.19233876499831148
   0.7577456441210488
   0.4144471706823838
   -0.12569971417129988
   0.2393655179954628
   0.38542286299824546
   -1.2017473743624023),
  :raw
  (-0.2288118051025041
   0.1646033370538591
   -0.6622152012023861
   1.5649028724865257
   -0.552955272824196
   -0.47795603930751085
   -1.519348370934047
   0.4044580010172647
   -1.1568679865547737
   -0.4508955749140302
   1.4375408809486063
   -1.711874286904794
   -0.9979149288597808
   0.3171121984869545
   0.4601400770850823
   -0.4342369380532425
   -0.01169480717427529
   0.4178230074820686
   0.8018558976196826
   1.1925890310562846
   -0.4838746090352446
   -0.09677068144628986
   0.5463642777595328
   0.586452744773883
   0.3062822283496889
   0.2512066842279186
   1.2104310042922668
   0.45191869576873545
   -0.016456330861003954
   -0.36356807598719776
   -0.8735918995507279
   -1.1517151634895462
   1.3798285110698778
   -2.0877086514697996
   -1.0954086114728767
   0.12509653583696978
   1.3934709984611722
   0.8757770089612711
   -1.2042344490388253
   -0.24828267020289374
   -0.05570488775390148
   0.7215098717731898
   0.6108505053678357
   0.19233876499831148
   0.7577456441210488
   0.4144471706823838
   -0.12569971417129988
   0.2393655179954628
   0.38542286299824546
   -1.2017473743624023),
  :loocv
  (-0.24098827709882234
   0.17322028788507254
   -0.719654492885639
   1.6102733473457083
   -0.5853212531086679
   -0.5068951123256392
   -1.6398367279797565
   0.4160467202382848
   -1.1935909053293294
   -0.48879335349640973
   1.5267572783446979
   -1.775274851128044
   -1.0188672629870412
   0.33109333920315603
   0.47958251902001103
   -0.44318631244669654
   -0.011949307663658936
   0.4265802707610643
   0.8290837215404698
   1.2170073405298063
   -0.49401541579191066
   -0.09911615420498436
   0.5751462317287435
   0.6255717544279965
   0.3135421103542509
   0.2563872866643102
   1.236594207069957
   0.4691212723227117
   -0.017030515037072876
   -0.3749784920507121
   -0.9014093326271503
   -1.191862852095386
   1.465419086442016
   -2.130699561858337
   -1.1200377009272628
   0.13038166332869497
   1.4881612853945143
   0.9271593964994245
   -1.278772593573741
   -0.25435896536690994
   -0.05686099826195035
   0.7789105194184449
   0.6641823318738636
   0.19762021548275813
   0.7760029347286751
   0.42297934249649477
   -0.13178578241256914
   0.24945396671597378
   0.41523068148711145
   -1.2327857273582778)},
 :fitted
 (11.908050594350051
  -3.3400704407330837
  14.97580715781681
  8.2056491555666
  12.485673699180147
  12.692851304308576
  14.390685179401153
  0.3422632607311842
  -0.32258217860733485
  14.768942403176307
  -4.375666361308374
  -1.2707397261512923
  5.277880785884123
  -2.3131923660713047
  -2.0598335460066135
  3.624011942469598
  5.817057617948156
  3.2244303334550155
  -0.7422057720045938
  4.586521267693986
  5.243067584790331
  1.5761546463246079
  -3.377819953494791
  -4.823056866227646
  1.7678981350800065
  4.86514991825207
  2.740369554576131
  -1.435857778449996
  -0.908816779261298
  8.719753119410889
  -0.34225932972284223
  -0.9031511332299784
  12.841600265567212
  4.818687805545399
  6.193504505668185
  -2.0591396474396464
  13.408103824701499
  -4.0309370763328705
  -4.359252441082539
  6.973286254445431
  3.4340477908642706
  14.411335219229109
  15.0189862480252
  0.6328837184826632
  6.842950642788581
  3.6592286062353176
  -2.8718033655258495
  10.513873408849534
  14.228947200018862
  1.0744168610655143),
 :df {:residual 48, :model 1, :intercept 1},
 :observations 50,
 :names ["Intercept" "x"],
 :cv 0.8873651015577338,
 :r-squared 0.9815646584786514,
 :adjusted-r-squared 0.9811805888636232,
 :sigma2 0.7547274907167489,
 :sigma 0.8687505342253028,
 :tss 1965.0799261002169,
 :rss 36.22691955440395,
 :regss 1928.8530065458128,
 :msreg 1928.8530065458128,
 :qt 2.010634757624228,
 :f-statistic 2555.694645114916,
 :p-value 0.0,
 :ll
 {:log-likelihood -62.89141343884106,
  :aic 131.7828268776821,
  :bic 137.51889589396654,
  :aic-rss -12.111026442785157,
  :bic-rss -8.286980431928871},
 :analysis #<Delay@181faa4b: :not-delivered>,
 :decomposition :cholesky,
 :augmentation nil,
 :hat
 (0.05052723785117988
  0.049745621234220404
  0.07981509495332315
  0.028175635480751787
  0.05529609613964076
  0.057090850383929775
  0.07347582536107028
  0.027854369851498452
  0.030766754849244807
  0.07753333450893146
  0.058435219966869636
  0.035713097711582
  0.020564341291949992
  0.042227188109099395
  0.040540347414367385
  0.02019325539195308
  0.021298346025322594
  0.020528992734173586
  0.0328408618012629
  0.020064225301132232
  0.020527308323790322
  0.0236638797934368
  0.050042845421588496
  0.06253320962338325
  0.023154408179365687
  0.020206159610303136
  0.02115746833367623
  0.036669785765209406
  0.03371502123212406
  0.030429521440315487
  0.030859934626313132
  0.033684822490488
  0.05840689272032688
  0.020176899248546237
  0.02198951824031991
  0.0405358188934993
  0.06362904872118051
  0.05541915201652734
  0.058288819223601365
  0.023888661267556323
  0.020332223200212614
  0.07369350678189833
  0.08029696658079152
  0.026725254152490552
  0.02352734737273872
  0.02017160404040698
  0.0461815237566082
  0.040442125869249425
  0.07178616566124629
  0.025177411051300294),
 :effective-dimension 1.9999999999999978}

Printing the model gives a tabular summary: We’ll capture the printed output and display it via Kindly for cleaner formatting.

(kind/code
 (with-out-str
   (println
    simple-linear-data-model)))

Residuals:

|      :min |       :q1 | :median |      :q3 |     :max |
|-----------+-----------+---------+----------+----------|
| -2.087709 | -0.501145 | 0.14485 | 0.556386 | 1.564903 |

Coefficients:

|     :name | :estimate |  :stderr |   :t-value | :p-value | :confidence-interval |
|-----------+-----------+----------+------------+----------+----------------------|
| Intercept | -4.964537 | 0.219559 | -22.611392 |      0.0 | [-5.40599 -4.523084] |
|         x |   2.00029 | 0.039567 |  50.553879 |      0.0 |  [1.920734 2.079845] |

F-statistic: 2555.694645114916 on degrees of freedom: {:residual 48, :model 1, :intercept 1}
p-value: 0.0

R2: 0.9815646584786514
Adjusted R2: 0.9811805888636232
Residual standard error: 0.8687505342253028 on 48 degrees of freedom
AIC: 131.7828268776821

As you can see, the estimated coefficients match our intercept \(b\) and slope \(a\) (the coefficient of \(x\)).

10.2.2 Dataset ergonomics

Below are a couple of helper functions that simplify how we use regression with datasets and display model summaries. We have similar ideas under development in the Tablemath library, but it is still in an experimental stage and not part of Noj yet.

(defn lm
  "Compute a linear regression model for `dataset`.
  The first column marked as target is the target.
  All the columns unmarked as target are the features.
  The resulting model is of type `fastmath.ml.regression.LMData`,
  created via [Fastmath](https://github.com/generateme/fastmath).
  
  See [fastmath.ml.regression.lm](https://generateme.github.io/fastmath/clay/ml.html#lm)
  for `options`."
  ([dataset]
   (lm dataset nil))
  ([dataset options]
   (let [inference-column-name (-> dataset
                                   ds-mod/inference-target-column-names
                                   first)
         ds-without-target (-> dataset
                               (tc/drop-columns [inference-column-name]))]
     (reg/lm
      ;; ys
      (get dataset inference-column-name)
      ;; xss
      (tc/rows ds-without-target)
      ;; options
      (merge {:names (-> ds-without-target
                         tc/column-names
                         vec)}
             options)))))

(defn summary
  "Generate a summary of a linear model."
  [lmdata]
  (kind/code
   (with-out-str
     (println
      lmdata))))

(-> simple-linear-data
    (ds-mod/set-inference-target :y)
    lm
    summary)

Residuals:

|      :min |       :q1 | :median |      :q3 |     :max |
|-----------+-----------+---------+----------+----------|
| -2.087709 | -0.501145 | 0.14485 | 0.556386 | 1.564903 |

Coefficients:

|     :name | :estimate |  :stderr |   :t-value | :p-value | :confidence-interval |
|-----------+-----------+----------+------------+----------+----------------------|
| Intercept | -4.964537 | 0.219559 | -22.611392 |      0.0 | [-5.40599 -4.523084] |
|        :x |   2.00029 | 0.039567 |  50.553879 |      0.0 |  [1.920734 2.079845] |

F-statistic: 2555.694645114916 on degrees of freedom: {:residual 48, :model 1, :intercept 1}
p-value: 0.0

R2: 0.9815646584786514
Adjusted R2: 0.9811805888636232
Residual standard error: 0.8687505342253028 on 48 degrees of freedom
AIC: 131.7828268776821

10.2.3 Prediction

Once we have a linear model, we can generate new predictions. For instance, let’s predict \(y\) when \(x=3\):

(simple-linear-data-model [3])

1.036332083222998

10.2.4 Displaying the regression line

We can visualize the fitted line by adding a smooth layer to our scatter plot. Tableplot makes this convenient:

(-> simple-linear-data
    (plotly/layer-point {:=name "data"})
    (plotly/layer-smooth {:=name "prediction"}))

Alternatively, we can build the regression line explicitly. We’ll obtain predictions and then plot them:

(-> simple-linear-data
    (tc/map-columns :prediction
                    [:x]
                    simple-linear-data-model)
    (plotly/layer-point {:=name "data"})
    (plotly/layer-smooth {:=y :prediction
                          :=name "prediction"}))

10.3 Multiple linear regression

We can easily extend these ideas to multiple linear predictors.

(def multiple-linear-data
  (let [rng (rand/rng 1234)
        n 50
        a0 2
        a1 -3
        b -5]
    (-> {:x0 (repeatedly n #(rand/frandom rng 0 10))
         :x1 (repeatedly n #(rand/frandom rng 0 10))}
        tc/dataset
        (tc/map-columns :y
                        [:x0 :x1]
                        (fn [x0 x1]
                          (+ (* a0 x0)
                             (* a1 x1)
                             b
                             (rand/grandom rng)))))))

(def multiple-linear-data-model
  (-> multiple-linear-data
      (ds-mod/set-inference-target :y)
      lm))

(summary multiple-linear-data-model)

Residuals:

|     :min |       :q1 |   :median |      :q3 |     :max |
|----------+-----------+-----------+----------+----------|
| -1.61903 | -0.794329 | -0.090853 | 0.432747 | 2.469702 |

Coefficients:

|     :name | :estimate |  :stderr |   :t-value | :p-value |  :confidence-interval |
|-----------+-----------+----------+------------+----------+-----------------------|
| Intercept | -4.817659 | 0.419325 | -11.489091 |      0.0 | [-5.661231 -3.974086] |
|       :x0 |  2.021816 | 0.057178 |  35.359944 |      0.0 |   [1.906789 2.136844] |
|       :x1 | -3.042843 | 0.050406 |  -60.36609 |      0.0 | [-3.144248 -2.941438] |

F-statistic: 2406.7857030885066 on degrees of freedom: {:residual 47, :model 2, :intercept 1}
p-value: 0.0

R2: 0.9903303549989472
Adjusted R2: 0.9899188807435833
Residual standard error: 1.0051695978293465 on 47 degrees of freedom
AIC: 147.31571126783544

Visualizing multiple dimensions is more involved. In the case of two features, we can use a 3D scatterplot and a 3D surface. Let us do that using Tableplot’s Plotly API.

(-> multiple-linear-data
    (plotly/layer-point {:=coordinates :3d
                         :=x :x0
                         :=y :x1
                         :=z :y})
    (plotly/layer-surface {:=dataset (let [xs (range 11)
                                           ys (range 11)]
                                       (tc/dataset
                                        {:x xs
                                         :y ys
                                         :z (for [y ys]
                                              (for [x xs]
                                                (multiple-linear-data-model
                                                 [x y])))}))
                           :=mark-opacity 0.5}))

10.4 Coming soon: Polynomial regression 🛠

10.5 Coming soon: One-hot encoding 🛠

10.6 Coming soon: Regularization 🛠

10.7 Example: Predicting Bicycle Traffic

As in the Python Data Science Handbook, we’ll try predicting the daily number of bicycle trips across the Fremont Bridge in Seattle. The features will include weather, season, day of week, and related factors.

10.7.1 Reading and parsing data

(def column-name-mapping
  {"Fremont Bridge Sidewalks, south of N 34th St" :total
   "Fremont Bridge Sidewalks, south of N 34th St Cyclist West Sidewalk" :west
   "Fremont Bridge Sidewalks, south of N 34th St Cyclist East Sidewalk" :east
   "Date" :datetime})

(column-name-mapping
 "Fremont Bridge Sidewalks, south of N 34th St")

:total

(def counts
  (tc/dataset "data/seattle-bikes-and-weather/Fremont_Bridge_Bicycle_Counter.csv.gz"
              {:key-fn column-name-mapping
               :parser-fn {"Date" [:local-date-time "MM/dd/yyyy hh:mm:ss a"]}}))

counts

data/seattle-bikes-and-weather/Fremont_Bridge_Bicycle_Counter.csv.gz [106608 4]:

:datetime	:total	:west	:east
2012-10-02T13:00	55	7	48
2012-10-02T14:00	130	55	75
2012-10-02T15:00	152	81	71
2012-10-02T16:00	278	167	111
2012-10-02T17:00	563	393	170
2012-10-02T18:00	381	236	145
2012-10-02T19:00	175	104	71
2012-10-02T20:00	86	51	35
2012-10-02T21:00	63	35	28
2012-10-02T22:00	42	27	15
…	…	…	…
2024-11-30T13:00	147	62	85
2024-11-30T14:00	154	73	81
2024-11-30T15:00	118	57	61
2024-11-30T16:00	88	51	37
2024-11-30T17:00	46	11	35
2024-11-30T18:00	46	18	28
2024-11-30T19:00	69	14	55
2024-11-30T20:00	18	8	10
2024-11-30T21:00	49	15	34
2024-11-30T22:00	14	4	10
2024-11-30T23:00	10	5	5

(def weather
  (tc/dataset "data/seattle-bikes-and-weather/BicycleWeather.csv.gz"
              {:key-fn keyword}))

weather

data/seattle-bikes-and-weather/BicycleWeather.csv.gz [1340 26]:

:STATION	:STATION_NAME	:DATE	:PRCP	:SNWD	:SNOW	:TMAX	:TMIN	:AWND	:WDF2	:WDF5	:WSF2	:WSF5	:FMTM	:WT14	:WT01	:WT17	:WT05	:WT02	:WT22	:WT04	:WT13	:WT16	:WT08	:WT18	:WT03
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120101	0	0	0	128	50	47	100	90	89	112	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120102	109	0	0	106	28	45	180	200	130	179	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	1	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120103	8	0	0	117	72	23	180	170	54	67	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120104	203	0	0	122	56	47	180	190	107	148	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	1	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120105	13	0	0	89	28	61	200	220	107	165	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120106	25	0	0	44	22	22	180	180	45	63	-9999	1	1	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120107	0	0	0	72	28	23	170	180	54	63	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	1	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120108	0	0	0	100	28	20	160	200	45	63	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120109	43	0	0	94	50	34	200	200	67	89	-9999	1	1	-9999	-9999	-9999	-9999	-9999	1	1	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20120110	10	0	0	61	6	34	20	30	89	107	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999	-9999
…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150822	0	0	0	267	122	25	20	20	63	76	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150823	0	0	0	278	139	18	10	10	67	81	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	1	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150824	0	0	0	239	122	23	190	190	54	67	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150825	0	0	0	256	122	34	350	360	63	76	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150826	0	0	0	283	139	17	30	40	58	67	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150827	0	0	0	294	144	21	230	200	45	63	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150828	5	0	0	233	156	26	230	240	81	103	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150829	325	0	0	222	133	58	210	210	157	206	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150830	102	0	0	200	128	47	200	200	89	112	-9999	-9999	1	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150831	0	0	0	189	161	58	210	210	112	134	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999
GHCND:USW00024233	SEATTLE TACOMA INTERNATIONAL AIRPORT WA US	20150901	58	0	0	194	139	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999	-9999

10.7.2 Preprocessing

The bike counts come in hourly data, but our weather information is daily. We’ll need to aggregate the hourly counts into daily totals before combining the datasets.

In the Python handbook, one does:

daily = counts.resample('d').sum()

Since Tablecloth’s time series features are still evolving, we’ll be a bit more explicit:

(def daily-totals
  (-> counts
      (tc/group-by (fn [{:keys [datetime]}]
                     {:date (datetime/local-date-time->local-date
                             datetime)}))
      (tc/aggregate-columns [:total :west :east]
                            tcc/sum)))

daily-totals

_unnamed [4443 4]:

:date	:total	:west	:east
2012-10-02	1938.0	1165.0	773.0
2012-10-03	3521.0	1761.0	1760.0
2012-10-04	3475.0	1767.0	1708.0
2012-10-05	3148.0	1590.0	1558.0
2012-10-06	2006.0	926.0	1080.0
2012-10-07	2142.0	951.0	1191.0
2012-10-08	3537.0	1708.0	1829.0
2012-10-09	3501.0	1742.0	1759.0
2012-10-10	3235.0	1587.0	1648.0
2012-10-11	3047.0	1468.0	1579.0
…	…	…	…
2024-11-20	2300.0	787.0	1513.0
2024-11-21	2382.0	775.0	1607.0
2024-11-22	1473.0	536.0	937.0
2024-11-23	1453.0	652.0	801.0
2024-11-24	727.0	311.0	416.0
2024-11-25	1483.0	486.0	997.0
2024-11-26	2173.0	727.0	1446.0
2024-11-27	1522.0	548.0	974.0
2024-11-28	631.0	261.0	370.0
2024-11-29	833.0	368.0	465.0
2024-11-30	1178.0	509.0	669.0

10.7.3 Prediction by day-of-week

Next, we’ll explore a simple regression by day of week.

(def days-of-week
  [:Mon :Tue :Wed :Thu :Fri :Sat :Sun])

We’ll convert the numeric day-of-week to the corresponding keyword:

(def idx->day-of-week
  (comp days-of-week dec))

For example,

(idx->day-of-week 1)

:Mon

(idx->day-of-week 7)

:Sun

Now, let’s build our dataset:

(def totals-with-day-of-week
  (-> daily-totals
      (tc/add-column :day-of-week
                     (fn [ds]
                       (map idx->day-of-week
                            (datetime/long-temporal-field
                             :day-of-week
                             (:date ds)))))
      (tc/select-columns [:total :day-of-week])))

totals-with-day-of-week

_unnamed [4443 2]:

:total	:day-of-week
1938.0	:Tue
3521.0	:Wed
3475.0	:Thu
3148.0	:Fri
2006.0	:Sat
2142.0	:Sun
3537.0	:Mon
3501.0	:Tue
3235.0	:Wed
3047.0	:Thu
…	…
2300.0	:Wed
2382.0	:Thu
1473.0	:Fri
1453.0	:Sat
727.0	:Sun
1483.0	:Mon
2173.0	:Tue
1522.0	:Wed
631.0	:Thu
833.0	:Fri
1178.0	:Sat

(def totals-with-one-hot-days-of-week
  (-> (reduce (fn [dataset day-of-week]
                (-> dataset
                    (tc/add-column day-of-week
                                   #(-> (:day-of-week %)
                                        (tcc/eq day-of-week)
                                        ;; convert booleans to 0/1
                                        (tcc/* 1)))))
              totals-with-day-of-week
              days-of-week)
      (tc/drop-columns [:day-of-week])
      (ds-mod/set-inference-target :total)))

(-> totals-with-one-hot-days-of-week
    (tc/select-columns ds-mod/inference-column?))

_unnamed [0 0]

Since the binary columns sum to 1, they’re collinear, and we won’t use an intercept. This way, each coefficient directly reflects the expected bike count for that day of week.

(def days-of-week-model
  (lm totals-with-one-hot-days-of-week
      {:intercept? false}))

Let’s take a look at the results:

(-> days-of-week-model
    println
    with-out-str
    kind/code)

Residuals:

|         :min |         :q1 |    :median |        :q3 |        :max |
|--------------+-------------+------------+------------+-------------|
| -3034.944882 | -841.755906 | -90.179528 | 828.820472 | 3745.244094 |

Coefficients:

| :name |   :estimate |   :stderr |  :t-value | :p-value |      :confidence-interval |
|-------+-------------+-----------+-----------+----------+---------------------------|
|  :Mon | 2785.741325 | 44.521658 | 62.570476 |      0.0 | [2698.456664 2873.025986] |
|  :Tue | 3115.527559 | 44.486587 | 70.032964 |      0.0 | [3028.311653 3202.743465] |
|  :Wed | 3109.944882 | 44.486587 | 69.907472 |      0.0 | [3022.728976 3197.160788] |
|  :Thu | 2961.179528 | 44.486587 | 66.563423 |      0.0 | [2873.963622 3048.395433] |
|  :Fri | 2623.755906 | 44.486587 | 58.978583 |      0.0 |     [2536.54 2710.971811] |
|  :Sat | 1688.029921 | 44.486587 | 37.944693 |      0.0 | [1600.814015 1775.245827] |
|  :Sun | 1576.178233 | 44.521658 | 35.402506 |      0.0 | [1488.893572 1663.462894] |

F-statistic: 3472.719277560978 on degrees of freedom: {:residual 4436, :model 7, :intercept 0}
p-value: 0.0

R2: 0.8456776967289252
Adjusted R2: 0.8454341764126724
Residual standard error: 1121.0266948292917 on 4436 degrees of freedom
AIC: 75015.17638480052

We can clearly see weekend versus weekday differences.

10.7.4 Coming soon: more predictors for the bike counts 🛠

source: notebooks/noj_book/linear_regression_intro.clj