Exercise 4: Analytical applications

Work in progress. These tutorial materials are still under active development and will continue to change until the tutorial takes place on October 20, 2026. Content, links, and exercises may be incomplete or shift without notice.

Goal

Design an analysis-ready exposure feature and identify threats to validity and federated execution requirements. Exposure records are longitudinal facts; models require index-aligned covariates, and turning one into the other requires deliberate choices that can be gotten wrong in specific, checkable ways.

Steps

  1. Choose a research question that uses the exposure metric you calculated in Exercise 3 (or a similar one).
  2. Specify the feature: lookback window, lag, aggregation method, missingness handling, unit harmonization, geography, and update frequency.
  3. Decide whether this feature belongs in cohort eligibility or as a post-index predictor — these must stay separate. If your feature involves a spatially informed cohort, check it against appropriate vs. inappropriate uses.
  4. Identify three specific threats to validity for your feature, drawn from: ecological fallacy, exposure misclassification, residential mobility, MAUP (modifiable areal unit problem), temporal mismatch, collider/selection bias, spatial autocorrelation/clustering, site effects, or region-level confounding. For each, state concretely how it could bias your result, not just name it.
  5. Note what would need to stay local vs. what could travel in a federated setting: standardized concepts and code travel; precise locations and patient-level geographies remain local; only diagnostics and aggregate estimates are shared.
  6. Write the feature definition as reproducible SQL (or pseudocode if SQL access isn’t available in your environment).

Suggested worked example

If you do not have your own question, use this one. It is a cohort study of whether higher PM2.5 in the two years before an index date predicts incident COPD, with the temporal alignment made explicit:

  • Index date: 2016-01-01. Exposure: day-weighted mean PM2.5 over the 730 days before the index date. Outcome: first COPD diagnosis on or after the index date (through 2019-12-31).
  • Eligibility: two years of observation before the index date, a full 730 days of exposure coverage, and no COPD before the index date (prevalent cases are excluded).
  • Covariates: age at index, sex, and county SES (the SDOH or county_reference value at the index date).

The feature table below implements exactly that. It builds on the exposure from Exercise 3 and keeps eligibility (pre-index) separate from the outcome (post-index):

CREATE TABLE demo.ex4_features AS
WITH params AS (SELECT DATE '2016-01-01' AS index_date, 730 AS lookback_days),
cohort AS (
    SELECT p.person_id, p.year_of_birth, p.gender_concept_id, pa.index_date, pa.lookback_days
    FROM omopgis.person p
    CROSS JOIN params pa
    JOIN omopgis.observation_period op
      ON op.person_id = p.person_id
     AND op.observation_period_start_date <= pa.index_date - pa.lookback_days
     AND op.observation_period_end_date >= pa.index_date
    WHERE NOT EXISTS (SELECT 1 FROM omopgis.condition_occurrence co
                      WHERE co.person_id = p.person_id AND co.condition_concept_id = 255573
                        AND co.condition_start_date < pa.index_date)
),
expo AS (
    SELECT c.person_id,
           SUM(ee.value_as_number * (LEAST(ee.exposure_end_date, c.index_date - 1)
                                     - GREATEST(ee.exposure_start_date, c.index_date - c.lookback_days) + 1))
             / SUM(LEAST(ee.exposure_end_date, c.index_date - 1)
                   - GREATEST(ee.exposure_start_date, c.index_date - c.lookback_days) + 1) AS pm25_24m,
           SUM(LEAST(ee.exposure_end_date, c.index_date - 1)
               - GREATEST(ee.exposure_start_date, c.index_date - c.lookback_days) + 1) AS days_covered
    FROM cohort c
    JOIN omopgis.external_exposure ee
      ON ee.person_id = c.person_id AND ee.exposure_concept_id = 2052499839
     AND ee.exposure_start_date <= c.index_date - 1
     AND ee.exposure_end_date >= c.index_date - c.lookback_days
    GROUP BY c.person_id
),
ses AS (
    SELECT c.person_id, cr.ses_index
    FROM cohort c
    JOIN omopgis.location_history lh
      ON lh.entity_id = c.person_id AND lh.start_date <= c.index_date AND lh.end_date >= c.index_date
    JOIN omopgis.location l ON l.location_id = lh.location_id
    JOIN omopgis.county_reference cr ON cr.county_ref_id = l.county_ref_id
)
SELECT c.person_id,
       (EXTRACT(YEAR FROM c.index_date) - c.year_of_birth)::int AS age_at_index,
       CASE WHEN c.gender_concept_id = 8532 THEN 1 ELSE 0 END AS female,
       s.ses_index,
       e.pm25_24m,
       EXISTS (SELECT 1 FROM omopgis.condition_occurrence co
               WHERE co.person_id = c.person_id AND co.condition_concept_id = 255573
                 AND co.condition_start_date >= c.index_date)::int AS copd_after_index
FROM cohort c
JOIN expo e USING (person_id)
JOIN ses s USING (person_id)
WHERE e.days_covered = c.lookback_days;

Extracting the exposure covariate with FeatureExtractionForExtensions

The SQL above is the reference implementation. In a real study you would let FeatureExtractionForExtensions build the exposure covariate from external_exposure so it plugs into FeatureExtraction, PatientLevelPrediction, and CohortMethod. First turn the eligible persons into a standard cohort table:

CREATE TABLE demo.ex4_cohort AS
SELECT 1 AS cohort_definition_id, person_id AS subject_id,
       DATE '2016-01-01' AS cohort_start_date, DATE '2019-12-31' AS cohort_end_date
FROM demo.ex4_features;

Then, in RStudio in the pre-built tutorial image (FeatureExtractionForExtensions is already installed):

library(DatabaseConnector)
library(FeatureExtraction)
library(FeatureExtractionForExtensions)

connectionDetails <- createConnectionDetails(dbms = "postgresql", server = "gaia-db/gaiacore",
                                             port = 5432, user = "postgres", password = "<provided>")
connection <- connect(connectionDetails)

pm25Settings <- createExtensionCovariateSettings(
  analysisId = 998,
  extensionDatabaseSchema = "omopgis",
  extensionTableName = "external_exposure",
  extensionFields = c("exposure_concept_id", "value_as_number", "exposure_start_date"),
  joinField = "person_id",
  covariateIdField = "exposure_concept_id",
  covariateValueField = "value_as_number",
  dateField = "exposure_start_date",
  endDateField = "exposure_end_date",
  startDay = -730,
  endDay = -1,
  conceptSet = list(Capr::cs(2052499839, name = "PM2.5, CDC monthly county mean")),
  valueAggregation = "mean"
)

covariateData <- getDbCovariateData(
  connection = connection,
  cdmDatabaseSchema = "omopgis",
  cohortDatabaseSchema = "demo",
  cohortTable = "ex4_cohort",
  cohortIds = 1,
  rowIdField = "subject_id",
  covariateSettings = pm25Settings,
  aggregated = FALSE
)

Read the SQL the package generates and compare it with your reference query. Three things differ, and you should be able to say what each does to the estimate:

  1. The covariate value is an aggregate of the monthly rows in the window (valueAggregation: the default is the maximum; here it is the mean). The mean treats every monthly row as one observation, so partial months (a move in the middle of a month) count as full months and the result is not the day-weighted mean of the SQL reference.
  2. A row is included when its interval overlaps the window (endDateField), but it is used whole, not clipped to the window. Without endDateField the row’s start date has to fall inside the window.
  3. The covariate ID is the concept ID times 1000 plus the analysis ID, so the same concept extracted by different analyses (for example the last 12 months and the 12 months before that) gets different IDs.

Join the covariate to the outcome by subject_id (rowId), fit the model with this covariate and with the weighted-mean covariate, and compare the two estimates. Which definition recovers the simulated effect better, and why? (Aggregated covariates are not supported by this package yet, so this is a person-level extraction.) Check the repository README for the exported function names in your installed version.

Fit a logistic regression of copd_after_index on pm25_24m, once crude and once adjusted for ses_index, age_at_index, and female (for example glm(copd_after_index ~ pm25_24m + ses_index + age_at_index + female, family = binomial, data = features) in R). Then repeat with a negative-control outcome (for example CKD, concept 46271022, replacing the COPD concept in both places in the query).

What to expect, and what the generator knows. The true simulated PM2.5 effect on COPD is in demo.generator_truth (0.15 log-odds per µg/m³), and it is deliberately exaggerated (see the data description). On the reference build the cohort has about 9,600 persons and about 725 incident cases, the crude estimate is about 0.14 and the SES-adjusted estimate is about 0.11 per µg/m³ (standard error about 0.026), and the negative-control outcome shows a crude association of about 0.09 that disappears (about 0.00) after adjusting for SES. Your numbers will vary slightly with the build. Why is the adjusted estimate below the truth? (Think about exposure measured at the county level, residential moves, and the pre-index window.) Why does the negative control have a crude association at all?

Comparing the Exercise 3 cohorts with CohortMethod

The exposure-stratified cohorts you built in Exercise 3 (1 = high PM2.5 at index, 2 = low PM2.5 at index, 3 = COPD, 4 = CKD, all in demo.cohort) are the inputs to a comparative cohort analysis with CohortMethod: high versus low exposure, with propensity-score adjustment for confounders, and incident COPD as the outcome.

Choose the confounders carefully. County SES drives both PM2.5 and disease, and you only observe it through the SDOH observations (poverty rate, education, and so on), which are noisy. Extract them with FeatureExtractionForExtensions from the standard observation table; the concept set limits the extraction to the indicators you list. Leave out AQI_CATEGORY: it is computed from the person’s PM2.5, so adjusting for it would adjust away the exposure. Do not include the PM2.5 covariate from the previous section either, since it defines the groups. Why is that?

library(CohortMethod)
library(FeatureExtraction)
library(FeatureExtractionForExtensions)

# SDOH observation concepts, without AQI_CATEGORY (2052499437)
sdohIds <- c(2052499478, 2052497092, 2052498758, 2052498464, 2052499459, 2052498073,
             2052497364, 2052497761, 2052498200, 2052498483, 2052497310)

sdohSettings <- createExtensionCovariateSettings(
  analysisId = 997,
  extensionDatabaseSchema = "omopgis",
  extensionTableName = "observation",
  extensionFields = c("observation_concept_id", "value_as_number"),
  joinField = "person_id",
  covariateIdField = "observation_concept_id",
  covariateValueField = "value_as_number",
  conceptSet = list(Capr::cs(sdohIds, name = "SDOH indicators")),
  isBinary = FALSE,
  missingMeansZero = FALSE
)
demographics <- createCovariateSettings(useDemographicsGender = TRUE, useDemographicsAge = TRUE)

cmData <- getDbCohortMethodData(
  connectionDetails = connectionDetails,
  cdmDatabaseSchema = "omopgis",
  targetId = 1, comparatorId = 2, outcomeIds = c(3, 4),
  exposureDatabaseSchema = "demo", exposureTable = "cohort",
  outcomeDatabaseSchema = "demo", outcomeTable = "cohort",
  removeDuplicateSubjects = "remove all",
  firstExposureOnly = TRUE,
  washoutPeriod = 0,
  restrictToCommonPeriod = FALSE,
  covariateSettings = list(demographics, sdohSettings)
)

studyPop <- createStudyPopulation(
  cohortMethodData = cmData,
  outcomeId = 3,
  removeSubjectsWithPriorOutcome = TRUE, priorOutcomeLookback = 99999,
  minDaysAtRisk = 1,
  riskWindowStart = 0, startAnchor = "cohort start",
  riskWindowEnd = 0, endAnchor = "cohort end"
)
getAttritionTable(studyPop)

# The default cross-validated prior shrinks everything to zero for these weakly predictive covariates,
# so use a fixed ridge prior.
ps <- createPs(
  cohortMethodData = cmData, population = studyPop,
  prior = Cyclops::createPrior("normal", variance = 1, useCrossValidation = FALSE),
  control = Cyclops::createControl(maxIterations = 5000)
)
computePsAuc(ps)
matchedPop <- matchOnPs(ps, maxRatio = 1)
balance <- computeCovariateBalance(matchedPop, cmData)

crude    <- fitOutcomeModel(population = studyPop,   modelType = "cox")
adjusted <- fitOutcomeModel(population = matchedPop, modelType = "cox")
hr <- function(m) with(m$outcomeModelTreatmentEstimate, c(HR = exp(logRr), lower = exp(logLb95), upper = exp(logUb95)))
rbind(crude = hr(crude), matched = hr(adjusted))

Before you model anything, characterize the two exposure groups. Aggregated covariates summarize a cohort without downloading person-level rows (the statistics are computed in the database), and computeStandardizedDifference() shows how different the groups are before any adjustment:

pm25Settings$aggregateOnServer   # TRUE: statistics are computed in the database
groupData <- lapply(1:2, function(cohortId) {
  getDbCovariateData(
    connection = connection, cdmDatabaseSchema = "omopgis",
    cohortDatabaseSchema = "demo", cohortTable = "cohort", cohortIds = cohortId,
    covariateSettings = list(demographics, sdohSettings, pm25Settings), aggregated = TRUE)
})
sdiff <- computeStandardizedDifference(groupData[[1]], groupData[[2]], cohortId1 = 1, cohortId2 = 2)
head(sdiff[order(-abs(sdiff$stdDiff)), c("covariateName", "mean1", "mean2", "stdDiff")], 6)

On the reference build the two groups differ by about 9.9 versus 7.5 µg/m³ in mean PM2.5 over the 24 months before the index date (a standardized difference of about 1.6; the sign follows FeatureExtraction’s convention, so a negative value means the first cohort is higher). The SDOH indicators differ by about 0.5 standard deviations, in the direction of a higher-SES profile in the low-exposure group: broadband access, physician access, food insecurity, and neighborhood disadvantage head the list. That imbalance is the confounding the propensity model has to remove, and you can compare it with the balance after matching below.

Run the same pipeline with outcome 4 (CKD, replacing outcomeId = 3 in createStudyPopulation) as a negative control. Answer these in your write-up:

  • Which persons does the attrition table remove, and why (prior COPD, no time at risk)?
  • How well does the propensity model separate the groups, and is covariate balance below 0.1 standardized difference after matching? Which covariates remain imbalanced?
  • How do the crude and adjusted hazard ratios compare, and how do they compare with the negative control?

What to expect. The image has CohortMethod 5.5; the code above is the version tested on the reference build. The attrition table shows 2,717 and 3,099 persons in the two groups, 2,608 and 3,009 after removing prior COPD, and about 2,100 matched pairs. The propensity model separates the groups only weakly (AUC about 0.65, because the SDOH indicators only weakly distinguish high from low exposure), and matching brings the largest standardized difference below 0.1. For COPD the crude hazard ratio is about 1.43 (95% CI 1.17 to 1.73) and the matched one about 1.26 (1.01 to 1.57). For the CKD negative control the crude hazard ratio is about 1.20 (1.01 to 1.43) and the matched one about 0.94 (0.77 to 1.15). Your numbers will vary slightly with the build. Two things go wrong if you change the setup: if you add AQI_CATEGORY (2052499437) to sdohIds, createPs stops with “High correlation between covariate(s) and treatment” (a useful safety net that catches the leak), and with the default cross-validated prior the model fits but every coefficient is zero (AUC 0.5), so matching does nothing. The matched estimate is pulled toward the null by the confounding adjustment, but it stays above 1 for COPD and about 1 for CKD, which is the pattern the simulation was built to show; it is also attenuated because the groups are defined by a single month’s exposure and SES is only measured through noisy indicators. You can compare it with the regression estimate from the SQL feature table: they answer slightly different questions (a two-group contrast versus a slope per µg/m³), so compare the direction, the confounding pattern, and the negative control rather than the numbers.

Deliverable

One feature specification, plus three anticipated biases (with the concrete mechanism for each, not just the name), and, from the worked example, the crude and adjusted estimates for COPD and for the negative control (the regression estimates and the CohortMethod hazard ratios) with a short interpretation.

Check Your Work


Closing synthesis

Dataset → metadata → source layer → location interval → exposure fact → covariate → evidence. Each session in this tutorial produced one link in that chain; the chain only holds if every link stays traceable back to its source.


Previous: Exercise 3