## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE
)
library(Desimml)
# 设置随机种子，保证结果的严格可重复性
set.seed(2026)

# 生成模拟数据: 样本量 n = 200, 协变量维度 p = 5
# family = "gaussian" 表示生成连续型高斯分布的响应变量
dat <- generate.data(
  n = 200, 
  p = 5, 
  family = "gaussian", 
  sigma = 0.4, 
  pi.1 = 0.5, 
  s = 2,      # 交互效应类型：非线性 (s=2 代表指数平方形式的非线性交互)
  delta = 1   # 主效应强度
)

# 检查生成数据的基本结构和维度
dim(dat$X)
table(dat$A)
head(dat$y)
# 拟合 SIMML 模型
# bs = "ps" 指定使用 P-splines
# k = 8 指定样条基底函数的自由度/维度
fit_res <- simml(
  y = dat$y, 
  A = dat$A, 
  X = dat$X, 
  family = "gaussian",
  bs = "ps",  
  k = 8       
)

# 输出最终收敛估计的单指数系数向量 beta
# 这些系数量化了 X 中的各个特征在决定最佳治疗方案时的权重
print(fit_res$beta.coef)
# 计算链接函数的一阶导数
# eps = 10^(-6) 控制数值偏导计算的微小步长
derivative_vals <- der.link(fit_res$g.fit, eps = 10^(-6))

# 查看部分观测样本对应的导数值
head(derivative_vals)
# 提取前 5 个观测样本的协变量作为新患者数据演示
new_X <- dat$X[1:5, ]

# 进行预测，假设目标是最大化响应变量值 (maximize = TRUE)
predictions <- pred.simml(
  simml.obj = fit_res, 
  newX = new_X, 
  type = "response", 
  maximize = TRUE
)

# 查看针对每个可能治疗方案的预测期望结果 
# 矩阵大小: n_new x L，其中每列代表一种治疗方案的预期 outcome
print(predictions$pred.new)

# 查看为这 5 个患者推荐的最优治疗方案 (即 pred.new 中每行最大值对应的治疗分配)
print(predictions$trt.rule)