# Which logit or probit model should I use for multiple response / dependent variables?

I have $300$ time series objects that constitute the $300$ columns of matrix $X$. This matrix has $5$ rows and represents $5$ days of time series information for each $300$ columns.

I set up a $300\textrm{x}5$ matrix of binary values. So the first row might look something like $(1 1 0 0 1)$, which would mean that column $1$ of $X$ had negative elements in the 1st, 2nd and 5th rows (modeled by category "$1$"), and positive elements in the 3rd and 4th rows (modeled by category "$0$").

I have $8$ predictor variables that explains/relates well to whether there are more negative elements early (i.e. if we observe something like $(1 1 0 0 0)$) versus whether there are more negative elements later (i.e. we observe something like $(0 0 0 1 1)$).

How do I model this? My knowledge is not expansive enough to know if there's any statistical framework where I can use this matrix as my response:

$Y = \left( \begin{array}{ccccc} 1 & 1 & 0 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 1 \\ . & . & . & . & . \\ \end{array} \right)$

Using this I want to do something like:

$\mathbf{Y} = a + \mathbf{B}\mathbf{Z}$ (I understand that this is not a correct logit/probit specification, but I think it's fine for communicating what I want).

Then the ideal interpretation of some coefficient $b_i \in \mathbf{B}$ would be; if it's negative, then we have increased odds of seeing something like $(11000)$ instead of something like $(00011)$ if the associated predictor $x_i \in \mathbf{Z}$ is larger.

• Please clarify: Your question uses "X", but it seems there are two distinct X matrices. The first X is a 5x300 matrix mentioned in the first paragraph. But the second X, in your regression pseudo-equation, is a matrix of the 8 predictor values, different from the first X, correct? (Otherwise, the question seems kind of screwy.) – pteetor Oct 13 '12 at 23:03

You have two decent suggestions, but I don't think either of them is optimal. If you turn your individual, five-element vectors into a single, ordinal scalar, you will lose information. This can be acceptable if it's necessary, but if a better way is possible, you may want to avoid it. Multivariate generalized linear models treat your response as a singular point in a multidimensional space, rather than five points ordered in time. Multivariate methods are typically used / understood for cases where you have five different kinds of measurements (here binary) that are all related to each other, but I gather you have a sequence of five instances of the same kind of measurement. It would be better to fit a model that is designed for that.

Fortunately, there are models that are designed exactly for this type of situation. You will want to use a Generalized Linear Mixed effects Model or use the Generalized Estimating Equations. Which you should choose depends on the question you want to ask, GLiMMs provide information on the effects of the covariates for the individual study units, whereas the GEE provides information on the effects of the covariates for the population average. There are several threads on CV that discuss these:

Regarding whether to use the logit link or the probit link, I discussed that fairly extensively here: Difference between logit and probit models. (Actually, the answer there is a little more fundamental in nature, so it may be worth reading that one first.)

Actually, I don't believe that either logit or probit regression is needed here. First, I would reduce the Y matrix to a simple 300x1 column vector of scores. This R code, for example, will reduce each each row of Y to a number between -3 and +3, where larger values correspond to "more negative, later":

f <- function(r) sum(r * c(-2, -1, 0, +1, +2))
Z <- apply(Y, 1, f)


Then, use linear regression to model those scores based on your predictors.

model <- lm.fit(X, Z)


(Here, X is your 300x8 matrix of predictor values, not the 5x300 matrix mentioned in your first paragraph.) The coefficients of the regression will have the interpretation you desire: larger values indicate stronger odds of "more negative, later".

If you really prefer the logistic model, the R code becomes

model <- glm.fix(X, Z, family=binomial())


The question for you is simply, which model works better for your application. The application does not strike me as intrinsically categorical; rather, you constructed the Y matrix to be categorical.

• Nice answer! Your latter suggestion is ordinal logit right? However I heard that there was some support for matrix response variables (clustered logit or something) which is why I asked what I did. – user14281 Oct 14 '12 at 10:00
• Thanks. Yes, ordinal logit regression. So perhaps the polr function of the MASS package is an alternative. But, again, I don't think you have a categorical problem here. – pteetor Oct 14 '12 at 19:46

You could try multivariate generalized linear models, if you wish to follow a regression approach. See SABRE package and http://www.amazon.com/Multivariate-Generalized-Linear-Mixed-Models/dp/1439813264.