r - 在 poLCA R 包中寻找最佳 LCA 模型

Question

我正在使用 PoLCA R 包应用 LCA 分析，但分析三天后没有结果(尚未找到最佳模型)，偶尔会出现以下错误:“警报:迭代完成，未找到最大似然” 。所以我在35潜类取消了这个过程。我正在分析 16 个变量(全部是分类变量)和 36036 行数据。当我测试 Boruta 包中 16 个变量的变量重要性时，所有 16 个变量的结果都是重要的，因此我在使用 poLCA 的 LCA 分析中使用了所有 16 个变量。我应该走哪条路？我是否应该使用其他聚类方法(例如 k 模式)来对该数据集中的分类变量进行聚类？我使用具有 500 次迭代和 nrep=10 模型估计数的参数。我用来寻找 LCA 中最佳模型的 R 脚本，输出之一如下:

    for(i in 2:50){
    lc <- poLCA(f, data, nclass=i, maxiter=500, 
                tol=1e-5, na.rm=FALSE,  
                nrep=10, verbose=TRUE, calc.se=TRUE)
    if(lc$bic < min_bic){
        min_bic <- lc$bic
        LCA_best_model<-lc
    }
}

========================================================= Fit for 35 latent classes: ========================================================= number of observations: 36036
number of estimated parameters: 2029 residual degrees of freedom: 34007
maximum log-likelihood: -482547.1
AIC(35): 969152.2 BIC(35): 986383 G^2(35): 233626.8 (Likelihood ratio/deviance statistic)
X^2(35): 906572555 (Chi-square goodness of fit)
ALERT: iterations finished, MAXIMUM LIKELIHOOD NOT FOUND

Answer 1

您使用的脚本会按顺序测试 2 到 50 个类别的每个模型，并保留 BIC 最低的模型。 BIC 并不是选择“最佳”模型的唯一一种或最佳方式，但足够公平。

问题是，您正在估计很多参数，尤其是在最后一步中。您适合的类(class)越多，该过程就越耗时。此外，在这种情况下，收敛问题是可以预料的，因为您要拟合如此多的类。这就是错误消息报告的内容，它无法找到具有 35 个类别的模型的最大可能性。

我不知道你想解决什么问题，但是超过 10 个类别的模型在 LCA 中并不常见。您执行 LCA 是为了尽可能降低数据的复杂性。如果您需要拟合具有超过 10 个类别的模型:

• 将它们一一适应，这样 RAM 消耗就不再是问题了。
• 增加调用中的 nrep= 参数，这样模型偶然找不到最大似然值(坏随机初始数)的概率就会降低。还增加了计算时间。

或者，您可以减少并行运行模型的计算时间。几乎每台现代 PC 都有 2 个或更多内核。下一个 block 中的函数 acl() 使用 foreach() 和 %dopar% 执行此操作，因此与操作系统无关。

library(poLCA)
library(foreach)
library(doParallel)
registerDoParallel(cores=2) #as many physical cores as available.

acl <- function(datos,   #a data.frame with your data 
                k,       #the maximum number of classes to fit
                formula) {  
  foreach(i=1:k, .packages="poLCA") %dopar% poLCA(formula, datos, nclass=i
  )
}

acm() 返回模型列表，您可以稍后选择“最佳”。下一个函数将从列表中检索兴趣的数量，并创建一个格式良好的 data.frame，其中包含有用的信息来选择正确的类别数量。

comparar_clases_acl <- function(modelo) {
  entropy<-function (p) sum(-p*log(p)) #to asses the quality of classification
    tabla_LCA <- data.frame(Modelo=0, BIC=0, Lik_ratio=0, Entropia=0, MenorClase=0)   #empty data.frame to prealocate memory. 
for(i in 1:length(modelo)){
  tabla_LCA [i,1] <- paste("Modelo", i)
  tabla_LCA [i,2] <- modelo[[i]]$bic
  tabla_LCA [i,3] <- modelo[[i]]$Gsq
  error_prior <- entropy(modelo[[i]]$P)
  error_post <- mean(apply(modelo[[i]]$posterior,1, entropy),na.rm = TRUE)
  tabla_LCA [i,4]<-round(((error_prior-error_post) / error_prior),3)
  tabla_LCA [i,5] <- min(modelo[[i]]$P)*100
  }
return(tabla_LCA)
}

它只需要一个参数:一个带有 LCA 模型列表的对象，正是 acl() 返回的内容。

这种并行方法应该可以减少计算时间。但 50 个类(class)仍然太多，您可能会在 50 个类(class)之前获得最小的 BIC。请记住，随着估计参数数量的增加，BIC 会对模型进行惩罚，帮助您找到模型中额外类别的 yield 递减点。