首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >用于多重线性回归的R中嵌套应用而不是双for循环

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

Stack Overflow用户
提问于 2019-03-09 03:34:29
回答 2查看 271关注 0票数 1

我有两个数据矩阵,a和b(具有多个列)和2个协变量矩阵(每个列1列)。我想要应用多元线性回归,分别得到a的每一列与b的因子之间的回归系数。

协变量是c1和c2。

我希望输出如下所示:

代码语言:javascript
复制
        Estimate Std. Error t value Pr(>|t|)
a1 b1    
a1 b2
...

a2 b1

a2 b2
...

a3 b1

a3 b2
...

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

我尝试了这种嵌套应用

代码语言:javascript
复制
apply(a, 2, function(x) apply(b, 2, function(y) summary(lm(y~x+c1+c2))$coefficients)[2,])

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

代码语言:javascript
复制
         a1  a2  a3

b1

b2

我还尝试了这个:

代码语言:javascript
复制
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]')"

这方面的任何帮助都将不胜感激。

EN

回答 2

Stack Overflow用户

发布于 2019-03-09 04:06:38

试试这个:

代码语言:javascript
复制
# 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)) } )
代码语言:javascript
复制
票数 1
EN

Stack Overflow用户

发布于 2019-03-09 03:58:17

我认为诀窍是在应用/映射命令期间将数据矩阵的列作为变量写入。

代码语言:javascript
复制
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
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/55069864

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档