How to convert a Gamma GLM to a Pearl's Structural Causal Model (SCM) with exogenous noise terms? Assume both X and Y are univariate for this. So in Judea Pearl's SCM, you have endogenous variables and exogenous noise. A regular causal Bayesian Network like X->Y where X and Y|X is normal can be represented as both a PGM and an SCM as follows:
PGM representation:
$X \sim ~ N(\mu_X,\sigma_X)$
$Y \sim ~ N(\beta_0 + \beta_1*X,\sigma_Y)$
SCM representation:
$\epsilon_X \sim N(0,\sigma_X)$
$ X \sim~ Dirac(\mu_X + \epsilon_X) $
$\epsilon_Y \sim N(0,\sigma_Y)$
$Y \sim ~ Dirac(\beta_0 + \beta_1*X + \epsilon_Y,\sigma_Y)$
Essentially, the mapping Y=f(X) is now deterministic, and the noise is separated. This enables one to answer counterfactual queries in Pearl's framework where by inferring the noise, it is possible to see what would have happened at an individual level if X were set to something else.
I'm wondering, how can one do the above if Y|X is say Gamma distributed or some other GLM distribution?
In the case of logistic regression, it seems possible following the latent variable formulation here https://en.wikipedia.org/wiki/Logistic_regression#As_a_latent-variable_model where it is possible to simulate the error exogenously. But I am wondering how do it it for a general case like the Gamma.
The hard part seems to be that with a Gamma, the noise is not independent of the mean, so I'm not sure how to separate it. Does that imply SCMs cannot capture other distributions?

With GLMs we usually presume that the distribution is determined by the mean. So $\mu = \beta_1 X + \beta_0$ determines some Gamma distribution $g_\mu(y)$.
Next, consider the inverse $G^{-1}_\mu(e)$ of the cumulative function $G_\mu(y)$ of $g_\mu(y)$. Now, use as exogenous variable $\epsilon_Y$ a uniformly over $[0, 1]$ distributed variable $\epsilon_Y\sim U(0, 1)$. Then you obtain as your SCM representation for $Y$:
$$
Y = G^{-1}_\mu(\epsilon_Y),
$$
which utilizes the endogenous variable $X$ through $\mu = \beta_1 X + \beta_0$ and an exogenous random variable $\epsilon_Y$.
