嵌套应用而不是R中的双重for循环用于多元线性回归

Question

我有两个数据矩阵，a和b（有多个cols）和2个协变量矩阵（每个1 col）。我想应用一个多元线性回归，并得到a的每一列之间的回归系数和b的因子。

协变量是c1和c2。

我希望输出看起来像这样：

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

a2 b1

a2 b2
...

a3 b1

a3 b2
...

线性回归的基本公式是lm（y~x + c1 + c2）

我试过这个嵌套的应用程序

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

但它只给出了以下格式的p值：

         a1  a2  a3

乙1

B2

我也试过这个：

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

这给出了一个错误"variable lengths differ (found for 'a[i]')"

任何有关这方面的帮助将非常感激。

Answer 1

试试这个 ：

# 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)) } )

```

Answer 2

我想诀窍是在apply / map命令期间将数据矩阵的列写为变量。

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