将含样条效应的SAS混合模型转换为R实现的技术咨询

1. 原SAS模型代码 #

proc hpmixed data = eplong  noclprint noitprint ;
    class ID;
    effect spl = spline(GA / details naturalcubic basis=tpf knotmethod=percentiles(3));
    model TAC = spl / s;  
    random int spl / type=un subject = ID s;  
run;

2. 测试数据集（前48行） #

post.df <- structure(list(ID = structure(c("8366", "8366", "8366", 
"8366", "8367", "8367", "8367", "8367", "8368", "8368", "8368", 
"8368", "8369", "8369", "8369", "8369", "8370", "8370", "8370", 
"8370", "8371", "8371", "8371", "8371", "8372", "8372", "8372", 
"8372", "8374", "8374", "8374", "8374", "8393", "8393", "8393", 
"8393", "8394", "8394", "8394", "8394", "8395", "8395", "8395", 
"8395", "8396", "8396", "8396", "8396"), format.sas = "$"), TAC = c(79.7, 
89.4, 82.6, 93, 79.7, 123.8, 100.4, 83.9, 96.3, 80, 81.8, 87.8, 
34.4, 91.2, 80.6, 76.8, 84.8, 94.1, 81.9, 82.4, 49.4, 46.7, 54, 
NA, 86.8, 76.6, 77.4, 79.2, 85.4, 93.1, 83.7, 84.3, 100.5, 82.3, 
86.5, 81.1, 95.3, 96.5, 98.5, 94.8, 87, 87.8, 82.9, 86.1, NA, 
35.8, 43.6, NA), GA = c(11.86, 19.14, 27.86, 38.71, 12.71, 17.86, 
25.71, 35.57, 13.29, 23.29, 28.29, 37.29, 11, 17.86, 25.29, 35.14, 
12, 15.86, 24.14, 35.86, 13.86, 22.86, 30, 39.43, 13.14, 19.14, 
25.14, 34.14, 11.71, 20, 27, 36.29, 13.86, 24.71, 27.29, 35.71, 
13.71, 25.57, 28, 36, 13.86, 18.43, 27.14, 35, 12.29, 21.29, 
31.29, 38.29), visit = c(0, 1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4, 
0, 1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4, 0, 
1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4, 0, 1, 2, 4)), row.names = c(NA, 
-48L), class = c("tbl_df", "tbl", "data.frame"))

3. 尝试的R代码及问题 #

• 首次尝试代码：

  lmer(TAC ~ 1 + (bs(GA, df=3) | ID), data = post.df)
  

  报错：随机效应数量(=48)多于观测值(45)，原因是直接在随机效应中使用样条函数会导致随机参数个数超过有效观测数。
• 二次尝试代码：

  lmer(TAC ~ 1 + bs(GA) + (1 | ID), data = post.df)
  

  可运行但仅包含随机截距，未实现SAS模型中每个ID的随机样条系数，与原模型不一致。

4. SAS模型输出参考 #

随机效应输出 #

effect      spl   estimate  Std Err DF t-value Pr>|t|
Intercept ID 8366 -6.7274 11.0589 1125 -0.61 0.5431
2 spl 1 ID 8366 -5.5371 12.3271 1125 -0.45 0.6534
3 spl 2 ID 8366 0.06818 0.1409 1125 0.48 0.6287
4 spl 3 ID 8366 -0.00007 0.001157 1125 -0.06 0.9544
5 Intercept ID 8367 -8.3470 11.4438 1125 -0.73 0.4659
6 spl 1 ID 8367 -23.3954 12.9557 1125 -1.81 0.0712
7 spl 2 ID 8367 0.3962 0.1527 1125 2.59 0.0096
8 spl 3 ID 8367 -0.00374 0.001314 1125 -2.85 0.0045
9 Intercept ID 8368 8.3946 11.7727 1125 0.71 0.4760
10 spl 1 ID 8368 10.9419 12.6986 1125 0.86 0.3891
11 spl 2 ID 8368 -0.1611 0.1445 1125 -1.11 0.2652
12 spl 3 ID 8368 0.001065 0.001188 1125 0.90 0.3705
13 Intercept ID 8369 -41.4729 10.6336 1125 -3.90 0.0001
14 spl 1 ID 8369 -42.0124 11.9225 1125 -3.52 0.0004
15 spl 2 ID 8369 0.5689 0.1344 1125 4.23 <.0001
16 spl 3 ID 8369 -0.00368 0.001243 1125 -2.96 0.0032
17 Intercept ID 8370 -0.6361 10.8564 1125 -0.06 0.9533
18 spl 1 ID 8370 -0.3350 12.4578 1125 -0.03 0.9786
19 spl 2 ID 8370 0.002890 0.1455 1125 0.02 0.9842
20 spl 3 ID 8370 -0.00024 0.001273 1125 -0.19 0.8501
21 Intercept ID 8371 -24.6320 12.4400 1125 -1.98 0.0479
22 spl 1 ID 8371 -2.3020 14.2970 1125 -0.16 0.8721
23 spl 2 ID 8371 -0.08930 0.1777 1125 -0.50 0.6153
24 spl 3 ID 8371 0.001215 0.001947 1125 0.62 0.5325
25 Intercept ID 8372 -0.2713 11.7812 1125 -0.02 0.9816
26 spl 1 ID 8372 4.9727 13.1356 1125 0.38 0.7051
27 spl 2 ID 8372 -0.09524 0.1543 1125 -0.62 0.5372
28 spl 3 ID 8372 0.000628 0.001358 1125 0.46 0.6441
29 Intercept ID 8374 -3.2900 11.0061 1125 -0.30 0.7651
30 spl 1 ID 8374 -4.7893 12.1409 1125 -0.39 0.6933
31 spl 2 ID 8374 0.07259 0.1368 1125 0.53 0.5959
32 spl 3 ID 8374 -0.00069 0.001200 1125 -0.58 0.5641
33 Intercept ID 8393 10.7260 12.0891 1125 0.89 0.3751
34 spl 1 ID 8393 9.5581 13.0069 1125 0.73 0.4626
35 spl 2 ID 8393 -0.1225 0.1497 1125 -0.82 0.4133
36 spl 3 ID 8393 0.000234 0.001277 1125 0.18 0.8546
37 Intercept ID 8394 2.2753 12.0000 1125 0.19 0.8497
38 spl 1  8394 -3.6781 13.0104 1125 -0.28 0.7775
39 spl 2 ID 8394 0.08162 0.1508 1125 0.54 0.5883
40 spl 3 ID 8394 -0.00053 0.001300 1125 -0.41 0.6810
41 Intercept ID 8395 -0.3884 12.2819 1125 -0.03 0.9748
42 spl 1 ID 8395 0.5401 14.0418 1125 0.04 0.9693
43 spl 2 ID 8395 -0.01228 0.1717 1125 -0.07 0.9430
44 spl 3 ID 8395 0.000087 0.001455 1125 0.06 0.9522
45 Intercept ID 8396 -16.6502 17.6935 1125 -0.94 0.3469
46 spl 1 ID 8396 10.1673 18.3253 1125 0.55 0.5791
47 spl 2 ID 8396 -0.2815 0.2314 1125 -1.22 0.2240
48 spl 3 ID 8396 0.001955 0.001850 1125 1.06 0.2910

固定效应输出 #

effect      spl   estimate  Std Err    DF     t-value  Pr>|t|
Intercept           93.99      2.33   1125      40.18   <.0001
spl          1       0          .       .        .        .
spl          2      -0.052     0.0197  1125     -2.64   0.0083
spl          3       0.000304  0.00004  1125    0.62   .5362

核心思路 #

SAS模型的核心是：

1. 对GA构建自然立方样条，采用TPF基，基于3分位数设置内部节点；
2. 固定效应包含截距+样条项；
3. 随机效应为每个ID的截距+样条系数，协方差矩阵为无结构（type=un）。

R中需先预计算样条基，再将其作为变量纳入模型，避免直接在随机效应中使用样条函数导致参数过多。

步骤1：匹配SAS的样条设置并预计算基函数 #

使用rms::rcs()（限制性立方样条，对应SAS的TPF基自然立方样条），或splines::ns()（B样条，拟合结果等价），先计算样条基：

library(lme4)
library(rms) # 用于匹配SAS的TPF基样条

# 清理数据集，移除NA
post.df_clean <- na.omit(post.df)

# 计算GA的3分位数节点（匹配SAS的knotmethod=percentiles(3)）
knots <- quantile(post.df_clean$GA, probs = c(1/3, 2/3))

# 生成限制性立方样条基（对应SAS的naturalcubic basis=tpf）
spl_matrix <- rcs(post.df_clean$GA, parms = knots)

# 将样条基转为数据集变量
post.df_clean$spl1 <- spl_matrix[,1]
post.df_clean$spl2 <- spl_matrix[,2]
post.df_clean$spl3 <- spl_matrix[,3]

步骤2：构建等价的混合 #