【原】power analysis

这几天很多外院的朋友，问我纵向数据如何做功效分析，其实很简单，要明白功效分析的基本原理，大概分成三步：1/模拟数据生成：该过程包括生成模拟的纵向数据和生存时间数据。这里包括了年龄和性别两种协变量，纵向数据响应变量 y 基于时间和年龄的线性函数生成，并添加正态随机误差。生存时间数据则基于年龄和性别的线性组合通过指数分布生成，同时处理了右删失。

2/模型拟合：使用 jmbayes2 包拟合一个联合模型，包含纵向子模型（响应变量 y 关于年龄和时间的关系）和生存子模型（生存时间关于性别的影响）。

3/统计功效评估：通过重复模拟数据生成和模型拟合过程多次（如1000次），计算在这些模拟研究中，能够以 p 值小于 0.05 检测到性别效应的比例。这个比例就是模型在当前设定下的统计功效。

# 定义模拟数据生成函数

generate_data <- function(n, beta_long, beta_surv) {

    # 随机生成一些协变量

    age <- rnorm(n, 50, 10)

    sex <- rbinom(n, 1, 0.5)

    # 生成纵向数据的响应

    times <- seq(0, 5, length.out = 10)

    y <- outer(age, times, function(age, t) 0.5 * age + beta_long * t) +

         matrix(rnorm(n * length(times), sd = 2), ncol = length(times))

    # 生成生存时间数据

    lambda <- exp(0.04 * age + 0.2 * sex)

    survival_time <- rexp(n, lambda)

    censoring_time <- 5

    event <- ifelse(survival_time <= censoring_time, 1, 0)

    observed_time <- pmin(survival_time, censoring_time)

    return(list(age = age, sex = sex, y = y, time = observed_time, status = event))

}

# 模拟一次数据

sim_data <- generate_data(100, beta_long = 0.1, beta_surv = 0.2)

# 拟合模型

fit <- jmbayes2(fixed = y ~ age + time, random = ~ 1 | id, data = sim_data,

                survival = Surv(time, status) ~ sex, data = sim_data)

# 检查模型结果

summary(fit)

# 设置模拟次数

n_sim <- 1000

power_count <- 0

n <- 100  # 初始样本大小

for (i in 1:n_sim) {

    sim_data <- generate_data(n, beta_long = 0.1, beta_surv = 0.2)

    fit <- jmbayes2(fixed = y ~ age + time, random = ~ 1 | id, data = sim_data,

                    survival = Surv(time, status) ~ sex, data = sim_data)

    if (summary(fit)$coefficients[2, "p-value"] < 0.05) {

        power_count <- power_count + 1

    }

}

power <- power_count / n_sim

print(paste("Estimated power:", power))

Estimated power：这个值表示在当前模型设定和样本大小下，检测到性别对生存时间有显著影响的概率。例如，如果得到的功效是 0.8，这意味着在重复相同的研究100次中，有80次能够显著检测到性别的效应。