Testing complex hypotheses involving glm/linear model coefficients I have one independent variable (X) and three dependent variables (W, Y and Z). I am fitting a generalized linear model to each:
$Y = g(\alpha X + \epsilon)$
$Z = g(\beta X + \epsilon)$
$W = g(\gamma X + \epsilon)$
I solve the above in R and get estimates for $\alpha, \beta\, \gamma$ along with p-values that they are not zero. 
However, what I am really interested in is a complex hypothesis such as the following:
$ H_0 : |\alpha - \gamma| < |\beta - \gamma|$ and $(\alpha - \gamma)(\beta - \gamma) > 0$
$H_1 : |\alpha - \gamma| > |\beta - \gamma|$ if $(\alpha - \gamma)(\beta - \gamma) > 0 $ or $(\alpha - \gamma)(\beta - \gamma) < 0 $


*

*I'd like to get a p-value for the above hypothesis, and if possible

*an "effect size" for something like the estimate of
$s = (\alpha - \beta)$ if $H_1$, 0 otherwise
How would I go about doing this? If this is too complicated to explain, how would I test any non-trivial hypothesis involving glm coefficients? Could we use the estimates and their distributions? Or could we transform it into an equivalent model whose coefficient has the same p-value?
ps. I am actually solving the glms in LIMMA.
 A: One simple way to do this is bootstrapping. I.e. draw samples with replacement, fit your regression equations for each bootstrap sample, calculate the quantities you are interested in for each bootstrap sample and then obtain hypothesis tests/confidence intervals from this bootstrap distribution for these quantities. Another standard approach for this type of situation is to derive standard errors for the quantities of interest using the delta method and some software would do this for you automatically (e.g. PROC NLMIXED in SAS). If your software does not do this for you automatically, then bootstrapping should require much less effort on your part and might be attractive for that reason alone. For either approach the right transformation of the parameters may matter - especially for small sample sizes. E.g. it would typically be better to work with a log-odds than a probability, even if you are interested in the probability (in which case you would derive a confidence interval/test for the log-odds and back-transform afterwards).
This should also be pretty straightforward to do, if you were willing to be Bayesian and to specify some reasonable vague priors on the individual model parameters. You would get any desired posterior probabilities simply from calculating how often any condition you are interested in holds for the set of MCMC samples you obtain.
A: I guess that this is difficult by using a frequentist method since $H_0$ is not a single value of the parameters $\alpha$, $\beta$ and $\gamma$. You can not compute, unambiguously, the likelihood function $\mathcal{L}(H_0 \vert X,Y,Z,W)$. But what you could do is maximize over all possible values that coincide with $H_0$ and use a likelihood ratio test:
$$\Lambda = \frac{\sup \lbrace \mathcal{L}(\alpha, \beta, \gamma \in H_0 \vert X,Y,Z,W) \rbrace}{\sup \lbrace \mathcal{L}(\alpha, \beta, \gamma \in H_1 \vert X,Y,Z,W) \rbrace}$$ 

For this purpose it may be useful to re-parametrize the variables (if $g$ is monotonic):
$$Y^\star = g^{-1}(Y) = \alpha X + \epsilon_Y $$
$$Z^\star = g^{-1}(Z) = \beta X + \epsilon_Z $$
$$W^\star = g^{-1}(W) = \gamma X + \epsilon_W $$
and further:
$$U = (Y^\star-W^\star) = (\alpha-\gamma)X+\epsilon_Y-\epsilon_W$$
$$V = (Z^\star-Z^\star) = (\beta-\gamma)X+\epsilon_Z-\epsilon_W$$
such that the residuals are distributed as multivariate normal:
$$\lbrace U-(\alpha-\gamma)X , V-(\beta-\gamma)X \rbrace \sim N(0,\Sigma)$$
where you could convert this into a multivariate distribution for $(\alpha-\gamma)$ and $(\beta-\gamma)$.
The hypotheses now coincide with simpler conditions: 


*

*$H_0$: $(\alpha-\gamma) > (\beta-\gamma)$ in the 1st and 3rd quadrant

*$H_1$: $(\alpha-\gamma) < (\beta-\gamma)$ in the 1st and 3rd quadrant plus  $(\alpha-\gamma) > (\beta-\gamma)$ in the 2nd and 4th quadrant

