「转」纵向数据的轨迹模型trajectory model—全攻略及R代码

统计碎碎念实战医学统计2023-08-30 20:48:32广东

上回，我们说到SAS如何实现组轨迹模型(GBTM)拟合、评价与正式结果报告，以及其结果解释和呈现方式,得到了各位看官好评。

今天为支持各路研究者的口味，总结各大轨迹模型分享给大家，并提供R实现代码！后续将会针对每个模型，从模型介绍、模型结果正式报告、模型实战和总结四个板块带领大家对每个模型进行系统学习，为您论文方法的选择多一条路径，让您科学合理的使用模型，规范的报告模型，顺利发得高分文章！

目录：

一、纵向数据分析常见研究目的

二、纵向数据特点及数据格式

三、针对纵向数据数据轨迹建模模型

（一）组轨迹模型(GBTM)

（二）潜变量增长混合模型(LGMM)

（三）潜在转移分析(LTA)

（四）时依性混合模型(TVEM)

（五）纵向k-均值聚类(LKMA)

（六）层次凝聚聚类(AHC)

四、全监督聚类—基于特征的聚类

（一）个体时间序列分析(ITSR)

（二）界标模型/标段模型(LME)

话说R虽然有function函数，但是本人还是更偏好SAS宏程序一点，主打一个稳定，现在用SAS运行10年前的代码，还是完美运行。R包更新较快，方法较新，但是SAS确实需要一定的数理统计基础，才能实现各种新奇的方法。本人R语言水平一般，各位大佬批评指正。

正题

今天，小编将使用R软件辅助介绍各路您见过或者没见过的轨迹模型，由于轨迹模型评价以及绘图时需要用到大量的function函数，如您需要，让和您一样从事科研工作的小伙伴点点关注，之后客服会给大家发送本教程所有代码。

后续R代码都是基于本次模拟数据，现将模型数据代码奉上！！

rm(list = ls());gc()generate_osa_data = function(patients=500,                             times=seq(1, 365, by=1),                             nAggr=14,                             props=c(GU=.24, SI=.13, SD=.14, VU=.17, OA=.08, ED=.13, NU=.11),                             # N=c(354, 344, 280, 299, 106, 55, 10),                             # sd.N=c(31, 49, 68, 55, 64, 30, 4),                             dropoutTimes=c(Inf, Inf, Inf, Inf, Inf, 80, 20),                             sd.dropoutTimes=c(0, 0, 0, 0, 0, 30, 10),                             attemptProbs=c(354, 344, 280, 299, 106, 55, 14) / pmin(dropoutTimes, 365), # based on median                             intercepts=c(GU=6.6, SI=5.8-1, SD=6.1, VU=4.9-.5, OA=3.2, ED=4.0, NU=2.5),                             sd.intercepts=c(GU=.81, SI=1.6-.1, SD=.95, VU=1.3, OA=1.6, ED=1.6, NU=1.4) * .667,                             slopes=c(GU=0, SI=.0058*3, SD=-.0038*5, VU=.0004*24, OA=-.003, ED=-.0014, NU=-.015),                             sd.slopes=c(GU=.0016, SI=.0031/2, SD=.0027/2, VU=.0032*0, OA=.0091, ED=.01, NU=.01),                             quads=c(GU=0, SI=-.00003, SD=.00003, VU=-.00003, OA=0, ED=-.0001, NU=-.0001),                             sd.quads=c(GU=0, SI=0, SD=0, VU=0, OA=0, ED=0, NU=0),                             vars=c(GU=2.0, SI=3.6, SD=3.2, VU=3.4, OA=3.6, ED=5.0, NU=3.0),                             sd.vars=c(.82, 1.3, .85, 1.2, 1.8, 2.6, 1.7),                             autocors=c(GU=.056, SI=.11, SD=.073, VU=.048, OA=.006, ED=-.044, NU=-.31),                             groupNames=c('Good users', 'Slow improvers', 'Slow decliners', 'Variable users', 'Occasional attempters', 'Early drop-outs', 'Non-users'),                             missing=FALSE,                             seed=1) {    set.seed(seed)    groupCounts = floor(patients * props)    incrIdx = order((patients * props) %% 1) %>%        head(patients - sum(groupCounts))    groupCounts[incrIdx] = groupCounts[incrIdx] + 1 # increment the groups that were closest to receiving another patient    assert_that(sum(groupCounts) == patients)    groupNames = factor(groupNames, levels=groupNames)
    # generate patient coefficients    groupCoefs = data.table(Group=groupNames, Patients=groupCounts,                            TDrop=dropoutTimes, Sd.TDrop=sd.dropoutTimes,                            AProb=attemptProbs,                            Int=intercepts, Sd.Int=sd.intercepts,                            Slope=slopes, Sd.Slope=sd.slopes,                            Quad=quads, Sd.Quad=sd.quads,                            Var=vars, Sd.Var=sd.vars,                            R=autocors)    assert_that(nrow(groupCoefs) == 7)    patCoefs = groupCoefs[, .(TDrop=rnorm(Patients, TDrop, Sd.TDrop) %>% round %>% pmax(7),                              AProb=AProb,                              Intercept=rnorm(Patients, Int, Sd.Int) %>% pmax(0),                              Slope=rnorm(Patients, Slope, Sd.Slope),                              Quad=rnorm(Patients, Quad, Sd.Quad),                           Variance=rnorm(Patients, Var, Sd.Var) %>% pmax(.75),                           R=rep(R, Patients)), by=Group]
    # generate patient measurements    genTs = function(N, Intercept, Slope, Quad, Variance, R, AProb, TDrop, ...) {        y = as.numeric(Intercept + times * Slope + times^2 * Quad + arima.sim(list(ar=R), n=length(times), sd=sqrt(Variance)))
        skipMask = !rbinom(length(times), size=1, prob=AProb)        y[skipMask] = 0
        if(missing) {            obsMask = times <= TDrop            list(Time=times[obsMask], Usage=pmax(y[obsMask], 0))        } else {            y[times > TDrop] = 0            list(Time=times, Usage=pmax(y, 0))        }    }
    patNames = paste0('P', 1:patients)    alldata = patCoefs[, do.call(genTs, .SD), by=.(Group, Id=factor(patNames, levels=patNames))] %>%        setkey(Id, Time)
    # generate group trajectories    groupTrajs = groupCoefs[, .(Time=times,                                Usage=(pmax(Int + times * Slope + times^2 * Quad, 0) * AProb) %>% ifelse(times > TDrop, 0, .)    ), by=Group]
    setattr(alldata, 'patCoefs', patCoefs)    setattr(alldata, 'groupCoefs', groupCoefs)    setattr(alldata, 'groupTrajs', groupTrajs)
    # Extra step: downsampling    if(is.finite(nAggr)) {        alldata = transformToAverage(alldata, binSize=nAggr)    }    return(alldata[])}
export_data = function(file, ..., seed = 1) {    data = generate_osa_data(..., seed = seed)    write.csv(data, file = file, quote = FALSE, row.names = FALSE)}

一、纵向数据分析研究目的

Puepose—任何研究，目的为王，围绕研究目的选择统计分析方法，纵向研究中笔者整理了以下常见研究目的：

（1）变化趋势与预测、影响因素的研究。纵向数据分析最根本的分析目的就是研究整个随访期间总体的平均变化趋势及个体增长趋势的差异。有的模型通过定义不同的方差-协方差结构来解决重复测量观测间具有的相关性，还有模型通过在模型中纳入随机截距和/或随机斜率来反映随机效应。进一步可以进行影响因素分析。可用模型有：多元方差分析(MANOVA)、重复测量方差分析(RM-ANOVA)、广义估计方程(GEE)、混合效应模型(MEM, 又叫“多水平模型”，根据因变量分为线性混合效应模型、广义线性混合效应模型、非线性混合效应模型)、潜变量增长曲线模型(LGCM)等。

（2）轨迹聚类分析的研究。研究目的①和②可以解决大部分纵向数据分析就问题，但不论哪种分析方法，前提都假设全部研究对象的发展趋势是相同的，然而这一假设往往并非总能满足，特别是以人为研究中心，其重点研究随时间推移个人行为、生物标志物或其他一些感兴趣现象的变化轨迹，总体往往具有较大的异质性(Heterogeneity)。

为了解决这一问题，不少学者开发了将变量分析和以人为中心分析相结合的方法，他们根据总体中个人的反应模式将个体划分为不同的轨迹组或类群体，一旦确定了子群体，轨迹成员可以作为因变量来确定健康轨迹的预测因素，或自变量来探索它们对未来健康结果的影响。它有助于探索健康状况的异质性，识别需要更好医疗服务的弱势群体，找到施加干预的最佳时机。可用模型有：组轨迹模型(GBTM，又叫潜类别增长模型)、潜变量增长混合模型(LGMM)、潜在转移分析(LTA)。

二、纵向数据特点及数据格式

（1）纵向数据：前面文章小编已经进行的详细的介绍！复习的同学点击：纵向数据分析—组轨迹模型(GBTM)：（一）模型介绍，纵向数据有哪些特点呢？

①时间序列性：随着时间的推移重复测量的反应变量可能发生变化。

②测量次数与测量时间间隔的非均衡性：在真正的纵向研究中，无论是实验性研究还是自然人群队列，个体在重复测量的次数上或者获得重复测量的时间上都有所不同。尽管这对于基于方差的传统分析方法是一个挑战（高度不平衡的设计，可能产生有偏差的参数估计和假设检验），但更现代的纵向数据分析方法几乎不受其影响。

③自相关性：是由于对受试者进行多次测量，同一观测单位的各测量之间存在自相关性，不同测量之间的相关系数不同，从而可以定义不同的相关矩阵。相邻数据点存在强自相关性，离得较远的数据随着时间间隔的增大，其相关性在减少。

④生态单位聚集性：纵向研究往往以几个生态单位（学校、教室、诊所、医院、县等）为研究对象进行抽样研究。生态单位内的个体聚集产生一个额外的相关性来源，尽管个体聚集产生相关的大小通常比在同一个体内重复测量产生的相关（自相关）小得多，但忽略这一重要的可变性来源会导致低估模型参数的标准误差，以及I型错误率（即假阳性率）升高。

⑤数据的不完整性：纵向研究中最显著的困难就是数据的不完整性。简单来说并非所有的受试者在整个研究期间都留在研究中，各种原因造成的失访，使得纵向数据完整性难以保证，故常常存在缺失。

⑥因果关联的时序性：对于横断面数据，其原因和结果都是同时存在的，往往无法判断其时序性。而纵向研究满足结果变量Y（效应）和协变量X（原因）之间有一个时间滞后，即原因必须发生在效应之前。

⑦变量类型与数据分布的复杂性：纵向研究往往要收集不同类型的变量，可以是连续型变量，也可以是离散型变量。具有维度高、变量多等特点，且往往不服从正态分布。

⑧线性与非线性：反应变量随自变量的变化趋势可能是线性的，也可能是非线性的。

⑨群体异质性：特别是在行为科学中，普遍存在个体发展趋势的差异。

（2）纵向数据数据整理格式：根据不用软件和模型的需要可整理成“宽形”或“长形”数据格式。例如：采用SAS构建组轨迹模型时需要“宽形”数据，而R则需要“长形”数据。

图 “宽形”数据

图 “长形”数据

三、针对纵向数据数据轨迹模型

（一）组轨迹模型(GBTM)

组轨迹模型(Group-based trajectory modelling, GBTM)，又叫潜类别增长模型(Latent Class Growth Modelling, LCGM)。Nagin于1999年对该方法进行了介绍，并随后将其定义为：有限混合模型的应用，使用轨迹组作为统计工具，用于近似人口成员的未知轨迹。轨迹组被定义为：随着时间的推移，在结果上遵循相似轨迹的个体组成的群。GBTM是一种半参数建模方法，其用于识别群体中具有相似发展轨迹的群/类。尽管每个个体都有一个独特的发展过程，但数据中个体间的异质性或分布可以通过一组有限且唯一的多项式函数来体现，每个多项式函数对应于一个离散的轨迹。具体来说通过多项式函数建立时间与潜变量之间的联系。GBTM可用于分析三种类型的资料:连续型数据(特别是心理测量数据)、计数数据和分类数据。适用于二项分布(Binomial)；删失正态分布(Censored normal)、泊松(Poisson)、零膨胀泊松(Zero-inflated Poisson)分布。

library(lcmm)library(splines)data = generate_osa_data() %>%    .[, Id := as.integer(Id)] %>%    .[, NormTime := (Time - min(Time)) / (max(Time) - min(Time))]
makeGbtmCall = function(k) {    substitute(        hlme(fixed=Usage ~ NormTime + I(NormTime^2),             mixture=~NormTime + I(NormTime^2),             random=~-1,             subject='Id', ng=k, data=data),        env=list(k=k)    )}
computeHlmeTrajectories = function(model) {    times = sort(unique(data$Time))    normTimes = sort(unique(data$NormTime))    predictY(model, newdata=data.frame(NormTime=normTimes))$pred %>%        data.table(Time=times) %>%        melt(id.vars='Time', value.name='Usage', variable.name='Group') %>%        .[, Group := factor(Group, levels=paste0('Ypred_class', 1:model$ng), labels=LETTERS[1:model$ng])] %>%        .[]}
# Single-group analysis ####mod00 = hlme(fixed=Usage ~ 1, random=~1, subject='Id', ng=1, data=data)summary(mod00)residuals(mod00) %T>% qqnorm %>% qqline
mod11 = hlme(fixed=Usage ~ NormTime, random=~NormTime, subject='Id', ng=1, data=data)summary(mod11)residuals(mod11) %T>% qqnorm %>% qqline
mod22 = hlme(fixed=Usage ~ poly(NormTime, 2, raw=TRUE), random=~poly(NormTime, 2, raw=TRUE), subject='Id', ng=1, data=data)summary(mod22)residuals(mod22) %T>% qqnorm %>% qqline
mod33 = hlme(fixed=Usage ~ poly(NormTime, 3, raw=TRUE), random=~poly(NormTime, 3, raw=TRUE), subject='Id', ng=1, data=data)summary(mod33)residuals(mod33) %T>% qqnorm %>% qqline
modbs = hlme(fixed=Usage ~ bs(NormTime), random=~bs(NormTime), subject='Id', ng=1, data=data)summary(modbs)residuals(modbs) %T>% qqnorm %>% qqline
# Estimation ####gbtms = list()# gbtms = readRDS('save/gbtm.rds')gbtms[['1']] = hlme(fixed=Usage ~ NormTime + I(NormTime^2), random=~-1, subject='Id', ng=1, data=data)
fitGbtm = function(k) {    start = Sys.time()    model = do.call(gridsearch, list(m=makeGbtmCall(k), rep=20, maxiter=1, minit=gbtms[['1']]))    model$runTime = Sys.time() - start    return(model)}
gbtms[['2']] = fitGbtm(2)gbtms[['3']] = fitGbtm(3)gbtms[['4']] = fitGbtm(4)gbtms[['5']] = fitGbtm(5)gbtms[['6']] = fitGbtm(6)gbtms[['7']] = fitGbtm(7)gbtms[['8']] = fitGbtm(8)
#! saveRDS(gbtms, file='save/gbtm.rds')
# Solutions ####plotMetric(sapply(gbtms, '[[', 'BIC'), as.integer(names(gbtms)), 'BIC')# ggsave('save/gbtm_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####k = 4bestGbtm = gbtms[[as.character(k)]]
pp = bestGbtm$pprob[paste0('prob', 1:k)] %>%    as.matrix() %>%    set_colnames(LETTERS[1:k])groupProps = colMeans(pp) %T>% print()
computeHlmeTrajectories(bestGbtm) %>% plotGroupTrajectories(groupProps)# ggsave('save/gbtm_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')
appa(pp)relativeEntropy(pp)

（二）潜变量增长混合模型(LGMM)

Muthen、Shedden等人于1999年提出潜变量增长混合模型(Latent growth mixture modelling ,LGMM)，也有学者称之为增长混合模型(Growth mixture modelling, GMM)。LGMM是有限混合模型(FMM)在潜变量增长曲线模型(latent growth curve modeling，LGCM)的扩展，LGCM在假设样本中所有个体具有相同的发展趋势下估计平均增长曲线，可以同时考虑个体和整体的变化情况，然而其假设总体具有同质性。LGMM通过在模型中引入分类潜变量来解决这一问题，分类潜变量将群体分成互斥的类来描述总体的异质性，每个类遵循传统的LGCM。不仅可以估计每个类的平均增长曲线，又可以估计每一类内个体的变化。它假设在任何给定的种群中，存在有限数量的未观测到的具有相似行为或经验的亚种群或类。

library(lcmm)data = generate_osa_data() %>%    .[, Id := as.integer(Id)] %>%    .[, NormTime := (Time - min(Time)) / (max(Time) - min(Time))]
makeGmmCall = function(k) {    substitute(        hlme(fixed=Usage ~ NormTime + I(NormTime^2),             mixture=~NormTime + I(NormTime^2),             random=~NormTime,             idiag=TRUE,             nwg=TRUE,             subject='Id', ng=k, data=data),        env=list(k=k)    )}
computeHlmeTrajectories = function(model) {    times = sort(unique(data$Time))    normTimes = sort(unique(data$NormTime))    predictY(model, newdata=data.frame(NormTime=normTimes))$pred %>%        data.table(Time=times) %>%        melt(id.vars='Time', value.name='Usage', variable.name='Group') %>%        .[, Group := factor(Group, levels=paste0('Ypred_class', 1:model$ng), labels=LETTERS[1:model$ng])] %>%        .[]}
# Estimation ####gmms = list()# gmms = readRDS('save/gmm.rds')gmms[['1']] = hlme(fixed=Usage ~ NormTime + I(NormTime^2), random=~NormTime, subject='Id', idiag=TRUE, ng=1, data=data)
fitGmm = function(k) {    set.seed(1)    start = Sys.time()    model = do.call(gridsearch, list(m=makeGmmCall(k), rep=20, maxiter=1, minit=gmms[['1']]))    model$runTime = Sys.time() - start    return(model)}
gmms[['2']] = fitGmm(2)gmms[['3']] = fitGmm(3)gmms[['4']] = fitGmm(4)gmms[['5']] = fitGmm(5)gmms[['6']] = fitGmm(6)gmms[['7']] = fitGmm(7)gmms[['8']] = fitGmm(8) #nwg=TRUE takes 60s per iter as opposed to 7s# saveRDS(gmms, file='save/gmm.rds')
# Solutions ####plotMetric(sapply(gmms, '[[', 'BIC'), as.integer(names(gmms)), 'BIC')# ggsave('save/gmm_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####k = 7bestGmm = gmms[[as.character(k)]]
pp = bestGmm$pprob[paste0('prob', 1:k)] %>%    as.matrix() %>%    set_colnames(LETTERS[1:k])groupProps = colMeans(pp) %T>% print()
computeHlmeTrajectories(bestGmm) %>% plotGroupTrajectories(groupProps)# ggsave('save/gmm_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')
appa(pp)relativeEntropy(pp)

（三）潜在转移分析(LTA)/潜在剖面分析(LPA)

首先，LTA和LPA都是用于研究个体在时间上的状态或类别变化的方法。它们都试图识别个体在不同时间点上可能存在的状态或类别，并分析个体在这些状态或类别之间的转移模式。不同之处在于，LTA关注的是个体在时间上的状态转移，即个体从一个状态转变到另一个状态的概率，而LPA关注的是个体在某个时间点上所处的状态的概率。

其次，LTA和LPA在模型估计和模型选择方面使用的统计方法也有些不同。LTA通常使用最大似然估计或贝叶斯估计来估计状态转移的概率，而LPA使用最大似然估计或贝叶斯估计来估计每个状态的概率和与指标之间的关系。在模型选择方面，LTA通常使用信息准则（如BIC、AIC）来比较不同模型，而LPA也常使用贝叶斯信息准则（BIC）来选择最优的模型。

LCA

学习潜在转移分析(LTA/LPA)之前，我们先学习潜在类别分析(Latent class analysis , LCA)。LCA通过潜在类别模型，用内在的潜在类别变量来解释外显的类别变量之间的关系，使得外显变量之间的关系经过潜在类别估计后，能够维持其局部独立性。LCA的基本假设是对各外显变量各种反应的概率分布可以由少数互斥的潜在类别变量变量来解释，每种类别对各外显变量的反应选择都有特定的倾向，LCA适用于多个分类变量。

LPA

潜在剖面分析,也被称为高斯混合模型或潜在类分析：是针对多个断面进行聚类的方法，其能够探索不同时间点之间观察对象的模式(组)变化。与先前LKMA不同的是，其聚类的准则是观察对象在各组的后验概率，哪一类更大即观察对象更可能属于哪个组，与LCA。它在拟合数据方面具有更大的灵活性。LPA在许多软件包中都可用，包括MPLUS、latent GOLD4和R包mclust5。此外SAS的PROC LCA也是可以做这个模型的，只是需要去github上加载SAS软件包。

library(mclust)data = generate_osa_data()dataMat = transformToRepeatedMeasures(data)
fitLlpa = function(k) {    cat('Fitting for k =', k, '\n')    start = Sys.time()    bestModel <<- NULL    for(n in 1:10) {        cat(sprintf('fit %d ', n))        model <<- Mclust(dataMat, G=k, modelNames='VVI',                         prior=priorControl(functionName='defaultPrior'),                         verbose=FALSE)
        if(!is.null(model)) {            if(is.null(bestModel) || model$bic > bestModel$bic) {                bestModel <<- model            }        } else {            cat('conv error')        }        cat('\n')    }    if(is.null(bestModel)) {        bestModel = list(bic=NA)    }    bestModel$runTime = Sys.time() - start    return(bestModel)}
computeLlpaTrajectories = function(model) {    times = sort(unique(data$Time))    model$parameters$mean %>%        set_rownames(NULL) %>%        set_colnames(LETTERS[seq_len(model$G)]) %>%        melt(varnames=c('Time', 'Group'), value.name='Usage') %>%        as.data.table() %>%        .[, Time := times[Time]] %>%        .[]}
# Estimation ####llpas = list()# llpas = readRDS('save/llpa.rds')llpas[['1']] = fitLlpa(1)llpas[['2']] = fitLlpa(2)llpas[['3']] = fitLlpa(3)llpas[['4']] = fitLlpa(4)llpas[['5']] = fitLlpa(5)llpas[['6']] = fitLlpa(6)llpas[['7']] = fitLlpa(7)llpas[['8']] = fitLlpa(8)
# saveRDS(llpas, file='save/llpa.rds')
# Solutions ####plotMetric(-1 * sapply(llpas, '[[', 'bic'), as.integer(names(llpas)), 'BIC')# ggsave('save/llpa_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####k = 5

LTA

潜在转移分析(Latent transition analysis, LTA)是LCA的扩展，也可以说是LCA的纵向版本，用于分析多个时点多个类别变量随时间的变化，如：是/否、轻度/中度/严重。LTA在每个时间点的数据经过LCA处理。该模型假设个人可以随着时间推移改变他们所在的类。例如：在三个疼痛强度组(轻度/中度/重度)中，LTA允许个体在一个时间点从重度组切换到下一个时间点的轻度或中度组，以此类推。故LTA分析的主要目的研究从一个时间点到下一个时间点的状态或行为的转移概率。

由于LTA和LCA模型针对于分析分类变量的纵向数据中的应用，不适合本次数据案例，这里仅仅展示常用的LPA，关于LCA/LPA/LTA的具体实现方式，将在后续本话题的推文中展示R代码和Mplus代码。

（四）时依性混合模型(TVEM)

TVEM方法最初由Rebecca A. Ferrer等人于2014年提出，该模型又被称为时变效应模型或动态广义线性模型，其发展和应用可以追溯到混合模型和时间序列模型在统计学和计量经济学领域的相关研究。通过结合这两个模型的优势，TVEM方法提供了一种强大而灵活的工具，用于研究长期数据中时间变化的效应。

TVEM方法结合了混合模型（mixture model）和时间序列模型（time series model）的思想。它假设观测数据是由若干个时间段组成的混合分布，其中每个时间段都具有特定的时间效应。是一种用于分析长期数据中时间变化效应的统计建模方法。其主要用于探索数据中的时间效应是否随着时间的推移而变化，并识别可能存在的变化点。需要注意的是TVEM方法还需要考虑数据的缺失、选择适当的时间单位和时间段的划分等问题。

注：目前install这个包可能有点问题需要去github自己下载包，然后Source函数引用，同样跑八个组的轨迹，这个模型跑了我四十分钟。

source('D:/DATA/研究PAPER/6月统计年鉴/demo-clustering-longitudinal-data-main/MixTVEM-master/MixTVEM.r')
# install.packages("MixTVEM")data = generate_osa_data() %>%    .[, NormTime := (Time - min(Time)) / (max(Time) - min(Time))]
fitMixTvem = function(k) {    start = Sys.time()    if(k == 1) {        # needed to prevent consistent estimation errors for the single-group solution        y = ifelse(data$Usage == 0, rnorm(nrow(data), sd=.1), data$Usage)    } else {        y = data$Usage    }    model = TVEMMixNormal(dep=y,                  id=data$Id,                  numInteriorKnots=6,                  deg=3,                  doPlot=FALSE,                  numClasses=k,                  numStarts=20,                  gridSize=365,                  maxVarianceRatio=NA,                  convergenceCriterion=1e-4,                  getSEs=FALSE,                  tcov=rep(1, nrow(data)),                  time=data$NormTime)    model$runTime = Sys.time() - start    return(model)}
computeMixTvemTrajectories = function(model) {    k = ncol(model$fittedValues)    times = model$timeGrid * (max(data$Time) - min(data$Time)) + min(data$Time)    data.table(model$betaByGrid[[1]], Time=times) %>%        setnames(paste0('V', 1:k), LETTERS[1:k]) %>%        melt(id.vars='Time', value.name='Usage', variable.name='Group')}
# mixtvems = readRDS(file='save/mixtvems.rds')mixtvems = list()mixtvems[['1']] = fitMixTvem(1)mixtvems[['2']] = fitMixTvem(2)mixtvems[['3']] = fitMixTvem(3)mixtvems[['4']] = fitMixTvem(4)mixtvems[['5']] = fitMixTvem(5)mixtvems[['6']] = fitMixTvem(6)mixtvems[['7']] = fitMixTvem(7)mixtvems[['8']] = fitMixTvem(8)# saveRDS(mixtvems, file='save/mixtvems.rds')
# Solutions ####plotMetric(sapply(mixtvems, function(m) m$bestFit$bic), as.integer(names(mixtvems)), 'BIC')# ggsave('save/mixtvem_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####k = 4bestMixTvem = mixtvems[[as.character(k)]]
pp = bestMixTvem$bestFit$postProbsBySubgroupProps = colMeans(pp) %T>% print()
computeMixTvemTrajectories(bestMixTvem) %>% plotGroupTrajectories(groupProps)# ggsave('save/mixtvem_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')
relativeEntropy(pp)confusionMatrix(pp) %>% {round(. * 100, .1)}

（五）纵向k-均值聚类(Lntiitudinal k-means clustering, LKMA)

纵向k-均值聚类（Longitudinal k-means clustering, LKMA）由斯坦福大学的Manda Wilson和Susan A. Murphy在2005年提出。他们将传统的k均值聚类算法扩展到了时间序列数据中，引入了时间维度的概念，考虑了样本点在不同时间点上的相互关联性。

传统k均值聚类算法主要用于对静态数据进行聚类，它假设数据中的样本点是相互独立的，没有时间上的关联。然而，在一些应用领域，例如医学、生物学和社会科学，数据往往是沿时间轴采集的，其具有时间上的相关性。LKMA的发展和应用使得在时间序列数据中进行聚类分析成为可能，为研究人员提供了一种新的工具来识别发展模式、预测趋势和理解变化的原因。它在医学中的应用尤为广泛，例如对疾病进展的监测、患者群体的分型等。

LKMA模型核心思想是将时间序列数据转换为具有固定维度的向量表示。该向量是由每个时间点上的特征值组成的。然后，使用传统的k-均值聚类算法对这些向量进行聚类。在聚类过程中，除了考虑特征值之间的空间距离，还考虑了时间维度上的距离。其优点为算法简单，收敛快，目前多用作其他轨迹建模前的探索性分析(半监督中聚类组数的探索)。每个簇的质心的选择是这个模型的关键，现目前提出了许多拓展模型，例如k-means++，Mini-Batch k-means，k-modoids等方法。

library(kml)data = generate_osa_data()
dataMat = transformToRepeatedMeasures(data)cld = clusterLongData(dataMat, idAll=rownames(dataMat), time=sort(unique(data$Time)), varNames='Usage')par = parALGO(startingCond='kmeans++')
computeKmlTrajectories = function(k) {    calculTrajMeanC(cld@traj, clust=getClusters(cld, bestK, asInteger=TRUE)) %>%        set_rownames(LETTERS[seq_len(nrow(.))]) %>%        melt(varnames=c('Group', 'Time'), value.name='Usage') %>%        as.data.table() %>%        .[, Time := cld@time[Time]]}
# Evaluate for all clusters ##### cld = readRDS('save/kml.rds')start = Sys.time(); kml(cld, nbClusters=1:8, nbRedrawing=20, parAlgo=par); runTime = Sys.time() - start# saveRDS(cld, file='save/kml.rds')
# Select best solution ####evalGroups = mapply(slot, list(cld), paste0('c', 1:26)) %>%    lengths() %>%    {which(. > 0)}kmls = mapply(slot, list(cld), paste0('c', evalGroups), SIMPLIFY=FALSE) %>%    lapply(first) %>%    set_names(evalGroups)
bics = -1 * sapply(kmls, function(part) part@criterionValues['BIC'])plotMetric(bics, evalGroups, 'BIC')+ expand_limits(y=6e4)# ggsave('save/kml_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####bestK = 7bestKml = kmls[[as.character(bestK)]]print(bestKml)
pp = bestKml@postProba %>%    set_colnames(levels(bestKml@clusters))groupProps =colMeans(pp) %T>% print()
computeKmlTrajectories(bestK) %>% plotGroupTrajectories(groupProps)# ggsave('save/kml_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')
relativeEntropy(pp)confusionMatrix(pp) %>% round(4)

（六）层次凝聚聚类(Agglomerative hierarchical clustering,AHC)

AHC的由来可以追溯到20世纪60年代。主要有两个独立的研究小组，分别由S. S. Rand和J. C. Dunn提出。S. S. Rand在1968年首次提出了AHC的概念，并且发表了名为《Objective Criteria for the Evaluation of Clustering Methods》的论文。J. C. Dunn则在1973年发表了《A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters》这篇论文。他对AHC方法进行了改进，提出了一种模糊划分的准则来解决一些难以分类的问题。

AHC与先前提到的LKMA们在原理、计算复杂度以及聚类结果等方面存在一些异同。

Item	AHC	K-means
原理	层次聚类方法	划分聚类方法
时间复杂度	较高	较低
簇个数确定	通过迭代合并样本簇，可自动确定不同数量的簇	需要事先确定簇的数量k，结果受初始簇中心的选择影响
数据分布适应性	适用于不同形状和分布的数据，可处理非球形簇和嵌套簇	对于分布不均匀、非球形簇或具有噪声的数据，聚类效果可能较差
聚类形式	可得到层次结构，提供了不同层次的聚类结果	每个样本只能分配到一个确定的簇中，无法提供层次结构

具体聚类过程大家看下图就很容易理解啦。

library(cluster)library(ggdendro)library(cluster)data = generate_osa_data()dataMat = transformToRepeatedMeasures(data)
computeMeanTrajectories = function(clusters) {    assert_that(is.factor(clusters), length(clusters) == uniqueN(data$Id))    data[, .(Usage=mean(Usage)), by=.(Group=clusters[as.integer(Id)], Time)]}
# Estimation ####start = Sys.time()D = dist(dataMat, method='euclidean')htree = hclust(D, method='average')
Ks = 1:8ahcsLabels = lapply(Ks, function(k) cutree(htree, k=k)) %>%    set_names(Ks)runTime = Sys.time() - start# Solutions ####dendro = dendro_data(htree, type='rectangle')ggplot(segment(dendro)) +    geom_segment(aes(x=x, y=y, xend=xend, yend=yend), size=.001) +    scale_x_continuous(breaks=pretty_breaks(4)) +    coord_flip() +    labs(x='Patient', y='Height') +    theme(panel.grid.major.y=element_blank(),          panel.grid.minor.y=element_blank())ggsave('save/ahc_tree.png', width=15, height=3, units='cm')
# Compute silhouette widthssils = lapply(ahcsLabels[as.character(setdiff(Ks, 1))], silhouette, dist=D) %>%    lapply(summary) %>%    sapply('[[', 'avg.width')
# Compute residual sum of squaresrsss = sapply(ahcsLabels, function(clusters) {    clusters = factor(clusters, labels=LETTERS[1:uniqueN(clusters)])    dtPatClus = data.table(Id=levels(data$Id), Group=clusters)    dtTraj = computeMeanTrajectories(clusters)    dtFitted = data[dtPatClus, on='Id'] %>%        merge(dtTraj, by.x=c('i.Group', 'Time'), by.y=c('Group', 'Time'))    dtFitted[, sum((Usage.x - Usage.y)^2)]})
plotMetric(c(NA, sils), Ks, 'Silhouette width')# ggsave('save/ahc_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
plotMetric(rsss, Ks, 'RSS')
# Assess result ####k = 4clusters = ahcsLabels[[as.character(k)]] %>%    factor(labels=LETTERS[1:k])groupProps = table(clusters) %>% prop.table() %>% as.numeric()
computeMeanTrajectories(clusters) %>% plotGroupTrajectories(groupProps)# ggsave('save/ahc_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')

先前研究均为半监督聚类的方法，下面我们介绍一下另一种轨迹聚类思路，即全监督聚类方法，基于轨迹某些参数的特征对轨迹分类，从而达到聚类的目的。

四、全监督聚类——基于特征的聚类

在基于特征的聚类方法中，使用参数模型来描述每个轨迹的关键特征。在其最简单的形式中，轨迹特征可以用汇总统计数据来表示，如平均值、标准差、偏度、范围、平稳度、周期性、自相关或熵。其优点是容易拟合，缺点是未考虑轨迹分组异质性问题，拟合效果可能不如前述模型。

（一）个体时间序列分析(Individual time series representations)

Bushway 等人在2009年提出ITS分析在轨迹聚类中的应用，当我们想要描述一个轨迹时，可以用许多不同的方式来表示它。其中一种直观的方法是使用时间相关的线性模型来描述每个轨迹，这个模型可以包含一些参数，比如截距和斜率。在多层次建模中，我们可以估计每个轨迹的随机效应，并使用聚类算法，比如k-means，对这些轨迹进行聚类。

这种方法特别适用于每个轨迹的观测数据较少的情况。通过独立地估计每个轨迹的表示，我们可以更好地适应人群的异质性，也可以得到更好的拟合效果。这一方法被称为个体时间序列分析，已经在过去的研究中广泛应用过。简单来说，我们可以用这种方法来更好地理解和描述不同轨迹的特征。

library(ggdendro)library(MASS)library(cluster)data = generate_osa_data() %>%    .[, NormTime := (Time - min(Time)) / (max(Time) - min(Time))]
computeMeanTrajectories = function(clusters) {    assert_that(is.factor(clusters), length(clusters) == uniqueN(data$Id))    data[, .(Usage=mean(Usage)), by=.(Group=clusters[as.integer(Id)], Time)]}
# Compute patient representation ####{    start = Sys.time()    patientFeatures = data[, {        mod = lm(Usage ~ poly(NormTime, 2), data=.SD)        as.list(coef(mod)) %>% c(logN=log(sum(Usage > 0)), sigma=sigma(mod))        }, keyby=Id]
    X = as.matrix(patientFeatures[, -c('Id', 'sigma')]) %>%        set_rownames(patientFeatures$Id) %>%        scale()
    # Hierarchical clustering    D = dist(X, method='euclidean')
    htree = hclust(D, method='average')
    Ks = 1:8    ahcs = lapply(Ks, function(k) pam(D, k=k)) %>%        set_names(Ks)    sils = lapply(ahcs, function(ahc) ahc$silinfo$avg.width) %>%        sapply(function(x) ifelse(is.null(x), NA, x))    runTime = Sys.time() - start}max(sils, na.rm=TRUE)which.max(sils)
# Solutions ####dendro = dendro_data(htree, type='rectangle')
ggplot(segment(dendro)) +    geom_segment(aes(x=x, y=y, xend=xend, yend=yend)) +    scale_x_continuous(breaks=pretty_breaks(10)) +    coord_flip() +    labs(x='Trajectories\n\n', y='Height') +    theme(panel.grid.major=element_blank(),          panel.grid.minor=element_blank())
plotMetric(sils, Ks, 'Silhouette width')# ggsave('save/featurebased_bic.pdf', width=bicPlotSize[1], height=bicPlotSize[2], units='cm')
# Assess the best solution ####k = 7bestAhc = ahcs[[as.character(k)]]clusters = factor(bestAhc$clustering, labels=LETTERS[1:k])groupProps = table(clusters) %>% prop.table() %>% as.numeric()
computeMeanTrajectories(clusters) %>% plotGroupTrajectories(groupProps)# ggsave('save/featurebased_groups.pdf', width=groupPlotSize[1], height=groupPlotSize[2], units='cm')

（二）界标模型/标段模型(Landmark Models,LME)

界标模型用于动态预测个体的事件风险，是一种时间依赖的轨迹模型。其最大的特点是对个体可以进行个性化风险预测。界标模型背后的关键思想是根据研究目的事先定义具有标志性的里程碑事件发生的时间。所以这个模型并不是严格意义上的聚类模型，但是可以实现聚类模型想达到的目的，即找到个体的轨迹的规律，后进行聚类分析。标志性时间是我们希望进行风险预测的时间点（即“landmark”）。在标志性时间点上分别建立模型。

例如下图我们想要预测45到50岁之间任何年龄的心血管疾病（CVD）风险，我们定义了六个分段年龄（45、…、50岁），并为每个年龄拟合了六个生存分析模型。

总结：

在纵向数据分析中都涉及到许多和专业相关的决策。对科学研究的结论需要的不仅是模型驱动的结果，还应该基于专业驱动，使我们的研成为可重复的，可指导实践的研究。

篇幅有限，后续我将持续更新。

文字编辑：想不出吸睛好名字、天涯二毛君

程序编辑：想不出吸睛好名字

审稿：老陈