Выполнение выбора t-теста вручную

Я пытаюсь написать код моделирования, который генерирует данные и выполняет выборку t-теста (отбрасывая те предикторы, p-значение t-критерия которых превышает 0,05, оставляя остальные). Моделирование в значительной степени является адаптацией прикладной эконометрики с R Клейбера и Цейлиса (2008, стр. 183–189).

При запуске кода обычно происходит сбой. Тем не менее, с некоторыми семенами (например, 1534 г.) он дает правдоподобный результат. Если он не производит вывод (например, 1911), он терпит неудачу из-за: "Error in x[, ii] : subscript out of bounds", который восходит к na.omit.data.frame(). Итак, по какой-то причине мой способ обращения с НА кажется неудачным, но я не могу понять, как это сделать.

  coef <- rep(coef[,3], length.out = pdim+1)
  err <- as.vector(rnorm(nobs, sd = sd))
  uX <- c(rep(1, times = nobs))
  pX <- matrix(scale(rnorm(nobs)), byrow = TRUE, ncol = pdim, nrow = nobs)
  X <- cbind(uX, pX)
  y <- coef %*% t(X) + err
  y <- matrix(y)

  tTp <- (summary(lm(y ~ pX)))$coefficients[,4]  
  tTp <- tTp[2:length(tTp)]
  TTT <- matrix(c(tTp, rep(.7, ncol(pX)-length(tTp))))

  tX <- matrix(NA, ncol = ncol(pX), nrow = nrow(pX))
  for(i in 1:ncol(pX)) {ifelse(TTT[i,] < ALPHA, tX[,i] <- pX[,i], NA)}
  tX <- matrix(Filter(function(x)!all(is.na(x)), tX), nrow = nobs)
  TTR <- lm(y ~ tX)

Первый блок вряд ли является причиной ошибки. Он просто генерирует данные и хорошо работает как сам по себе, так и с другими методами, такими как PCA. Второй блок извлекает p-значения из выходных данных регрессии; удаляет p-значение точки пересечения (beta_0); и заполняет вектор таким количеством семерок, которое необходимо, чтобы иметь ту же длину, что и количество переменных, чтобы обеспечить такую ​​же размерность для матричных вычислений. Семь является произвольной и может быть любым числом больше 0,05, чтобы цикл не прошел проверку. Я считаю, что это становится необходимым, если R отбрасывает предикторы из-за мультиколлинеарности.

Последний блок создает пустую матрицу исходных размеров; вставляет исходные данные, если p-значение t-критерия меньше 0,05, иначе сохраняет NA; в то время как предпоследняя строка удаляет все столбцы, содержащие NA ((здесь исключительно NA или одна NA), взятые из ответа mnel на Удалите столбцы из фрейма данных, где ВСЕ значения НЕТ); наконец, измененные данные снова принимают форму линейной регрессии.

Кто-нибудь знает, что вызывает такое поведение или как оно будет работать по назначению? Я ожидал, что это либо сработает, либо нет, но не то и другое одновременно. В идеале первое.

Рабочая версия кода:

set.seed(1534)
Sim_TTS  <- function(nobs = c(1000, 15000), pdim = pdims, coef = coef100, 
    model = c("MLC", "MHC"), ...){
 DGP_TTS <- function(nobs = 1000, model = c("MLC", "MHC"), coef = coef100, 
     sd = 1, pdim = pdims, ALPHA = 0.05)
 {
  model <- match.arg(model)
  if (model == "MLC") {
   coef <- rep(coef[,1], length.out = pdim+1)
   err <- as.vector(rnorm(nobs, sd = sd))
   uX <- c(rep(1, times = nobs))
   pX <- matrix(scale(rnorm(nobs)), byrow = TRUE, ncol = pdim, nrow = nobs)
   X <- cbind(uX, pX)
   y <- coef %*% t(X) + err
   y <- matrix(y)

   tTp <- (summary(lm(y ~ pX)))$coefficients[,4]  
   tTp <- tTp[2:length(tTp)]
   TTT <- matrix(c(tTp, rep(.7, ncol(pX)-length(tTp)))) 

   tX <- matrix(NA, ncol = ncol(pX), nrow = nrow(pX)) 
   for(i in 1:ncol(pX)) {ifelse(TTT[i,] < ALPHA, tX[,i] <- pX[,i], NA)}
   tX <- matrix(Filter(function(x)!all(is.na(x)), tX), nrow = nobs) 
   TTR <- lm(y ~ tX) 
   } else {
   coef <- rep(coef[,2], length.out = pdim+1)
   err <- as.vector(rnorm(nobs, sd = sd))
   uX <- c(rep(1, times = nobs))
   pX <- matrix(scale(rnorm(nobs)), byrow = TRUE, ncol = pdim, nrow = nobs)
   X <- cbind(uX, pX)
   y <- coef %*% t(X) + err
   y <- matrix(y)

   tTp <- (summary(lm(y ~ pX)))$coefficients[,4]  
   tTp <- tTp[2:length(tTp)]
   TTT <- matrix(c(tTp, rep(.7, ncol(pX)-length(tTp))))

   tX <- matrix(NA, ncol = ncol(pX), nrow = nrow(pX))
   for(i in 1:ncol(pX)) {ifelse(TTT[i,] < ALPHA, tX[,i] <- pX[,i], NA)}
   tX <- matrix(Filter(function(x)!all(is.na(x)), tX), nrow = nobs)
   TTR <- lm(y ~ tX)
   }
  return(TTR)
  }
  PG_TTS <- function(nrep = 1, ...)
  {
   rsq <- matrix(rep(NA, nrep), ncol = 1)
   rsqad <- matrix(rep(NA, nrep), ncol = 1)
   pastr <- matrix(rep(NA, nrep), ncol = 1)
   vmat <- cbind(rsq, rsqad, pastr)
   colnames(vmat) <- c("R sq.", "adj. R sq.", "p*")
   for(i in 1:nrep) {
     vmat[i,1] <- summary(DGP_TTS(...))$r.squared
     vmat[i,2] <- summary(DGP_TTS(...))$adj.r.squared
     vmat[i,3] <- length(DGP_TTS(...)$coefficients)-1
     }
   return(c(mean(vmat[,1]), mean(vmat[,2]), round(mean(vmat[,3]))))
  }
  SIM_TTS <- function(...)
  {
   prs <- expand.grid(pdim = pdim, nobs = nobs, model = model)
   nprs <- nrow(prs)

   pow <- matrix(rep(NA, 3 * nprs), ncol = 3)
   for(i in 1:nprs) pow[i,] <- PG_TTS(pdim = prs[i,1],
       nobs = prs[i,2], model = as.character(prs[i,3]), ...)

   rval <- rbind(prs, prs, prs)
   rval$stat <- factor(rep(1:3, c(nprs, nprs, nprs)),
       labels = c("R sq.", "adj. R sq.", "p*"))
   rval$power <- c(pow[,1], pow[,2], pow[,3])
   rval$nobs <- factor(rval$nobs)
   return(rval)
  }

 psim_TTS <- SIM_TTS()
 tab_TTS <- xtabs(power ~ pdim + stat + model + nobs, data = psim_TTS)
 ftable(tab_TTS, row.vars = c("model", "nobs", "stat"), col.vars = "pdim")}

 FO_TTS <- Sim_TTS()
 FO_TTS
}

Предшествовал:

pdims <- seq(12, 100, 4)
coefLC12 <- c(0, rep(0.2, 4), rep(0.1, 4), rep(0, 4))/1.3
rtL <- c(0.2, rep(0, 3))/1.3
coefLC100 <- c(coefLC12, rep(rtL, 22))
coefHC12 <- c(0, rep(0.8, 4), rep(0.4, 4), rep(0, 4))/1.1
rtH <- c(0.8, rep(0, 3))/1.1
coefHC100 <- c(coefHC12, rep(rtH, 22))
coef100 <- cbind(coefLC100, coefHC100)

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