21

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?

conjugateprior
  • 20,883
Juan
  • 399
  • 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
  • See the VGAM package, especially for multiple correlated binomial outcomes (loglinb2 and loglinb3 functions) where running separate glms as in the accepted answer may be erroneous. – dcsuka Dec 29 '22 at 20:15

2 Answers2

15

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.

conjugateprior
  • 20,883
8

I was told Multivariate Generalized Linear (Mixed) Models exists that address your problem. I'm not an expert about it, but I would have a look to SABRE documentation and this book on multivariate GLMs. Maybe they help...

conjugateprior
  • 20,883
Giorgio Spedicato
  • 3,692
  • 2
    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
  • @conjugateprior How is sphericity of errors inferred? Are you referring the result of the Mauchly's Sphericity Test? https://en.wikipedia.org/wiki/Mauchly%27s_sphericity_test – Galen Feb 05 '20 at 17:48
  • @conjugateprior Can you unpack the equivalence between spherical errors and separate regressions a little more? I'm slightly uncertain on what is being proposed, and especially unclear about why the equivalence holds. – Galen Feb 05 '20 at 17:52