I appreciate the help and resources from this community so I am glad to share my codes for this project. Here they are:
One thing I am not sure is the variance of the independent of "albumin_used" (a binary variable)which is the intervention of interest. I used the IPTW technique to control the known probable confounders and it produced a "pseudo" cohort. I used the PROC UNIVARIATE procedure with WEIGHT statement to compute the variance of "albumin_used", and it turned out to be 0.70+, which surprised me. With the IPTW technique, the PROC FREQ procedure showed the percentage of albumin_used = 1 is nearly 0.50, therefore, I think the variance of it should be ~0.25 for the "pseudo" cohort.
The other consideration of the variance's value is that it direct affect the width of the 95%CI for both the survival and hazard function. If it's too big, the interval estimation would be too obscure.
Okay, please view the codes.
proc logistic data = albumin;
model albumin_used(event = '1') = sex
age
sofa
mv
rrt
hypertension
diabetes
copd
coronary_disease;
output out = albumin_pred p = p1;
run;
proc sql noprint;
create table albumin_iptw as
select *,
case
when albumin_used = 1 then 1/p1
when albumin_used = 0 then 1/(1 - p1)
end as wt
from albumin_pred;
quit;
PROC SQL noprint;
UPDATE albumin_iptw
SET time = los_hospital
WHERE time = .; /* Optional: to update specific rows */
QUIT;
proc print data = albumin_iptw (obs=5);
run;
proc lifetest data = albumin_iptw plots = survival(cb test atrisk) nelson
outsurv = surv1 notable;
strata albumin_used;
time time*mortality(0);
weight wt;
run;
proc lifereg data = albumin_iptw order = data;
class
albumin_used;
model time*mortality(0) = albumin_used
/ distribution = weibull;
weight wt;
run;
proc phreg data = albumin_iptw;
class albumin_used(ref = '0');
model time*mortality(0) = albumin_used / rl=wald;
weight wt;
run;
proc univariate data = albumin_iptw;
var albumin_used;
weight wt;
run;
proc freq data = albumin_iptw;
tables albumin_used;
weight wt;
run;
data hazard1;
do t = 0 to 365 by 0.1; /* Adjust the range and increment as appropriate for your data */
lambda = exp(-(4.1239 + 0.3456*1)/1.0101); /* Replace `intercept` with your actual intercept value */
rho = 0.9900; /* Shape parameter from PROC LIFEREG output */
h_t = lambda * rho * (t**(rho - 1)); /* Calculate hazard */
vht = (lambda*t**(rho-1)*rho*(-0.3456/1.0101))**2*0.25;
lower_h_t = h_t - 1.96*vht**(0.5);
upper_h_t = h_t + 1.96*vht**(0.5);
suriv = exp(-(lambda)*(t**(rho)));
vsuriv = (suriv*lambda*(t**rho)*0.3456/1.0101)**2*0.25;
suriv_l = suriv - 1.96*vsuriv**(0.5);
suriv_u = suriv + 1.96*vsuriv**(0.5);
output;
end;
*drop lambda rho rho_l rho_u;
run;
data hazard2;
do t = 0 to 365 by 0.1; /* Adjust the range and increment as appropriate for your data */
lambda = exp(-(4.1239 + 0.3456*0)/1.0101); /* Replace `intercept` with your actual intercept value */
rho = 0.9900; /* Shape parameter from PROC LIFEREG output */
h_t = lambda * rho * (t**(rho - 1)); /* Calculate hazard */
vht = (lambda*t**(rho-1)*rho*(-0.3456/1.0101))**2*0.25;
lower_h_t = h_t - 1.96*vht**(0.5);
upper_h_t = h_t + 1.96*vht**(0.5);
suriv = exp(-(lambda)*(t**(rho)));
vsuriv = (suriv*lambda*(t**rho)*0.3456/1.0101)**2*0.25;
suriv_l = suriv - 1.96*vsuriv**(0.5);
suriv_u = suriv + 1.96*vsuriv**(0.5);
output;
end;
*drop lambda rho rho_l rho_u;
run;
proc sql noprint;
create table newhazard1 as
select *, /* Use * to select all existing columns */
rho*0+1 as subgroup /* Define the new column using an expression */
from hazard1;
quit;
proc sql noprint;
create table newhazard2 as
select *, /* Use * to select all existing columns */
rho*0 as subgroup /* Define the new column using an expression */
from hazard2;
quit;
proc sql noprint;
create table hzf as
select one.*
from newhazard1 as one
union
select two.*
from newhazard2 as two;
quit;
proc sgplot data = hzf;
*styleattrs datacolors=('CX115F97' 'CX9D7729');
series x=t y=suriv / markers markerattrs=(symbol=circlefilled size=4) group=subgroup;
band x=t lower=suriv_l upper=suriv_u / transparency = 0.7 group=subgroup;
xaxis label="Time (day)" values = (0, 7, 14, 21, 30, 90 180);
yaxis label="Survival (%)" values = (0, 0.25, 0.5, 0.75, 1.0);
title "Weibull Survival Function";
refline 0.50 / axis=y lineattrs=(color = black);
inset 'IPTW Model';
run;
proc sgplot data = hzf;
*styleattrs datacolors=('CX115F97' 'CX9D7729');
series x=t y=h_t / markers markerattrs=(symbol=circlefilled size=4) group=subgroup;
band x=t lower=lower_h_t upper=upper_h_t / transparency = 0.7 group=subgroup;
xaxis label="Time (day)" values = (0, 7, 14, 21, 30, 90 180);
yaxis label="Hazard" values = (0, 0.05);
title "Weibull Hazard Function";
inset 'IPTW Model';
run;
