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

Question

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$?

Answer 1

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

guys=c("fred","bob","andrew","joe");
## [1] "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
))
## [1] 145
sum(abs(residuals(fit)))
## [1] 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
))
## [1] 2.842171e-14
sum(abs(residuals(fit_int)))
## [1] 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!

Answer 2

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.