Monotonic constraints<\/h2>\n\n\n\n
On ML competition platforms like Kaggle<\/a>, complex and unintuitively behaving models dominate. In this respect, reality is completely different. There, the majority of models do not serve as pure prediction machines but rather as fruitful source of information. Furthermore, even if used as prediction machine, the users of the models might expect a certain degree of consistency when “playing” with input values. <\/p>\n\n\n\n A classic example are statistical house appraisal models. An additional bathroom or an additional square foot of ground area is expected to raise the appraisal, everything else being fixed (ceteris paribus). The user might lose trust in the model if the opposite happens.<\/p>\n\n\n\n One way to enforce such consistency is to monitor the signs of coefficients of a linear regression model. Another useful strategy is to impose monotonicity constraints<\/strong> on selected model effects.<\/p>\n\n\n\n Monotonicity constraints are especially simple to implement for decision trees. The rule is basically as follows: Tree ensembles like boosted trees or random forests will automatically inherit this property.<\/p>\n\n\n\n Most implementations of boosted trees offer monotonicity constraints. Here is a selection:<\/p>\n\n\n\n Unfortunately, the picture is completely different for random forests. At the time of writing, I am not aware of any random forest implementation in R or Python offering this useful feature.<\/p>\n\n\n\n Some options<\/strong><\/p>\n\n\n\n For the moment, let’s stick to option 3. In our last R <-> Python blog post<\/a>, we demonstrated that XGBoost’s random forest mode works essentially as good<\/strong> as standard random forest implementations, at least in regression settings and using sensible defaults.<\/p>\n\n\n\n Ask yourself: does the constraint really make sense for all possible values of other features? You will see that the answer is often “no”.<\/p>\n\n\n\n An example: If your house price model uses the features “number of rooms” and “living area”, then a monotonic constraint on “living area” might make sense (given any number of rooms), while such constraint would be non-sensical for the number of rooms. Why? Because having six rooms in a 1200 square feet home is not necessarily better than having just five rooms in an equally sized home<\/em>.<\/p>\n\n\n\n We use a nice dataset containing information on over 20,000 sold houses in Kings County. Along with the sale price, different features describe the size and location of the properties. The dataset is available on OpenML.org<\/a> with ID 42092.<\/p>\n\n\n\n The following R and Python codes <\/p>\n\n\n\n An ICE curve for variable X shows how the prediction of one specific observation changes if the value of X changes. Repeating this for multiple observations gives an idea of the effect of X. The average over multiple ICE curves produces the famous partial dependent plot.<\/em><\/p>\n\n\n <\/p>\n\n\n\n We clearly see many non-monotonic (and in this case counterintuitive) ICE curves.<\/p>\n\n\n\n What would a model give with monotonically increasing constraint on the ground area?<\/p>\n\n\n We see:<\/strong><\/p>\n\n\n\n The Python notebook and R code can be found at:<\/p>\n\n\n\n “R <-> Python” continued… Random forests with monotonic constraints<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[16],"tags":[10,6,5],"class_list":["post-575","post","type-post","status-publish","format-standard","hentry","category-machine-learning","tag-lost-in-translation","tag-python","tag-r"],"featured_image_src":null,"author_info":{"display_name":"Michael Mayer","author_link":"https:\/\/lorentzen.ch\/index.php\/author\/michael\/"},"_links":{"self":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/575","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/comments?post=575"}],"version-history":[{"count":34,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/575\/revisions"}],"predecessor-version":[{"id":688,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/575\/revisions\/688"}],"wp:attachment":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/media?parent=575"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/categories?post=575"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/tags?post=575"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}Trees and monotonic constraints<\/h2>\n\n\n\n
If a monotonicity constraint would be violated by a split on feature X<\/em>, it is rejected.<\/strong> (Or a large penalty is subtracted from the corresponding split gain.) This will imply monotonic behavior of predictions in X<\/em>, keeping all other features fixed. <\/p>\n\n\n\nBoosted trees<\/h3>\n\n\n\n
What about random forests?<\/h3>\n\n\n\n
Warning: Be careful with imposing monotonicity constraints<\/strong> <\/h3>\n\n\n\n
Let’s try it out<\/h2>\n\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/05\/image-1024x224.png\")
library(farff)\nlibrary(OpenML)\nlibrary(dplyr)\nlibrary(xgboost)\n\nset.seed(83454)\n\nrmse <- function(y, pred) {\n sqrt(mean((y-pred)^2))\n}\n\n# Load King Country house prices dataset on OpenML\n# ID 42092, https:\/\/www.openml.org\/d\/42092\ndf <- getOMLDataSet(data.id = 42092)$data\nhead(df)\n\n# Prepare\ndf <- df %>%\n mutate(\n log_price = log(price),\n log_sqft_lot = log(sqft_lot),\n year = as.numeric(substr(date, 1, 4)),\n building_age = year - yr_built,\n zipcode = as.integer(as.character(zipcode))\n )\n\n# Define response and features\ny <- \"log_price\"\nx <- c(\"grade\", \"year\", \"building_age\", \"sqft_living\",\n \"log_sqft_lot\", \"bedrooms\", \"bathrooms\", \"floors\", \"zipcode\",\n \"lat\", \"long\", \"condition\", \"waterfront\")\n\n# random split\nix <- sample(nrow(df), 0.8 * nrow(df))\ny_test <- df[[y]][-ix]\n\n# Fit untuned, but good(!) XGBoost random forest\ndtrain <- xgb.DMatrix(data.matrix(df[ix, x]),\n label = df[ix, y])\n\nparams <- list(\n objective = \"reg:squarederror\",\n learning_rate = 1,\n num_parallel_tree = 500,\n subsample = 0.63,\n colsample_bynode = 1\/3,\n reg_lambda = 0,\n max_depth = 20,\n min_child_weight = 2\n)\n\nsystem.time( # 25 s\n unconstrained <- xgb.train(\n params,\n data = dtrain,\n nrounds = 1,\n verbose = 0\n )\n)\n\npred <- predict(unconstrained, data.matrix(df[-ix, x]))\n\n# Test RMSE: 0.172\nrmse(y_test, pred)\n\n# ICE curves via our flashlight package\nlibrary(flashlight)\n\npred_xgb <- function(m, X) predict(m, data.matrix(X[, x]))\n\nfl <- flashlight(\n model = unconstrained,\n label = \"unconstrained\",\n data = df[ix, ],\n predict_function = pred_xgb\n)\n\nlight_ice(fl, v = \"log_sqft_lot\", indices = 1:9,\n evaluate_at = seq(7, 11, by = 0.1)) %>%\n plot()<\/pre><\/div>\n\n<\/div># Imports\nimport numpy as np\nfrom sklearn.datasets import fetch_openml\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.inspection import PartialDependenceDisplay\nfrom xgboost import XGBRFRegressor\n\n# Fetch data from OpenML\ndf = fetch_openml(data_id=42092, as_frame=True)[\"frame\"]\n\n# Prepare data\ndf = df.assign(\n year=lambda x: x.date.str[0:4].astype(int),\n zipcode=lambda x: x.zipcode.astype(int),\n log_sqft_lot=lambda x: np.log(x.sqft_lot),\n building_age=lambda x: x.year - x.yr_built,\n)\n\n# Feature list\nxvars = [\n \"grade\",\n \"year\",\n \"building_age\",\n \"sqft_living\",\n \"log_sqft_lot\",\n \"bedrooms\",\n \"bathrooms\",\n \"floors\",\n \"zipcode\",\n \"lat\",\n \"long\",\n \"condition\",\n \"waterfront\",\n]\n\n# Data split\ny_train, y_test, X_train, X_test = train_test_split(\n np.log(df[\"price\"]), df[xvars], train_size=0.8, random_state=766\n)\n\n# Modeling - wall time: 39 seconds\nparam_dict = dict(\n n_estimators=500,\n max_depth=20,\n learning_rate=1,\n subsample=0.63,\n colsample_bynode=1 \/ 3,\n reg_lambda=0,\n objective=\"reg:squarederror\",\n min_child_weight=2,\n)\n\nunconstrained = XGBRFRegressor(**param_dict).fit(X_train, y_train)\n\n# Test RMSE 0.176\npred = unconstrained.predict(X_test)\nprint(f\"RMSE: {mean_squared_error(y_test, pred, squared=False):.03f}\")\n\n# ICE and PDP - wall time: 47 seconds\nPartialDependenceDisplay.from_estimator(\n unconstrained,\n X=X_train,\n features=[\"log_sqft_lot\"],\n kind=\"both\",\n subsample=20,\n random_state=1,\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/11\/ice-1.png\")
# Monotonic increasing constraint\n(params$monotone_constraints <- 1 * (x == \"log_sqft_lot\"))\n\nsystem.time( # 179s\n monotonic <- xgb.train(\n params,\n data = dtrain,\n nrounds = 1,\n verbose = 0\n )\n)\n\npred <- predict(monotonic, data.matrix(df[-ix, x]))\n\n# Test RMSE: 0.176\nrmse(y_test, pred)\n\nfl_m <- flashlight(\n model = monotonic,\n label = \"monotonic\",\n data = df[ix, ],\n predict_function = pred_xgb\n)\n\nlight_ice(fl_m, v = \"log_sqft_lot\", indices = 1:9,\n evaluate_at = seq(7, 11, by = 0.1)) %>%\n plot()<\/pre><\/div>\n\n<\/div>
# One needs to pass the constraints as single string, which is rather ugly\nmc = \"(\" + \",\".join([str(int(x == \"log_sqft_lot\")) for x in xvars]) + \")\"\nprint(mc)\n\n# Modeling - wall time 49 seconds\nconstrained = XGBRFRegressor(monotone_constraints=mc, **param_dict)\nconstrained.fit(X_train, y_train)\n\n# Test RMSE: 0.178\npred = constrained.predict(X_test)\nprint(f\"RMSE: {mean_squared_error(y_test, pred, squared=False):.03f}\")\n\n# ICE and PDP - wall time 39 seconds\nPartialDependenceDisplay.from_estimator(\n constrained,\n X=X_train,\n features=[\"log_sqft_lot\"],\n kind=\"both\",\n subsample=20,\n random_state=1,\n)<\/pre><\/div>\n\n<\/div><\/div>\n\t\t<\/div>\n\n\n![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2021\/11\/ice_m.png\")
Summary<\/h2>\n\n\n\n