Regression with multiple dependent variables?

Question

Is it possible to have a (multiple) regression equation with two or more dependent variables? Sure, you could run two separate regression equations, one for each DV, but that doesn't seem like it would capture any relationship between the two DVs?

Answer 1

Yes, it is possible. What you're interested is is called "Multivariate Multiple Regression" or just "Multivariate Regression". I don't know what software you are using, but you can do this in R.

Here's a link that provides examples.

Answer 2

@Brett's response is fine.

If you are interested in describing your two-block structure, you could also use PLS regression. Basically, it is a regression framework which relies on the idea of building successive (orthogonal) linear combinations of the variables belonging to each block such that their covariance is maximal. Here we consider that one block $X$ contains explanatory variables, and the other block $Y$ responses variables, as shown below:

We seek "latent variables" who account for a maximum of information (in a linear fashion) included in the $X$ block while allowing to predict the $Y$ block with minimal error. The $u_j$ and $v_j$ are the loadings (i.e., linear combinations) associated to each dimension. The optimization criteria reads

$$ \max_{\mid u_h\mid =1,\mid v_h\mid =1}\text{cov}(X_{h-1}u_h,Yv_h)\quad \big(\equiv \max\text{cov}(\xi_h,\omega_h)\big) $$

where $X_{h-1}$ stands for the deflated (i.e., residualized) $X$ block, after the $h^\text{th}$ regression.

The correlation between factorial scores on the first dimension ($\xi_1$ and $\omega_1$) reflects the magnitude of the $X$-$Y$ link.

Answer 3

Multivariate regression is done in SPSS using the GLM-multivariate option.

Put all your outcomes (DVs) into the outcomes box, but all your continuous predictors into the covariates box. You don't need anything in the factors box. Look at the multivariate tests. The univariate tests will be the same as separate multiple regressions.

As someone else said, you can also specify this as a structural equation model, but the tests are the same.

(Interestingly, well, I think it's interesting, there's a bit of a UK-US difference on this. In the UK, multiple regression is not usually considered a multivariate technique, hence multivariate regression is only multivariate when you have multiple outcomes / DVs.)

Answer 4

I would do this by first transforming the regression variables to PCA calculated variables, and then I would to the regression with the PCA calculated variables. Of course I would store the eigenvectors to be able to calculate the corresponding pca values when I have a new instance I wanna classify.

Answer 5

As mentionned by caracal, you can use mvtnorm package in R. Assuming you made a lm model (named "model") of one of the response in your model, and called it "model", here is how to obtain the multivariate predictive distribution of several response "resp1", "resp2", "resp3" stored in a matrix form Y:

library(mvtnorm)
model = lm(resp1~1+x+x1+x2,datas) 
         # this is only a fake model to get
                                  #the X matrix out of it
Y = as.matrix(datas[, c("resp1", "resp2", "resp3")])
X =  model.matrix(delete.response(terms(model)), 
           data, model$contrasts)
XprimeX  = t(X) %*% X
XprimeXinv = solve(xprimex)
hatB =  xprimexinv %*% t(X) %*% Y
A = t(Y - X%*%hatB)%*% (Y-X%*%hatB)
F = ncol(X)
M = ncol(Y)
N = nrow(Y)
nu= N-(M+F)+1 #nu must be positive
C_1 =  c(1  + x0 %*% xprimexinv %*% t(x0)) 
     # for a prediction of the factor setting x0 
     # (a vector of size F=ncol(X))
varY = A/(nu) 
postmean = x0 %*% hatB
nsim = 2000
ysim = rmvt(n=nsim, delta=postmux0, C_1*varY, df=nu) 

Now, quantiles of ysim are beta-expectation tolerance intervals from the predictive distribution, you can of course directly use the sampled distribution to do whatever you want.

To answer Andrew F., degrees of freedom are hence nu=N-(M+F) + 1 ... N being the # of observations, M the # of responses and F the # of parameters per equation model. nu must be positive.

(You may want to read my work on in this document :-) )

Answer 6

Did you already come across the term "canonical correlation"? There you have sets of variables on the independent as well as on the dependent side. But maybe there are more modern concepts available, the descriptions I have are all of the eighties/nineties...

Answer 7

For Bayesian multivariate regression, one can use R package BNSP. For example, the dataset ami that comes with the package includes 3 responses and 3 covariates.

# First load the package and dataset
require(BNSP)
data(ami)

# Second, centre and scale variables - 
# this is specific to the dataset
sc <- function(x){return((x-mean(x))/sd(x))}
ami$ratio <- sc(log(ami$ami)-log(ami$tot))
ami$tot <- sc(log(ami$tot))
ami$amt <- sc(log(ami$amt)) 
ami$pr <- sc(ami$pr)
ami$qrs <- sc(ami$qrs)
ami$bp <- sc(ami$bp)

##Third, define the mode: on the left of ~ 
# are the 3 responses, separated by |. 
# On the right of ~ is the model for the mean that 
# includes smooth functions sm
# of the 3 covariates: amt, tot, and ratio  

model <- pr | qrs | bp ~ sm(amt, k = 5) + sm(tot, k = 5) + 
           sm(ratio, k = 5) 

#Fourth, fit the model
multiv <- mvrm(formula = model, data = ami, sweeps = 10000, 
            burn = 5000, thin = 2, seed = 1, 
             StorageDir = getwd())

# And last, plot the fitted curves and estimated 
# correlation matrix
plot(multiv, nrow = 3)
plotCorr(multiv)

Results are shown below. For the correlation matrix, the plot of the left shows posterior means and the one on the right posterior credible intervals.

Answer 8

It's called structural equation model or simultaneous equation model.