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 $

  1. I'd like to get a p-value for the above hypothesis, and if possible
  2. 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.

  • $\begingroup$ THis might not help but what you need to do is create a 'test statistic' that has an observable or expected behaviour when H0 is true and that is rarely trivial. You could use Monte Carlo simulation and induce errors to observe if your null hypothesis holds true in large number of simulations but that isn't Hypothesis Testing in the traditional sense $\endgroup$ Aug 6 '18 at 16:37
  • $\begingroup$ @GauravTaneja I can construct the test statistic manually but I don't know it's distribution under H_0. You're saying I could use Monte Carlo to arrive at an approximate distribution (histogram). But then I'd need to know the variance etc. for the data to simulate. H_0 doesn't provide information about the variance. GLM/LIMMA etc. estimate variance from the data to get p-values. How do I close the loop? $\endgroup$ Aug 6 '18 at 16:53

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

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.

  • $\begingroup$ I can probably try to research this, but would I sample under H0 (i.e., by randomizing the labels) or under H1? For the first, I get the sample the distribution under H0 and I understand this will give me a p-value for the statistic. But then why not just take N permuted samples (without replacement). If I resample under H1, I'll get the variance of the estimator - and I can see how often it is more than zero. $\endgroup$ Aug 27 '18 at 1:29
  • $\begingroup$ For a resampling test you would resample under $H_0$ (easy when type null hypothesis is a point null hypothesis like two groups don't differ - when you just resample ignoring the group label), but given the complex shape of the hypothesis that seems non-trivial. Bootstrapping (i.e. resampling the observations with replacementis to create multiple datasets that you then analyse) is a lot easier to implement. $\endgroup$
    – Björn
    Aug 27 '18 at 5:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.