I would like to perform a regression analysis on a dataset comprising one independent variable (X) and two dependent variables (Y1 and Y2) which may be affected by correlated errors. R's stats::lm function handles this situation nicely, producing a 4x4 covariance matrix of the regression parameters (slopes and intercepts). I noticed that the correlation coefficient of the intercepts equals the correlation coefficient of the residuals. Is this true for all other fitting models, including Generalised Linear Models? Am I correct in thinking that:
- The regression coefficients for the dependent variable may be obtained by fitting each of them separately?
- The covariance/correlation matrix of the intercepts may always be obtained from the residuals?
Here's an example:
library(MASS) a1 <- 5 # intercept 1 a2 <- 10 # intercept 2 b1 <- 1 # slope 1 b2 <- 2 # slope 2 x <- 1:10 # independent variable e <- mvrnorm(n=length(x),mu=c(0,0),Sigma=matrix(c(1,0.5,0.5,1),nrow=2)) # residuals # generating some synthetic data: y1 <- a1 + b1 * x + e[,1] # dependent variable 1 y2 <- a2 + b2 * x + e[,2] # dependent variable 2 y <- cbind(y1,y2) # linear regression fit <- lm(y ~ x)
Which produces the following output:
> fit Call: lm(formula = y ~ x) Coefficients: y1 y2 (Intercept) 4.603 9.591 x 1.095 2.066 > vcov(fit) y1:(Intercept) y1:x y2:(Intercept) y2:x y1:(Intercept) 0.39622492 -0.056603560 0.35146659 -0.050209513 y1:x -0.05660356 0.010291556 -0.05020951 0.009129002 y2:(Intercept) 0.35146659 -0.050209513 0.44698258 -0.063854654 y2:x -0.05020951 0.009129002 -0.06385465 0.011609937
You can see that
stats::lm has correctly retrieved both the regression coefficients and the covariance matrix of the regression coefficients. However,
stats::glm does not accept multiple dependent variables. So I was wondering if I could calculate their covariance structure from the residuals as well.