Yes, I do know that the underlying data is non-parametric. I don't understand why that make medians derived from the data non-parametric, as well. Consider the classic `flights`

database in the `nycflights13`

library. Let `arr_delay`

represent the difference in medians and arbitrarily model it on dep_time + distance.

```
library(nycflights13)
library(tidyverse)
#> Registered S3 methods overwritten by 'ggplot2':
#> method from
#> [.quosures rlang
#> c.quosures rlang
#> print.quosures rlang
data(flights)
delta <- flights %>% select(arr_delay, dep_time, distance)
delta.fit <- lm(arr_delay ~ dep_time + distance, data = delta)
summary(delta.fit)
#>
#> Call:
#> lm(formula = arr_delay ~ dep_time + distance, data = delta)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -95.25 -23.31 -8.88 9.21 1294.06
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -1.792e+01 2.486e-01 -72.07 <2e-16 ***
#> dep_time 2.116e-02 1.551e-04 136.40 <2e-16 ***
#> distance -3.554e-03 1.029e-04 -34.53 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 43.33 on 327343 degrees of freedom
#> (9430 observations deleted due to missingness)
#> Multiple R-squared: 0.0574, Adjusted R-squared: 0.05739
#> F-statistic: 9967 on 2 and 327343 DF, p-value: < 2.2e-16
confint(delta.fit)
#> 2.5 % 97.5 %
#> (Intercept) -18.402990518 -17.42857004
#> dep_time 0.020853423 0.02146147
#> distance -0.003755629 -0.00335215
```

^{Created on 2019-06-13 by the reprex package (v0.3.0)}

If you `plot(delta.fit)`

, you'll see a Q-Q plot that's half normal and then climbs to its cruising altitude like a jet.

My point is that to get confidence intervals on the differences, you have to have some idea of how far your point departs from your test statistic.

For what that test statistic should be, I have to invoke Wittgenstein

Whereof one cannot speak, thereof one must be silent.