时间戳

首页 > 解决方案 > “if (!all(pars < 1e-06)) pars[pars < 1e-06] <- 0”错误在带有 depmixS4 的模型中

r - “if (!all(pars < 1e-06)) pars[pars < 1e-06] <- 0”错误在带有 depmixS4 的模型中

问题描述

问题

我想在 r 中使用 depmixS4 包运行潜在类分析。尝试拟合只有一个类（或 depmixS4 包中的状态）的模型时会出现问题。当我尝试使用 6000 个案例的数据集调整模型时，出现以下错误。但是，当病例数为 5000 时，就没有问题了。

Error in if (!all(pars < 1e-06)) pars[pars < 1e-06] <- 0 : 
  missing value where TRUE/FALSE needed

问题出在哪里？有人可以帮我理解为什么会发生这个错误吗？

一个可重现的例子

案例 A (n = 6000) 同样的案例也发生在偶然变量上。为了有一个可重现的例子，首先我生成一个数据集（n = 6000），其中包含两个随机变量（a 和 b）和两个可能的值（0 和 1）。

library(depmixS4)
#> Loading required package: nnet
#> Loading required package: MASS
#> Loading required package: Rsolnp
a <- sample(0:1, size = 6000, replace = T)
b <- sample(0:1, size = 6000, replace = T)

foo_large <- data.frame(a,b)

set.seed(123)
mod1 <- mix(response = list(a~1, b~1),
            data=foo_large, # the dataset to use
            nstates=1, # the number of latent classes
            family=list(multinomial("identity"),multinomial("identity")))

fmod1 <- fit(mod1, verbose=TRUE)
#> Error in if (!all(pars < 1e-06)) pars[pars < 1e-06] <- 0: missing value where TRUE/FALSE needed

案例 B (n = 5000)但是，对于一个数据集 (n = 5000 )，其中两个随机变量的特征与前面的相同，没有错误。

library(depmixS4)
#> Loading required package: nnet
#> Loading required package: MASS
#> Loading required package: Rsolnp
c <- sample(0:1, size = 5000, replace = T)
d <- sample(0:1, size = 5000, replace = T)

foo_short <- data.frame(c,d)

set.seed(123)
mod1 <- mix(response = list(c~1, d~1),
            data=foo_short, # the dataset to use
            nstates=1, # the number of latent classes
            family=list(multinomial("identity"),multinomial("identity")))

fmod1 <- depmixS4::fit(mod1, verbose=TRUE)
#> iteration 0 logLik: -6928.943 
#> converged at iteration 1 with logLik: -6928.943

标签： rglmmodel-fitting

解决方案

我做了一些挖掘，错误似乎是由于 depmixS4 提供随机起始值来初始化 EM 算法的方式（它为具有 Dirichlet 分布的类成员资格生成随机概率，并且我们用来从该分布中提取的代码不起作用很好的一维狄利克雷）。我们将在即将发布的版本中解决此问题。现在，您可以使用以下命令在没有随机起始值的情况下运行 EM：

fmod1 <- fit(mod1, emcontrol=em.control(random.start=FALSE), verbose=TRUE)

这在您的两个示例中都有效。

请注意，问题不是由于观察次数(n=5000)或(n=6000)的差异。代码收敛于(n=6000)是使用set.seed(123). 删除此行后，您很可能会遇到与(n=6000)相同的错误。如果你设置了后者，你可以巧合地开始工作set.seed(1234)。

推荐阅读