I was trying to fit a 2-level "hierarchical model" all in one go, in MATLAB. But then realised it might be better to do the lower level first, then the higher level.

Simply, I have 80 subjects, from 2 populations, each of which performed 100 trials, with a continuous predictor (dependent variable) over the 100 trials. The predictors are matched across subjects. SUBJECT is therefore a random effect nested within GROUP.

I am interested in the effect of group on the outcome, and the interaction of group with the continuous predictor -- does one group have a higher slope than the other?

So I have the design matrix

 0     1 0 0 0 0 ...   0.2    0
 0     1 0 0 0 0 ...   0.6    1 
 0     1 0 0 0 0 ...   0.3    0
 0     1 0 0 0 0 ...   0.7    1
 0     0 1 0 0 0 ...   0.4    1
 0     0 1 0 0 0 ...   0.5    0
 0     ...
 1     0 0 1 0 0 ...   0.3    0


Then I used nlmefit as follows

nlmefit(X,y,subj, [], @(psi, x) 1./(1+exp(x*psi'))  , zeros(1,size(X,2)) )

Because it's trying to fit >100 parameters (each subject's interecpt), it didn't finish after running for 24 hours!

So I reasoned that a simple approach is to to a logistic regression for each subject, to get an intercept and slope. I used glmfit for each subject:

glmfit( X, y, 'binomial')

That gives me

  0      1     ...     ...        ...             ...
  0      2     ...     ...
  1      3     ...     ...


And by dividing the two fitted parameters by their estimated standard errors, I get a t-statistic for the intercepts and slopes.

So the question is, is it valid to compare the t-statistics between the two groups as follows?

ttest2( T_SLOPE(GROUP==0), T_SLOPE(GROUP==1) )

Your approach looks very odd to me. I would go at it with a standard hierarchical linear model, where you have a random intercept and slope. I do not know Matlab, but the equations would be:




Where i is the index of the observation in j individual. I'm treating "group" here as an individual factor, though if you considered group random, you could have 3 levels, with group as the third, though with only 2 levels I don't think it will get you very far.

You can use a random effects regression model (also called mixed effect model or hierarchical linear model), where the intercept $\beta_0$ is assumed to follow a random normal distribution around the expected value of the intercept $\gamma_{01}$, the continuous predictor $X_j$, and the grouping variable, with an error term $\mu$. Similarly, the slope is dependent on the same variables, and has it's own error term. You can specify the error terms $\mu_0$ and $\mu_1$ to be correlated or independent. You can also test for whether the variables in both the intercept and the slope are significant, and whether a random slope is significant and substantive. If group is significant in the slope equation you can take the coefficient $\gamma_{32}$ as the mean difference between the two groups in slope, after accounting for the predictor and the repeated measures structure of the data.

This very common approach does not require using up all the degrees of freedom of the interaction terms, and also does not give you a million interaction terms to interpret.


This looks valid in the sense that the samples are independent. But have you tried the fitglme function? It looks like OUTC is a binary outcome, PRED and GROUP are fixed effects, and SUBJ is a nested random effect.


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.