机器算法验证 - 混合模型模拟的问题 - 吾爱随笔录

混合模型模拟的问题

机器算法验证 r 混合模式 lme4-nlme

2022-03-30 09:43:18

所以我只是度过了一个愉快的夜晚，试图自己模拟一些层次模型并估计它们的参数。我尝试模拟和估计的第一个模型是其中。所以实际上我估计的是。我已经生成了一些数据和估计参数以及误差项的方差，一切都很好：

y_{i j} = β_{0} + ε_{i j},

$y_{ij}=\beta_0+\varepsilon_{ij},$

b_{0} = γ_{00} + u_{0 j}

$b_0=\gamma_{00}+u_{0j}$

y_{i j} = γ_{00} + u_{0 j} + ε_{i j}

$y_{ij}=\gamma_{00}+u_{0j}+\varepsilon_{ij}$

γ_{00}

$\gamma_{00}$

set.seed(1)
N <- 100
nj <- 100
g00 <- 10
e <- rnorm(N*nj)
j <- c(sapply(1:N, function(x) rep(x, nj)))
uj <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))
d <- data.frame(j, y=g00+uj+e)
library(nlme)
lme(y~1, data=d, random=~1|j)

Linear mixed-effects model fit by REML
  Data: d 
  Log-restricted-likelihood: -14520.94
  Fixed: y ~ 1 
(Intercept) 
   10.00215 

Random effects:
 Formula: ~1 | j
        (Intercept) Residual
StdDev:   0.7752422 1.012683

Number of Observations: 10000
Number of Groups: 100

然后我尝试了不同的模型：其中和。所以我不得不估计方程我做了和以前一样的事情，但这一次没有收敛。我试过这个，但它似乎没有用。

y_{i j} = β_{0} + β_{1} x_{i j} + ε_{i j},

$y_{ij}=\beta_0+\beta_1 x_{ij}+\varepsilon_{ij},$

β_{0} = γ_{00} + u_{0 j}

$\beta_0=\gamma_{00}+u_{0j}$

β_{1} = γ_{10} + u_{1 j}

$\beta_1=\gamma_{10}+u_{1j}$

y_{i j} = γ_{00} + u_{0 j} + (γ_{10} + u_{1 j}) \cdot x_{i j} + ε_{i j} = γ_{00} + γ_{10} x_{i} j + u_{0 j} + u_{1 j} x_{i j} + ε_{i j} .

$y_{ij}=\gamma_{00}+u_{0j}+(\gamma_{10}+u_{1j}) \cdot x_{ij}+\varepsilon_{ij}=\gamma_{00}+\gamma_{10}x_ij+u_{0j}+u_{1j}x_{ij}+\varepsilon_{ij}.$ lme

g10 <- 10
u0j <- uj
u1j <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))
x1 <- rnorm(N*nj)
d1 <- data.frame(j, y=g00+u0j+(g10+u1j)*x1, x1)
lme(y~1+x1, data=d1, random=~1+x1|j)

这是最后一次lme吐出的电话：

Error in lme.formula(y ~ 1 + x1, data = d1, random = ~1 + x1 | j) : 
  nlminb problem, convergence error code = 1
  message = false convergence (8)

你能建议我做什么？也许问题出在我的模型规范、冗余参数、奇异矩阵或其他方面？

4个回答

如前所述?lmeControl，默认优化器函数是nlminb而不是optim. ?nlminb，但是，说这optim是首选，实际上，以下工作。

lme(y~1+x1, data=d1, random=~1+x1|j, control = lmeControl(opt = "optim"))
Linear mixed-effects model fit by REML
  Data: d1 
  Log-restricted-likelihood: 320824.3
  Fixed: y ~ 1 + x1 
(Intercept)          x1 
   8.302459    8.183053 

Random effects:
 Formula: ~1 + x1 | j
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev       Corr  
(Intercept) 1.699207e+00 (Intr)
x1          1.952189e+00 0.752 
Residual    1.347943e-15       

Number of Observations: 10000
Number of Groups: 100

很难说为什么会这样。一般来说，可以看出这false convergence (8)意味着

梯度可能计算不正确，其他停止公差可能太紧，或者或附近可能不连续 $\nabla f(x)$ $f$ $\nabla f$ $x$

d1 <- data.frame(j, y=g00+u0j+(g10+u1j)*x1+e, x1) # 需要添加误差项

我的直觉是，当我们模拟没有误差项的混合模型数据集时，有时可能难以收敛。这很可能是一个零残差问题。

在此处输入图像描述

我可能错过了一个更简单的答案（即模型中的错误规范），但是 /if/ 我确信混合模型应该收敛，我会继续使用openbugs，这会让你拖出更长的 mcmc样本。

这对我有用。我没有使用 nlme 或 lmer，而是使用 mgcv::gamm 并指定一个高斯链接，如下所示：

GAMM(y~1+x1, data=d1, random=~1+x1|j,family = gaussian)

估计结果接近nlme/lmer混合模型，适用于大部分情况。事实上，HLM/混合模型是 GAMM 的一个特例。

其它你可能感兴趣的问题

上一篇带有代码示例的机器学习书下一篇估计谱密度