1

我正在尝试编写模拟代码，生成数据并在其上运行 t 检验选择（丢弃那些 t 检验 p 值超过 0.05 的预测变量，保留其余部分）。该模拟很大程度上是对 Kleiber 和 Zeileis (2008, pp. 183–189) 的 Applied Econometrics with R 的改编。

运行代码时，通常会失败。然而，对于某些种子（例如 1534），它会产生合理的输出。如果它不产生输出（例如 1911），它会由于以下原因而失败："Error in x[, ii] : subscript out of bounds"，它追溯到na.omit.data.frame()。因此，出于某种原因，我尝试处理 NA 的方式似乎失败了，但我无法弄清楚是怎么回事。

  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)；并根据需要用尽可能多的 7 填充向量，使其长度与变量数相同，以确保矩阵计算的维度相同。七是任意的，可以是大于 0.05 的任何数字，以不通过循环测试。如果 R 由于多重共线性而丢弃预测变量，这将成为 - 我相信 - 必要的。

最后一个块创建一个原始维度的空矩阵；插入原始数据，如果t检验p值低于0.05，否则保留NA；而倒数第二行删除了所有包含 NA 的列（（此处仅 NA 或一个 NA 相同）取自 mnel 对Remove columns from dataframe where ALL values are NA的回答）；最后，修改后的数据再次以线性回归的形式放置。

有谁知道导致这种行为的原因或它如何按预期工作？我希望它要么工作要么不工作，但不是两者兼而有之。理想情况下，前者。

代码的工作版本是：

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)

我知道不建议通过单个预测变量的重要性来选择模型，但这就是重点——它旨在与更复杂的方法进行比较。

4

0 回答 0