10

我正在尝试模拟一个三变量数据集，以便可以在其上运行线性回归模型。“X1”和“X2”将是连续自变量（均值=0，sd=1），“Y”将是连续因变量。

变量将是回归模型将产生如下系数：Y = 5 + 3(X1) - 2(X2)

我想模拟这个数据集，使得得到的回归模型的 R 平方值为 0.2。如何确定“sd.value”的值，以便回归模型具有这个 R 平方？

n <- 200 
set.seed(101) 
sd.value <- 1

X1 <- rnorm(n, 0, 1)
X2 <- rnorm(n, 0, 1)
Y <- rnorm(n, (5 + 3*X1 - 2*X2), sd.value)

simdata <- data.frame(X1, X2, Y)

summary(lm(Y ~ X1 + X2, data=simdata))
4

4 回答 4

8

看看这段代码，它应该足够接近你想要的：

simulate <- function(n.obs=10^4, beta=c(5, 3, -2), R.sq=0.8) {
    stopifnot(length(beta) == 3)
    df <- data.frame(x1=rnorm(n.obs), x2=rnorm(n.obs))  # x1 and x2 are independent
    var.epsilon <- (beta[2]^2 + beta[3]^2) * (1 - R.sq) / R.sq
    stopifnot(var.epsilon > 0)
    df$epsilon <- rnorm(n.obs, sd=sqrt(var.epsilon))
    df$y <- with(df, beta[1] + beta[2]*x1 + beta[3]*x2 + epsilon)
    return(df)
}
get.R.sq <- function(desired) {
    model <- lm(y ~ x1 + x2, data=simulate(R.sq=desired))
    return(summary(model)$r.squared)
}
df <- data.frame(desired.R.sq=seq(from=0.05, to=0.95, by=0.05))
df$actual.R.sq <- sapply(df$desired.R.sq, FUN=get.R.sq)
plot(df)
abline(a=0, b=1, col="red", lty=2)

基本上，您的问题归结为找出 var.epsilon 的表达式。因为我们有 y = b1 + b2*x1 + b3*x2 + epsilon，并且 Xs 和 epsilon 都是独立的，所以我们有 var[y] = b2^2 * var[x1] + b3^2 * var[x2] + var[eps]，假设 var[Xs]=1。然后，您可以将 var[eps] 作为 R 平方的函数求解。

于 2013-09-30T15:01:03.163 回答
3

所以 R^2 的公式是 1-var(residual)/var(total)

在这种情况下， 的方差Y将是3^2+2^2+sd.value^2，因为我们要添加三个独立的随机变量。而且，渐近地，剩余方差将是简单的sd.value^2。

因此，您可以使用此函数显式计算 rsquared：

rsq<-function(x){1-x^2/(9+ 4+x^2)}

用一点代数，你可以计算这个函数的逆：

rsqi<-function(x){sqrt(13)*sqrt((1-x)/x)}

所以设置sd.value<-rsqi(rsquared)应该给你你想要的。

我们可以这样测试：

simrsq<-function(x){
  Y <- rnorm(n, (5 + 3*X1 - 2*X2), rsqi(x))
  simdata <- data.frame(X1, X2, Y)
  summary(lm(Y ~ X1 + X2, data=simdata))$r.squared
}

> meanrsq<-rep(0,9)
> for(i in 1:50)
+   meanrsq<-meanrsq+Vectorize(simrsq)((1:9)/10)
> meanrsq/50
[1] 0.1031827 0.2075984 0.3063701 0.3977051 0.5052408 0.6024988 0.6947790
[8] 0.7999349 0.8977187

所以看起来是正确的。

于 2013-09-30T15:03:20.400 回答
2

这就是我的做法（盲迭代算法，假设不知道，因为当你纯粹对“如何模拟这个”感兴趣时）：

simulate.sd <- function(nsim=10, n=200, seed=101, tol=0.01) {
  set.seed(seed)
  sd.value <- 1
  rsquare <- 1:nsim
  results <- 1:nsim
  for (i in 1:nsim) {
    # tracking iteration: if we miss the value, abort at sd.value > 7.
    iter <- 0 
    while (rsquare[i] > (0.20 + tol) | rsquare[i] < (0.2 - tol)) {
      sd.value <- sd.value + 0.01
      rsquare[i] <- simulate.sd.iter(sd.value, n)
      iter <- iter + 1
      if (iter > 3000) { break }
    }
    results[i] <- sd.value  # store the current sd.value that is OK!
    sd.value <- 1
  }
  cbind(results, rsquare)
}

simulate.sd.iter <- function(sd.value, n=200) {  # helper function
  # Takes the sd.value, creates data, and returns the r-squared
  X1 <- rnorm(n, 0, 1)
  X2 <- rnorm(n, 0, 1)
  Y <- rnorm(n, (5 + 3*X1 - 2*X2), sd.value)
  simdata <- data.frame(X1, X2, Y)
  return(summary(lm(Y ~ X1 + X2, data=simdata))$r.squared)
}

simulate.sd()

需要注意的几点：

  • 我让 X1 和 X2 变化，因为这会影响这个寻求sd.value的 .
  • 容差是您希望这个估计值有多精确。r-squared 为 ~0.19 或 ~0.21 是否合适？公差为 0.01。
  • 请注意，过于精确的容差可能无法让您找到结果。
  • 1 的值是一个非常糟糕的起始值，使得这个迭代算法非常慢。

10 个结果的结果向量是：

[1] 5.64 5.35 5.46 5.42 5.79 5.39 5.64 5.62 4.70 5.55,

在我的机器上大约需要 13 秒。

我的下一步将从 4.5 开始，在迭代中添加 0.001 而不是 0.01，并可能降低容差。祝你好运！

好的，nsim=100 的一些汇总统计，耗时 150 秒，步长增加了 0.001，容差仍为 0.01：

  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 4.513   4.913   5.036   5.018   5.157   5.393

你为什么对这个感兴趣？

于 2013-09-30T15:30:02.143 回答
-1

这是另一个生成多元线性回归的代码，错误遵循正态分布：OPS sorry this code just generated multiple regression 

sim.regression<-function(n.obs=10,coefficients=runif(10,-5,5),s.deviation=.1){
  
  n.var=length(coefficients)  
  M=matrix(0,ncol=n.var,nrow=n.obs)
    
  beta=as.matrix(coefficients)
  
  for (i in 1:n.var){
    M[,i]=rnorm(n.obs,0,1)
  }
  
  y=M %*% beta + rnorm(n.obs,0,s.deviation)
  
  return (list(x=M,y=y,coeff=coefficients))
  
}

于 2014-09-26T11:59:05.383 回答