机器算法验证 - 如何模拟随机效应模型？ - 吾爱随笔录

如何模拟随机效应模型？

机器算法验证混合模式重复测量面板数据多层次分析随机效应模型

2022-04-11 12:00:06

模拟线性模型非常简单：

set.seed(42)    
years <- rnorm(100, 12, 8)
work_hours <- rnorm(100, 8, 2)
income <- 2*years + 0.5*work_hours + 2000 + rnorm(100, 0, 10)
plot(work_hours, income2)
lmmodel <- lm(income ~ years + work_hours)
summary(lmmodel)

或逻辑模型：

set.seed(42)
x1 <-  rnorm(100) 
x2 <-  rnorm(100)
z <- 1 + 2*x1 + 3*x2   
pr <- 1/(1+exp(-z)) 

y = rbinom(100,1,pr) 

df <- data.frame(y=y,x1=x1,x2=x2)
logitmodel <- glm( y~x1+x2,data=df,family="binomial")
summary(logitmodel)

那么，如何模拟随机效应模型呢？我的意思是，这类模型有很多“味道”。看看 Faraway 的 [book][1] 有：

块作为随机效应
裂区
嵌套效果
交叉效应
多级模型
重复措施
纵向/面板数据
非正态响应的混合效应模型

我将如何模拟它们以便我可以玩弄它们？

[1]：用 R 扩展线性模型 - John Faraway

2个回答

只需写下模型的（代数）公式，并根据该描述进行模拟。我将举一个非常简单的例子，一个对同一主题有多个观察的模型，具有可交换的协方差结构。这种结构可以用每个受试者的随机截距来表示。也是一个主题级别的协变量： for and在每个主题内. 所以这是一个平衡的模型。相同的原理用于不平衡的情况，但这提供了更多的编程。然后我们必须指定固定参数的值和随机效应的分布。一些简单的 R 代码是：

y_{i j} = μ + α x_{i} + ϵ_{i} + ϵ_{i j}

$y_{ij}=\mu + \alpha x_i + \epsilon_i + \epsilon_{ij}$

i = 1, 2, \dots, n

$i=1,2,\dotsc,n$

j = 1, \dots, k

$j=1,\dotsc,k$

ϵ_{i}, ϵ_{i j}

$\epsilon_i, \epsilon_{ij}$

N <- 20 # Number subjects
k <- 4  # Number obs within subject
set.seed(7*11*13) # My public seed

id <- as.factor(1:N)
x <-  runif(N, 1, 5)
idran <- rnorm(N, 0, 1)
obsran <- rnorm(N*k, 0, 2)
mu <- 10.
alpha <- 1.

X <- rep(x, each=k)
Y <- mu + alpha*X + rep(idran, each=k) + obsran

该模拟数据的图是：

对于更复杂的情况，它会帮助一些预编程的包，simstudyCRAN 上有一个包可以提供帮助。另请参阅混合效应模型的模型矩阵和https://stackoverflow.com/questions/30896540/extract-raw-model-matrix-of-random-effects-from-lmer-objects-lme4-r，https://stackoverflow .com/questions/55199251/how-to-create-a-simulation-of-a-small-data-set-in-r。

这是我模拟随机效果的方法。我将演示线性回归，但将其扩展到不同的 GLM 应该是直截了当的。

让我们从一个随机截距模型开始。模型通常写成

y = X B + Z γ

$y = XB + Z\gamma$

其中是该组的指标，均值为 0 和一些方差的正态分布。该模型的仿真如下... $Z$ $\gamma_i$

groups<-1:5
N <- 250
g <- factor(sample(groups, replace = T, size = N), levels = groups)
x <- rnorm(N)
X <- model.matrix(~x)
Z <- model.matrix(~g-1)

beta <- c(10, 2)
gamma <- rnorm(length(groups), 0, 0.25)

y = X %*% beta + Z%*% gamma + rnorm(N, 0, 0.3)

让我们拟合一个混合模型，看看我们是否能恢复其中的一些估计


library(lme4)

model = lmer(y~x + (1|g), data = d)

summary(model)
inear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | g)
   Data: d

REML criterion at convergence: 136.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.85114 -0.65429 -0.00888  0.65268  2.63459 

Random effects:
 Groups   Name        Variance Std.Dev.
 g        (Intercept) 0.05771  0.2402  
 Residual             0.09173  0.3029  
Number of obs: 250, groups:  g, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  9.95696    0.10914   91.23
x            2.00198    0.01993  100.45

Correlation of Fixed Effects:
  (Intr)
x -0.008

固定效应看起来不错，组标准差 (0.25) 估计得非常准确，残差标准差也是如此。

随机斜率模型类似。假设每个斜率都来自正态分布，那么我们可以将斜率写为

y = B x + β_{i} x

$y = Bx + \beta_i x$

这里是总体均值，是第 i 组的影响。这是一个模拟 $B$ $\beta_i$

library(tidyverse)

groups<-1:5
N <- 250
g <- sample(groups, replace = T, size = N)
x <- rnorm(N)
X <- model.matrix(~x)

B <- c(10, 2)
beta <- rnorm(length(groups), 0, 2)

y = X %*% B + x*beta[g] + rnorm(N, 0, 0.3)

和一个模特……



library(lme4)

d = tibble(y, x, g)

model = lmer(y ~ x + (x|g), data = d)

summary(model)

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (x | g)
   Data: d

REML criterion at convergence: 158.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.95141 -0.65904  0.02218  0.61932  2.66614 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 g        (Intercept) 2.021e-05 0.004496     
          x           3.416e+00 1.848314 1.00
 Residual             9.416e-02 0.306856     
Number of obs: 250, groups:  g, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept) 10.00883    0.01984  504.47
x            2.05913    0.82682    2.49

Correlation of Fixed Effects:
  (Intr)
x 0.099

这是5组的系数

coef(model)
$g
  (Intercept)         x
1    10.00135 -1.015180
2    10.01335  3.919787
3    10.00934  2.270760
4    10.01081  2.873636
5    10.00928  2.246626

并将它们与真实值进行比较

B[2] + beta

-0.9406479  3.9195119  2.2976457  2.8536623  2.3539863

```

其它你可能感兴趣的问题

上一篇我如何需要从 diff(log(data),1) 进行反转换？下一篇应用具有外生变量的 ARIMA 模型进行预测