# How to do a generalized linear model with multiple dependent variables in R?

I have six dependent variables (count data) and several independent variables, I see that in a MMR the script goes like this:

my.model <- lm(cbind(DV1,DV2,DV3,DV4,DV5,DV6) ~ IV1 + IV2 + ... + IVn)


But, since my data are counts, I want to use a generalized linear model and I tried this:

my.model <- glm(cbind(DV1,DV2,DV3,DV4,DV5,DV6) ~ IV1 + IV2 + ... + IVn, family="poisson")


And appears this error message:

Error in glm.fit(x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  :
(subscript) logical subscript too long


Can anyone explain me this error message or a way to solve my problem?

-
Following up on @Giorgio Spedicato's answer: are we to suppose that you do want a set of models that treat each dependent variable separately, as lm does when you give it a matrix? –  conjugateprior Apr 17 '12 at 22:02
I miss the second part of the analysis. In a MMR (Multivariate Multiple Regression) after: lm(cbind(DV1,DV2,DV3,DV4,DV5,DV6) ~ IV1 + IV2 + ... + IVn) I must –  Juan Apr 18 '12 at 15:22
I might have just adjusted my answer to answer this question. Also, remember not to press return in the comments :-) –  conjugateprior Apr 18 '12 at 15:32
I miss the second part of the analysis. In a MMR (Multivariate Multiple Regression) after: lm(cbind(DV1,DV2,DV3,DV4,DV5,DV6) ~ IV1 + IV2 + ... + IVn) I must use the manova() command like this: summary(manova(my.model)) to do a multivariate analysis of variance and see the significance of each independent variable. That is the final target. –  Juan Apr 18 '12 at 16:27
Neither manova nor anova are defined for this kind of data, hence it's not offered. But if you want to see the impact of each IV then the regression table provided by summary will give them to you for each DV. –  conjugateprior Apr 18 '12 at 17:20

The short answer is that glm doesn't work like that. The lm will create mlm objects if you give it a matrix, but this is not widely supported in the generics and anyway couldn't easily generalize to glm because users need to be able to specify dual column dependent variables for logistic regression models.

The solution is to fit the models separately. Assume your IVs and DVs live in a data.frame called dd and are labelled the way they are in your question. The following code makes a list of fitted models indexed by the name of the dependent variable they use:

models <- list()
dvnames <- paste("DV", 1:6, sep='')
ivnames <- paste("IV", 1:n, sep='') ## for some value of n

for (y in dvnames){
form <- formula(paste(y, "~", ivnames))
models[[y]] <- glm(form, data=dd, family='poisson')
}


To examine the results, just wrap your usual functions in a lapply, like this:

lapply(models, summary) ## summarize each model


There are no doubt more elegant ways to do this in R, but that should work.

-
You bring up an interesting point (+1). Multivariate GLMs certainly exist. On the other hand, giving lm` a matrix for a dependent variable should probably be seen more as syntactic sugar, than as the expression of a multivariate model: if it were a multivariate (normal) model it'd be the one where the errors are 'spherical', i.e. one where you could have run separate regressions on each element of the dependent variable and gotten the same answer. –  conjugateprior Apr 17 '12 at 21:57