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
GROUP SUBJ PRED OUTC
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
etc.
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
GROUP SUBJ INTERCEPT SLOPE SE(INTERCEPT) SE(SLOPE)
0 1 ... ... ... ...
0 2 ... ...
...
1 3 ... ...
etc
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?
T_INTERCEPT = INTERCEPT / SE_INTERCEPT
T_SLOPE = SLOPE / SE_SLOPE
ttest2( T_INTERCEPT(GROUP==0), T_INTERCEPT(GROUP==1) )
ttest2( T_SLOPE(GROUP==0), T_SLOPE(GROUP==1) )