r - 在 R 模拟中估计线性回归和线性混合模型中的偏差
问题描述
我想运行模拟来估计线性模型和线性混合模型中的偏差。偏差是 E(beta)-beta,其中 beta 是我的 X 和 Y 之间的关联。
我从正态分布生成 X 变量,从多元正态分布生成 Y。
我了解如何从模拟中计算 E(beta),即所有模拟的 beta 估计值之和除以模拟总数,但我不确定如何估计真正的 beta。
meanY <- meanY + X*betaV
这就是我如何生成均值 Y(betaV 是效应大小),然后用于生成多元 Y 结果,如下所示。
Y[jj,] <- rnorm(nRep, mean=meanY[jj], sd=sqrt(varY))
我了解如何从模拟中计算 E(beta),即所有模拟的 beta 估计值之和除以模拟总数,但我不确定如何估计真实的 beta。
根据我有限的理解,真正的 beta 不是从数据中获得的,而是从我设置固定 beta 值的设置中获得的。
根据我生成数据的方式,我如何估计真正的 beta?
解决方案
有几种模拟偏差的方法。我将使用线性模型举一个简单的例子。线性混合模型可能会使用类似的方法,但是我不确定它是否适用于广义线性混合模型(我只是不确定)。
当使用简单的线性模型时,估计偏差的一种简单方法是“选择”从哪个模型来估计偏差。比方说Y = 3 + 4 * X + e
。我选择beta <- c(3,4)
了,因此我只需要模拟我的数据。对于线性模型,模型假设是
- 观察是独立的
- 观测值是正态分布的
- 均值可以用线性预测器来描述
使用这 3 个假设,模拟固定设计很简单。
set.seed(1)
xseq <- seq(-10,10)
xlen <- length(xseq)
nrep <- 100
#Simulate X given a flat prior (uniformly distributed. A normal distribution would likely work fine as well)
X <- sample(xseq, size = xlen * nrep, replace = TRUE)
beta <- c(3, 4)
esd = 1
emu <- 0
e <- rnorm(xlen * nrep, emu, esd)
Y <- cbind(1, X) %*% beta + e
fit <- lm(Y ~ X)
bias <- coef(fit) -beta
>bias
(Intercept) X
0.0121017239 0.0001369908
这表明偏差很小。为了测试这种偏差是否显着,我们可以执行 wald 检验或 t 检验(或重复该过程 1000 次,并检查结果的分布)。
#Simulate linear model many times
model_frame <- cbind(1,X)
emany <- matrix(rnorm(xlen * nrep * 1000, emu, esd),ncol = 1000)
#add simulated noise. Sweep adds X %*% beta across all columns of emany
Ymany <- sweep(emany, 1, model_frame %*% beta, "+")
#fit many models simulationiously (lm is awesome!)
manyFits <- lm(Y~X)
#Plot density of fitted parameters
par(mfrow=c(1,2))
plot(density(coef(manyFits)[1,]), main = "Density of intercept")
plot(density(coef(manyFits)[2,]), main = "Density of beta")
#Calculate bias, here i use sweep to substract beta across all rows of my coefficients
biasOfMany <- rowMeans(sweep(coef(manyFits), 1, beta, "-"))
>biasOfMany
(Intercept) X
5.896473e-06 -1.710337e-04
在这里,我们看到偏差减少了很多,并且改变了 betaX 的符号,让我们有理由相信偏差是微不足道的。
改变设计将允许人们使用相同的方法研究交互、异常值和其他东西的偏差。
对于线性混合模型,可以执行相同的方法,但是在这里您必须设计随机变量,这需要更多的工作,并且lmer
据我所知,实现并不适合所有列的模型Y
。
然而b
(随机效应)可以被模拟,任何噪声参数也可以。但是请注意,由于b
是包含单个模拟结果(通常为多元正态分布)的单个向量,因此必须为每个模拟重新运行模型b
。基本上,这将增加人们必须重新运行模型拟合过程的次数,以便获得对偏差的良好估计。
推荐阅读
- php - update.php 不返回错误,但不更新数据库
- macos - 当前的 Apple 沙盒文档在哪里
- clojure - Clojure + Seesaw:为什么这个 GUI 不起作用?
- python - 一切似乎都被定义为一个字符串?Python
- html - Flexbox div 隐藏在小屏幕上
- automated-tests - 用于并行执行的上下文注入
- python - Python 3.7,MySql-Python 的构建轮失败
- java - 如何存储从mysql检索到的数据并显示为用户在java中搜索的结果
- image - NativeScript 在动态网格中显示图像
- mysql - 续集模型关系