潜类别分析相对熵计算：两种方法结果存异的原因探究

我用R语言的poLCA包运行包含4个分类指标（3水平、3水平、9水平、5水平）的潜类别模型，由于poLCA不计算相对熵，我从相关资料中找到两种手动计算的公式，但在自有数据中得到了有差异的结果：

• 示例数据中，两种公式结果分别为0.7379364和0.7254486；
• 全量数据（N>6000）中差异更大，分别为0.72和0.68；
  但在poLCA包自带的癌（carcinoma）示例数据中，两种公式结果完全一致。我已剔除所有含缺失值的观测，想知道这两种公式是否存在差异，或是我的应用方式有误。

附可复现代码：

var1<-c(1,1,1,1,1,1,1,3,1,1,3,1,2,2,1,1,1,1,1,1,1,3,2,1,1,1,1,1,1,1,1,1,1,1,3,2,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,3,2,2,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,2,3,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,2,1,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,3,3,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,2,1,1,3,1,1,1,3,1,1,3,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
var2<-c(3,2,2,1,2,1,2,2,1,1,1,3,2,3,1,3,2,1,1,3,1,2,1,3,1,1,3,3,3,3,1,2,1,3,2,3,1,1,1,3,3,3,1,3,3,2,3,3,3,3,2,2,2,3,3,1,3,1,1,1,2,3,1,3,2,1,1,1,1,1,3,1,3,2,1,1,1,1,1,1,2,1,1,1,1,2,1,2,3,3,3,2,1,1,2,3,3,2,1,3,3,3,3,3,3,1,2,3,3,1,3,3,3,3,3,2,1,3,2,3,1,1,1,2,2,2,1,2,1,2,2,1,1,3,1,3,1,2,3,2,1,2,1,3,1,1,1,1,2,2,2,2,1,2,3,1,3,1,1,1,2,2,1,2,2,3,2,3,1,2,3,3,3,3,3,3,3,3,2,3,3,3,3,3,1,3,1,3,3,1,1,2,1,1,1,3,2,3,3,1,3)
var3<-c(3,8,2,3,1,8,1,1,8,8,1,8,2,8,6,6,8,9,8,4,2,2,8,6,6,6,5,6,2,6,8,2,2,9,2,9,2,8,8,4,4,2,5,8,6,2,2,2,3,2,8,8,2,4,5,9,1,1,1,8,5,3,8,3,4,3,6,1,1,2,8,1,6,5,8,4,8,8,8,8,9,8,4,3,4,1,9,1,4,3,1,2,1,2,5,8,8,4,9,4,8,8,8,4,8,8,2,8,5,2,3,6,4,9,8,2,2,1,1,3,8,1,1,4,2,5,8,1,2,8,4,1,8,8,8,4,9,4,8,5,8,4,8,4,3,8,9,8,4,9,4,4,9,9,3,8,8,8,8,8,4,3,8,4,9,4,4,4,8,4,9,4,5,8,6,8,4,4,1,2,3,3,8,4,3,3,2,6,9,2,8,4,4,8,9,8,9,2,4,1,6)
var4<-c(1,2,1,1,1,1,1,2,1,1,2,2,1,2,1,1,1,2,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,4,1,1,1,1,1,1,1,1,1,1,1,1,2,1,4,1,1,2,2,1,1,2,1,3,2,1,1,1,1,1,1,1,4,1,1,4,1,1,1,1,1,2,1,1,3,1,1,1,2,1,1,1,1,1,1,1,3,2,1,1,2,2,1,1,1,1,1,1,1,1,1,2,2,1,3,1,1,1,1,1,1,2,1,1,2,1,2,1,1,2,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,2,1,1,3,1,2,1,1,1,1,1,1,2,2,1,1,1,4,1,1,2,1,2,1,1,1,1,3,1,2,1,1,1,1,1,1,1,1,1,1,1,1,4,2,1,1,1,1)

ex.data<-data.frame(var1, var2, var3, var4)

f <- cbind(var1, var2, var3, var4)~1

lc.ex<-poLCA(f, ex.data, nclass=3) # 运行3类别模型

# 第一种：Israel Souza的公式
nume.E<- -sum(lc.ex$posterior * log(lc.ex$posterior), na.rm=T)
deno.E<-201*log(3)
ent.ex<-1-(nume.E/deno.E)
ent.ex
# [1] 0.7379364

# 第二种：Daniel Oberski的公式
entropy<-function (p) sum(-p*log(p))

error_prior <- entropy(lc.ex$P)
error_post <- mean(apply(lc.ex$posterior, 1, entropy), na.rm=T)
ent.ex2 <- (error_prior - error_post) / error_prior
ent.ex2
# [1] 0.7254486

这两种公式本质是两种不同定义的相对熵，你的应用方式都是正确的，结果差异源于计算逻辑的不同：

1. 公式差异解析 #

• Israel Souza公式：
  分母用的是完全随机分类的最大熵（即假设每个观测被均匀分配到K个类别，总熵为N*log(K)），相对熵的含义是：实际后验分布的熵，相对于完全随机分类的熵减少了多少比例。公式为：
  相对熵 = 1 - (实际后验总熵 / 最大可能熵)

• Daniel Oberski公式：
  分母用的是模型估计的类别先验熵（基于lc.ex$P，即各潜类别的先验概率计算的熵），相对熵的含义是：后验分布的不确定性，相对于先验分布的不确定性减少了多少比例。公式为：
  相对熵 = (先验熵 - 平均后验熵) / 先验熵

2. 为何癌示例数据结果一致？ #

当模型估计的潜类别先验概率完全均匀（即每个类别的P=1/K）时，先验熵等于log(K)，此时N*log(K)（Souza公式分母）与N*先验熵（Oberski公式分母的N倍，因error_prior是单样本先验熵，平均后验熵也是单样本的，整体比例一致）等价，因此两个公式的计算结果会完全相同。而你的自有数据中，潜类别的先验概率并非均匀分布，所以两个公式的分母逻辑不同，结果自然出现差异。

3. 选择建议 #

两种相对熵都有合理的统计解释，选择哪一种取决于你的研究需求：

• 如果想衡量分类结果相对于完全随机猜测的提升，选Souza的公式；
• 如果想衡量模型后验分布相对于先验分布的不确定性降低程度，选Oberski的公式。

内容的提问来源于stack exchange，提问作者Sointu