1.方法介绍
目前在教育学和心理学的研究中出现了两种不同的观点，以变量为中心的观点(Variable-centred)提倡探究变量与变量之间的关系，例如很多横断面研究测量某个时刻参与率与成绩之间的关系。而以事件为中心的观点(event-centred)则认为学习不是一个静态的过程，参与率和学习成绩会随着时间发生动态变化，探索学习在特定时间点的输入和输出之间的因果关系是没有意义的。因此研究者需要在动态的时间变化中捕捉变量之间的联系。在这种情况下，KML(k-means designed to work specifically on longitudinal data)和KML3d(k-means designed to work specifically on longitudinal joint trajectories)应运而生。
KML和KML3d都是一种纵向聚类技术。即记录下某一个时间段内变量的纵向变化过程，然后基于一个(KML)或者多个(KML3d)变量的纵向轨迹进行聚类，最终得到了多种不同模式的纵向变化的分组。本期内容我们只介绍KML技术，有关KML3d的内容下期会介绍。如图1所示，研究者记录下一个学期(15次授课)中每节课的参与率，绘制了所有班级的参与率的纵向变化曲线，然后基于纵向轨迹进行聚类，区分出了两大类班级，一种是参与率持续走低的班级，一种是参与率普遍更高的班级。区分出班级之后呢，研究者就可以基于组别去进行一些其他的处理，例如比较两组的学习成绩等。
之前的研究一致认为，KmL技术是一种成熟、可靠、更不耗时的纵向聚类方法。它首先应用于流行病学研究中的队列研究，然后逐步扩展到教育和心理学领域。KmL技术是一种使用k-means聚类纵向数据的方法，它不关注单一变量，而是涉及同质个体的轨迹。K-Means是一种EM（期望最大化）爬山算法。由于该算法不需要簇内的任何正态性或参数假设，它在数值收敛方面更鲁棒(翻译成人话就是更稳定)。更重要的是，在纵向数据的背景下，它不需要对轨迹的形状进行任何假设。自从Genolini & Falissard (2010)开发了基于R语言的KmL软件包以来，KmL技术得到了广泛的应用。
2.统计原理
Kml其实就是K-means聚类方法在纵向数据上的使用。K-means算法是一种非参数爬山算法，将n个具有p个变量的对象划分为k个同质子群(这个定义包含了KML3d，KML3d主要是基于多个变量的纵向轨迹进行联合聚类，由于当变量等于2个时，再加上时间轴，这就是一个3d图，如果再增加变量就是4维，我们很难通过图形表示出来)。通过利用KmL包，我们可以通过多次迭代获得最优的簇数及其轨迹。一个典型的KML的步骤包括以下步骤：
1)定义了集群（种子）的初始中心。
2)根据初始种子，得到对象与种子向量之间的距离。然后，每个观测结果都被分配到最近的集群中心。KmL可以使用欧几里德距离法(Euclidean distance method)和曼哈顿距离法(Manhattan distance method)。对于给定的n个对象的集合S，下图两个公式分别代表欧式距离（左）和曼哈顿距离（右）的计算方法：
公式中的yi和yj分别表示在t个不同的时刻被试i和被试j所对应的变量y的值。DE(yi,yj)和DM(yi,yj)分别表示i被试和j被试y轨迹之间的的欧式距离和曼哈顿距离。yil表示在l时刻被试i的y变量的数值，而yjl则表示在l时刻被试j的y变量的数值。
3)根据新的分配情况计算每个集群的中心，并替换初始种子。
4)计算对象和新种子之间的距离，并将每个观测值重新分配到距离最小的集群中。
5)使用更新后的集群成员关系重新计算每个集群的中心。
6)重复步骤4和步骤5，直到集群中没有发生进一步的变化为止。
此外，为了降低算法收敛到局部最大值的风险，KML多次执行该过程，并保留最优解。KML软件包的默认合并将运行算法20次。不仅如此，KML使用Calinski and Harabasz准则来确定集群的最佳数量，默认情况下测试2、3、4、5和6个集群。在实际运用过程中，通过比较多个准则的最优聚类数，我们能够得到最终的类别数。
3.如何实现
我们先简单介绍一下kml包自带的epipageShort数据集，方便我们了解后续的处理。epipageShort是一项法国对严重早产儿的多地区队列调查。其中包括4000多名出生在33周以下的儿童，以及两个对照样本，分别在33-34周出生和足月出生。对约2600名严重早产儿童和400名和600名对照组进行随访，直至5岁，然后至8岁。图中的gender表示性别，sdq是Strengths and Difficulties Questionnaire（优势与困难问卷）的得分，3，4，5，8分别表示在3岁，4岁，5岁，8岁时的得分。
没有安装kml包的小伙伴需要先安装这个包，然后将这个包导入到目前的工作环境中。
>>> install.packages("kml") #安装kml包
>>> library(kml) #导入kml包
>>> data("epipageShort", package = "kml") #这一步是导入kml中自带的epipageShort数据集
>>> head(epipageShort) #简单看一下epipageShort数据集长什么样子
从数据中，我们发现在追访过程中出现了缺失值，如果数据量过大我们可能很难肉眼看出来，我们可以通过下列代码进行数据缺失的简单查看。
>>> sum(is.na(epipageShort)) 计算出数据集中的所有缺失值个数，输出结果405
然而kml不允许出现缺失，否则就会报错，所以研究者在正式运行kml程序前需要解决缺失值问题。KML包提供了自带的缺失值处理方法，方便研究者更好地处理数据缺失问题。缺失值处理这块内容过于庞大，我们后期打算单独开一个章节讲述在横断面研究和纵向研究中常用的缺失值处理技术。在这里我们仅仅简单进行介绍。
在横断面研究中，被试i在j时刻在变量Y上出现缺失，我们往往会用其他所有被试在j时刻Y变量上的均值等一系列指标进行插值，而在纵向研究中，被试i在j时刻在变量Y上出现了缺失，那我们就需要用被试i在其他所有没有缺失时刻的Y数据进行插补操作，例如均值、中位数、Hot Deck、LOCF、线性插值以及样条插值等。
而KML中则使用了相对更新的结合横断面和纵向构成的一种新的插值方法：Copy Mean（关于这种方法这里不做过多介绍，总的来说就是将纵向特征和人群特征相结合）。在KML包中可以通过下列代码实现这一过程。
>>> imputation(as.matrix(epipageShort[, 3:6])) #kml的默认插值方法就是copy mean，由于imputation函数的对象是矩阵，因此我们需要想将其转换成矩阵。
当然，如果这种方法结果不好，你也可以通过修改method参数调整其他的方法。
>>> new_data <- imputation(as.matrix(epipageShort[, 3:6]), method = "locf") #也可以换成linearInterpol
>>> head(new_data)#此时再看数据可以发现已经插补完成了。
完成了对数据的简单处理，我们可以开始构建纵向轨迹了，sdq的轨迹在数据框的第3到第6列中。所以我们可以构建类“cld”的对象
>>> cldSDQ <- cld(new_data, timeInData = 1:4) #1:4表示的是数据框第1列至4列的数据是我需要用来纵向聚类的时间序列数据然后使用函数kml()对数据进行分区。首先，我们希望“看到”集群过程。请注意，由于我们程序运行的尽可能快所以将nbRedrawing 设置为2，但“nbRedrawing”参数代表每个类别数K-means必须运行的次数，一般设置为20次。“toplot”参数有4个取值，=traj意味着只画出轨迹曲线，=criterion意味着只画出分类标准图，=both意味着两张图都画出来，=none表示什么都不画。
>>> kml(cldSDQ, nbRedrawing = 20, toPlot = "none") # ~fast kml~
>>> kml(cldSDQ, nbRedrawing = 20, toPlot = "both") # ~slow kml~
需要注意的是，在kml背后，有两种不同的程序，一种是fast程序，一种是slow程序。当将参数distance设置为"euclidean"时或者将toPlot 设置为 'none' 或者 'criterion' 时，kml运行的时fast程序；而当用户打算自定义距离或者想要去看聚类的结构(将toPlot 设置为 'traj' 或者 'both' 时)运行的则是slow程序。Fast程序比slow程序快25倍，但保持其他参数一致的情况下，结果上并无差异。
>>> kml(cldSDQ) #仅仅输入cldSDQ默认运行fast程序
>>> choice(cldSDQ) #使用函数choice()来探索该算法找到的所有分区上述两行代码一起运行可以生成下图。这幅图左边表示根据Calinski.Harabatz准则在多次迭代过程中所确定的每个类别所对应的Calinski.Harabatz指数。程序一般会选择对应Calinski.Harabatz指标最大的类别作为最佳类别数。因此在右图中按照这一标准分成了两大类，其中A类表示量SDQ较低组，B类表示SDQ较高组。
>>> kml(cldSDQ)
>>> plotAllCriterion(cldSDQ) # 绘制出选择最佳类别数的5个不同的质量标准
然而需要注意的是，质量标准并不总是有效的。使用其中的几种可能会加强结果的可靠性。函数plotAllCriterion显示了算法估计的5个准则所对应的最佳聚类数。为了将它们绘制在同一个图上，它们映射到[0,1]的尺度上。图中共5个准则，每个准则对应不同的颜色，其中1，3，5对应的准则推荐分为2类，而2，3准则推荐分为6类，这意味着没有明确的证据表明有任何“最佳集群数”。因此，在实际的研究中，最终选择聚成多少类别不仅仅需要参考质量标准，还需要结合研究领域和研究问题进行确定。
例如，就上述问题，研究者根据专家的建议最终分成了4类。
>>> kml(cldSDQ, 4, parAlgo = parALGO(distance = function(x, y) + cor(x, y), saveFreq = 10))#这里的4意味着分为4类， parAlgo参数这里的设置意味着以相关距离聚类，保存频率为10
一个需要注意的细节时，当你决定聚为4类时，请清空工作空间中的所有历史记录，否则你在运行时程序仍然为给你输出两类，因为它数据驱动地系统比较了各个类别的质量标准，认为2类最好。
>>> epipageShort$clusters <- getClusters(cldSDQ, 4)#提取每个被试对应的类别到epipageShort中的clusters列，用作后续处理
>>> epipageGroupAD <- epipageShort[epipageShort$clusters %in% c("A", "D"), ]#只提取出A类和D类进行比较
>>> summary(glm(clusters ~ gender, data = epipageGroupAD, family = "binomial"))# 进行逻辑回归
由上图可知，性别似乎并不能解释A类和D类之间的差异，或许存在其他因素的影响。
如果你需要更详细的代码和进一步了解，请关注公众号（Psych统计自习室）并在后台留言，我们会及时进行回复。
参考文献：
Wang, F., Zhu, X., Pi, L. et al. Patterns of participation and performance at the class level in English online education: A longitudinal cluster analysis of online K-12 after-school education in China. Education and Information Technologies (2024). https://doi.org/10.1007/s10639-024-12451-2.
Genolini, C., Alacoque, X., Sentenac, M., & Arnaud, C. (2015). kml and kml3d: R packages to cluster longitudinal data. Journal of statistical software, 65, 1-34.
Genolini, C., & Falissard, B. (2011). KmL: A package to cluster longitudinal data. Computer methods and programs in biomedicine, 104(3), e112-e121.
Dullo, T. T., Kalyanapu, A. J., & Teegavarapu, R. S. (2017). Evaluation of changing characteristics of temporal rainfall distribution within 24-hour duration storms and their influences on peak discharges: case study of Asheville, North Carolina. Journal of Hydrologic Engineering, 22(11), 05017022.
5284
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
