Development of a high-resolution indoor radon map using a new machine learning-based probabilistic model and German radon survey data

In plain words

Radon is a naturally occurring gas that leaks out of the ground and can build up inside houses. Germany ran a survey in which volunteers left detectors in their homes for a year. Those homes are not a random picture of the country: some floors and some soil types are over- or under-represented, so the plain average of the survey can be off.

The authors instead learn how radon depends on things that can be looked up for every building: the soil, the climate, the terrain, the age and type of the building, and the floor. A forest of decision trees is trained on the survey, but instead of one average it predicts a whole range of likely values for each floor of each building. Random draws from those ranges are then combined and weighted by how many people live on each floor, and the results are summarised for regions from federal states down to single municipalities. The output is a map of typical levels, of how often homes exceed a reference value, and of how many people that affects, together with an honest sense of how uncertain each number is.

TL;DR

A quantile regression forest, trained on a national survey and on mapped environmental and building predictors, gives a predicted radon distribution for each floor of about 21 million residential buildings. Monte Carlo sampling combines these distributions, weighted by residents per floor, into national and regional statistics. Everything is author-reported. The code is public, but the training data and the fitted model are not, so the predictions cannot be reproduced from the repository alone.

Problem

National indoor radon averages usually come from survey campaigns. Survey participation is voluntary and the sample does not match the population on floor level or on radon in soil p. 2 fig. 4 fig. 5. Descriptive statistics therefore carry bias, and they give no resolution below the survey's own regions p. 2. A model that uses predictors available across the whole building stock can correct for sampling bias and predict below survey level p. 2. The model must also carry prediction uncertainty to the aggregate level, because a single predicted value at a building is much less reliable than a predicted distribution. Quantities such as exceedance probabilities need that distribution, not a conditional mean p. 2.

Method

Overview: a survey of indoor radon gives the training rows. A quantile regression forest learns conditional quantiles from environmental and building predictors. Each building's quantiles are turned into a lognormal distribution per floor. The distributions are sampled, weighted by residents per floor, and aggregated p. 6 fig. 3.

1. Survey. Two habitable rooms per household were measured with solid-state nuclear track detectors for one year, from autumn 2019 to spring 2021. Each household was asked about building and occupancy characteristics p. 2. Sampling density was made roughly proportional to population at district level fig. 1 p. 3. The paper reports 7,479 households returning detectors p. 3 and 14,053 measurements from 7,061 households in the analysis p. 7.
2. Environmental predictors. Eight continuous predictors came from maps at a 500 m grid: radon in soil, outdoor temperature, slope, soil gas permeability, soil moisture, precipitation, wind exposure and tectonic fault density p. 7 fig. 2. Soil radon and soil gas permeability were themselves modelled from field data p. 4.
3. Building predictors. Four building-related predictors: floor level, building age (five classes), building type and number of living units p. 7. The paper says it started with 15 candidate predictors p. 6.
4. Building stock. The 2021 building dataset of the federal mapping agency has 21,866,120 buildings; 21,343,204 with at least one inhabitant were used p. 4. Missing floor levels (8%), building years (1%) and building types (1%) were imputed with linear models, and building type was filled from the number of households p. 5.
5. Floor levels. Each building is expanded to one row per level. Residents are spread evenly over the floors, and the basement carries a share. The share is 30% of an upper floor for single- and two-family houses and townhouses, and 5% for other types. The paper calls this an educated guess p. 5 p. 11.
6. Forest. Predictor selection is forward selection by RMSE over 10 spatial cross-validation folds, with blocks of 40 km p. 7. The loss is RMSE on the untransformed radon concentration p. 6. The final model is a conditional random forest with 500 trees and mtry 4 p. 7. Its cross-validated RMSE is 110.5 Bq/m³ and r² is 0.24 p. 8.
7. Quantile prediction. For each floor row, nine quantiles are predicted, at the 10th, 25th, 50th, 75th, 80th, 85th, 90th, 95th and 98th percentile p. 7.
8. Lognormal per floor. The local outdoor radon estimate is used as an offset. It is subtracted from the predicted quantiles, and a lognormal is fitted to the remainder. The fit uses a three-parameter lognormal: mean log, SD log and the offset p. 7. The paper says higher percentiles are weighted more heavily p. 7.
9. Monte Carlo sampling. For each floor, the number of samples is proportional to the expected residents on that floor. The paper uses a factor of 10 and reports 840,224,604 samples in total p. 7. Each sample is drawn from the floor's lognormal, the offset is added back, and the result is used as the indoor value p. 7.
10. Aggregation. Samples are grouped by municipality, district, federal state and Germany. Summary statistics are the arithmetic mean and SD, geometric mean and GSD, the 50th, 90th and 95th percentiles, and exceedance probabilities for 100, 300, 600 and 1,000 Bq/m³ p. 7 aggregate predictions.R.

# Indoor radon samples for one residential building  [p. 5] [p. 7]
# [code: IRC prediction individual buildings.R#L136-L223]
# Input:  building (residents, floors, type, outdoor radon), forest
# Output: integer indoor radon samples (Bq/m^3)

def sample_building(building, forest):
    # Basement share: 0.3 for single/two-family and row houses, else 0.05
    # [p. 5] [code#L136-L137]
    if building.is_small_house:
        basement_share = 0.3
    else:
        basement_share = 0.05
    # Residents per upper floor; the basement carries its share  [p. 5]
    # [code#L146-L153]
    per_floor = building.residents / (building.n_floors + basement_share)
    samples = []
    for level in [BASEMENT] + floors(building):
        # Expected samples = 10 x residents on this level, rounded  [p. 7]
        # [code#L172]
        weight = per_floor
        if level == BASEMENT:
            weight = per_floor * basement_share
        n_samples = round(10 * weight)
        if n_samples == 0:
            continue
        # Nine quantiles for this level  [p. 7] [code#L182-L186]
        quantiles = forest.predict_quantiles(building, level, PERCENTILES)
        # Offset = outdoor radon, held 1 Bq/m^3 below the 10th percentile
        # when it is at or above it  [code#L191-L195] (not in paper)
        offset = building.outdoor_radon
        if offset >= quantiles[P10]:
            offset = quantiles[P10] - 1
        # Lognormal fitted to quantiles minus offset  [p. 7] [code#L198-L205]
        # FIT_WEIGHTS is in the code only: [1, 2, 3, 4, 5, 6, 6, 6, 6]
        shifted = [q - offset for q in quantiles]
        mu, sigma = fit_lognormal(PERCENTILES, shifted, FIT_WEIGHTS)
        # Draw, add the offset back, round to whole Bq/m^3  [p. 7]
        # [code#L213-L220]
        draws = lognormal(mu, sigma, n_samples) + offset
        samples.extend(round(draws))
    return samples

Training, as the authors' code does it p. 7 model building.R:

# Forest training  [p. 6-7] [code: model building.R#L66-L105]
# Input:  survey rows (radon, 15 candidate predictors), spatial folds
# Output: conditional random forest for quantile prediction

def train_forest(survey, candidates, spatial_folds):
    # Forward feature selection, RMSE over spatial folds  [p. 7]
    # [code#L66-L75]
    selected = forward_select(survey, candidates, spatial_folds)
    # Target is untransformed radon; loss is RMSE  [p. 6] [code#L68-L69]
    # The code also tunes mtry over 2 to 12  [code#L82], but the final
    # fit hard-codes mtry = 4 and ntree = 500  [code#L105]
    return fit_forest(survey[selected], ntree=500, mtry=4)

Summary statistics per unit p. 7 aggregate predictions.R:

# Summary of one administrative unit  [p. 7]
# [code: aggregate predictions.R#L109-L122]
# Input:  integer indoor radon samples (Bq/m^3) inside the unit
# Output: arithmetic and geometric statistics, percentiles, exceedances

def summarise_unit(samples):
    # Arithmetic mean and SD of the raw samples  [code#L110-L111]
    arithmetic_mean = mean(samples)
    arithmetic_sd = std(samples)
    # Geometric mean and GSD from the log samples  [code#L113-L114]
    geometric_mean = exp(mean(log(samples)))
    geometric_gsd = exp(std(log(samples)))
    # Percentiles 50, 90, 95 and 99  [code#L115-L118]
    percentiles = quantiles(samples, [50, 90, 95, 99])
    # Exceedance shares use >= the threshold  [code#L119-L122]
    exceed_300 = share(samples >= 300)
    return arithmetic_mean, arithmetic_sd, geometric_mean, geometric_gsd

Figures are described in words only: fig. 1 shows survey sampling density, fig. 2 the predictor maps, fig. 3 the workflow, fig. 4 and fig. 5 the representativeness checks, fig. 6 the predictor ranking, fig. 7 and fig. 8 the cross-validated predictions and prediction intervals, and fig. 9 the output maps. The paper has no numbered equations p. 7. Supplementary Figures S1 and S2 and Excel tables S1 to S10 are referred to but are not in the source we hold, so they were not checked.

Hyperparameters

Results

Author-reported. Setup: national survey 2019-2021 p. 2 and model predictions for all residential buildings in Germany p. 7. Measured figures are the survey; predicted figures are the model's output.

Quantile coverage in 10-fold spatial CV (nominal → actual):

Source: tab. 1 p. 11.

Predictor importance, from highest to lowest, by permutation fig. 6 p. 9: floor level, radon in soil, building age, outdoor temperature, slope, soil gas permeability, soil moisture, precipitation, number of living units, wind exposure, building type, tectonic faults. Importance values are read from the plot and are not stated numerically.

Paper vs code

The paper links the code at https://github.com/EricPetermann/scripts_for_IRC_prediction_Germany p. 7. The owner is an author, so this is the official code, at commit aa075a332698534069ed18aa5b768c67f738e09f. The repository has no LICENSE file. It has three R scripts and a README. The training data and the fitted model are not in the repository.

• Training data and model are absent. The model script loads a data file that the authors say cannot be shared model building.R. The final forest is never saved; the save line is commented out model building.R. The prediction script loads a model file that the repository does not produce IRC prediction individual buildings.R. Predictions of the published model are available at model.radonmap.info p. 7.
• Predictor list. The code's 15 candidate predictors match the paper's 15 model building.R p. 6. The final eight-plus-four set is selected at runtime and not written to disk, so it cannot be checked against the paper's list.
• Final forest. ntree 500 and mtry 4 match p. 7 model building.R. The tuning step's result is not used; mtry is hard-coded model building.R model building.R.
• Outdoor radon. The paper subtracts the outdoor estimate from the predicted quantiles p. 7. The code also clamps it: when outdoor radon is at or above the 10th percentile, it is set to one below that percentile IRC prediction individual buildings.R. The paper does not mention this step.
• Fit weights. The lognormal fit weights are 1 to 6, rising with percentile IRC prediction individual buildings.R. The paper gives no weights p. 7.
• Floor categories. The paper lists basement, souterrain, ground, first, second and third-or-higher floors p. 7. The code's floor mapping has five labels and no souterrain IRC prediction individual buildings.R.
• Basement share. The 30% share is applied to two building-type labels, "EFH/ZFH" and "RH/DHH" IRC prediction individual buildings.R. The paper says "single- and two-family houses and townhouses" p. 5. Mapping "RH/DHH" to row or semi-detached houses is our reading of the code label; the paper's own building list includes such types p. 7.
• Exceedance threshold. The code counts a sample as exceeding a threshold when it is at or above it aggregate predictions.R. Parts of the paper write "> 300 Bq/m³" p. 13.
• Sample rounding. Samples are rounded to integers, and a floor with a weight that rounds to zero draws no samples IRC prediction individual buildings.R IRC prediction individual buildings.R IRC prediction individual buildings.R. The paper does not state this.
• Total sample count. The README gives "~800 million" random samples README.md; the paper gives 840,224,604 p. 7.
• Chunk coverage. Prediction loops over chunks 1 to 4400 IRC prediction individual buildings.R. The aggregation script reads only chunks up to 4374 aggregate predictions.R. If chunks 4375 to 4400 hold data, they are dropped.
• Aggregation as shipped. The aggregation script loops over seven files from a placeholder directory, ".../Predictions" aggregate predictions.R. It does not run unchanged.
• Model-fit script as shipped. The fold-metric function ends with a trailing comma inside c() model building.R. R rejects that, so the cross-validation step does not run unchanged.
• Spatial folds. The block size of 40 km p. 7 cannot be checked: the folds are loaded from a binary file model building.R.
• Age and type codes. The code maps 12 age classes to five labels IRC prediction individual buildings.R and nine building characteristics to labels IRC prediction individual buildings.R. The paper refers to a supplement for these mappings p. 5, which is not in our source.

Pitfalls

• Basement occupancy is the largest unknown. The authors call the 30% and 5% shares an educated guess p. 11. The survey has no basement occupancy data p. 5, and basements carry the highest radon levels in a building p. 5. Test the sensitivity of any downstream result to these two values.
• Predicted mean differs from the survey mean. The survey gives AM 78 Bq/m³ p. 7; the national prediction gives 63 Bq/m³ p. 9. The authors attribute this to the survey's over-sampling of basements and of high-radon areas p. 11. Do not treat the prediction as a corrected survey mean without that caveat.
• Smoothing at both ends. High observations are predicted too low and low observations too high p. 9 fig. 7. Tails of individual distributions are not reliable, so exceedance results for single buildings are weaker than the national figures.
• Prediction intervals are wide, and the 80% coverage is near nominal. The 80% interval covers 78.2% of test observations p. 9; per-building intervals are wide fig. 8. Report intervals, not single predictions.
• Offset guard. The outdoor-radon clamp keeps the indoor distribution above the outdoor value IRC prediction individual buildings.R. It is not in the paper, so an implementation that omits it will differ on buildings near the threshold.
• Floors of zero. A building with a floor count of zero gets a single basement row, and all residents go to that row IRC prediction individual buildings.R. This follows from the code's arithmetic and was not tested. Check your data for zero floor counts before relying on the output.
• Rounding drops buildings. A floor's samples are round(10 × weight), so floors whose weight is below 0.05 draw none IRC prediction individual buildings.R. Small buildings are underrepresented in the sampled totals.
• Sample volume. The run produces about 840 million samples. The authors split the output into several CSV files to save memory and aggregated them with arrow p. 7. Plan storage and memory accordingly.
• Lognormal fit has no error handling. The fit is run for every floor, with no check for a failed or degenerate fit IRC prediction individual buildings.R. The paper does not discuss failure cases.
• Exceedance wording. The paper writes thresholds both as "exceedance of 300" p. 9 and as "> 300 Bq/m³" p. 13, and never says whether a value equal to the threshold counts. The code counts equal values aggregate predictions.R. Pick one and state it.
• Household count. The paper gives 7,479 households returning detectors p. 3 and 7,061 households in the analysis p. 7. It does not explain the difference.
• Survey years. The survey ran from 2019 to 2021 p. 2. Lockdowns and remote work may have changed occupancy and ventilation during that time p. 12.
• Reproduction. The training data are not public model building.R. Without the survey, an implementer can reproduce the method but not the fitted forest.

Lineage

• Builds on: the authors' earlier geogenic-radon map of Germany p. 6 [p. 15, ref. 71], which used a similar predictor set. The current paper adds building data and floor-level modelling p. 13.
• Compared against: Menzler et al. 2008, an earlier national estimate p. 13 p. 16, ref. 88. The authors also compare with their earlier building-level estimate p. 13 p. 15, ref. 33.
• Superseded by: not known.