Stata has dtable and etable procedure, respectively, to produce the demographic and regression outcome table. See the screenshot.   I wonder if SAS has the similar syntaxes.   ... View more

Hello, pal   Thanks for your code about the spline in SAS. However, I found the result derived from the code in SAS is quite different from the result I got by Stata procedures, given the exactly same dataset. I have not yet identified the causes. In spite of this dispute, I would like to share my SAS code, which handles the parameters got by Stata spline regression and then do the prediction through matrix calculation.   The code use the covariance-variance matrix to compute the variance of the outcome variable. Additionally, the mean of the outcome variable could derive from the matrix calculation by regression coefficient matrix.   I hope this update will facilitate the handle of the data to make a restricted cubic spline. Have a great day!   Tom   /*Matrix calculation*/ proc iml; use mx; read all var {bnp} into S; close mx; use mx; /* 1. Read all three variables into a single matrix */ read all var {bnps1 bnps2 bnps3 interaction} into X; close mx; use cov; read all var {bnps1 bnps2 bnps3 interaction} into V; close cov; use coefficients; read all var {bnps1 bnps2 bnps3 interaction} into M; close coefficients; *print X; *print V; Pred = X * V * X`; PV = vecdiag(Pred); *print PV; SE = sqrt(PV); print SE; Mean = X * M`; print X; print Mean; Upper = Mean + SE*1.96; Lower = Mean - SE*1.96; HR = exp(Mean); lo = exp(Lower); hi = exp(Upper); print HR lo hi; /* Combine vectors into one matrix */ outmatrix = S || HR || lo || hi; /* Extract the columns from a matrix */ /* a. The CREATE statement defines the new dataset and its column names. */ create HR_Results from outmatrix [colname={'bnp' 'hr' 'lower' 'upper'}]; /* b. The APPEND statement writes the data from the matrix into the new dataset. */ append from outmatrix; quit; ... View more

Thanks for the additional information and I will try your method later. Have a great day! ... View more

So far, I have to get a dataset of predicted values from Stata, then make the spline diagram by proc sgplot. This combined procedure makes sure the spline is continuous.   On the other hand, I learned the proc iml a little and by the matrix algebra, I can get the variance (actually the covariance-variance matrix) of the outcome variable. After extraction the variance from the diagonal line from the matrix, standard error for each observation is known. However, this approach requires Stata two, as I don't now how to transfer a continuous variable to a few spline variables in SAS. These new variables transferred are those taken to construct the regression model.   ... View more

Thanks, pal. Hmmm, looks like we should avoid the IPTW technique for time-dependent COX model. ... View more

Thanks for your feedback, pal.   I have never heard of the time-dependent IPTW. I have read thousands of clinical study publications but not yet encountered a case with the time-based IPTW technique. It must be very cool. ... View more

Hi, Jacob Thank you for your feedback and I agree with you. The likelihood function on the textbook is a good way for parameter estimation. I looked up in the SAS document and found SAS can provide two distributions for U, including gamma distribution and the other which I forget the name. Much appreciated! ... View more

The source of the likelihood function of a cure model above also lists the PSPMCM macro as the solution for the mixture survival model.   ... View more

@ChrislHemedinger and @SAS_Rob Hey, guys Thanks a bunch for your response. I know there are some MACRO but I am not sure if the follow the maximum likelihood method for the cure model. I read from a book that the likelihood function of a cure model include a mixture of a logistic and a Cox model. I really wish the SAS can have a procedure to handle this type of survival analysis because I see its academic potential in evidence-based medicine, especially cancer, etc.   The likelihood function is below and it simply follows the partial likelihood function of a regular survival analysis model.   ... View more

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; ... View more

After some coding effort, I got the survival and hazard function graph, respectively. It makes me satisfied to solve this issue and I'm grateful we have this community platform and you guys' kind help.   ... View more