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
(:as ds]
[tech.v3.dataset :as tc]
[tablecloth.api :as tcc]
[tablecloth.column.api :as datetime]
[tech.v3.datatype.datetime :as ds-mod]
[tech.v3.dataset.modelling :as reg]
[fastmath.ml.regression :as kind]
[scicloj.kindly.v4.kind :as rand]
[fastmath.random :as plotly])) [scicloj.tableplot.v1.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)
(50
n 2
a 5]
b --> {:x (repeatedly n #(rand/frandom rng 0 10))}
(
tc/dataset:y
(tc/map-columns :x]
[fn [x]
(+ (* a x)
(
b (rand/grandom rng)))))))
simple-linear-data
_unnamed [50 2]:
:x | :y |
---|---|
9.98543358 | 15.33216867 |
0.94265521 | -3.23259698 |
2.10963964 | -0.23138979 |
0.07797003 | -4.86687826 |
5.61702967 | 6.89369896 |
2.12225914 | -0.36947672 |
0.53472459 | -2.63802288 |
9.60548687 | 14.11914213 |
6.92310095 | 10.24232095 |
0.66484571 | -4.05731139 |
… | … |
8.57214832 | 12.60793533 |
5.31104183 | 6.17687199 |
9.83097076 | 14.41649007 |
2.51579762 | -0.96940266 |
2.72616148 | 1.14498384 |
2.24961567 | -0.75149865 |
1.75070643 | -0.57715399 |
1.93591714 | -0.69747216 |
1.78568423 | -2.90318172 |
0.24327755 | -4.62719145 |
4.01331377 | 3.29549263 |
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:
:y)
(simple-linear-data ;; xss - a sequence of "rows", each containing `x` values:
;; (one `x` per row, in our case):
-> simple-linear-data
(:x])
(tc/select-columns [
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
0x33faf5e0 "Array2DRowRealMatrix{{0.0688868419,-0.0103696919},{-0.0103696919,0.0021995798}}"],
#object[org.apache.commons.math3.linear.Array2DRowRealMatrix :intercept -5.013691312721078,
:beta [2.0279785456833213],
:coefficients
:estimate -5.013691312721078,
[{:stderr 0.23017108956642138,
:t-value -21.782454617413006,
:p-value 1.5679059554058387E-26,
:confidence-interval [-5.476481305603564 -4.550901319838592]}
:estimate 2.0279785456833213,
{:stderr 0.04112941171364424,
:t-value 49.30725875202735,
:p-value 0.0,
:confidence-interval [1.945282320931231 2.1106747704354114]}],
: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.09561491406262412
[0.13059019839997577
-0.5039975853154866
0.011308490560706375
-0.5161746210728131
0.3403185880161257
1.2912584289391873
0.34688784942445494
-1.2161120704661776
0.3919128986879459
-0.7920056275160121
0.963733138905635
1.2331670135017454
-0.5585882329243337
1.7439676496032401
0.7003299579882221
0.48251564399468294
-0.6982067692995866
0.23172087897930638
0.3435287269719609
-0.31072606471314224
-2.2886379050289407
-0.37060561035305284
1.0284416640751104
-0.5473002704406538
-0.005626130663584528
0.011323227129404145
0.3258914770877763
-0.9431498075141876
1.97553149785786
-1.1840625525684816
-0.6296930043452711
-0.007943662251562955
-1.2852233160963165
1.1566753185716268
0.6098461718373587
-0.19710386338865238
-0.8107753475359667
1.4422347711696872
-0.23749375069142964
0.419884416517597
0.5068164126284813
-1.0576949329974001
-0.6300781586103774
0.2999796508909305
-0.8861422311729502
0.39022073268990254
1.5108197143284197
-0.10686179353496605
-0.17026972005166519],
:raw
0.09561491406262412
[0.13059019839997577
-0.5039975853154866
0.011308490560706375
-0.5161746210728131
0.3403185880161257
1.2912584289391873
0.34688784942445494
-1.2161120704661776
0.3919128986879459
-0.7920056275160121
0.963733138905635
1.2331670135017454
-0.5585882329243337
1.7439676496032401
0.7003299579882221
0.48251564399468294
-0.6982067692995866
0.23172087897930638
0.3435287269719609
-0.31072606471314224
-2.2886379050289407
-0.37060561035305284
1.0284416640751104
-0.5473002704406538
-0.005626130663584528
0.011323227129404145
0.3258914770877763
-0.9431498075141876
1.97553149785786
-1.1840625525684816
-0.6296930043452711
-0.007943662251562955
-1.2852233160963165
1.1566753185716268
0.6098461718373587
-0.19710386338865238
-0.8107753475359667
1.4422347711696872
-0.23749375069142964
0.419884416517597
0.5068164126284813
-1.0576949329974001
-0.6300781586103774
0.2999796508909305
-0.8861422311729502
0.39022073268990254
1.5108197143284197
-0.10686179353496605
-0.17026972005166519]},
:fitted
15.23655375380493
(3.1020067795307877
-0.7353873743030475
-4.8555697688857204
-6.377524342220142
0.7097953083910351
-3.9292813099382116
-14.466029980100522
9.026208880099603
3.665398486725773
-10.243808107148851
9.199005851566955
5.947772005019181
0.20548875561741387
9.441362183895173
10.409974253150164
0.8518102416150013
-5.439799715900937
9.467220179077007
13.369869529899809
5.381493514723352
0.0713849764167298
8.919131250836909
8.917903141118103
7.2460406343146415
1.3345054923661142
9.510273644353205
2.3878065834205127
-1.3703757023196097
-13.991254366995634
3.5293193587951714
0.8533732026969867
14.226324237330768
7.605567336956892
5.618936441160812
3.629547594845725
-2.0913450071162285
-6.127322612924449
14.651128356082012
12.370441576858
5.756987577741161
14.923306480235684
0.08829227588416355
0.5149056806857146
0.4515189994480613
-1.4632962242725776
-1.0876928885002601
-1.3923620090867361
-4.520329661194544
-3.1252229105309786),
:df {:residual 48, :model 1, :intercept 1},
:observations 50,
:names ["Intercept" "x"],
:r-squared 0.980638960805409,
:adjusted-r-squared 0.9802356058221884,
:sigma2 0.7690689396103255,
:sigma 0.8769657573761507,
:tss 1906.6801492560794,
:rss 36.91530910129563,
:regss 1869.7648401547838,
:msreg 1869.7648401547838,
:qt 2.010634757624228,
:f-statistic 2431.205765639375,
:p-value 0.0,
:ll
:log-likelihood -63.36201016955901,
{:aic 132.724020339118,
:bic 138.46008935540246,
:aic-rss -11.169832981349257,
:bic-rss -7.3457869704929655},
:analysis #<Delay@3f3539c7: :not-delivered>}
Printing the model gives a tabular summary: We’ll capture the printed output and display it via Kindly for cleaner formatting.
(kind/codewith-out-str
(println
( simple-linear-data-model)))
Residuals:
:min | :q1 | :median | :q3 | :max |
|
|-----------+-----------+----------+---------+----------|2.288638 | -0.488591 | 0.008475 | 0.64711 | 1.743968 |
| -
Coefficients:
:name | :estimate | :stderr | :t-value | :p-value | :confidence-interval |
|
|-----------+-----------+----------+------------+----------+-----------------------|5.013691 | 0.230171 | -21.782455 | 0.0 | [-5.476481 -4.550901] |
| Intercept | -2.027979 | 0.041129 | 49.307259 | 0.0 | [1.945282 2.110675] |
| x |
2431.205765639375 on degrees of freedom: {:residual 48, :model 1, :intercept 1}
F-statistic: 0.0
p-value:
0.980638960805409
R2: 0.9802356058221884
Adjusted R2: 0.8769657573761507 on 48 degrees of freedom
Residual standard error: 132.724020339118
AIC:
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]nil))
(lm dataset
([dataset options]let [inference-column-name (-> dataset
(
ds-mod/inference-target-column-namesfirst)
-> dataset
ds-without-target (
(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-namesvec)}
options)))))
defn summary
("Generate a summary of a linear model."
[lmdata]
(kind/codewith-out-str
(println
( lmdata))))
-> simple-linear-data
(:y)
(ds-mod/set-inference-target
lm summary)
Residuals:
:min | :q1 | :median | :q3 | :max |
|
|-----------+-----------+----------+---------+----------|2.288638 | -0.488591 | 0.008475 | 0.64711 | 1.743968 |
| -
Coefficients:
:name | :estimate | :stderr | :t-value | :p-value | :confidence-interval |
|
|-----------+-----------+----------+------------+----------+-----------------------|5.013691 | 0.230171 | -21.782455 | 0.0 | [-5.476481 -4.550901] |
| Intercept | -:x | 2.027979 | 0.041129 | 49.307259 | 0.0 | [1.945282 2.110675] |
|
2431.205765639375 on degrees of freedom: {:residual 48, :model 1, :intercept 1}
F-statistic: 0.0
p-value:
0.980638960805409
R2: 0.9802356058221884
Adjusted R2: 0.8769657573761507 on 48 degrees of freedom
Residual standard error: 132.724020339118
AIC:
10.2.3 Prediction
Once we have a linear model, we can generate new predictions. For instance, let’s predict \(y\) when \(x=3\):
3]) (simple-linear-data-model [
1.0702443243288853
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
("data"})
(plotly/layer-point {:=name "prediction"})) (plotly/layer-smooth {:=name
Alternatively, we can build the regression line explicitly. We’ll obtain predictions and then plot them:
-> simple-linear-data
(:prediction
(tc/map-columns :x]
[
simple-linear-data-model)"data"})
(plotly/layer-point {:=name :prediction
(plotly/layer-smooth {:=y "prediction"})) :=name
10.3 Multiple linear regression
We can easily extend these ideas to multiple linear predictors.
def multiple-linear-data
(let [rng (rand/rng 1234)
(50
n 2
a0 3
a1 -5]
b --> {:x0 (repeatedly n #(rand/frandom rng 0 10))
(:x1 (repeatedly n #(rand/frandom rng 0 10))}
tc/dataset:y
(tc/map-columns :x0 :x1]
[fn [x0 x1]
(+ (* a0 x0)
(* a1 x1)
(
b (rand/grandom rng)))))))
def multiple-linear-data-model
(-> multiple-linear-data
(:y)
(ds-mod/set-inference-target lm))
(summary multiple-linear-data-model)
Residuals:
:min | :q1 | :median | :q3 | :max |
|
|-----------+-----------+-----------+----------+----------|1.928904 | -0.761525 | -0.101411 | 0.485421 | 3.106313 |
| -
Coefficients:
:name | :estimate | :stderr | :t-value | :p-value | :confidence-interval |
|
|-----------+-----------+----------+------------+----------+-----------------------|4.499123 | 0.490708 | -9.168643 | 0.0 | [-5.486299 -3.511947] |
| Intercept | -:x0 | 1.982262 | 0.061884 | 32.031664 | 0.0 | [1.857767 2.106758] |
| :x1 | -3.051798 | 0.055615 | -54.873541 | 0.0 | [-3.163681 -2.939915] |
|
2047.7721187786567 on degrees of freedom: {:residual 47, :model 2, :intercept 1}
F-statistic: 0.0
p-value:
0.9886543154871138
R2: 0.9881715204014591
Adjusted R2: 1.1345224500096749 on 47 degrees of freedom
Residual standard error: 159.42126448491905
AIC:
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
(:3d
(plotly/layer-point {:=coordinates :x0
:=x :x1
:=y :y})
:=z let [xs (range 11)
(plotly/layer-surface {:=dataset (range 11)]
ys (
(tc/dataset:x xs
{:y ys
:z (for [y ys]
for [x xs]
(
(multiple-linear-data-model
[x y])))}))0.5})) :=mark-opacity
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
("data/seattle-bikes-and-weather/Fremont_Bridge_Bicycle_Counter.csv.gz"
(tc/dataset :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
("data/seattle-bikes-and-weather/BicycleWeather.csv.gz"
(tc/dataset :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:
= counts.resample('d').sum() daily
Since Tablecloth’s time series features are still evolving, we’ll be a bit more explicit:
def daily-totals
(-> counts
(fn [{:keys [datetime]}]
(tc/group-by (:date (datetime/local-date-time->local-date
{
datetime)})):total :west :east]
(tc/aggregate-columns [ 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,
1) (idx->day-of-week
:Mon
7) (idx->day-of-week
:Sun
Now, let’s build our dataset:
def totals-with-day-of-week
(-> daily-totals
(:day-of-week
(tc/add-column fn [ds]
(map idx->day-of-week
(
(datetime/long-temporal-field:day-of-week
:date ds)))))
(:total :day-of-week]))) (tc/select-columns [
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
1)))))
(tcc/*
totals-with-day-of-week
days-of-week):day-of-week])
(tc/drop-columns [:total))) (ds-mod/set-inference-target
-> 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] |
|
3472.719277560978 on degrees of freedom: {:residual 4436, :model 7, :intercept 0}
F-statistic: 0.0
p-value:
0.8456776967289252
R2: 0.8454341764126724
Adjusted R2: 1121.0266948292917 on 4436 degrees of freedom
Residual standard error: 75015.17638480052
AIC:
We can clearly see weekend versus weekday differences.