/******************************************************************
* 两个独立 Poisson 发生率比的条件精确区间与检验
* 数据: example.csv（原创模拟汇总数据）
* 方向: Drug A / Placebo
* 方法: 条件于总事件数的二项分布 + Clopper-Pearson 中心区间
* 环境: Base SAS 9.4
******************************************************************/

proc import datafile="example.csv" out=rate_input dbms=csv replace;
  guessingrows=max;
run;

proc sql noprint;
  select events, exposure_py into :x_a, :t_a
  from rate_input where upcase(trt01p)='DRUG A';
  select events, exposure_py into :x_p, :t_p
  from rate_input where upcase(trt01p)='PLACEBO';
quit;

data exact_rate_ratio;
  alpha=0.05;
  x_a=&x_a; t_a=&t_a;
  x_p=&x_p; t_p=&t_p;
  n=x_a+x_p;
  rate_a=x_a/t_a;
  rate_p=x_p/t_p;
  if rate_p>0 then rate_ratio=rate_a/rate_p;
  else if rate_a>0 then rate_ratio=constant('BIG'); /* shell 显示为 NE/∞ */
  else call missing(rate_ratio);                    /* 两组均 0 事件：点估计未定义 */

  /* 条件分布 X_A | (X_A+X_P=n) ~ Binomial(n,p) */
  if x_a=0 then p_lower=0;
  else p_lower=quantile('BETA',alpha/2,x_a,n-x_a+1);
  if x_a=n then p_upper=1;
  else p_upper=quantile('BETA',1-alpha/2,x_a+1,n-x_a);

  if p_lower=0 then rr_lower=0;
  else rr_lower=(p_lower/(1-p_lower))*(t_p/t_a);
  if p_upper=1 then rr_upper=constant('BIG');
  else rr_upper=(p_upper/(1-p_upper))*(t_p/t_a);

  /* H0: rate ratio=1；中心双侧条件精确 p 值 */
  p0=t_a/(t_a+t_p);
  lower_tail=cdf('BINOMIAL',x_a,p0,n);
  upper_tail=1-cdf('BINOMIAL',x_a-1,p0,n);
  p_exact=min(1,2*min(lower_tail,upper_tail));
  label rate_ratio='Rate ratio (Drug A / Placebo)'
        rr_lower='Exact 95% lower confidence limit'
        rr_upper='Exact 95% upper confidence limit'
        p_exact='Two-sided conditional exact p-value';
  format rate_a rate_p rate_ratio rr_lower rr_upper 10.4 p_exact pvalue8.4;
run;

proc print data=exact_rate_ratio noobs label;
  var x_a t_a rate_a x_p t_p rate_p rate_ratio rr_lower rr_upper p_exact;
run;