# What is the proper way to do vector based linear regression in R

I want to do linear regression between vector inputs and vector output. That is each $y$ is a vector with $M$ components, and each $x$ is a vector with $N$ components and the answer should look like $y \sim Ax + b$ where $A$ is an $M \times N$ matrix and $b$ is a vector with $M$ components.

I have a very clear understanding of the concept and what I want R to do, but it is the proper syntax I am lacking.

Trying to google around to find this has been quite difficult because terms like multivariable seem to always point me to answers of the form

$$y \sim x_1 + x_2 + x_3 + \dots + x_n$$

where there are multiple input sources (or rather, a multidimensional input), but never with multidimensional outputs.

If I just feed in matrices for $y$ and $x$ that MIGHT give what I want, but it might also just treat each $y$ component as directly related to each $x$ component and give an answer based on that ($M = N$ for the important instance I have). So I have to be sure that I am doing it correctly.

What is the correct means for using R to do linear regression of the sort

$$y \sim A x + b$$

where the solution $A$ is an $M \times N$ matrix, and $b$ is a vector of length $M$, and each datum $x$ is a vector of length $N$ and each corresponding datum $y$ is a vector of length $M$?

• I don't think you need help choosing an R function, I think you need assistance choosing a statistical method. If you have multiple response per individual, there are many ways you can model that, but you need to decide what model is right for you. A simple linear regression is probably not the right choice. If you need help choosing a statistical model, consider posting to Cross Validated instead as such matters are off-topic for Stack Overflow. Jan 30, 2015 at 0:59
• That doesn't make sense. Simple linear regression assuming a single response value and essentially does a least-squares fit on the regressors. Do you want to independently model each value in the response vector? What's shared here? What's the benefit of modeling jointly? Jan 30, 2015 at 1:06
• There is strong physical cause to believe the real relationship in my case is y = Ax + b plus some noise. The matrix A is important as it is related to a rotation in the system. I really want the matrix A. It is crucial information for my purpose. Jan 30, 2015 at 1:22
• This is nearly closed but I fear not on a statistically correct basis. In 2007 Peter Dalgaard wrote an extremely useful summary of the R functions that support exactly what you request. Here is the first google hit on a search for "Dalgaard multivariate rnews" but It's not really the Ur-repository for all the Rnews articles: lib.stat.cmu.edu/R/CRAN/doc/Rnews/Rnews_2007-2.pdf
– DWin
Jan 30, 2015 at 3:24
• This may be relevant: stats.stackexchange.com/questions/11127/… Nov 17, 2016 at 6:28

Set up the problem. We have 4 "data" items

guys=c("fred","bob","andrew","joe");
##  "fred"   "bob"    "andrew" "joe"


Each row $$x$$ of $$X$$ is a vector with $$N=3$$ components

# By default, matrices are filled by column
X=matrix(c(1,0,1,0,  1,0,0,1,  1,1,1,0), ncol=3,
dimnames=list(guys,c("tea","coffee","yoga")))
##        tea coffee yoga
## fred     1      1    1
## bob      0      0    1
## andrew   1      0    1
## joe      0      1    0


Each row $$y$$ of $$Y$$ is a vector with $$M=2$$ components

Y=matrix(c(185,209,138,167,  36, 32, 30, 31), ncol=2,
dimnames=list(guys, c("weight","girth")))
##        weight girth
## fred      185    36
## bob       209    32
## andrew    138    30
## joe       167    31


This is how we can fit a linear model $$Y=X\times A$$ with no intercept

fit = lm(Y~X+0)
## ...
## Coefficients:
##          weight   girth
## Xtea     -111.00   -10.33
## Xcoffee   127.00    22.67
## Xyoga     209.00    32.00


It is the same as calling lm.fit(X,Y)

fit_ = lm.fit(X,Y)
sum(abs(coef(fit)-coef(fit_)))
## 0


except that you get a decent "summary()" etc.

summary(fit)
## ...
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)
## Xtea      -10.33      18.63  -0.555    0.678
## Xcoffee    22.67      11.79   1.923    0.305
## Xyoga      32.00      14.43   2.217    0.270


Check that the coefficient matrix has the intended meaning of $$A$$ in $$Y=X\times A$$

sum(abs(
X %*% coef(fit) - Y
))
##  145
sum(abs(residuals(fit)))
##  145


Here $$A$$ is an $$N \times M$$ matrix, while you asked for $$M \times N$$. If you really wanted $$Y=A \times X$$ instead of $$Y=X \times A$$, then take the transpose of everything. Each "guy" will then be a column rather than a row. But usually in R we think of observations being stored in rows, that's why I did it this way.

Here is the version with an intercept

fit_int = lm(Y~X+1)     # same as default: lm(Y~X)
##              weight  girth
## (Intercept)  120      25
## Xtea         -71      -2
## Xcoffee       47       6
## Xyoga         89       7


The intercept is part of the "coef()" matrix, so instead of $$Y=X\times A+b$$ we have to modify $$X$$, it's like $$Y=[1,X]\times A$$:

sum(abs(
cbind(1,X) %*% coef(fit_int) - Y
))
##  2.842171e-14
sum(abs(residuals(fit_int)))
##  0   # same as above, more or less


In my setup, after introducing an intercept term, the number of unknowns now matches the number of constraints and the equation has an exact solution as you can see.

I added row and column names for clarity, unfortunately R does not check that these match when doing matrix arithmetic or when creating a formula. So be careful you don't get rows and columns mixed up, or columns out of order. The following produces no error message for example:

t(cbind(X,1)) %*% coef(fit_int)  # wrong!


I suggest that your question is somewhat statistically naive in not actually describing the uses to which this operation might be employed. Matrix outcomes can be of various sorts. You ask for "correct" version of code and then do not offer descriptions of the measurements or their expected correlations. You do hint in your comment that you expect the outcome to be a "rotation" so an example using a within subjects design would might not be directly applicable to your (highly underspecified) problem.

Nonetheless, one can offer examples of code that does not throw errors which might be one "correctness criterion". Using Dalgaard's example where subjects respond to stimuli presented at different angles, the code below

 reacttime <- matrix(c(
420, 420, 480, 480, 600, 780,
420, 480, 480, 360, 480, 600,
480, 480, 540, 660, 780, 780,
420, 540, 540, 480, 780, 900,
540, 660, 540, 480, 660, 720,
360, 420, 360, 360, 480, 540,
480, 480, 600, 540, 720, 840,
480, 600, 660, 540, 720, 900,
540, 600, 540, 480, 720, 780,
480, 420, 540, 540, 660, 780),
ncol = 6, byrow = TRUE,
dimnames=list(subj=1:10,
cond=c("deg0NA", "deg4NA", "deg8NA",
"deg0NP", "deg4NP", "deg8NP")))

mlmfit <- lm(reacttime~1)
estVar(mlmfit)
# produces the variance covariance matrix
mlmfit0 <- update(mlmfit, ~0)
anova(mlmfit, mlmfit0, X=~1)

#----------
Analysis of Variance Table

Model 1: reacttime ~ 1
Model 2: reacttime ~ 1 - 1

Contrasts orthogonal to
~1

Res.Df Df Gen.var. Pillai approx F num Df den Df   Pr(>F)
1      9      1249.6
2     10  1   2013.2 0.9456   17.381      5      5 0.003534 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


This is a test for independence of the response within subjects. Dalgaard then goes into considerable detail in refinements and alternatives to such analysis.