Ggtcf

I have two data matrices, a and b (with multiple cols) and 2 covariate matrices (1 col each). I want to apply a multiple linear regression and get the coefficients for the regression between each column of a with the factors of b, repsectively.

Covariates are c1 and c2.

I want the output to look like this:

 Estimate Std. Error t value Pr(>|t|)
a1 b1 
a1 b2
...

a2 b1

a2 b2
...

a3 b1

a3 b2
...

The basic formula for linear regression is lm(y~x+c1+c2)

I tried this nested apply

apply(a, 2, function(x) apply(b, 2, function(y) summary(lm(y~x+c1+c2))$coefficients)[2,])

but it only gives me the p-values in the following format:

 a1 a2 a3

b1

b2

I also tried this:

for (i in dim(a)[2])
 pvals= apply(b, 2, function(y) summary(lm(y~a[i]+c1+c2))$coefficients)[2,]

This gives an error "variable lengths differ (found for 'a[i]')"

Any help with this would be much appreciated.

Try this :

# transform your data matrices into data.frames
a <- as.data.frame( matrix(rnorm(1:(250*4)), ncol = 4) )
colnames(a) <- paste0("A", 1:ncol(a))
b <- as.data.frame( matrix(rnorm(1:(250*6)), ncol = 6) )
colnames(b) <- paste0("B", 1:ncol(b))
c1 <- rnorm(1:250)
c2 <- rnorm(1:250)

# get the explanatory variables, RHS of the formula
X <- paste(c(colnames(b), "c1", "c2"), collapse = "+")

# get the dependent variables, LHS of the formula
Y <- colnames(a)

# Create a single data.frame
dat <- data.frame(a, b, c1, c2)

# Do the regressions
results <- lapply(Y, function(y) 
 coefficients( lm(
 as.formula( paste0(y, " ~ ", X) ), data=dat)) )

```

I suppose the trick is to write the column of the data matrix as variable during your apply / map command.

library(broom) # to clean the regression output
library(tidyverse)

a <- matrix(rnorm(1:1000), ncol = 4)
head(a)
 [,1] [,2] [,3] [,4]
[1,] 0.9214791 0.3273086 -0.456702485 1.504571891
[2,] -0.6705181 1.3443408 1.496302280 0.516068092
[3,] -0.9122278 0.2392211 -0.163004516 -0.041937414
[4,] -0.6614763 1.1596926 2.004846224 -0.001818212
[5,] -0.7902421 0.3022333 -0.002848944 0.265987941
[6,] 0.3451988 0.3187038 -0.149836811 0.122283166

b <- matrix(rnorm(1:500), ncol = 2)
head(b)
 [,1] [,2]
[1,] 1.6100023 0.4861797
[2,] 0.2128886 -1.0762123
[3,] -0.7645170 -0.4972273
[4,] -0.4084541 0.8930468
[5,] -0.1471686 -1.3193856
[6,] 0.4331506 -0.4044583

c <- matrix(rnorm(1:500), ncol = 2)
head(c)
 [,1] [,2]
[1,] -0.9476932 0.1292495
[2,] -0.8653959 -1.3278809
[3,] -1.5162128 0.2765994
[4,] -0.5140617 1.8684472
[5,] 0.8104582 1.7564293
[6,] 1.4162302 -1.5383332

(col_a <- seq(dim(a)[2])) # to map to the columns of matrix a
[1] 1 2 3 4

(col_b <- seq(dim(b)[2])) # to map to the columns of matrix b
[1] 1 2

map_df(col_a, ~ map2_df(.x, col_b, ~ lm(b[,.y] ~ a[,.x] + c) %>% # the first ".x" uses the mapping output from the first "map_df" in the second "map2_df"
 tidy() %>% # clean regression output
 mutate(y = str_c("b", .y, sep = "_"), # add variable y with indicator for matrix b
 x = str_c("a", .x, sep = "_")))) %>% # add variable x with indicator for matrix a
 select(y, x, 1:5) # rearrange columns
# A tibble: 32 x 7
 y x term estimate std.error statistic p.value
 <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
 1 b_1 a_1 (Intercept) -0.0747 0.0645 -1.16 0.248 
 2 b_1 a_1 a[, .x] 0.0653 0.0638 1.02 0.307 
 3 b_1 a_1 c1 -0.117 0.0672 -1.74 0.0834
 4 b_1 a_1 c2 0.0219 0.0617 0.355 0.723 
 5 b_2 a_1 (Intercept) 0.0145 0.0618 0.234 0.815 
 6 b_2 a_1 a[, .x] -0.142 0.0612 -2.33 0.0208
 7 b_2 a_1 c1 0.0458 0.0644 0.711 0.478 
 8 b_2 a_1 c2 0.0450 0.0591 0.761 0.447 
 9 b_1 a_2 (Intercept) -0.0779 0.0645 -1.21 0.229 
10 b_1 a_2 a[, .x] -0.0502 0.0678 -0.741 0.459 
# ... with 22 more rows