Chapter 9 Modelling: a federated logistic GLM

Because a dsOMOP extract is just a server-side data frame, the whole dsBaseClient modelling toolchain applies to it. We fit a federated logistic regression of gastrointestinal haemorrhage on demographics, NSAID exposure and prior peptic ulcer — pooling only the model’s sufficient statistics across the three sites. No patient row ever leaves a server.

9.1 Assembling the modelling frame

ds.glm needs numeric responses, so we cast the boolean features to 0/1 and bind them — together with the harmonised sex factor, the derived age, and the NSAID burden — into one frame M:

ds.asNumeric("D$gi_bleed",     "gi_bleed")
ds.asNumeric("D$diclofenac",   "diclofenac")
ds.asNumeric("D$peptic_ulcer", "peptic_ulcer")

ds.dataFrame(
  c("gi_bleed", "D$sex", "age", "diclofenac", "D$nsaid_n", "peptic_ulcer"),
  newobj = "M")
ds.colnames("M")[[1]]
#> [1] "gi_bleed"     "sex"          "age"          "diclofenac"   "nsaid_n"     
#> [6] "peptic_ulcer"

9.2 Fitting the model (with its live training trace)

ds.glm fits by iteratively re-weighted least squares, exchanging only aggregated sufficient statistics each round. The federation reports its progress every iteration — this is the genuine convergence trace from the live run:

fit <- ds.glm(
  formula = "M$gi_bleed ~ M$sex + M$age + M$diclofenac + M$nsaid_n + M$peptic_ulcer",
  family  = "binomial", datasources = conns)
#> Iteration 1...
#> CURRENT DEVIANCE:      3734.67700885699
#> Iteration 2...
#> CURRENT DEVIANCE:      2402.84773061472
#> Iteration 3...
#> CURRENT DEVIANCE:      2344.91106777003
#> Iteration 4...
#> CURRENT DEVIANCE:      2343.48010937631
#> Iteration 5...
#> CURRENT DEVIANCE:      2343.47859539296
#> Iteration 6...
#> CURRENT DEVIANCE:      2343.47859539105
#> SUMMARY OF MODEL STATE after iteration 6
#> Current deviance 2343.47859539105 on 2688 degrees of freedom
#> Convergence criterion TRUE (8.15548972641514e-13)
#> 
#> beta: -2.17298353130696 0.00052672527275702 0.00674726691578815 -0.00110142982525188 -0.0859587376367709 1.37487667184113
#> 
#> Information matrix overall:
#>                (Intercept)   M$sexmale       M$age M$diclofenac  M$nsaid_n M$peptic_ulcer
#> (Intercept)       365.6699   180.46686   22734.622    100.45630  1322.2588      177.00793
#> M$sexmale         180.4669   180.46686   11142.069     48.54682   642.6945       88.60275
#> M$age           22734.6217 11142.06909 1508111.822   6286.98725 88291.8068    11028.42073
#> M$diclofenac      100.4563    48.54682    6286.987    100.45630   371.1896       29.09972
#> M$nsaid_n        1322.2588   642.69454   88291.807    371.18961  6102.6082      657.89954
#> M$peptic_ulcer    177.0079    88.60275   11028.421     29.09972   657.8995      177.00793
#> 
#> Score vector overall:
#>                         [,1]
#> (Intercept)    -7.331602e-10
#> M$sexmale      -3.561182e-10
#> M$age          -4.539100e-08
#> M$diclofenac   -2.680922e-10
#> M$nsaid_n      -3.194190e-09
#> M$peptic_ulcer -1.896350e-11
#> 
#> Current deviance: 2343.47859539105

The model converged in 6 iterations over the full 2694-person cohort. Here are the pooled coefficients, with standard errors, p-values and odds ratios:

round(fit$coefficients, 4)
#>                Estimate Std. Error z-value p-value low0.95CI.LP high0.95CI.LP   P_OR
#> (Intercept)     -2.1730     0.2263 -9.6025  0.0000      -2.6165       -1.7295 0.1022
#> M$sexmale        0.0005     0.1047  0.0050  0.9960      -0.2046        0.2057 1.0005
#> M$age            0.0067     0.0039  1.7410  0.0817      -0.0008        0.0143 1.0068
#> M$diclofenac    -0.0011     0.1208 -0.0091  0.9927      -0.2378        0.2356 0.9989
#> M$nsaid_n       -0.0860     0.0329 -2.6157  0.0089      -0.1504       -0.0215 0.9176
#> M$peptic_ulcer   1.3749     0.1080 12.7298  0.0000       1.1632        1.5866 3.9546
#>                low0.95CI.P_OR high0.95CI.P_OR
#> (Intercept)            0.0681          0.1507
#> M$sexmale              0.8150          1.2283
#> M$age                  0.9992          1.0144
#> M$diclofenac           0.7884          1.2656
#> M$nsaid_n              0.8604          0.9787
#> M$peptic_ulcer         3.2001          4.8869

9.3 Reading the result

The headline is unambiguous and clinically coherent: a prior history of peptic ulcer multiplies the odds of gastrointestinal haemorrhage by about 3.95 (p ≈ 4.0e-37), adjusting for everything else. NSAID burden (nsaid_n) carries a smaller but significant signal (p ≈ 0.0089), while age, sex and the diclofenac flag are not individually associated with the outcome in this cohort — exactly the kind of mixed, honest result a real analysis produces.

That diclofenac and the any-NSAID flag add little is itself the story from Chapter 8: this is an NSAID-saturated cohort, so it is the amount of exposure and the prior ulcer, not NSAID use per se, that move the outcome.

9.4 The same model, study-level

ds.glm above pools individual-level sufficient statistics — one model over the combined federation. ds.glmSLMA instead fits the model separately on each site and meta-combines the per-site estimates (study-level meta-analysis), sharing no individual-level statistic across sites. The two routes agree closely, which is the reassurance that the pooled fit is not driven by any single site:

slma <- ds.glmSLMA(
  formula = "M$gi_bleed ~ M$sex + M$age + M$diclofenac + M$nsaid_n + M$peptic_ulcer",
  family  = "binomial", datasources = conns)
round(slma$SLMA.pooled.ests.matrix, 4)
#>                pooled.ML  se.ML pooled.REML se.REML pooled.FE  se.FE
#> (Intercept)      -2.1766 0.2277     -2.1766  0.2277   -2.1766 0.2277
#> M$sexmale         0.0048 0.1050      0.0048  0.1050    0.0048 0.1050
#> M$age             0.0069 0.0046      0.0068  0.0057    0.0071 0.0039
#> M$diclofenac     -0.0004 0.1212     -0.0004  0.1212   -0.0004 0.1212
#> M$nsaid_n        -0.0887 0.0332     -0.0887  0.0362   -0.0887 0.0332
#> M$peptic_ulcer    1.3815 0.1085      1.3815  0.1085    1.3815 0.1085

9.5 Disconnect

The analysis is complete, so we release the OMOP handles and log out of the federation:

ds.omop.disconnect(symbol = "omop", conns = conns)
DSI::datashield.logout(conns)

9.6 In summary

Starting from nothing but three independent OMOP databases, we explored their catalogues, fused three tables into one person-level frame on each server, described it with federated summaries and plots, and fitted a federated logistic regression — recovering a clinically sensible result over the full 2 694-person cohort while no patient-level data ever left a site. The OMOP layer produced and harmonised the analysis dataset; the ordinary DataSHIELD toolchain took it from there.