如何使用真实模型编写一个循环来模拟 t 统计量在 null 下的抽样分布？

Posted 2023-03-12

【中文标题】如何使用真实模型编写一个循环来模拟 t 统计量在 null 下的抽样分布？

【英文标题】：How to write a loop to simulate sampling distribution of t-statistic under null using a true model?

【发布时间】：2019-09-24 09:45:35

【问题描述】：

我目前遇到的问题是了解如何模拟 10,000 次绘制并修复协变量。

Y
<int>
X1
<dbl>
X2
<dbl>
X3
<int>
1   4264    305.657 7.17    0
2   4496    328.476 6.20    0
3   4317    317.164 4.61    0
4   4292    366.745 7.02    0
5   4945    265.518 8.61    1
6   4325    301.995 6.88    0
6 rows

那是杂货店代码的头部。

到目前为止，我针对其他相关问题所做的工作：

#5.
#using beta_hat
#create a matrix with all the Xs and numbers from 1-52
X <- cbind(rep(1,52), grocery$X1, grocery$X2, grocery$X3)
beta_hat <- solve((t(X) %*% X)) %*% t(X) %*% grocery$Y
round(t(beta_hat), 2)

#using lm formula and residuals
#lm formula
lm0 <- lm(formula = Y ~ X1 + X2 + X3, data = grocery)

#6.
residuals(lm0)[1:5]

以下是原始函数中的 lm() ：

Call:
lm(formula = Y ~ X1 + X2 + X3, data = grocery)

Residuals:
    Min      1Q  Median      3Q     Max 
-264.05 -110.73  -22.52   79.29  295.75 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4149.8872   195.5654  21.220  < 2e-16 ***
X1             0.7871     0.3646   2.159   0.0359 *  
X2           -13.1660    23.0917  -0.570   0.5712    
X3           623.5545    62.6409   9.954 2.94e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 143.3 on 48 degrees of freedom
Multiple R-squared:  0.6883,    Adjusted R-squared:  0.6689 
F-statistic: 35.34 on 3 and 48 DF,  p-value: 3.316e-12

结果应该是一个循环，可以做t检验中的抽样分布。现在我要解决的另一个问题是基于数据拟合模型。

这里给出了真实模型（对于真实假设），但不确定从哪里开始循环。

【问题讨论】：

我已经进行了编辑，因此代码更具可读性并包含问题的图像，因为符号很奇怪。如果有人可以提供第一步的指导，非常感谢。

所以，您有 X1、X2 和 X3。您可以使用 beta1、beta2 和 beta3 的真实参数来生成一些 Y*。您将随机噪声 epsilon 添加到 Y* 以获得 Y。然后您运行回归。您重复 10000 次的整个过程（并且在每次迭代中，您都会得到新的随机噪声）。你到底在哪里卡住了？

@coffeinjunky 这很有帮助，但是您从哪里获得随机噪声 epsilon？我知道 X1 X2 和 X3 是从数据集中提取的。

根据您的等式 5，随机噪声正态分布与给定的均值和方差。只需从正态分布中提取即可。

如果以下内容回答了您的问题，请考虑通过单击答案开头旁边的相应按钮来接受答案。这是为了让其他人可以看到这不再是一个悬而未决的问题。如果是，请澄清尚不清楚的地方。

【参考方案1】：

好的，请看以下内容：

# get some sample data:
set.seed(42)
df <- data.frame(X1 = rnorm(10), X2 = rnorm(10), X3 = rbinom(10, 1, 0.5))
# note how X1 gets multiplied with 0, to highlight that the null is imposed.
df$y_star <- with(df, 4200 + 0*X1 - 15*X2 + 620 * X3)
head(df)
            X1         X2 X3   y_star
1   1.37095845  1.3048697  0 4180.427
2  -0.56469817  2.2866454  0 4165.700
3   0.36312841 -1.3888607  0 4220.833
4   0.63286260 -0.2787888  1 4824.182
5   0.40426832 -0.1333213  0 4202.000

# define function to get the t statistic
get_tstat <- function()
  # declare the outcome, with random noise added:
  # The added random noise here will be different in each draw
  df$y <- with(df, y_star + rnorm(10, mean = 0, sd = sqrt(20500)))
  # run linear model
  mod <- lm(y ~ X1 + X2 + X3, data = df)
  return(summary(mod)$coefficients["X1", "t value"])


# get 10 values from the t-statistic:
replicate(10, get_tstat())
 [1] -0.8337737 -1.2567709 -1.2303073  0.3629552 -0.1203216 -0.1150734  0.3533095  1.6261360
 [9]  0.8259006 -1.3979176