scitable包复现一篇charl数据库4分马尔科夫链文章

原创于 2026-06-22 09:07:50 发布 · 190 阅读

3 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#数据库

R语言同时被 2 个专栏收录

401 篇文章

订阅专栏

CHARLS数据库挖掘系列教程

41 篇文章

订阅专栏

马尔科夫链既往我在系列文章《手把手带你复现一篇一区9.1分肌少症和马尔科夫链》的文章已经有过介绍，它的特点是：“未来只取决于现在，与过去无关”。主要研究的是状态和概率的变化。比如我在手机打字：你nih，会出现很多选项：你好，你和，你还等等，这就是软件统计了咱们打各个字符的概率，这个就是马尔科夫链的一个应用。自然语言处理：用马尔科夫链生成文本（比如根据前一个词预测下一个词）

在这里插入图片描述
本期带大家复现一篇4分马尔科夫链文章，就是下篇这篇，为什么想复现这篇呢？因为作者用了一个多模型结构来分析问题，和咱们的回归模型有点相似，引起了我的兴趣

在这里插入图片描述

我们先来看看这篇文章说什么,作者指出：肌少症和虚弱是密切相关的老年综合征，但肌少症如何与虚弱的动态进展相关仍不清楚。这项研究调查了肌少症与中国老年人多状态虚弱转变之间的关系。因此，作者研究的是肌少症和虚弱之前的联系。

这里先说明一下我是模拟一个数据，复现的主要是作者的方法，最主要的就是图3，其他的既往文章都有说过，我不是按步骤收集数据来复现，为什么不收集数据复现呢？这个问题文章末尾再说

先导入数据

library(scitable)
library(ggscitable)
library(tidyverse)
setwd("E:/公众号文章2026年/中国老年人肌少症与虚弱转变的关联：一项多状态马尔可夫模型研究")
markov_data<-read.csv("markov_data.csv",sep=',',header=TRUE)
names(markov_data)

在这里插入图片描述

“ID” 就是患者ID，，“time”，是时间，"state"就是患者的状态，"X1-X4"是协变量，其中X4是分类变量，X4刚好是ABC，咱们用来表示肌少症。

下面咱们正式开始,基线表我就不做了，先做图2，图二是啥，就是不同年份每个状态的转移概率，注意这里不是转移强度。
咱们先生成基础的马尔科夫链

out<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,
              username=username,token=token)
##转移强度
out[["tb2"]]

在这里插入图片描述
咱们注意一下out[[“tb2”]]这里是它的转移强度，不是概率，概率需要从这个基础的马尔科夫链方程生成，下面这个注意一下，time = c(1,2,3,4)表示生成1.2.3.4年的转移概率

probability.msm(out,time = c(1,2,3,4))

在这里插入图片描述
这样图2的结果就生成了，不过图需要自己画了，使用PPT等软件画起来还是很简单。下面咱们来画图3，这个是整个文章的核心图，咱们先要看懂它表达的是什么意思。

在这里插入图片描述

作者协变量肌少症作为关联，建立了3个模型，一个未调整模型（模型1），一个调整了人口学和社会经济因素的模型（模型2：年龄、性别、教育程度、婚姻状况、居住地），以及一个进一步调整了生活方式因素的完全调整模型（模型3：吸烟、饮酒、BMI）。

这个怎么做呢，下面咱们来演示一下，假设X4这里就表示肌少症，
先做未调整模型

out1<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,cov="X4",
            username=username,token=token)
out1<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,cov="X4",
            username=username,token=token)
bc2fit<-out1[["fit"]]
summary(bc2fit)

在这里插入图片描述
咱们可以看下X4B，X4C就是模型1的结果了，接下来咱们来做模型2.模型是：微调整模型

out2<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,cov=c("X4","X1"),
              username=username,token=token)
bc2fit<-out2[["fit"]]
summary(bc2fit)

在这里插入图片描述
接下来就是全模型

out3<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,cov=c("X4","X1","X2"),
              username=username,token=token)
bc2fit<-out3[["fit"]]
summary(bc2fit)

在这里插入图片描述

这个图的结果就全部出来啦，下面这个图其实不是真正意义的森林图

在这里插入图片描述

我们看到上图没有P值，也没有交互的P值，这个不算真正意义上的森林图，其实就是把数据对变量取亚组，然后再重复一下图2的操作，这里我就不做了，

但是咱们还可以扩展一下思维，比如在图1中咱们已经算出每个年份的转移概率

###########扩展
out4<-sci_msm(data=markov_data,ID="ID",time = "time",state = "state",digits = 7,cov="X4",
              username=username,token=token)

out5<-probability.msm(out4,time = c(1,2,3,4),state = 1,cov="X4")

然后把概率和年份提取出来

#提取数据
data5<-out5[["data"]]
data5$cov<-as.factor(data5$cov)
data5$cov<-factor(data5$cov,labels = c("n", "m","l"))
names(data5)<-c("State11","State12","State13","State14","t","cov")
colors <- c("n" = "#FFD1DC", "m" = "#FF4D4D","l"="blue")
data5$cov<-as.factor(data5$cov)

在这里插入图片描述
比如我想了解不同年份状态1到状态2的概率情况

library(ggplot2)

ggplot(data5, aes(x = t, y = State12, group = cov, color = cov)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  # 关键代码：按 cov 变量分面
  facet_wrap(~ cov) +
  labs(title = "按协变量(cov)分面展示", x = "时间(t)", y = "状态值(State12)") +
  theme_minimal()