捐助
  1. 贪吃的夜猫子Home
  2. 信息技术
  3. 统计分析

组别轨迹模型分析（Group-based Trajectory Model, GBTM）(SUN Bin)

xujunzju • • 统计分析 • 5 views 
library(dplyr)
library(gbmt)
library(tidyr)


####模拟1000组轨迹数据，5个时间点，轨迹分成三类：300例一次曲线,300例2次曲线，
####400例3次曲线。函数分别是y=-0.25t+20、y=(5/16)*(t-4)^2+15、y=-(5/32)*(t-4)^3+10
####
time_points <- seq(0, 4, length.out = 5)
# 生成符合时间的1次幂的轨迹数据
linear_trajectories <- matrix(NA, nrow = 300, ncol = 6)
for (i in 1:300) {
  linear_trajectories[i, 1] <- "Linear"
  linear_trajectories[i, 2:6] <- -0.25*time_points +20+ 1.2*rnorm(5)
}
# 生成符合时间的2次幂的轨迹数据
quadratic_trajectories <- matrix(NA, nrow = 300, ncol = 6)
for (i in 1:300) {
  quadratic_trajectories[i, 1] <- "Quadratic"
  quadratic_trajectories[i, 2:6] <- (5/16)*(time_points-4)^2 +15+ 1.2*rnorm(5)
}
# 生成符合时间的3次幂的轨迹数据
cubic_trajectories <- matrix(NA, nrow = 400, ncol = 6)
for (i in 1:400) {
  cubic_trajectories[i, 1] <- "Cubic"
  cubic_trajectories[i, 2:6] <- -(5/32)*(time_points-4)^3 +10+ 1.2*rnorm(5)
}
# 将所有数据合并
trajectory_data <- rbind(linear_trajectories, quadratic_trajectories, cubic_trajectories)
# 将数据转换成数据框
trajectory_df <- as.data.frame(trajectory_data)
# 设置列名
colnames(trajectory_df) <- c("Trajectory", paste0("Time", 1:5))
data_with_id <- trajectory_df %>% 
  mutate(ID = seq(1, n()))



###绘制模拟数据的轨迹线图
data_with_id<-as.data.frame(data_with_id)

vars1 <- c("Time1", "Time2", "Time3", "Time4","Time5")
data_with_id[vars1] <- lapply(data_with_id[vars1], as.numeric)
average_scores <- colMeans(data_with_id[,c("Time1", "Time2", "Time3", "Time4","Time5")])
average_scores_group_A <- data_with_id %>% 
  filter(Trajectory == "Linear") %>% 
  summarize(
    Time1_mean = mean(Time1, na.rm = TRUE),
    Time2_mean = mean(Time2, na.rm = TRUE),
    Time3_mean = mean(Time3, na.rm = TRUE),
    Time4_mean = mean(Time4, na.rm = TRUE),
    Time5_mean = mean(Time5, na.rm = TRUE)
  )
average_scores_group_B <- data_with_id %>% 
  filter(Trajectory == "Quadratic") %>% 
  summarize(
    Time1_mean = mean(Time1, na.rm = TRUE),
    Time2_mean = mean(Time2, na.rm = TRUE),
    Time3_mean = mean(Time3, na.rm = TRUE),
    Time4_mean = mean(Time4, na.rm = TRUE),
    Time5_mean = mean(Time5, na.rm = TRUE)
  )
average_scores_group_C<- data_with_id %>% 
  filter(Trajectory == "Cubic") %>% 
  summarize(
    Time1_mean = mean(Time1, na.rm = TRUE),
    Time2_mean = mean(Time2, na.rm = TRUE),
    Time3_mean = mean(Time3, na.rm = TRUE),
    Time4_mean = mean(Time4, na.rm = TRUE),
    Time5_mean = mean(Time5, na.rm = TRUE)
  )
##不区分分组的线图
plot(1, type="n", xlim=c(1, 5), ylim=range(data_with_id[,c("Time1", "Time2", "Time3", "Time4","Time5")]), 
     xlab="Time", ylab="Symptom Score", main="Symptom Scores Over Time")
# 添加每个样本的症状分线
for(i in 1:nrow(data_with_id)) {
  lines(1:5, data_with_id[i, c("Time1", "Time2", "Time3", "Time4","Time5")], col="grey")
}
lines(1:5, average_scores, col="black", lwd=2, lty=1)

##区分分组的线图
plot(1, type="n", xlim=c(1, 5), ylim=range(data_with_id[,c("Time1", "Time2", "Time3", "Time4","Time5")]), 
     xlab="Time", ylab="Symptom Score", main="Symptom Scores Over Time")
# 添加每个样本的症状分线
for(i in 1:nrow(data_with_id)) {
  lines(1:5, data_with_id[i, c("Time1", "Time2", "Time3", "Time4","Time5")], col="grey")
}
lines(1:5, average_scores_group_A, col="#CCCCFF", lwd=2, lty=1)
lines(1:5, average_scores_group_B, col="#CCCCCC", lwd=2, lty=1)
lines(1:5, average_scores_group_C, col="#CCCC99", lwd=2, lty=1)


plot(1, type="n", xlim=c(1, 5), ylim=range(data_with_id[,c("Time1", "Time2", "Time3", "Time4","Time5")]), 
     xlab="Time", ylab="Symptom Score", main="Symptom Scores Over Time")
# 添加每个样本的症状分线
for(i in 1:nrow(data_with_id)) {
  lines(1:5, data_with_id[i, c("Time1", "Time2", "Time3", "Time4","Time5")], 
        col=ifelse(data_with_id$Trajectory[i]=="Linear", "#CCCCFF", ifelse(data_with_id$Trajectory[i]=="Quadratic", "#CCCCCC", "#CCCC99")))
}
lines(1:5, average_scores_group_A, col="#333366", lwd=2, lty=1)
lines(1:5, average_scores_group_B, col="#333366", lwd=2, lty=1)
lines(1:5, average_scores_group_C, col="#333366", lwd=2, lty=1)


###将宽型数据转成长型数据，以便使用gbmt包进行分析。
data_long_simu <- tidyr::gather(data_with_id, key = "Time", value = "Value", -ID,,-Trajectory)
data_processed_simu <- data_long_simu %>% 
  mutate(Time = recode(Time, "Time1" = 0, "Time2" = 1, "Time3" = 2,"Time4"=3,"Time5"=4) %>% as.numeric())
df_simu<-data_processed_simu
df_simu$Value <- as.numeric(df_simu$Value)
###使用gbmt包，分别建立1组、2组、3组、4组的3阶模型，scaling=0表示Y不做任何转换。
varNames<-c("Value")
m1_3_simu <- gbmt(x.names=varNames, unit="ID", time="Time", d=3, ng=1, data=df_simu,scaling=0)
m2_3_simu <- gbmt(x.names=varNames, unit="ID", time="Time", d=3, ng=2, data=df_simu,scaling=0)
m3_3_simu <- gbmt(x.names=varNames, unit="ID", time="Time", d=3, ng=3, data=df_simu,scaling=0)
m4_3_simu <- gbmt(x.names=varNames, unit="ID", time="Time", d=3, ng=4, data=df_simu,scaling=0)
summary(m1_3_simu)
summary(m2_3_simu)
summary(m3_3_simu)
summary(m4_3_simu)
#分别显示四组的AIC、BIC、CIC等指标
ICs<-rbind(m1_3_simu$ic,m2_3_simu$ic, m3_3_simu$ic,m4_3_simu$ic)
ICs<-as.data.frame(ICs)
#分别显示四组的average posterior probability.
data1 <- m1_3_simu$appa
data2 <- m2_3_simu$appa
data3 <- m3_3_simu$appa
data4 <- m4_3_simu$appa

max_length <- max(length(data1), length(data2), length(data3), length(data4))
traj_1 <- c(data1, rep(NA, max_length - length(data1)))
traj_2 <- c(data2, rep(NA, max_length - length(data2)))
traj_3 <- c(data3, rep(NA, max_length - length(data3)))
traj_4 <- c(data4, rep(NA, max_length - length(data4)))


avePP <- data.frame(traj_1,traj_2, traj_3, traj_4)
avePP_t <- t(avePP)
colnames(avePP_t) <- c("G1", "G2", "G3","G4")
avePP_t<-as.data.frame(avePP_t)
avePP_t$model <- row.names(avePP_t)
row.names(avePP_t) <- NULL
avePP_t<-avePP_t[,c(5,1,2,3,4)]

##把信息准则和平均后验概率合并
table_simu<-cbind(avePP_t,ICs)
print(table_simu)
##绘制四种轨迹模型的曲线
layout(matrix(1:4, nrow=2, byrow=TRUE))
plot(m1_3_simu,bands = F)
plot(m2_3_simu,bands = F)
plot(m3_3_simu,bands = F)
plot(m4_3_simu,bands = F)
##计算OCC
APP_3<-m3_3_simu$appa
posteri_3<-m3_3_simu$prior
l_ap_3<-(APP_3)/(1-APP_3)
l_pos_3<-posteri_3/(1-posteri_3)
OCC_3<-l_ap_3/l_pos_3



library(haven)
trajdata <- read_sav("trajdata.sav")
library(gbmt)
library(tidyr)
library(trajeR)
df<-trajdata
df<-as.data.frame(df)
varNames <- c("HAMD")
#分别建立1组、2组、3组、4组的3阶模型，scaling=0表示Y不做任何转换。
m1_3 <- gbmt(x.names=varNames, unit="id", time="time", d=3, ng=1, data=df,scaling=0)
m2_3 <- gbmt(x.names=varNames, unit="id", time="time", d=3, ng=2, data=df,scaling=0)
m3_3 <- gbmt(x.names=varNames, unit="id", time="time", d=3, ng=3, data=df,scaling=0)
m4_3 <- gbmt(x.names=varNames, unit="id", time="time", d=3, ng=4, data=df,scaling=0)
m5_3<- gbmt(x.names=varNames, unit="id", time="time", d=3, ng=5, data=df,scaling=0)

#分别显示四种分组的结果描述
summary(m1_3)
summary(m2_3)
summary(m3_3)
summary(m4_3)
summary(m5_3)
#分别显示三组的AIC、BIC、CIC等指标
ICs<-rbind(m1_3$ic,m2_3$ic, m3_3$ic,m4_3$ic,m5_3$ic)
ICs<-as.data.frame(ICs)
#分别显示三组的average posterior probability.
data1 <- m1_3$appa
data2 <- m2_3$appa
data3 <- m3_3$appa
data4 <- m4_3$appa
data5 <- m5_3$appa
max_length <- max(length(data1), length(data2), length(data3), length(data4), length(data5))
traj_1 <- c(data1, rep(NA, max_length - length(data1)))
traj_2 <- c(data2, rep(NA, max_length - length(data2)))
traj_3 <- c(data3, rep(NA, max_length - length(data3)))
traj_4 <- c(data4, rep(NA, max_length - length(data4)))
traj_5 <- c(data5, rep(NA, max_length - length(data5)))

avePP <- data.frame(traj_1,traj_2, traj_3, traj_4,traj_5 )
avePP_t <- t(avePP)
colnames(avePP_t) <- c("G1", "G2", "G3","G4","G5")
avePP_t<-as.data.frame(avePP_t)
avePP_t$model <- row.names(avePP_t)
row.names(avePP_t) <- NULL
avePP_t<-avePP_t[,c(6,1,2,3,4,5)]

##把信息准则和平均后验概率合并
table<-cbind(avePP_t,ICs)
print(table)

##计算OCC
#三轨迹的OCC
APP3<-m3_3$appa
posteri3<-m3_3$prior
l_ap3<-(APP3)/(1-APP3)
l_pos3<-posteri3/(1-posteri3)
OCC3<-l_ap3/l_pos3
#四轨迹的OCC
APP4<-m4_3$appa
posteri4<-m4_3$prior
l_ap4<-(APP4)/(1-APP4)
l_pos4<-posteri4/(1-posteri4)
OCC4<-l_ap4/l_pos4

##绘制轨迹图——三轨迹
layout(1)
plot(m3_3,bands=F,xlab="time",ylab = "HAMD total score",titles = "HAMD scale",transparency=90,
     add.grid = TRUE,add.legend = FALSE,xaxt = "n",col=c("#CCCCFF", "#CCCCCC","#CCCC99"))
axis(side = 1, at = c(1,2,3,4,5), labels = c("baseline", "week1","week2","week3","week6"), las = 1)
legend("topright", legend = c("trajectory 1", "trajectory 2","trajectory 3"), 
       col = c("#CCCCFF", "#CCCCCC","#CCCC99"), lty = 1, lwd = 2, bty = "n")
##绘制轨迹图——四轨迹
layout(1)
plot(m4_3,bands=F,xlab="time",ylab = "HAMD total score",titles = "HAMD scale",transparency=90,
     add.grid = TRUE,add.legend = FALSE,xaxt = "n",col=c("#CCCCFF", "#CCCCCC","#CCCC99","purple"))
axis(side = 1, at = c(1,2,3,4,5), labels = c("baseline", "week1","week2","week3","week6"), las = 1)
legend("topright", legend = c("trajectory 1", "trajectory 2","trajectory 3","trajectory 4"), 
       col = c("#CCCCFF", "#CCCCCC","#CCCC99","purple"), lty = 1, lwd = 2, bty = "n")

# 使用pivot_wider函数将长型数据转换为宽型数据
df_HAMD<-df
df_wide<-pivot_wider(df_HAMD,names_from=time,values_from=c("SDS","SAS","PSQI","HAMD","BeckD" ))
##将轨迹分析后的轨迹分组合并到数据中
df_wide$trajetory<-m4_3$assign
###导出数据为excel格式
library(xlsx)
write.xlsx(df_wide, file = "df_wide_traj1.xlsx")

原创文章（本站视频密码：66668888），作者：xujunzju，如若转载，请注明出处：https://zyicu.cn/?p=21550

Like (0)
Donate 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
xujunzju的头像xujunzju管理者
0 0
Previous 2026年2月2日 15:00
Next 2022年4月8日 21:13

相关推荐

发表回复

Please Login to Comment
联系我们
邮箱：
xujunzju@gmail.com
公众号：
xujunzju6174
捐赠本站
捐赠本站
SHARE
TOP