Author: BlueGreenCity-KI Research Team
Date: January 2026
Purpose: Investigate BART as an alternative or addon to static risk scoring via literature-based weights
Bayesian Additive Regression Trees (BART) is a sophisticated machine learning method that could complement or replace the current literature-based weighted risk scoring system in BlueGreenCity-KI.
Key Findings:
Recommendation: Start with a hybrid approach using BART predictions alongside static risk scores, then transition gradually as confidence in the model grows.
Imagine you need to predict which urban trees are at risk of drought stress. Instead of manually deciding how important each factor is (age, soil type, temperature, etc.), BART lets the data tell you.
Simple Analogy:
Think of BART as a committee of decision tree experts:
| For example: “Is the tree older than 50 years? Yes → more risk | No → check soil type” |
INPUT: Tree characteristics
├── Age: 45 years
├── Soil: Sealed
├── Temperature: 32°C
├── Species: Tilia cordata
└── Group size: 1 (isolated)
↓
┌─────────────────────────────────────────┐
│ BART MODEL (Committee of Trees) │
│ │
│ Tree 1: "Old + sealed → high risk" │
│ Tree 2: "Hot + isolated → high risk" │
│ Tree 3: "Linden + age > 40 → risk" │
│ ... (50-200 trees voting) │
│ │
│ Average vote: 0.78 │
│ Uncertainty: ±0.12 │
└─────────────────────────────────────────┘
↓
OUTPUT: Risk prediction
├── Risk Score: 0.78 (high risk)
└── Confidence Interval: [0.66, 0.90]
Key Insight: BART doesn’t just give you a prediction—it tells you how confident it is. This is crucial for decision-making!
The “Bayesian” part means BART uses probability theory to:
This is like a scientist who:
“Additive” means the model combines many small predictions:
For those with technical background, here’s the formal definition:
Model Structure:
Y = f(x) + ε
where:
f(x) = Σ(m=1 to M) g(x; T_m, M_m)
Y = outcome (e.g., tree vitality deterioration)
x = predictors (age, temperature, soil type, etc.)
g(x; T_m, M_m) = m-th regression tree
T_m = tree structure (splits)
M_m = terminal node parameters (predictions)
ε ~ N(0, σ²) = error term
M = number of trees (typically 50-200)
Bayesian Prior:
BART uses a regularizing prior that:
Prior on tree structure:
- P(tree depth = d) = α(1 + d)^(-β) [prefer shallow trees]
- Each terminal node: μ ~ N(0, σ²_μ) [small contributions]
- Variance: σ² ~ Inv-Gamma(ν/2, νλ/2)
Sampling Process (Markov Chain Monte Carlo):
Prediction:
For a new tree observation:
ŷ = (1/S) Σ(s=1 to S) f_s(x)
where:
S = number of posterior samples
f_s(x) = prediction from sample s
This gives both a point estimate (mean) and credible intervals (uncertainty).
How it works:
criteria:
- name: age
weight: 0.2
- name: vitality
weight: 0.3
- name: street_type
weight: 0.2
- name: air_temperature
weight: 0.15
- name: planting_structure
weight: 0.15
overall_risk = Σ (criterion_score × weight)
Strengths:
Limitations:
How it works:
1. Learn from historical data (2010-2024 trees)
2. Discover patterns and interactions automatically
3. Predict risk for current trees with uncertainty
Strengths:
Limitations:
| Aspect | Literature Weights | BART |
|---|---|---|
| Interpretability | ⭐⭐⭐⭐⭐ High | ⭐⭐⭐ Medium |
| Accuracy | ⭐⭐⭐ Good | ⭐⭐⭐⭐⭐ Excellent |
| Uncertainty quantification | ❌ None | ✅ Full posterior |
| Captures interactions | ❌ No | ✅ Yes |
| Training data required | ❌ No | ✅ Yes (1000+ trees) |
| Computational cost | ⭐⭐⭐⭐⭐ Instant | ⭐⭐⭐ Minutes |
| Adaptability | ❌ Manual updates | ✅ Learns from data |
| Stakeholder trust | ⭐⭐⭐⭐⭐ High | ⭐⭐⭐ Medium |
| Edge case handling | ⭐⭐⭐ Predictable | ⭐⭐⭐⭐ Robust |
Static approach:
# Manual assignment based on literature
age: 0.2
vitality: 0.3
temperature: 0.15
BART approach:
Example: Static scoring misses this:
"Young trees on sealed soil under high temperature
are MORE at risk than the sum of individual factors"
BART discovers: risk = f(age, soil, temp, age×soil, age×temp, soil×temp, ...)
Tree A: Risk = 0.75 ± 0.05 [very confident]
Tree B: Risk = 0.75 ± 0.25 [uncertain - maybe new conditions]
This helps prioritize:
Can rigorously evaluate:
Minimum requirements:
Challenge for BlueGreenCity-KI:
Training time:
Prediction time:
Stakeholder questions:
BART answer:
Implementation needs:
Python libraries:
pyBART (limited functionality)bartpy (experimental)BART package)Option 1: R with rpy2 (Recommended)
# R packages (in R console)
install.packages("BART") # Core BART implementation
install.packages("dbarts") # Alternative, faster implementation
install.packages("bartMachine") # Java-based, very fast
# Python integration
pip install rpy2 # Call R from Python
Option 2: Pure Python (Limited)
pip install pymc-bart # PyMC-based BART
pip install bartpy # Pure Python (experimental)
Training data structure:
# Features (X): DataFrame with predictors
X_train = pd.DataFrame({
'age': [45, 23, 67, ...],
'vitality_2010': [2, 1, 3, ...],
'soil_sealed': [1, 0, 1, ...],
'air_temp_mean_summer': [29.5, 28.2, 31.1, ...],
'species_group': [1, 2, 1, ...], # categorical
# ... all criteria scores
})
# Labels (y): Binary or continuous outcomes
y_train = pd.Series([1, 0, 1, ...]) # 1 = deteriorated, 0 = stable
# Test data (current trees to predict)
X_test = pd.DataFrame({...}) # same structure as X_train
"""
Step 1: Extract training data from historical tree records
"""
import pandas as pd
import numpy as np
from bluegreencity_ki.ml.labels import LabelExtractor
from sqlalchemy import create_engine
# Setup database connection
# engine = create_engine('postgresql://user:password@host:port/database')
# Load historical data (2010 baseline, 2024 outcomes)
extractor = LabelExtractor()
historical_data = extractor.extract_labels(
baseline_year=2010,
outcome_year=2024,
labels=['vitality_deterioration', 'maintenance_stress']
)
# Join with risk criteria (from pipeline)
# Note: Use the engine from above (uncomment and configure first)
risk_features = pd.read_sql(
"SELECT * FROM bgcki.analysis_tree_risk_scores_2010",
engine
)
# Combine
training_data = historical_data.merge(
risk_features,
on='gisid',
how='inner'
)
# Define features and labels
feature_cols = [
'age_risk', 'vitality_risk', 'street_type_risk',
'air_temperature_risk', 'planting_structure_risk',
'age_value', 'st_umfang', 'kr_durchm', # raw attributes
'species_drought_tolerance' # if available
]
X_train = training_data[feature_cols]
y_train = training_data['vitality_deterioration'] # binary outcome
print(f"Training samples: {len(X_train)}")
print(f"Positive rate: {y_train.mean():.2%}")
"""
Step 2: Train BART model using R's dbarts package
"""
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr
# Activate pandas conversion
pandas2ri.activate()
# Import R packages
dbarts = importr('dbarts')
base = importr('base')
# Convert data to R
X_train_r = pandas2ri.py2rpy(X_train)
y_train_r = ro.FloatVector(y_train)
# Train BART model
print("Training BART model...")
bart_model = dbarts.bart(
x_train=X_train_r,
y_train=y_train_r,
ntree=200, # number of trees
ndpost=1000, # posterior samples after burn-in
nskip=1000, # burn-in iterations
keepevery=10, # thinning (keep every 10th sample)
printevery=100, # progress reporting
nchain=4, # parallel chains for convergence check
)
print("Training complete!")
# Check convergence
# (examine trace plots, Gelman-Rubin statistic)
"""
Step 3: Predict risk for current trees
"""
import numpy as np
# Load current tree inventory
# Note: Use the engine configured in Step 1
current_trees = pd.read_sql(
"SELECT * FROM bgcki.analysis_tree_risk_scores",
engine
)
X_test = current_trees[feature_cols] # using feature_cols from Step 1
# Convert to R
X_test_r = pandas2ri.py2rpy(X_test)
# Predict (returns posterior samples)
predictions = dbarts.predict(bart_model, newdata=X_test_r)
# Convert back to Python
pred_samples = pandas2ri.rpy2py(predictions)
# Compute statistics
current_trees['bart_risk_mean'] = pred_samples.mean(axis=0)
current_trees['bart_risk_sd'] = pred_samples.std(axis=0)
current_trees['bart_risk_q025'] = np.percentile(pred_samples, 2.5, axis=0)
current_trees['bart_risk_q975'] = np.percentile(pred_samples, 97.5, axis=0)
print("Predictions complete!")
print(f"Mean risk: {current_trees['bart_risk_mean'].mean():.3f}")
print(f"High risk (>0.7): {(current_trees['bart_risk_mean'] > 0.7).sum()}")
"""
Step 4: Model evaluation and validation
"""
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score
from sklearn.metrics import roc_auc_score, classification_report
# Cross-validation (if implementing in scikit-learn compatible way)
# Note: True Bayesian BART doesn't fit sklearn API directly
# For validation, use temporal split
# Note: training_data and feature_cols from Step 1
X_train_2010 = training_data[training_data['year'] == 2010][feature_cols]
y_train_2010 = training_data[training_data['year'] == 2010]['vitality_deterioration']
X_val_2015 = training_data[training_data['year'] == 2015][feature_cols]
y_val_2015 = training_data[training_data['year'] == 2015]['vitality_deterioration']
# Train on 2010, predict 2015
# ... (similar to Step 2 and 3)
# Evaluate
# Note: pred_samples would come from predictions on validation data
y_pred = pred_samples.mean(axis=0)
y_pred_binary = (y_pred > 0.5).astype(int)
print("Validation Results:")
print(f"ROC-AUC: {roc_auc_score(y_val_2015, y_pred):.3f}")
print("\nClassification Report:")
print(classification_report(y_val_2015, y_pred_binary))
# Calibration plot
import matplotlib.pyplot as plt
def plot_calibration(y_true, y_pred_proba, n_bins=10):
"""Plot calibration curve"""
bins = np.linspace(0, 1, n_bins + 1)
bin_centers = (bins[:-1] + bins[1:]) / 2
bin_true_rates = []
for i in range(n_bins):
mask = (y_pred_proba >= bins[i]) & (y_pred_proba < bins[i+1])
if mask.sum() > 0:
bin_true_rates.append(y_true[mask].mean())
else:
bin_true_rates.append(np.nan)
plt.figure(figsize=(8, 6))
plt.plot([0, 1], [0, 1], 'k--', label='Perfect calibration')
plt.plot(bin_centers, bin_true_rates, 'o-', label='BART')
plt.xlabel('Predicted probability')
plt.ylabel('Observed frequency')
plt.title('Calibration Plot')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('bart_calibration.png', dpi=150, bbox_inches='tight')
plot_calibration(y_val_2015, y_pred)
"""
Step 5: Extract variable importance
"""
# Variable importance from BART
# (Based on how often each variable is used for splits)
# In R:
ro.r('''
variable_importance <- function(bart_model) {
# Extract variable selection counts
var_counts <- colMeans(bart_model$varcount)
var_props <- var_counts / sum(var_counts)
result <- data.frame(
variable = names(var_props),
importance = var_props
)
result <- result[order(-result$importance), ]
return(result)
}
''')
var_imp = ro.r['variable_importance'](bart_model)
var_imp_df = pandas2ri.rpy2py(var_imp)
print("Variable Importance:")
print(var_imp_df)
# Visualize
import matplotlib.pyplot as plt
import seaborn as sns
plt.figure(figsize=(10, 6))
sns.barplot(data=var_imp_df.head(15), x='importance', y='variable')
plt.xlabel('Importance (proportion of splits)')
plt.title('BART Variable Importance')
plt.tight_layout()
plt.savefig('bart_variable_importance.png', dpi=150)
"""
Step 6: Partial dependence plots (understand feature effects)
"""
import numpy as np
import matplotlib.pyplot as plt
# Note: pandas2ri and dbarts should be imported from Step 2 context
# If running standalone, uncomment:
# from rpy2.robjects import pandas2ri
# from rpy2.robjects.packages import importr
# dbarts = importr('dbarts')
def partial_dependence_bart(model, X, feature_idx, n_points=50):
"""
Compute partial dependence for a feature
"""
# Get range of feature values
feature_values = np.linspace(
X.iloc[:, feature_idx].min(),
X.iloc[:, feature_idx].max(),
n_points
)
# For each value, average predictions over dataset
pd_values = []
for val in feature_values:
X_temp = X.copy()
X_temp.iloc[:, feature_idx] = val
X_temp_r = pandas2ri.py2rpy(X_temp)
preds = dbarts.predict(model, newdata=X_temp_r)
pd_values.append(pandas2ri.rpy2py(preds).mean())
return feature_values, np.array(pd_values)
# Example: age effect
# Note: bart_model and X_test from earlier steps
age_values, age_pd = partial_dependence_bart(
bart_model, X_test, feature_idx=0, n_points=50
)
plt.figure(figsize=(8, 5))
plt.plot(age_values, age_pd, linewidth=2)
plt.xlabel('Age Risk Score')
plt.ylabel('Predicted Risk')
plt.title('Partial Dependence: Age Effect')
plt.grid(alpha=0.3)
plt.savefig('bart_pd_age.png', dpi=150, bbox_inches='tight')
# bluegreencity_ki/ml/models/bart_risk_model.py
from dataclasses import dataclass
import pandas as pd
import numpy as np
from typing import Optional
@dataclass
class BARTRiskResult:
"""Results from BART risk model"""
risk_mean: pd.Series
risk_std: pd.Series
risk_q025: pd.Series
risk_q975: pd.Series
variable_importance: pd.DataFrame
class BARTRiskModel:
"""BART-based risk prediction model"""
def __init__(self, config: dict):
self.config = config
self.model = None
self.feature_names = None
def train(self, X: pd.DataFrame, y: pd.Series):
"""Train BART model on historical data"""
# Implementation using rpy2 + dbarts
# (See Step 2 above)
pass
def predict(self, X: pd.DataFrame) -> BARTRiskResult:
"""Predict risk for new trees"""
# Implementation using rpy2 + dbarts
# (See Step 3 above)
pass
def save_model(self, path: str):
"""Save trained model to disk"""
pass
def load_model(self, path: str):
"""Load trained model from disk"""
pass
# bluegreencity_ki/pipelines/hybrid_risk_scoring.py
def run_hybrid(ctx):
"""
Compute both static and BART-based risk scores
"""
# 1. Run standard weighted scoring
static_results, stats, profiler, dists = risk_scoring.run(ctx)
# 2. If BART model available, add predictions
if ctx.config.get('bart', {}).get('enabled', False):
bart_model = BARTRiskModel(ctx.config['bart'])
bart_model.load_model(ctx.config['bart']['model_path'])
# Prepare features
feature_cols = [f"{c['name']}_risk" for c in ctx.config['criteria']]
X = static_results[feature_cols]
# Predict
bart_results = bart_model.predict(X)
# Add to results
static_results['bart_risk_mean'] = bart_results.risk_mean
static_results['bart_risk_std'] = bart_results.risk_std
static_results['bart_risk_lower'] = bart_results.risk_q025
static_results['bart_risk_upper'] = bart_results.risk_q975
# Ensemble prediction (average)
static_results['ensemble_risk'] = (
0.5 * static_results['overall_drought_risk'] +
0.5 * static_results['bart_risk_mean']
)
return static_results, stats, profiler, dists
# config/config.yaml
bart:
enabled: true
model_path: "./models/bart_drought_risk_v1.rds"
training:
baseline_year: 2010
outcome_year: 2024
label_type: "vitality_deterioration"
hyperparameters:
ntree: 200
ndpost: 1000
nskip: 1000
ensemble:
# How to combine static and BART scores
static_weight: 0.5
bart_weight: 0.5
Essential:
Recommended minimum:
From existing criteria:
features = [
# Risk scores (normalized)
'age_risk', 'vitality_risk', 'street_type_risk',
'air_temperature_risk', 'planting_structure_risk',
# Raw values (additional information)
'age_value', 'st_umfang', 'kr_durchm',
# Environmental
'soil_sealed', 'soil_water_capacity',
# Climate (if available)
'summer_temp_mean', 'precip_growing_season',
]
Derived features:
derived_features = [
# Interactions
'age_x_temperature': age_risk * air_temperature_risk,
'age_x_soil': age_risk * soil_sealed,
# Ratios
'crown_per_stem': kr_durchm / st_durchm,
# Spatial
'neighbor_count_10m': spatial_count(trees, radius=10),
'avg_neighbor_vitality': spatial_avg(trees, 'vitality', radius=10),
]
Handle missing data:
# BART can handle missing values via surrogate splits
# But better to impute if possible
from sklearn.impute import SimpleImputer
imputer = SimpleImputer(strategy='median')
X_imputed = pd.DataFrame(
imputer.fit_transform(X),
columns=X.columns,
index=X.index
)
Handle categorical variables:
# One-hot encoding for species, districts, etc.
X_encoded = pd.get_dummies(
X,
columns=['species_group', 'district'],
drop_first=True
)
Class imbalance:
# If positive cases are rare (e.g., 5% deterioration rate)
# Consider:
# 1. Oversampling minorities (SMOTE)
from imblearn.over_sampling import SMOTE
X_resampled, y_resampled = SMOTE().fit_resample(X, y)
# 2. Weighted loss (if available in BART implementation)
# 3. Stratified sampling
Use both systems in parallel:
┌─────────────────────────────────────────────┐
│ Current Tree Inventory │
└────────────┬────────────────────────────────┘
│
┌─────┴─────┐
│ │
▼ ▼
┌─────────────┐ ┌──────────────┐
│ Static │ │ BART │
│ Scoring │ │ Scoring │
│ (weighted) │ │ (ML-based) │
└──────┬──────┘ └──────┬───────┘
│ │
│ ▼
│ Risk: 0.73 ± 0.12
▼
Risk: 0.68
│
└───────┬───────┘
▼
┌──────────────────────┐
│ Combined Output │
│ ├─ Static: 0.68 │
│ ├─ BART: 0.73 │
│ ├─ Agreement: Good │
│ └─ Ensemble: 0.71 │
└──────────────────────┘
Advantages:
Use static scoring first, then adjust with BART:
# Step 1: Static risk
static_risk = compute_static_risk(tree)
# Step 2: BART adjustment
bart_adjustment = bart_model.predict_adjustment(tree, static_risk)
# Step 3: Final risk
final_risk = static_risk + bart_adjustment # or multiply
Advantages:
Replace static scoring once validated:
# Train comprehensive BART model
bart_model.train(X_historical, y_outcomes)
# Predict directly
risk = bart_model.predict(X_current)
# Static scores become features, not the answer
Requirements:
Combine multiple models:
models = {
'static': StaticRiskModel(),
'bart': BARTRiskModel(),
'random_forest': RandomForestModel(), # alternative
'xgboost': XGBoostModel(), # alternative
}
# Weighted ensemble
ensemble_risk = (
0.3 * models['static'].predict(X) +
0.4 * models['bart'].predict(X).mean +
0.2 * models['random_forest'].predict_proba(X)[:, 1] +
0.1 * models['xgboost'].predict_proba(X)[:, 1]
)
Advantages:
Problem: Limited resources, must prioritize which trees to water first.
Static approach:
Top 100 high-risk trees (risk > 0.8)
→ Water all 100
BART approach:
Top 100 by risk, but consider uncertainty:
- 50 trees: risk = 0.85 ± 0.05 [HIGH CONFIDENCE]
- 30 trees: risk = 0.82 ± 0.15 [MEDIUM CONFIDENCE]
- 20 trees: risk = 0.80 ± 0.30 [LOW CONFIDENCE]
→ Water 50 high-confidence trees first
→ Monitor 30 medium-confidence
→ Inspect 20 low-confidence (may be data issues)
Problem: Risk thresholds may vary by species.
BART discovers:
Linden (Tilia cordata):
High risk threshold: 0.75
Key factors: age, soil sealing
Oak (Quercus robur):
High risk threshold: 0.85 [more drought tolerant]
Key factors: soil, temperature
Application:
risk_threshold = bart_model.get_species_threshold(tree.species)
if tree.bart_risk > risk_threshold:
priority = "HIGH"
Problem: Predict which trees will deteriorate next year.
BART trained on:
X: Tree conditions in year t
y: Deterioration in year t+1
Application:
# Current conditions (2025)
current_conditions = get_tree_features(year=2025)
# Predict deterioration in 2026
risk_2026 = bart_model.predict(current_conditions)
# Proactive intervention
high_risk_trees = trees[risk_2026 > 0.7]
schedule_preventive_care(high_risk_trees)
Problem: What happens if we improve soil conditions?
BART allows:
# Current scenario
X_current = get_features()
risk_current = bart_model.predict(X_current).mean()
# Hypothetical: unsealed soil
X_unsealed = X_current.copy()
X_unsealed['soil_sealed'] = 0
risk_unsealed = bart_model.predict(X_unsealed).mean()
# Impact
risk_reduction = risk_current - risk_unsealed
print(f"Unsealing reduces risk by {risk_reduction:.2%}")
Extract knowledge from BART:
Example analysis:
# Partial dependence: age effect
age_effect = bart_model.partial_dependence('age_risk')
# Result: Risk increases slowly until age 40,
# then accelerates rapidly
# → Focus interventions on trees > 40 years
When to use BART:
When to stick with static scoring:
Recommended strategy:
Year 1 (starting from project approval):
├─ Implement BART in parallel
├─ Compare with static scoring
├─ Build trust and understanding
└─ No operational changes
Year 2:
├─ Use BART for prioritization refinement
├─ Flag high-uncertainty predictions
├─ Validate with field observations
└─ Gradual operational integration
Year 3+:
├─ BART as primary risk model
├─ Static scoring as interpretability aid
├─ Continuous model updates
└─ Full integration with decision-making
BART Methodology:
Applications in Ecology:
Urban Tree Risk:
R Packages:
Python Integration:
Tutorials:
Why BART over:
Urban tree datasets for testing:
Example BART workflow:
# R code for BART model
library(dbarts)
# Load data
data <- read.csv("tree_data.csv")
# Define features and outcome
X <- data[, c("age", "soil", "temp", "species")]
y <- data$deteriorated
# Train BART
bart_model <- bart(
x.train = X,
y.train = y,
ntree = 200,
ndpost = 1000,
nskip = 1000
)
# Predict
X_new <- new_tree_data[, c("age", "soil", "temp", "species")]
predictions <- predict(bart_model, X_new)
# Analyze
mean_risk <- colMeans(predictions)
sd_risk <- apply(predictions, 2, sd)
Where to get help:
Bayesian Additive Regression Trees represent a powerful and promising approach for urban tree risk assessment in the BlueGreenCity-KI system. While more complex than the current literature-based weighted scoring, BART offers significant advantages in predictive accuracy, uncertainty quantification, and automatic discovery of complex patterns.
Key Takeaways:
Next Steps:
→ Review this report with team and stakeholders
→ Assess data availability and quality
→ Decide on timeline for pilot implementation
→ Allocate resources (time, expertise, compute)
→ Begin Phase 1 (data preparation)
For questions or to discuss implementation details, contact the ML team or refer to the BlueGreenCity-KI ML module documentation.
Document Version: 1.0
Last Updated: January 2026
Authors: BlueGreenCity-KI Research Team
Status: Investigation Complete - Awaiting Decision on Implementation