# Using one linear model to test multiple hypotheses

Let’s assume the following:

1. I have data of 100 subjects each performing 100 trials of a reaction time task
2. At each trial, I measured the reaction time (RT) and a specific brain signal (BS)
3. The RT and S data of each subject are Gaussian distributed

Both RTs and BSs vary from trial to trial, but also from subject to subject (i.e. some subjects are on average slower than others; some subjects have on average larger brain signals than others)

I have two hypotheses:

H1: For each subject, RT and S are correlated on a trial-by-trial level.

H2: The average RT of each subject (mRT) and the average BS of each subject (mBS) are correlated.

In other words, H1 tests for an intra-individual effect of S (e.g. “a large signal S is accompanied by small RTs, independent of subject-to-subject differences”), while H2 tests for an inter-individual effect (e.g. “slower subjects have a larger average BS, independent of trial-by-trial variations”).

Ideally, I would like to test both hypotheses simultaneously in one model. However, the problem is that the RT and BS data have 10000 DF, whereas the mRT and mBS data have only 100 DF.

Thus, one “unproblematic” way of testing the hypotheses would be to fit two separate models, e.g.

M1: Linear mixed model with BS as response variable and RT as predictor, allowing for a random intercept and slope for each subject (alternatively: z-transform RTs for each subject and fit a linear model on all trials)

M2: Linear model with mBS as response variable and mRT as predictor

In my data, the predictors in both models are significant. Now I would like to rule out an interaction between inter-individual and intra-individual effects. Or possibly one effect (correlation between RT and BS) is explained by the other. Thus, I would ideally like to fit one single model.

One option would be a linear model without random effects, where BS is the response variable and two predictors: RT (z-transformed for each subject), and mRT (the same repeated value for each trial of one subject). The problem here is that because for each subject the predictor mRT is the same value repeated 100 times, the error measure for the variance of that predictor is presumably wrong, which results in extremely small (and likely wrong) SEs and thus p-values.

Another option would be to fit the same data in a mixed model allowing for random intercepts and slopes for subjects. The model results very well reflect those of the separate models M1 and M2, so I guess this model is correct despite the data having inherently different DFs.

My question is: Is this last model correct?

I think the following model might be appropriate. Let $$R_{it}$$ and $$B_{it}$$ be the reaction time and brain signal, respectively, of the $$i^{th}$$ subject on the $$t^{th}$$ trial. Assuming that these responses are approximately normally distributed (or can be transformed so as to be approximately normally distributed), you could assume a conditional bivariate normal distribution for $$(R_{it},B_{it})$$ with $$\left[ R_{it} \atop B_{it} \right] \sim BVN\left( \left[ \mu_R +U_i \atop \mu_B +V_i \right], \Lambda \right) ,$$ where $$\Lambda=\left[ \begin{array}{cc} \sigma_R^2 & c_{RB} \\ c_{RB} & \sigma_B^2 \end{array} \right]$$ and $$(U_i,V_i)$$ are iid bivariate normal random effects with mean $$\left[ 0 \atop 0 \right]$$ and variance-covariance matrix $$\Sigma=\left[ \begin{array}{cc} \sigma_U^2 & c_{UV} \\ c_{UV} & \sigma_V^2 \end{array} \right]$$.
H1 is then equivalent to the test of $$H_0: c_{RB}=0$$, while H2 is equivalent to the test of $$H_0: c_{UV}=0$$.
I don't know whether standard linear mixed model software will allow you to specify this model. But coding the likelihood yourself using the mvtnorm::mvnorm function and then maximizing it using nlm should be straightforward, I think. You can then do likelihood ratio tests to test your hypotheses.