МНК регрессия в R

Question

Я действительно борюсь со следующей проблемой, установленной с помощью R , Я хочу смоделировать набор данных с одной зависимой и 20 независимыми переменными, которые обычно являются i.i.d .. Каждая переменная должна иметь 100 наблюдений. (Мне удалось сделать эту часть)

(Теперь часть, с которой я борюсь): Мой план состоит в том, чтобы провести автоматические регрессии для всех возможных комбинаций до 5 регрессоров, используя собственную кодированную регрессионную функцию, которая имитирует вывод итогов (lm), который использует вектор y и матрицу или вектор x в качестве входных данных (поэтому my.lm ( у, х)). И затем приведение результатов в подходящую структуру данных.

Буду благодарен за каждый намек!

Answer 1

Я сомневаюсь в обоснованности того, что вы пытаетесь сделать, но здесь все идет.
Я составлю набор данных, поскольку вы его еще не опубликовали.

my.lm <- function(x, y, n = 5){
  f <- function(inx){
    inx_cols <- Combn[inx, ]
    inx_cols <- inx_cols[inx_cols != 0]
    X <- as.data.frame(x[, inx_cols])
    names(X) <- paste0("X", inx_cols)
    X <- cbind(X, y)
    name_y <- names(X)[length(names(X))]
    fmla <- as.formula(paste(name_y, ".", sep = "~"))
    tryCatch(lm(fmla, data = X), error = function(e) e)
  }

  nc_x <- ncol(x)
  nr <- sum(choose(nc_x, seq_len(n)))
  Combn <- matrix(0, nrow = nr, ncol = n)
  first <- 1
  for(i in seq_len(n)){
    last <- first + choose(nc_x, i) - 1
    Combn[first:last, seq_len(i)] <- t(combn(nc_x, i))
    first <- last + 1
  }

  apply(Combn, 1, f)
}

set.seed(6876)

regr <- replicate(20, rnorm(100))
coefs <- sample(-5:5, 20, TRUE)
resp <- regr %*% coefs + rnorm(100)

lm_list <- my.lm(regr, resp)
length(lm_list)
#[1] 21699

Таким образом, вышеприведенная функция создала столько объектов, сколько ожидалось.
Прежде чем продолжить, давайте посмотрим, сколько ошибок (например, единичная матрица).

err_list <- lapply(lm_list, function(x){
  if(inherits(x, "error")) x else NULL
})
err_list <- err_list[!sapply(err_list, is.null)]
length(err_list)
#[1] 0

Нет ошибок.
Так что получайте аннотации объектов класса "lm".

good_list <- lapply(lm_list, function(x){
  if(inherits(x, "lm")) x else NULL
})
good_list <- good_list[!sapply(good_list, is.null)]

smry_list <- lapply(good_list, summary)
smry_list[[1]]
#
#Call:
#  lm(formula = fmla, data = X)

#Residuals:
#    Min      1Q  Median      3Q     Max 
#-34.654  -9.487  -1.985   9.486  50.213 

#Coefficients:
#                Estimate Std. Error t value Pr(>|t|)    
#(Intercept)       0.6449     1.5237   0.423    0.673    
#X1               -7.3969     1.5074  -4.907 3.68e-06 ***
#  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 15.02 on 98 degrees of freedom
#Multiple R-squared:  0.1972,   Adjusted R-squared:  0.189 
#F-statistic: 24.08 on 1 and 98 DF,  p-value: 3.684e-06