From the help page for lm
:
If ‘response’ is a matrix a linear model is fitted separately by least-squares to each column of the matrix.
So it would seem that a simple approach would be to combine all the different y vectors into a matrix and pass that as the response in a single call to lm
. For example:
(fit <- lm( cbind(Sepal.Width,Sepal.Length) ~ Petal.Width+Petal.Length+Species, data=iris))
summary(fit)
summary(fit)[2]
coef(summary(fit)[2])
coef(summary(fit))[2]
sapply( summary(fit), function(x) x$r.squared )