R语言基于poLCA包进行潜类别分析

潜在类别分析是一种分析多元分类数据的统计技术。当观测数据以一系列分类响应的形式出现时- -例如,在民意调查、个人层面的投票数据、人与人之间可靠性的研究或消费者行为和决策中- -通常感兴趣的是调查观测变量之间的混淆来源,识别和表征相似案例的集群,并在许多感兴趣的变量中近似观测值的分布。潜在类别模型是实现这些目标的有用工具。
下面我们通过R语言poLCA包来演示一下,我们先导入R包和数据

library(poLCA)
election<-read.csv("E:/r/test/xuanju.csv",sep=',',header=TRUE)

在这里插入图片描述
这是一个调查问卷数据,调查数据来自2000年美国全国选举研究,公众号回复:选举数据,可以获得数据。两套共6个问题,每个问题4个答案,询问受访者对(道德性、关怀性、知识性、好领导、不诚实)描述总统候选人阿尔·戈尔和乔治· W ·布什的各种特质的看法。分为4个等级( 1 )极好;( 2 )相当好;( 3 )不太好;( 4 )一点也不好。该数据集还包括潜在的协变量VOTE3、被调查者的2000次投票选择(询问时);年龄,受访者年龄;EDUC,受访者的受教育程度;性别,受访者的性别;和第四部分,受访者的民主-共和两党认同。
poLCA包进行分析变量的formula(公司)基本格式为:

f <- cbind(Y1,Y2,Y3)~1

如果含有协变量的话格式为:

f <- cbind(Y1,Y2,Y3)~X1+X2*X3

Y1,Y2,Y3需要是等级变量,在本数据中就是这些问卷的条目,如果不指定协变量,就是

f <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
           MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~1

生成公式后就可以使用poLCA函数分析了,它的基本格式为

poLCA(formula, data, nclass=2, maxiter=1000, graphs=FALSE, tol=1e-10,
na.rm=TRUE, probs.start=NULL, nrep=1, verbose=TRUE, calc.se=TRUE)

data就是你的数据,nclass代表你要分成几个类型,maxiter表示算法最大迭代次数,默认1000次,graphs等于T的话会生成图形,tol:一个容差值用于判断何时达到收敛,一般默认值就可以了。nrep表示对模型估计的次数,nrep > 1可以自动搜索对数似然函数的全局最大值而不是局部最大值。poLCA只返回产生最大对数似然函数的模型对应的参数估计。
我们先分成一类看看

nes1 <- poLCA(f,election,nclass=1)  # log-likelihood: -18647.31

在这里插入图片描述
模型返回了很多数据,除了最大似然比还有AIC、bic、Estimated等,我们把它分成两类和三类看看

nes2 <- poLCA(f,election,nclass=2)  # log-likelihood: -17344.92

在这里插入图片描述

nes3 <- poLCA(f,election,nclass=3)  # log-likelihood: -16714.66

在这里插入图片描述
我们可以看到分成3类后最大似然比较前缩小,说明3类比前两类更加适合,模型还返回了每组变量的相应概率,图比较大,我截图一部分
在这里插入图片描述
这里表示先验和后验概率
在这里插入图片描述
如果把graphs=T,可以图示概率响应情况

nes3 <- poLCA(f,election,nclass=3,graphs=T)  # log-likelihood: -16714.66

在这里插入图片描述
参照别的博主做法,可以把数据提出来,使用ggplot绘图,这样美观一点

library(reshape2)
library(ggplot2)

dat1<- melt(nes3$probs, level = 3)

ggplot(dat1,aes(x = L3, y = value, fill = Var2)) +
  geom_bar(stat = 'identity', position = 'stack', width = 0.5) +
  facet_grid(Var1~.) +
  scale_fill_brewer(type = 'seq', palette = 'Greys') +
  theme_bw() +
  labs(x = '', fill = 'probabilities') +
  guides(fill = guide_legend(reverse = TRUE))

在这里插入图片描述
如果我们不知道这个数据应该分几类,可以写个循环来确定,最大似然比、AIC、bic都是很好的分类指标,我们可以写一个简单的小循环来确定

stat<-NULL
for(i in 2:10){
  out<- poLCA(f,election,nclass=i,graphs=F) 
  likelihood<-out[["aic"]]
  stat<-append(stat,likelihood)
}
stat
min(stat)

在这里插入图片描述
最后分成几类不能单看指标,还要它能不能解释,比如本例数据是个两个总统戈尔和布什的选举数据。一个合理的假设调查对象存在三个潜在类别:戈尔支持者、布什支持者和或多或少中立者。戈尔的支持者会倾向于对戈尔做出有利的反应,而对布什做出不利的反应,布什的支持者则相反。中立组的人不会对任何一个候选人有强烈的意见。民主党更倾向于支持戈尔,坚定的共和党更倾向于支持布什,而不太激烈的党派倾向于无动于衷。
本例数据中的PARTY为协变量,1是支持民族党,7支持共和党,3 - 4 – 5是属于中立,我们把PARTY加入模型中,并建立模型

f.party <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
                 MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~PARTY
nes2a <- poLCA(f2a,election,nclass=3,nrep=5)    # log-likelihood: -16222.32 graphs

这里我们把中立者参考对比对象

probs.start <- poLCA.reorder(nes.party$probs.start,
                             order(nes.party$P,decreasing=T))
nes.party <- poLCA(f.party,election,nclass=3,probs.start=probs.start)

在这里插入图片描述
我们可以看到模型把投票者分成了3个潜在人群,三个群体的分离,比例分别为27 %、34 %和39 %。
在模型输出结果列表中(下图),我们可以看到
在这里插入图片描述
因为模型的公式为
在这里插入图片描述
所以我们可以得知中立者和布什群体的对数比先验概率为ln ( p2i / p1i) = -3.82 + 0.79 × PARTY。同样中立者和戈尔群体的对数比先验概率为ln( p3i / p1i) = 1.16-0.57 × PARTY。
我们可以通过如下公式转换
在这里插入图片描述
R的操作为

pidmat <- cbind(1,c(1:7))
exb <- exp(pidmat %*% nes2a$coeff)

绘图

matplot(c(1:7),(cbind(1,exb)/(1+rowSums(exb))),ylim=c(0,1),type="l",
        main="Party ID as a predictor of candidate affinity class",
        xlab="Party ID: strong Democratic (1) to strong Republican (7)",
        ylab="Probability of latent class membership",lwd=2,col=1)
text(5.9,0.35,"Other")
text(5.4,0.7,"Bush affinity")
text(1.8,0.6,"Gore affinity")

在这里插入图片描述
亲民主党人有超过60 %的先验概率属于戈尔亲和集团,而亲共和党人有超过80 %的先验概率属于布什亲和集团。
假如我们想进一步了解年龄的影响是否会改变党派偏见对候选人的影响,可以进一步设立交互模型

f.3cov <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
                MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB)~PARTY*AGE
nes.3cov <- poLCA(f.3cov,election,nclass=3,verbose=F)

继续把中立者设为参考对象,可以看到多了一个交互选项值

probs.start <- poLCA.reorder(nes.3cov$probs.start,
                             order(nes.3cov$P,decreasing=T))
nes.3cov <- poLCA(f.3cov,election,nclass=3,probs.start=probs.start)

在这里插入图片描述
为了考察年龄对强势党派候选人亲和力的影响,设定协变量的假设值矩阵:民主党的strdems和共和党的strreps。然后计算并绘制每个假设值对应的潜在类别成员的预测先验概率。
民主党的先验概率

strdems <- cbind(1,1,c(18:80),(c(18:80)*1))
exb.strdems <- exp(strdems %*% nes.3cov$coeff)
matplot(c(18:80),(cbind(1,exb.strdems)/(1+rowSums(exb.strdems))),
          main="Age and candidate affinity for strong Democrats",
          xlab="Age",ylab="Probability of latent class membership",
          ylim=c(0,1),type="l",col=1,lwd=3)
text(50,0.6,"Gore affinity")
text(50,0.39,"Other")
text(50,0.1,"Bush affinity")

在这里插入图片描述
共和党的先验概率

strreps <- cbind(1,7,c(18:80),(c(18:80)*7))
exb.strreps <- exp(strreps %*% nes.3cov$coeff)
matplot(c(18:80),(cbind(1,exb.strreps)/(1+rowSums(exb.strreps))),
          main="Age and candidate affinity for strong Republicans",
          xlab="Age",ylab="Probability of latent class membership",
          ylim=c(0,1),type="l",col=1,lwd=3)
text(50,0.8,"Gore affinity")
text(50,0.2,"Other")
text(50,0.05,"Bush affinity")

在这里插入图片描述
由下两图可以得出民族党无论年龄大小,都不会亲布什,共和党人也不会倾向戈尔。,第一张图看到,共和党的成员在30岁一下更多属于中立派。
在这里插入图片描述
参考文献:

  1. poLCA包说明
  2. 张洁婷, 焦璨, 张敏强. 潜在类别分析技术在心理学研究中的应用[J]. 心理科学进展, 2010(12):8.
  3. 房立艳, 张大均, 武丽丽,等. 中学生心理素质的类别特征:基于个体中心的潜在类别分析[J]. 心理与行为研究, 2017, 15(1):6.
  4. Mirza S S , Wolters F J , Swanson S A , et al. 10-year trajectories of depressive symptoms and risk of dementia: a population-based study[J]. Lancet Psychiatry, 2016:628-635.
  5. https://mp.weixin.qq.com/s/Bvob1r2g6vr7gcavMOydAA
### R语言中`poLCA`计算熵值的方法 在R语言中的`poLCA`用于执行潜类别分析(Latent Class Analysis),而熵(Entropy)是一个衡量分类不确定性的重要指标。当利用`poLCA`函数完成模型拟合之后,可以通过访问返回对象中的特定属性来获取熵值。 对于由`poLCA()`创建的对象而言,熵通常不是直接作为输出的一部分显示出来,但是可以借助于该函数产生的结果间接地计算得到。具体来说,在调用了`poLCA()`并保存其结果到某个变量(比如命名为`result`)后,熵可通过如下方式获得: ```r entropy <- result$BIC - logLik(result)[1]*2 ``` 不过上述表达式实际上给出了一个简化版的信息准则差异,并不严格等于传统意义上的熵;为了更精确地反映潜在类别的纯度或者说成员归属的确信程度,则应该采用下面这种方法来估算样本级别的平均熵[^1]: ```r posterior_probs <- result$post.prob[, 1:result$nclass] sample_entropy <- apply(posterior_probs, 1, function(x){ sum(-x * log(x)) }) mean(sample_entropy) ``` 这段代码首先提取了每个观测属于各个隐含类的概率矩阵,接着针对每一个个体应用香农熵公式\[ H(X)=-\sum_{i} p(x_i)\log{p(x_i)} \] 来量化不确定性的大小,最后求取所有观察单位上这种测度的均值以代表整个数据集内的总体混乱状况。 值得注意的是,较低的熵意味着较高的置信水平——即大多数情况下能够清楚地区分不同组别之间的区别;相反,如果熵较高则表明存在较多模棱两可的情况使得划分变得困难起来。 ### 示例代码展示如何计算熵值 假设已经通过`poLCA`得到了名为`fit`的结果对象,那么可以根据上面提到的方式编写一段完整的脚本来实现熵的评估过程: ```r library(poLCA) # 构建测试用的数据框... data(carcinoma) f <- cbind(A,B,C,D,E,F,G)~1 fit <- poLCA(f, carcinoma, nclass=3) # 提取出后验概率表单 post_probabilities <- fit$post.prob[, 1:fit$nclass] # 应用Shannon Entropy Formula逐行处理 individual_entropies <- apply(post_probabilities, 1, function(prob_vector){ sum(-prob_vector * log(prob_vector)) }) # 输出整体平均熵 average_entropy <- mean(individual_entropies) print(paste("Average Sample-Level Entropy:", round(average_entropy, digits = 4))) ```
评论 13
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

天桥下的卖艺者

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值