# California Housing: From Exploration to Linear Regression
This notebook walks the full arc on the classic California housing dataset: we
load it, take stock of its structure and distributions, engineer a handful of
features that track house prices better than the raw columns, confirm those
features correlate with the target, and finally train a linear-regression model
with **Hasktorch** to predict the median house value of a district.
The exploration follows chapter 2 of *Hands-on Machine Learning with
Scikit-Learn, Keras & TensorFlow* (the prose is heavily derived from that book
and included here purely for didactic purposes). The modeling section is a port
of the [`dataframe` CaliforniaHousing example](https://github.com/mchav/dataframe/blob/main/examples/CaliforniaHousing.hs).
## Setup
This notebook depends on `dataframe` for the data work, its plotting and
Hasktorch bridges for visualisation and tensors, `granite` for the grammar of
graphics, and `Torch` for the model. We declare the packages here once; each
import is introduced where it is first used, so every cell carries its own
context.
```haskell
-- cabal: build-depends: dataframe ==2.1.0.3, dataframe-hasktorch ==0.2.0.1, hasktorch, text, random, granite
-- cabal: default-extensions: BangPatterns, NumericUnderscores, OverloadedStrings, ScopedTypeVariables, TemplateHaskell, TypeApplications
```
## Loading the data
The dataset lives alongside the notebook. We read it into a dataframe with
`readCsv`.
```haskell
import qualified DataFrame as D
df <- D.readCsv "./examples/data/housing.csv"
```
## A first look at the data
Let's take a look at the top five rows using the dataframe's `take` function.
```haskell
import DataFrame.Operators
df |> D.take 5
|> D.toMarkdown'
|> displayMarkdown
```
>
> | longitude
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
## Typed column references
Before preparing the data for a model, it's worth trying out various attribute
combinations. The total number of rooms in a district is not very useful if you
don't know how many households there are; what you really want is the number of
rooms *per household*. Similarly the total number of bedrooms is more telling
relative to the number of rooms, and the population per household is an
interesting combination too.
To reference existing attributes in type-safe expressions, we generate global
references to each column with `declareColumns`. This reads the loaded dataframe
and creates one typed handle per column (`total_rooms`, `households`,
`median_house_value`, and so on), so the compiler catches typos and type
mismatches in the expressions that follow.
```haskell
$(D.declareColumns df)
```
## Spatial structure
Let's bucket districts by price and plot them by location. We use `F.name` to
turn a typed column reference back into its text name, which keeps us from
mistyping column names in the plotting call.
```haskell
import qualified DataFrame.Functions as F
priceBucket p
| p > 500000 = (">500'000" :: Text)
| p > 450000 = ">450'000"
| p > 400000 = ">400'000"
| p > 300000 = ">300'000"
| p > 200000 = ">200'000"
| otherwise = "<=200'000"
withPriceBuckets = df
|> D.derive "price_bucket" (F.lift priceBucket median_house_value)
|> D.sortBy [D.Desc median_house_value]
scatter <- Plt.scatter ((Plt.mkScatter (F.name longitude) (F.name latitude)){ Plt.color = Just "price_bucket" }) withPriceBuckets
displayHtml scatter
```
>
>
>
>
>
>
### Rebuilding it with granite
The scatter above is drawn with Chart.js. Below we build the same
longitude/latitude view with [granite](https://github.com/mchav/granite)'s
grammar of graphics, one step at a time: prepare the data, then layer on the
cities and the price emphasis.
```haskell
import Granite.Data.Frame
import Granite.Render.Pipeline
import Granite.Spec
housingFrame =
fromColumns
[ ("lng", ColNum (D.columnAsList longitude withPriceBuckets))
, ("lat", ColNum (D.columnAsList latitude withPriceBuckets))
, ("bucket", ColCat (D.columnAsList (F.col @Text "price_bucket") withPriceBuckets))
, ("value", ColNum (D.columnAsList median_house_value withPriceBuckets))
]
(frameLength housingFrame, columnNames housingFrame)
```
>
> (20640,["lng","lat","bucket","value"])
That is a granite frame of all 20,640 districts with their coordinates and
price bucket. Next, the major population centres we want to label:
```haskell
cityName = ["Los Angeles", "San Francisco", "San Diego", "San Jose", "Sacramento", "Fresno", "Long Beach", "Bakersfield"] :: [Text]
cityLng = [-118.24, -122.42, -117.16, -121.89, -121.49, -119.77, -118.19, -119.02] :: [Double]
cityLat = [34.05, 37.77, 32.72, 37.34, 38.58, 36.75, 33.77, 35.37] :: [Double]
citiesFrame =
fromColumns
[ ("lng", ColNum cityLng)
, ("lat", ColNum cityLat)
, ("name", ColCat cityName)
]
zip3 cityName cityLng cityLat
```
>
> [("Los Angeles",-118.24,34.05),("San Francisco",-122.42,37.77),("San Diego",-117.16,32.72),("San Jose",-121.89,37.34),("Sacramento",-121.49,38.58),("Fresno",-119.77,36.75),("Long Beach",-118.19,33.77),("Bakersfield",-119.02,35.37)]
With both frames ready, the base scatter draws one point per district,
coloured by its price bucket:
```haskell
import qualified Data.Text as T
housingLayer =
(defLayer GeomPoint)
{ layerMapping =
emptyMapping
{ aesX = Just (ColumnRef "lng")
, aesY = Just (ColumnRef "lat")
, aesColor = Just (ColumnRef "bucket")
}
, layerAesDef = emptyAesDefaults { defSize = Just 2 }
}
baseChart =
emptyChart
{ chartData = housingFrame
, chartLayers = [housingLayer]
, chartTitle = Just "Median house value by location"
, chartSize = SizePixels 760 620
}
displayHtml (T.unpack (renderChartSvg baseChart))
```
>
>
> -120.0
> -115.0
> 35.0
> 40.0
> Median house value by location
> >500'000
> >450'000
> >400'000
> >300'000
> >200'000
> <=200'000
>
Cheap districts are by far the most numerous, so they crowd the plot. Mapping
point opacity to the raw price with `aesAlpha` fades the cheap districts and
keeps the expensive ones solid; overlaying the city names then shows which
cities sit in the pricey clusters:
```haskell
fadedLayer =
(defLayer GeomPoint)
{ layerMapping =
emptyMapping
{ aesX = Just (ColumnRef "lng")
, aesY = Just (ColumnRef "lat")
, aesColor = Just (ColumnRef "bucket")
, aesAlpha = Just (ColumnRef "value")
}
, layerAesDef = emptyAesDefaults { defSize = Just 2 }
}
cityLabels =
(defLayer GeomText)
{ layerData = Just citiesFrame
, layerMapping =
emptyMapping
{ aesX = Just (ColumnRef "lng")
, aesY = Just (ColumnRef "lat")
, aesLabel = Just (ColumnRef "name")
}
, layerAesDef = emptyAesDefaults { defColor = Just (Hex "#000000") }
}
spatialChart =
emptyChart
{ chartData = housingFrame
, chartLayers = [fadedLayer, cityLabels]
, chartTitle = Just "Median house value by location (faded by price), with major cities"
, chartSize = SizePixels 760 620
}
displayHtml (T.unpack (renderChartSvg spatialChart))
```
>
>
> -120.0
> -115.0
> 35.0
> 40.0
> Los Angeles
> San Francisco
> San Diego
> San Jose
> Sacramento
> Fresno
> Long Beach
> Bakersfield
> Median house value by location (faded by price), with major cities
> >500'000
> >450'000
> >400'000
> >300'000
> >200'000
> <=200'000
>
From the plot above we can see that there are specific population centres where
house prices are high. A really good model should be able to account for this
non-linearity. The expensive population centres are coastal, which suggests that
`ocean_proximity` will be useful for predicting house prices.
## Engineering features
Now let's build those derived attributes. We encode `ocean_proximity` as an
ordinal value. 5 means closest to the ocean, 1 means farthest inland.
```haskell
ordinalOceanProximity :: Text -> Double
ordinalOceanProximity op = case op of
"ISLAND" -> 5
"NEAR OCEAN" -> 4
"NEAR BAY" -> 3
"<1H OCEAN" -> 2
"INLAND" -> 1
_ -> error ("Unknown ocean proximity value: " ++ show op)
meanTotalBedrooms = D.meanMaybe total_bedrooms df
imputedTotalBedrooms = F.fromMaybe meanTotalBedrooms total_bedrooms
augmented = df |> D.deriveMany
[ "rooms_per_household" .= total_rooms / households
, "total_bedrooms" .= imputedTotalBedrooms
, "bedrooms_per_household" .= imputedTotalBedrooms / total_rooms
, "population_per_household".= population / households
, "ocean_proximity" .= F.lift ordinalOceanProximity ocean_proximity
]
```
## How well do the features correlate with price?
Let's see how each attribute correlates with our
target, `median_house_value`.
```haskell
import Data.Function (on)
import qualified Data.List as L
import Data.Maybe (fromMaybe)
correlationWithHouseValue columnName = (columnName, fromMaybe 0 (D.correlation columnName "median_house_value" augmented))
correlations = map correlationWithHouseValue (D.columnNames augmented)
mapM_ print (L.sortBy (flip compare `on` snd) correlations)
```
>
> ("median_house_value",1.0)
> ("median_income",0.6880752079585185)
> ("ocean_proximity",0.38324444779865596)
> ("rooms_per_household",0.1519482897414869)
> ("total_rooms",0.13415311380656322)
> ("housing_median_age",0.10562341249320993)
> ("households",6.584265057005642e-2)
> ("total_bedrooms",4.945354544321848e-2)
> ("population_per_household",-2.373741295613557e-2)
> ("population",-2.4649678888894886e-2)
> ("longitude",-4.5966615119112986e-2)
> ("latitude",-0.14416027687448565)
> ("bedrooms_per_household",-0.22004911702819213)
```haskell
corrAbbrev n = case n of
"median_house_value" -> "value"
"median_income" -> "income"
"ocean_proximity" -> "ocean"
"rooms_per_household" -> "rooms/hh"
"total_rooms" -> "rooms"
"housing_median_age" -> "age"
"households" -> "h.holds"
"total_bedrooms" -> "bedrms"
"population_per_household" -> "pop/hh"
"population" -> "pop"
"longitude" -> "lng"
"latitude" -> "lat"
"bedrooms_per_household" -> "bed/hh"
other -> other
corrColumns = D.columnNames augmented
corrCells =
[ (corrAbbrev a, corrAbbrev b, fromMaybe 0 (D.correlation a b augmented))
| a <- corrColumns, b <- corrColumns
]
corrFrame =
fromColumns
[ ("col", ColCat [a | (a, _, _) <- corrCells])
, ("row", ColCat [b | (_, b, _) <- corrCells])
, ("corr", ColNum [c | (_, _, c) <- corrCells])
, ("label", ColCat [T.pack (show (fromIntegral (round (c * 100)) / 100 :: Double)) | (_, _, c) <- corrCells])
]
```
```haskell
corrTile =
(defLayer GeomTile)
{ layerMapping =
emptyMapping
{ aesX = Just (ColumnRef "col")
, aesY = Just (ColumnRef "row")
, aesFill = Just (ColumnRef "corr")
}
, layerStat = StatIdentity
}
corrText =
(defLayer GeomText)
{ layerMapping =
emptyMapping
{ aesX = Just (ColumnRef "col")
, aesY = Just (ColumnRef "row")
, aesLabel = Just (ColumnRef "label")
}
, layerStat = StatIdentity
}
```
```haskell
corrChart =
emptyChart
{ chartData = corrFrame
, chartLayers = [corrTile, corrText]
, chartScales = defScales { scaleFill = Just (SColorContinuous [Hex "#2166ac", Hex "#f7f7f7", Hex "#b2182b"]) }
, chartTitle = Just "Feature correlation matrix"
, chartSize = SizePixels 820 760
}
displayHtml (T.unpack (renderChartSvg corrChart))
```
>
>
> lng
> lat
> age
> rooms
> bedrms
> pop
> h.holds
> income
> value
> ocean
> rooms/hh
> bed/hh
> pop/hh
> lng
> lat
> age
> rooms
> bedrms
> pop
> h.holds
> income
> value
> ocean
> rooms/hh
> bed/hh
> pop/hh
> 1.0
> -0.92
> -0.11
> 4.0e-2
> 7.0e-2
> 0.1
> 6.0e-2
> -2.0e-2
> -5.0e-2
> -9.0e-2
> -3.0e-2
> 8.0e-2
> 0.0
> -0.92
> 1.0
> 1.0e-2
> -4.0e-2
> -7.0e-2
> -0.11
> -7.0e-2
> -8.0e-2
> -0.14
> -0.16
> 0.11
> -9.0e-2
> 0.0
> -0.11
> 1.0e-2
> 1.0
> -0.36
> -0.32
> -0.3
> -0.3
> -0.12
> 0.11
> 0.21
> -0.15
> 0.13
> 1.0e-2
> 4.0e-2
> -4.0e-2
> -0.36
> 1.0
> 0.93
> 0.86
> 0.92
> 0.2
> 0.13
> -3.0e-2
> 0.13
> -0.18
> -2.0e-2
> 7.0e-2
> -7.0e-2
> -0.32
> 0.93
> 1.0
> 0.87
> 0.97
> -1.0e-2
> 5.0e-2
> 0.0
> 0.0
> 7.0e-2
> -3.0e-2
> 0.1
> -0.11
> -0.3
> 0.86
> 0.87
> 1.0
> 0.91
> 0.0
> -2.0e-2
> -3.0e-2
> -7.0e-2
> 1.0e-2
> 7.0e-2
> 6.0e-2
> -7.0e-2
> -0.3
> 0.92
> 0.97
> 0.91
> 1.0
> 1.0e-2
> 7.0e-2
> 2.0e-2
> -8.0e-2
> 3.0e-2
> -3.0e-2
> -2.0e-2
> -8.0e-2
> -0.12
> 0.2
> -1.0e-2
> 0.0
> 1.0e-2
> 1.0
> 0.69
> 0.15
> 0.33
> -0.51
> 2.0e-2
> -5.0e-2
> -0.14
> 0.11
> 0.13
> 5.0e-2
> -2.0e-2
> 7.0e-2
> 0.69
> 1.0
> 0.38
> 0.15
> -0.22
> -2.0e-2
> -9.0e-2
> -0.16
> 0.21
> -3.0e-2
> 0.0
> -3.0e-2
> 2.0e-2
> 0.15
> 0.38
> 1.0
> -0.11
> 7.0e-2
> -2.0e-2
> -3.0e-2
> 0.11
> -0.15
> 0.13
> 0.0
> -7.0e-2
> -8.0e-2
> 0.33
> 0.15
> -0.11
> 1.0
> -0.35
> 0.0
> 8.0e-2
> -9.0e-2
> 0.13
> -0.18
> 7.0e-2
> 1.0e-2
> 3.0e-2
> -0.51
> -0.22
> 7.0e-2
> -0.35
> 1.0
> 0.0
> 0.0
> 0.0
> 1.0e-2
> -2.0e-2
> -3.0e-2
> 7.0e-2
> -3.0e-2
> 2.0e-2
> -2.0e-2
> -2.0e-2
> 0.0
> 0.0
> 1.0
> Feature correlation matrix
>
Not bad! The engineered ratios are much more informative than the raw totals:
the number of rooms per household tells you more than the total number of rooms
in a district; obviously, the larger the houses, the more expensive they are.
And once we translate ocean proximity to an ordinal value (5 = close to the
ocean, 1 = far), we can see how predictive that feature is too. These are
exactly the features we'll feed the model.
## Modeling: linear regression with Hasktorch
Exploration done. We'll now train a linear-regression model on the `augmented`
features to predict `median_house_value`, using Hasktorch for the tensors and
gradient descent.
### Train/test split
We hold out 20% of the districts for testing, with a fixed seed for
reproducibility.
```haskell
import System.Random (mkStdGen)
let (train, test) = D.randomSplit (mkStdGen 42) 0.8 augmented
(D.dimensions train, D.dimensions test)
```
>
> ((16408,13),(4232,13))
### Building the tensors
Min-max normalisation rescales every numeric column to the `[0, 1]` range, which
keeps the different feature scales from dominating the gradient updates. We use
`D.fold` to apply the rescaling across every `Double` column.
```haskell
normalizeFeatures :: D.DataFrame -> D.DataFrame
normalizeFeatures df =
df
|> D.fold
( \name d ->
let
-- Convenience reference to the column.
col = F.col @Double name
in
D.derive name ((col - F.minimum col) / (F.maximum col - F.minimum col)) d
)
(D.columnNames (df |> D.selectBy [D.byProperty (D.hasElemType @Double)]))
```
We drop the target column (`median_house_value`) from the features, normalise
them, and convert both features and labels to Hasktorch tensors with
`toTensor`.
```haskell
import DataFrame.Hasktorch (toTensor)
let trainFeatures =
toTensor (train |> D.exclude [F.name median_house_value] |> normalizeFeatures)
let testFeatures =
toTensor (test |> D.exclude [F.name median_house_value] |> normalizeFeatures)
let trainLabels = toTensor (D.select [F.name median_house_value] train)
```
### The model
The model is one linear layer; `squeezeAll` drops the trailing size-1 dimension
so the output is a flat vector of predictions.
```haskell
import Torch
model :: Linear -> Tensor -> Tensor
model state input = squeezeAll $ linear state input
```
### Training
We initialise a linear layer with one output, then run 1,000 gradient-descent
steps minimising mean-squared error, logging the loss every 100 iterations.
(The original program runs 100,000 steps; 1,000 is plenty to watch the loss
fall and keeps the notebook responsive; bump `foldLoop init 1_000` higher for
a better fit.)
```haskell
import Control.Monad (when)
putStrLn "Training linear regression model..."
init <-
Torch.sample $
LinearSpec{in_features = snd (D.dimensions train) - 1, out_features = 1}
trained <- foldLoop init 1_000 $ \state i -> do
let labels' = model state trainFeatures
loss = mseLoss trainLabels labels'
when (i `mod` 100 == 0) $
putStrLn $ "Iteration: " ++ show i ++ " | Loss: " ++ show loss
(state', _) <- runStep state GD loss 0.1
pure state'
```
>
> Training linear regression model...
> Iteration: 100 | Loss: Tensor Float [] 8.3162e9
> Iteration: 200 | Loss: Tensor Float [] 6.9104e9
> Iteration: 300 | Loss: Tensor Float [] 6.2390e9
> Iteration: 400 | Loss: Tensor Float [] 5.8835e9
> Iteration: 500 | Loss: Tensor Float [] 5.6869e9
> Iteration: 600 | Loss: Tensor Float [] 5.5741e9
> Iteration: 700 | Loss: Tensor Float [] 5.5064e9
> Iteration: 800 | Loss: Tensor Float [] 5.4632e9
> Iteration: 900 | Loss: Tensor Float [] 5.4334e9
> Iteration: 1000 | Loss: Tensor Float [] 5.4111e9
### Predictions
Finally we run the trained model on the held-out test set and place the
predictions next to the true median house values.
```haskell
let predictionColumn = "predicted_house_value"
let predictions =
D.insert predictionColumn (asValue @[Float] (model trained testFeatures)) test
predictions |> D.select [F.name median_house_value, predictionColumn]
|> D.take 10
|> D.toMarkdown'
|> displayMarkdown
```
>
> | median_house_value