# Regression for multiple scores, obtained from multiple subjects

I have preference ratings (1-to-7 scale) for k=80 stimuli, obtained from N=30 subjects. I would like to do a multiple regression for these scores, against a small number of weakly-correlated predictors. Some of the predictors are continuous, others are categorical. Also, some of the predictors characterise the subject (e.g. age), while others characterise the stimulus (e.g. loudness).

I am confused whether the regression should take in cross-stimulus or cross-subject variability. I will be using the 'regress' function in Matlab:

b = regress(y,X) returns a p-by-1 vector b of coefficient estimates for a multilinear regression of the responses in y on the predictors in X. X is an n-by-p matrix of p predictors at each of n observations. y is an n-by-1 vector of observed responses.

I know the regression line is normally fit between data points that represent different subjects, not different stimuli (repeated measures). In this assumption, the Y matrix (i.e. the different observations) would contain all subjects' scores for a single stimulus. But this would mean I have to do k multiple regressions (one for each stimulus), which seems to me like a mass-univariate apporach, whereas I'd have thought this is a multivariate problem.

My questions are:

1) How should the ratings for the different subjects, and for the different stimuli, be arranged to define the Y and X matrices in the regression formula?

2) Do I need to do a separate regression for each stimulus, and if so, is a multiple comparisons correction necessary?

3) Should multiple regression not in fact be the most optimal tool in this case (e.g. due to the rather low sample size), would mixed-effects modelling be more suitable? If so, the same questions regarding how the cross-subject vs cross-stimulus variability should be arranged in the input and output matrices.

The most obvious approach here is indeed a single mixed model for all the data. Each subject–stimulus pair should get a row in $X$ and a value in $y$, so you'll have $30 \cdot 80 = 2,400$ cases. You can then include a random intercept for each stimulus and a random intercept for each subject. A vignette for the R package lme4, "Fitting Linear Mixed-Effects Models using lme4", says "Such models are common in item response theory, where subject and item factors are fully crossed." (p. 7)