diff --git a/README.md b/README.md index 2e6cfb6..8713033 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # ae-summary -| 🧩 程序名称 | ✨ 描述 | 📚 文档 | -| ----------- | ---------------------------------------------- | ------------------ | -| `%ass2` | 按系统器官分类、首选术语汇总不良事件 | [↗️](docs/ass2.md) | -| `%ass3` | 按系统器官分类、首选术语、严重程度汇总不良事件 | [↗️](docs/ass3.md) | +| 🧩 程序名称 | ✨ 描述 | 📚 文档 | +| ----------- | ---------------------------------------------- | -------------------- | +| `%ass2` | 按系统器官分类、首选术语汇总不良事件 | [↗️](docs/ass2.md) | +| `%ass3` | 按系统器官分类、首选术语、严重程度汇总不良事件 | [↗️](docs/ass3.md) | +| `%reg_pb` | Psssing-Bablok 回归 | [↗️](docs/reg_pb.md) | diff --git a/docs/reg_pb.md b/docs/reg_pb.md new file mode 100644 index 0000000..e6b6ccf --- /dev/null +++ b/docs/reg_pb.md @@ -0,0 +1,234 @@ +# %reg_pb + +## 简介 + +Passing-Bablok 回归，计算斜率和截距，及其 95% 置信区间。 + +## 语法 + +### 参数 + +#### 必选参数 + +- [indata](#indata) +- [outdata](#outdata) +- [x](#x) +- [y](#y) + +#### 可选参数 + +- [alpha](#alpha) +- [format](#format) + +#### 调试参数 + +- [debug](#debug) + +### 参数说明 + +#### indata + +**Syntax** : _data-set-name_<(_data-set-option_)> + +指定待分析的数据集，可使用数据集选项。 + +**Usage** : + +```sas +indata = adeff +``` + +--- + +#### outdata + +**Syntax** : _data-set-name_<(_data-set-option_)> + +指定保存汇总结果的数据集，可使用数据集选项。 + +汇总结果的数据集包含以下变量： + +| 变量名 | 含义 | +| -------------- | ----------- | +| _param_ | 参数 | +| _desc_ | 参数描述 | +| _estimate_ | 估计值 | +| _lower_ | 下限 | +| _upper_ | 上限 | +| _estimate_fmt_ | 估计值（C） | +| _lower_fmt_ | 下限（C） | +| _upper_fmt_ | 上限（C） | + +**Usage** : + +```sas +outdata = t_5_4_1_2 +``` + +--- + +#### x + +**Syntax** : _variable_ + +指定 x 轴变量，在体外诊断试剂临床试验中，通常是对比试剂检测结果。 + +**Default** : `crcd1` + +**Usage** : + +```sas +x = crcd1 +``` + +--- + +#### y + +**Syntax** : _variable_ + +指定 y 轴变量，在体外诊断试剂临床试验中，通常是考核试剂检测结果。 + +**Default** : `trcd1` + +**Usage** : + +```sas +y = trcd1 +``` + +--- + +#### alpha + +**Syntax**: _numeric_ + +指定显著性水平。 + +**Default** : `0.05` + +**Usage** : + +```sas +alpha = 0.05 +``` + +--- + +#### format + +**Syntax** : _format_ + +指定估计值和置信区间上下限的输出格式。 + +**Default** : `8.4` + +**Usage** : + +```sas +format_freq = percentn9.2 +``` + +--- + +#### debug + +**Syntax** : `true` | `false` + +指定是否启用调试模式。 + +**Default** : `false` + +> [!NOTE] +> +> 这是一个用于开发者调试的参数，通常不需要关注。 + +## Passing-Bablok 回归算法 + +### 1. 构建点对子，计算斜率 + +设样本数量为 $n$ ，任取两个样本 $(x_i, y_i)$ 和 $(x_j, y_j)$ ，构成一个点对子，每个点对子的连线都有一个斜率 $S_{ij}$ + +$$ +S_{ij} = \frac{y_i - y_j}{x_i - x_j} \ \ , 1 \le i < j \le n +$$ + +这样的点对子一共有 $\mathrm{C}_n^2$ 个。 + +### 2. 排除不做贡献的点对子 + +以下两种点对子需要排除，它们对 Passing-Bablok 回归的参数估计不作贡献。 + +- 重合的点对子，即 $x_i = x_j$ 且 $y_i = y_j$ ； +- 斜率为 -1 的点对子，即 $S_{ij} = -1$ ； + +记剩余样本数量为 $M$ 。 + +### 3. 构建斜率序列 + +对 $M$ 个点对子按照斜率从小到大排序，得到序列 $S_{(1)} \le S_{(2)} \le \cdots \le S_{(M)}$ 。 + +对于斜率不存在的点对子，处理方式如下： + +- 若 $x_i = x_j$ 且 $y_i \lt y_j$ ，则 $S_{ij} = -\infty$ ，将该点对子排在序列最前； +- 若 $x_i = x_j$ 且 $y_i \gt y_j$ ，则 $S_{ij} = +\infty$ ，将该点对子排在序列最后。 + +### 4. 估计斜率 + +计 $S_{ij} \lt -1$ 的点对子数量为 $K$ ，称作“偏移量”，Passing-Bablok 回归的斜率 $b$ 可以通过序列 $S_{(i)}$ 的偏移中位数来估计： + +$$ +b = +\begin{cases} +S_{\left(\frac{M + 1}{2} + K\right)} &,\text{ if } M \text{ is odd} \\ +\frac{1}{2} \cdot \left( S_{\left(\frac{M}{2} + K\right)} + S_{\left(\frac{M}{2} + 1 + K\right)}\right) &,\text{ if } M \text{ is even} \\ +\end{cases} +$$ + +### 5. 估计斜率的置信区间 + +$$ +C = Z_{1 - \frac{\alpha}{2}} \sqrt{\frac{n(n - 1)(2n + 5)}{18}} +$$ + +其中 $Z_{1 - \frac{\alpha}{2}}$ 为标准正态分布的 $1 - \frac{\alpha}{2}$ 分位数。 + +$$ +M_1 = \frac{M - C}{2}, \ \ M_2 = M - M_1 + 1 +$$ + +其中， $M_1$ 需要舍入到最接近的整数值。 + +$b$ 的置信区间为： + +$$ +b_L = S_{(M_1 + K)} +$$ + +$$ +b_U = S_{(M_2 + K)} +$$ + +### 6. 估计截距及其置信区间 + +截距 $a$ 的估计原理是：至少有一半的点在线 $y = a + bx$ 的上方，另一半则在线 $y = a + bx$ 的下方。 + +假设每个点都能落在线 $y = a + bx$ 上，进而对每个点，都能计算一个截距 $a_i = y_i - bx_i$ ，所有 $a_i$ 从小到大排序，取中位数，即可得到截距 $a$ 的估计值： + +$$ +a = \text{median} \lbrace y_i - bx_i \rbrace +$$ + +截距 $a$ 的置信区间为： + +$$ +a_L = \text{median} \lbrace y_i - b_Ux_i \rbrace +$$ + +$$ +a_U = \text{median} \lbrace y_i - b_Lx_i \rbrace +$$ + +## References + +1. Passing H, Bablok W. A new biometrical procedure for testing the equality of measurements from two different analytical methods. Application of linear regression procedures for method comparison studies in clinical chemistry, Part I[J]. Clinical Chemistry and Laboratory Medicine, 1983, 21(11): 709-720. diff --git a/src/gb18030/reg_pb.sas b/src/gb18030/reg_pb.sas new file mode 100644 index 0000000..fac4857 --- /dev/null +++ b/src/gb18030/reg_pb.sas @@ -0,0 +1,167 @@ +/* + * Macro Name: reg_pb + * Macro Purpose: Passing-Bablok ع + * Author: wtwang + * Version Date: 2026-03-24 +*/ + +%macro reg_pb(indata, outdata, x, y, alpha = 0.05, format = 8.4, debug = false) / parmbuff; + /* indata: ݼ + * outdata: عݼ + * x: x + * y: y + * alpha: ˫ˮƽĬΪ 0.05 + * format: ͳʽĬΪ 8.4 + * debug: ģʽ + */ + + /*ͳһСд*/ + %let indata = %sysfunc(strip(%superq(indata))); + %let outdata = %sysfunc(strip(%superq(outdata))); + %let x = %upcase(%sysfunc(strip(%bquote(&x)))); + %let y = %upcase(%sysfunc(strip(%bquote(&y)))); + %let alpha = %sysfunc(strip(%bquote(&alpha))); + %let format = %upcase(%sysfunc(strip(%bquote(&format)))); + %let debug = %upcase(%sysfunc(strip(%bquote(&debug)))); + + + /* indata*/ + data tmp_indata; + set %superq(indata) end = end; + _seq = _n_; + if end then call symputx("n", _n_); + run; + + /*͵*/ + proc sql noprint; + create table tmp_indata_paired as + select + p1._seq as p1_seq label = "һ", + p1.&x as p1_x label = "һ x ֵ", + p1.&y as p1_y label = "һ y ֵ", + p2._seq as p2_seq label = "ڶ", + p2.&x as p2_x label = "ڶ x ֵ", + p2.&y as p2_y label = "ڶ y ֵ", + ifc(p1_x = p2_x and p1_y = p2_y, "Y", "") as flag_overlap label = "ʶ-غϵĵ", + ifc(p1_x = p2_x and p1_y < p2_y, "Y", "") as flag_neg_inf label = "ʶ-", + ifc(p1_x = p2_x and p1_y > p2_y, "Y", "") as flag_pos_inf label = "ʶ-", + ifn(calculated flag_overlap ^= "Y" and + calculated flag_pos_inf ^= "Y" and + calculated flag_neg_inf ^= "Y", (p2_y - p1_y)/(p2_x - p1_x), .) + as slope label = "б", + ifc(not missing(calculated slope) and calculated slope = -1, "Y", "") + as flag_slope_eq_minus_one label = "ʶ-бʵڸһ", + (case when calculated flag_neg_inf = "Y" then 1 + when calculated flag_pos_inf = "Y" then 3 + else 2 + end) as flag_sort label = "ʶ-˳" + from tmp_indata as p1, tmp_indata as p2 + where (p1_seq < p2_seq) and calculated flag_overlap ^= "Y" and calculated flag_slope_eq_minus_one ^= "Y" + order by flag_sort, slope; + + /**/ + select count(*) into :M trimmed from tmp_indata_paired; + + /*¼еλ*/ + alter table tmp_indata_paired add seq num(8) label = ""; + update tmp_indata_paired set seq = monotonic(); + quit; + + /*бС -1 ĵ*/ + proc sql noprint; + select sum((not missing(slope) and slope < -1) or flag_neg_inf = "Y") into :K trimmed from tmp_indata_paired; + quit; + + /*洢ֵݼ*/ + proc sql noprint; + create table tmp_parameter_estimate + (param char(20) label = "", + desc char(20) label = "", + estimate num(8) label = "ֵ", + lower num(8) label = "", + upper num(8) label = "", + estimate_fmt char(20) label = "ֵC", + lower_fmt char(20) label = "ޣC", + upper_fmt char(20) label = "ޣC"); + quit; + + /*-б*/ + proc sql noprint; + insert into tmp_parameter_estimate + set param = "slope", desc = "б", estimate = (case when mod(&M, 2) = 1 then (select slope from tmp_indata_paired where seq = (&M + 1)/2 + &K) + when mod(&M, 2) = 0 then ((select slope from tmp_indata_paired where seq = (&M/2 + &K)) + + (select slope from tmp_indata_paired where seq = (&M/2 + 1 + &K))) / 2 + end); + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_est as + select + &x, + &y, + &y - (select estimate from tmp_parameter_estimate where param = "slope") * &x as intercept_individual label = "ڸϼĽؾ" + from tmp_indata; + insert into tmp_parameter_estimate + set param = "intercept", desc = "ؾ", estimate = (select median(intercept_individual) from tmp_indata_intercept_est); + quit; + + /*-б*/ + data _null_; + C = probit(1 - &alpha / 2) * sqrt(&n * (&n - 1) * (2 * &n + 5) / 18); + M1 = round((&M - C) / 2); + M2 = &M - M1 + 1; + call symputx("M1", M1); + call symputx("M2", M2); + run; + proc sql noprint; + update tmp_parameter_estimate + set lower = (select slope from tmp_indata_paired where seq = &M1 + &K), + upper = (select slope from tmp_indata_paired where seq = &M2 + &K) + where param = "slope"; + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_ci_est as + select + &x, + &y, + &y - (select upper from tmp_parameter_estimate where param = "slope") * &x as intercept_lower_individual label = "ڸϼĽؾ", + &y - (select lower from tmp_parameter_estimate where param = "slope") * &x as intercept_upper_individual label = "ڸϼĽؾ" + from tmp_indata; + update tmp_parameter_estimate + set lower = (select median(intercept_lower_individual) from tmp_indata_intercept_ci_est), + upper = (select median(intercept_upper_individual) from tmp_indata_intercept_ci_est) + where param = "intercept"; + quit; + + /*ʽ*/ + proc sql noprint; + update tmp_parameter_estimate + set estimate_fmt = put(estimate, &format -L), + lower_fmt = put(lower, &format -L), + upper_fmt = put(upper, &format -L); + quit; + + + /*ݼ*/ + data &outdata; + set tmp_parameter_estimate; + run; + + /*ɾмݼ*/ + %if %bquote(&debug) = %upcase(false) %then %do; + proc datasets library = work nowarn noprint; + delete tmp_indata + tmp_indata_paired + tmp_parameter_estimate + tmp_indata_intercept_est + tmp_indata_intercept_ci_est + ; + quit; + %end; + + %exit: + %put NOTE: reg_pb ѽУ; +%mend; diff --git a/src/gb2312/reg_pb.sas b/src/gb2312/reg_pb.sas new file mode 100644 index 0000000..fac4857 --- /dev/null +++ b/src/gb2312/reg_pb.sas @@ -0,0 +1,167 @@ +/* + * Macro Name: reg_pb + * Macro Purpose: Passing-Bablok ع + * Author: wtwang + * Version Date: 2026-03-24 +*/ + +%macro reg_pb(indata, outdata, x, y, alpha = 0.05, format = 8.4, debug = false) / parmbuff; + /* indata: ݼ + * outdata: عݼ + * x: x + * y: y + * alpha: ˫ˮƽĬΪ 0.05 + * format: ͳʽĬΪ 8.4 + * debug: ģʽ + */ + + /*ͳһСд*/ + %let indata = %sysfunc(strip(%superq(indata))); + %let outdata = %sysfunc(strip(%superq(outdata))); + %let x = %upcase(%sysfunc(strip(%bquote(&x)))); + %let y = %upcase(%sysfunc(strip(%bquote(&y)))); + %let alpha = %sysfunc(strip(%bquote(&alpha))); + %let format = %upcase(%sysfunc(strip(%bquote(&format)))); + %let debug = %upcase(%sysfunc(strip(%bquote(&debug)))); + + + /* indata*/ + data tmp_indata; + set %superq(indata) end = end; + _seq = _n_; + if end then call symputx("n", _n_); + run; + + /*͵*/ + proc sql noprint; + create table tmp_indata_paired as + select + p1._seq as p1_seq label = "һ", + p1.&x as p1_x label = "һ x ֵ", + p1.&y as p1_y label = "һ y ֵ", + p2._seq as p2_seq label = "ڶ", + p2.&x as p2_x label = "ڶ x ֵ", + p2.&y as p2_y label = "ڶ y ֵ", + ifc(p1_x = p2_x and p1_y = p2_y, "Y", "") as flag_overlap label = "ʶ-غϵĵ", + ifc(p1_x = p2_x and p1_y < p2_y, "Y", "") as flag_neg_inf label = "ʶ-", + ifc(p1_x = p2_x and p1_y > p2_y, "Y", "") as flag_pos_inf label = "ʶ-", + ifn(calculated flag_overlap ^= "Y" and + calculated flag_pos_inf ^= "Y" and + calculated flag_neg_inf ^= "Y", (p2_y - p1_y)/(p2_x - p1_x), .) + as slope label = "б", + ifc(not missing(calculated slope) and calculated slope = -1, "Y", "") + as flag_slope_eq_minus_one label = "ʶ-бʵڸһ", + (case when calculated flag_neg_inf = "Y" then 1 + when calculated flag_pos_inf = "Y" then 3 + else 2 + end) as flag_sort label = "ʶ-˳" + from tmp_indata as p1, tmp_indata as p2 + where (p1_seq < p2_seq) and calculated flag_overlap ^= "Y" and calculated flag_slope_eq_minus_one ^= "Y" + order by flag_sort, slope; + + /**/ + select count(*) into :M trimmed from tmp_indata_paired; + + /*¼еλ*/ + alter table tmp_indata_paired add seq num(8) label = ""; + update tmp_indata_paired set seq = monotonic(); + quit; + + /*бС -1 ĵ*/ + proc sql noprint; + select sum((not missing(slope) and slope < -1) or flag_neg_inf = "Y") into :K trimmed from tmp_indata_paired; + quit; + + /*洢ֵݼ*/ + proc sql noprint; + create table tmp_parameter_estimate + (param char(20) label = "", + desc char(20) label = "", + estimate num(8) label = "ֵ", + lower num(8) label = "", + upper num(8) label = "", + estimate_fmt char(20) label = "ֵC", + lower_fmt char(20) label = "ޣC", + upper_fmt char(20) label = "ޣC"); + quit; + + /*-б*/ + proc sql noprint; + insert into tmp_parameter_estimate + set param = "slope", desc = "б", estimate = (case when mod(&M, 2) = 1 then (select slope from tmp_indata_paired where seq = (&M + 1)/2 + &K) + when mod(&M, 2) = 0 then ((select slope from tmp_indata_paired where seq = (&M/2 + &K)) + + (select slope from tmp_indata_paired where seq = (&M/2 + 1 + &K))) / 2 + end); + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_est as + select + &x, + &y, + &y - (select estimate from tmp_parameter_estimate where param = "slope") * &x as intercept_individual label = "ڸϼĽؾ" + from tmp_indata; + insert into tmp_parameter_estimate + set param = "intercept", desc = "ؾ", estimate = (select median(intercept_individual) from tmp_indata_intercept_est); + quit; + + /*-б*/ + data _null_; + C = probit(1 - &alpha / 2) * sqrt(&n * (&n - 1) * (2 * &n + 5) / 18); + M1 = round((&M - C) / 2); + M2 = &M - M1 + 1; + call symputx("M1", M1); + call symputx("M2", M2); + run; + proc sql noprint; + update tmp_parameter_estimate + set lower = (select slope from tmp_indata_paired where seq = &M1 + &K), + upper = (select slope from tmp_indata_paired where seq = &M2 + &K) + where param = "slope"; + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_ci_est as + select + &x, + &y, + &y - (select upper from tmp_parameter_estimate where param = "slope") * &x as intercept_lower_individual label = "ڸϼĽؾ", + &y - (select lower from tmp_parameter_estimate where param = "slope") * &x as intercept_upper_individual label = "ڸϼĽؾ" + from tmp_indata; + update tmp_parameter_estimate + set lower = (select median(intercept_lower_individual) from tmp_indata_intercept_ci_est), + upper = (select median(intercept_upper_individual) from tmp_indata_intercept_ci_est) + where param = "intercept"; + quit; + + /*ʽ*/ + proc sql noprint; + update tmp_parameter_estimate + set estimate_fmt = put(estimate, &format -L), + lower_fmt = put(lower, &format -L), + upper_fmt = put(upper, &format -L); + quit; + + + /*ݼ*/ + data &outdata; + set tmp_parameter_estimate; + run; + + /*ɾмݼ*/ + %if %bquote(&debug) = %upcase(false) %then %do; + proc datasets library = work nowarn noprint; + delete tmp_indata + tmp_indata_paired + tmp_parameter_estimate + tmp_indata_intercept_est + tmp_indata_intercept_ci_est + ; + quit; + %end; + + %exit: + %put NOTE: reg_pb ѽУ; +%mend; diff --git a/src/gbk/reg_pb.sas b/src/gbk/reg_pb.sas new file mode 100644 index 0000000..fac4857 --- /dev/null +++ b/src/gbk/reg_pb.sas @@ -0,0 +1,167 @@ +/* + * Macro Name: reg_pb + * Macro Purpose: Passing-Bablok ع + * Author: wtwang + * Version Date: 2026-03-24 +*/ + +%macro reg_pb(indata, outdata, x, y, alpha = 0.05, format = 8.4, debug = false) / parmbuff; + /* indata: ݼ + * outdata: عݼ + * x: x + * y: y + * alpha: ˫ˮƽĬΪ 0.05 + * format: ͳʽĬΪ 8.4 + * debug: ģʽ + */ + + /*ͳһСд*/ + %let indata = %sysfunc(strip(%superq(indata))); + %let outdata = %sysfunc(strip(%superq(outdata))); + %let x = %upcase(%sysfunc(strip(%bquote(&x)))); + %let y = %upcase(%sysfunc(strip(%bquote(&y)))); + %let alpha = %sysfunc(strip(%bquote(&alpha))); + %let format = %upcase(%sysfunc(strip(%bquote(&format)))); + %let debug = %upcase(%sysfunc(strip(%bquote(&debug)))); + + + /* indata*/ + data tmp_indata; + set %superq(indata) end = end; + _seq = _n_; + if end then call symputx("n", _n_); + run; + + /*͵*/ + proc sql noprint; + create table tmp_indata_paired as + select + p1._seq as p1_seq label = "һ", + p1.&x as p1_x label = "һ x ֵ", + p1.&y as p1_y label = "һ y ֵ", + p2._seq as p2_seq label = "ڶ", + p2.&x as p2_x label = "ڶ x ֵ", + p2.&y as p2_y label = "ڶ y ֵ", + ifc(p1_x = p2_x and p1_y = p2_y, "Y", "") as flag_overlap label = "ʶ-غϵĵ", + ifc(p1_x = p2_x and p1_y < p2_y, "Y", "") as flag_neg_inf label = "ʶ-", + ifc(p1_x = p2_x and p1_y > p2_y, "Y", "") as flag_pos_inf label = "ʶ-", + ifn(calculated flag_overlap ^= "Y" and + calculated flag_pos_inf ^= "Y" and + calculated flag_neg_inf ^= "Y", (p2_y - p1_y)/(p2_x - p1_x), .) + as slope label = "б", + ifc(not missing(calculated slope) and calculated slope = -1, "Y", "") + as flag_slope_eq_minus_one label = "ʶ-бʵڸһ", + (case when calculated flag_neg_inf = "Y" then 1 + when calculated flag_pos_inf = "Y" then 3 + else 2 + end) as flag_sort label = "ʶ-˳" + from tmp_indata as p1, tmp_indata as p2 + where (p1_seq < p2_seq) and calculated flag_overlap ^= "Y" and calculated flag_slope_eq_minus_one ^= "Y" + order by flag_sort, slope; + + /**/ + select count(*) into :M trimmed from tmp_indata_paired; + + /*¼еλ*/ + alter table tmp_indata_paired add seq num(8) label = ""; + update tmp_indata_paired set seq = monotonic(); + quit; + + /*бС -1 ĵ*/ + proc sql noprint; + select sum((not missing(slope) and slope < -1) or flag_neg_inf = "Y") into :K trimmed from tmp_indata_paired; + quit; + + /*洢ֵݼ*/ + proc sql noprint; + create table tmp_parameter_estimate + (param char(20) label = "", + desc char(20) label = "", + estimate num(8) label = "ֵ", + lower num(8) label = "", + upper num(8) label = "", + estimate_fmt char(20) label = "ֵC", + lower_fmt char(20) label = "ޣC", + upper_fmt char(20) label = "ޣC"); + quit; + + /*-б*/ + proc sql noprint; + insert into tmp_parameter_estimate + set param = "slope", desc = "б", estimate = (case when mod(&M, 2) = 1 then (select slope from tmp_indata_paired where seq = (&M + 1)/2 + &K) + when mod(&M, 2) = 0 then ((select slope from tmp_indata_paired where seq = (&M/2 + &K)) + + (select slope from tmp_indata_paired where seq = (&M/2 + 1 + &K))) / 2 + end); + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_est as + select + &x, + &y, + &y - (select estimate from tmp_parameter_estimate where param = "slope") * &x as intercept_individual label = "ڸϼĽؾ" + from tmp_indata; + insert into tmp_parameter_estimate + set param = "intercept", desc = "ؾ", estimate = (select median(intercept_individual) from tmp_indata_intercept_est); + quit; + + /*-б*/ + data _null_; + C = probit(1 - &alpha / 2) * sqrt(&n * (&n - 1) * (2 * &n + 5) / 18); + M1 = round((&M - C) / 2); + M2 = &M - M1 + 1; + call symputx("M1", M1); + call symputx("M2", M2); + run; + proc sql noprint; + update tmp_parameter_estimate + set lower = (select slope from tmp_indata_paired where seq = &M1 + &K), + upper = (select slope from tmp_indata_paired where seq = &M2 + &K) + where param = "slope"; + quit; + + /*-ؾ*/ + proc sql noprint; + create table tmp_indata_intercept_ci_est as + select + &x, + &y, + &y - (select upper from tmp_parameter_estimate where param = "slope") * &x as intercept_lower_individual label = "ڸϼĽؾ", + &y - (select lower from tmp_parameter_estimate where param = "slope") * &x as intercept_upper_individual label = "ڸϼĽؾ" + from tmp_indata; + update tmp_parameter_estimate + set lower = (select median(intercept_lower_individual) from tmp_indata_intercept_ci_est), + upper = (select median(intercept_upper_individual) from tmp_indata_intercept_ci_est) + where param = "intercept"; + quit; + + /*ʽ*/ + proc sql noprint; + update tmp_parameter_estimate + set estimate_fmt = put(estimate, &format -L), + lower_fmt = put(lower, &format -L), + upper_fmt = put(upper, &format -L); + quit; + + + /*ݼ*/ + data &outdata; + set tmp_parameter_estimate; + run; + + /*ɾмݼ*/ + %if %bquote(&debug) = %upcase(false) %then %do; + proc datasets library = work nowarn noprint; + delete tmp_indata + tmp_indata_paired + tmp_parameter_estimate + tmp_indata_intercept_est + tmp_indata_intercept_ci_est + ; + quit; + %end; + + %exit: + %put NOTE: reg_pb ѽУ; +%mend; diff --git a/src/utf8/reg_pb.sas b/src/utf8/reg_pb.sas new file mode 100644 index 0000000..8119a65 --- /dev/null +++ b/src/utf8/reg_pb.sas @@ -0,0 +1,167 @@ +/* + * Macro Name: reg_pb + * Macro Purpose: Passing-Bablok 回归 + * Author: wtwang + * Version Date: 2026-03-24 +*/ + +%macro reg_pb(indata, outdata, x, y, alpha = 0.05, format = 8.4, debug = false) / parmbuff; + /* indata: 数据集名称 + * outdata: 保存回归分析结果的数据集名称 + * x: x 轴变量 + * y: y 轴变量 + * alpha: 双侧显著性水平，默认为 0.05 + * format: 统计量输出格式，默认为 8.4 + * debug: 调试模式 + */ + + /*统一参数大小写*/ + %let indata = %sysfunc(strip(%superq(indata))); + %let outdata = %sysfunc(strip(%superq(outdata))); + %let x = %upcase(%sysfunc(strip(%bquote(&x)))); + %let y = %upcase(%sysfunc(strip(%bquote(&y)))); + %let alpha = %sysfunc(strip(%bquote(&alpha))); + %let format = %upcase(%sysfunc(strip(%bquote(&format)))); + %let debug = %upcase(%sysfunc(strip(%bquote(&debug)))); + + + /*复制 indata*/ + data tmp_indata; + set %superq(indata) end = end; + _seq = _n_; + if end then call symputx("n", _n_); + run; + + /*点和点进行配对*/ + proc sql noprint; + create table tmp_indata_paired as + select + p1._seq as p1_seq label = "第一个点的序号", + p1.&x as p1_x label = "第一个点的 x 值", + p1.&y as p1_y label = "第一个点的 y 值", + p2._seq as p2_seq label = "第二个点的序号", + p2.&x as p2_x label = "第二个点的 x 值", + p2.&y as p2_y label = "第二个点的 y 值", + ifc(p1_x = p2_x and p1_y = p2_y, "Y", "") as flag_overlap label = "标识变量-重合的点", + ifc(p1_x = p2_x and p1_y < p2_y, "Y", "") as flag_neg_inf label = "标识变量-负无穷大", + ifc(p1_x = p2_x and p1_y > p2_y, "Y", "") as flag_pos_inf label = "标识变量-正无穷大", + ifn(calculated flag_overlap ^= "Y" and + calculated flag_pos_inf ^= "Y" and + calculated flag_neg_inf ^= "Y", (p2_y - p1_y)/(p2_x - p1_x), .) + as slope label = "斜率", + ifc(not missing(calculated slope) and calculated slope = -1, "Y", "") + as flag_slope_eq_minus_one label = "标识变量-斜率等于负一", + (case when calculated flag_neg_inf = "Y" then 1 + when calculated flag_pos_inf = "Y" then 3 + else 2 + end) as flag_sort label = "标识变量-排列顺序" + from tmp_indata as p1, tmp_indata as p2 + where (p1_seq < p2_seq) and calculated flag_overlap ^= "Y" and calculated flag_slope_eq_minus_one ^= "Y" + order by flag_sort, slope; + + /*计算对子数*/ + select count(*) into :M trimmed from tmp_indata_paired; + + /*记录对子在序列中的位置*/ + alter table tmp_indata_paired add seq num(8) label = "点对子序号"; + update tmp_indata_paired set seq = monotonic(); + quit; + + /*计算斜率小于 -1 的点对数*/ + proc sql noprint; + select sum((not missing(slope) and slope < -1) or flag_neg_inf = "Y") into :K trimmed from tmp_indata_paired; + quit; + + /*创建存储参数估计值的数据集*/ + proc sql noprint; + create table tmp_parameter_estimate + (param char(20) label = "参数", + desc char(20) label = "参数描述", + estimate num(8) label = "估计值", + lower num(8) label = "下限", + upper num(8) label = "上限", + estimate_fmt char(20) label = "估计值（C）", + lower_fmt char(20) label = "下限（C）", + upper_fmt char(20) label = "上限（C）"); + quit; + + /*点估计-斜率*/ + proc sql noprint; + insert into tmp_parameter_estimate + set param = "slope", desc = "斜率", estimate = (case when mod(&M, 2) = 1 then (select slope from tmp_indata_paired where seq = (&M + 1)/2 + &K) + when mod(&M, 2) = 0 then ((select slope from tmp_indata_paired where seq = (&M/2 + &K)) + + (select slope from tmp_indata_paired where seq = (&M/2 + 1 + &K))) / 2 + end); + quit; + + /*点估计-截距*/ + proc sql noprint; + create table tmp_indata_intercept_est as + select + &x, + &y, + &y - (select estimate from tmp_parameter_estimate where param = "slope") * &x as intercept_individual label = "在个体样本上计算的截距" + from tmp_indata; + insert into tmp_parameter_estimate + set param = "intercept", desc = "截距", estimate = (select median(intercept_individual) from tmp_indata_intercept_est); + quit; + + /*区间估计-斜率*/ + data _null_; + C = probit(1 - &alpha / 2) * sqrt(&n * (&n - 1) * (2 * &n + 5) / 18); + M1 = round((&M - C) / 2); + M2 = &M - M1 + 1; + call symputx("M1", M1); + call symputx("M2", M2); + run; + proc sql noprint; + update tmp_parameter_estimate + set lower = (select slope from tmp_indata_paired where seq = &M1 + &K), + upper = (select slope from tmp_indata_paired where seq = &M2 + &K) + where param = "slope"; + quit; + + /*区间估计-截距*/ + proc sql noprint; + create table tmp_indata_intercept_ci_est as + select + &x, + &y, + &y - (select upper from tmp_parameter_estimate where param = "slope") * &x as intercept_lower_individual label = "在个体样本上计算的截距下限", + &y - (select lower from tmp_parameter_estimate where param = "slope") * &x as intercept_upper_individual label = "在个体样本上计算的截距上限" + from tmp_indata; + update tmp_parameter_estimate + set lower = (select median(intercept_lower_individual) from tmp_indata_intercept_ci_est), + upper = (select median(intercept_upper_individual) from tmp_indata_intercept_ci_est) + where param = "intercept"; + quit; + + /*格式化*/ + proc sql noprint; + update tmp_parameter_estimate + set estimate_fmt = put(estimate, &format -L), + lower_fmt = put(lower, &format -L), + upper_fmt = put(upper, &format -L); + quit; + + + /*输出数据集*/ + data &outdata; + set tmp_parameter_estimate; + run; + + /*删除中间数据集*/ + %if %bquote(&debug) = %upcase(false) %then %do; + proc datasets library = work nowarn noprint; + delete tmp_indata + tmp_indata_paired + tmp_parameter_estimate + tmp_indata_intercept_est + tmp_indata_intercept_ci_est + ; + quit; + %end; + + %exit: + %put NOTE: 宏程序 reg_pb 已结束运行！; +%mend;