cox ph simulation

[Ben12345]

Hi there, I am trying to calculate the power of a cox ph regression model by running a simulation of 1000 samples and calculating the % of pvalues<=0.05. However I cant seem to recall the stderr's of the parameter estimates or associated pvalues or 95% CI. The only thing I have been successful doing is the mean parameter estimate to calculate bias. Having trouble with what I can get from the outest command. Can anyone please help with the code to generate a table that shows percentage of pvalues <=0.05. Progress of my code shown below. Thanks heaps

run simulation, then:

/*Sort dataset */
proc sort data=simcoxlog out=simcoxlog2; by sampleno; run;
/*run cox regression */
proc phreg data=simcoxlog2 outest=coxest noprint;
by sampleno;
model time*censored(1) = treatment;
run;
/*calculate means of sampling dist for cox regression */
proc means data=coxest mean stddev;
var treatment;
run;

[FreelanceReinh · Accepted solution]

Hi @Ben12345 and welcome to the SAS Support Communities!

Select the appropriate items from the list of ODS tables and create datasets containing the desired information by means of an ODS OUTPUT statement. In your example the ODS table ParameterEstimates includes the parameter estimates, standard errors, (Wald Chi-Square) p-values and more. Instead of the NOPRINT option use an ODS EXCLUDE (or ODS SELECT) statement to suppress listing or HTML output. The ODS NORESULTS statement can further improve performance by not sending output to the Results window. You may also want to suppress the notes about convergence (for each BY group) in the log.

ods exclude all;
ods noresults;
options nonotes;

ods output ParameterEstimates=est;
        
proc phreg data=simcoxlog2;
by sampleno;
model time*censored(1) = treatment;
run;

ods exclude none;
ods results;
options notes;

proc format;
value signif
low-0.05   = 'significant'
0.05<-high = 'not significant';
run;

proc freq data=est;
format ProbChiSq signif.;
tables ProbChiSq;
run;

proc means data=est mean stddev;
var estimate;
run;

Add more items of the form ODS_table = dataset to the ODS OUTPUT statement if you need more statistics (e.g. global p-values from ODS table GlobalTests). In general, check the structure of the ODS output datasets (here: est) and add WHERE (or other) statements to the evaluating PROC FREQ, PROC MEANS etc. steps as appropriate, e.g., when you add more effects to the MODEL statement or when you use a CLASS variable with more than two levels.

[Ben12345]

Man thank you so much I've been so stuck. This worked perfectly! Thankyou again

[Ben12345]

Man thank you so much I've been so stuck. This worked perfectly! Thankyou again

[Ben12345]

Thanks for being so helpful already, hoping I can stretch it with two more questions on this topic, sorry I'm new to SAS and have a long way to go.

Question 1

I'm trying to compare power of cox and logistic regression. I want to run a logistic regression on the same data however the intercept values are being included in the pvalues. I can't seem to use the 'by treatment' statement it says 'variable treatment not found'.

/*Create dataset with ODS*/
ods output ParameterEstimates=logest;
/*Run logistic regression specify descending to find P(event=1)*/
proc logistic data=subset2 descending;
by sampleno;
model event = treatment;
run;
/*Unsuppress ODS listing*/
ods exclude none;
ods results;
options notes;
/*Create Power criteria*/
proc format;
value signif
low-0.05 = 'significant'
0.05<-high = 'not significant';
run;
/*Display Power results*/
proc freq data=logest;
format ProbChiSq signif.;
tables ProbChiSq;
run;
/*calculate means of sampling dist for log regression */
proc means data=logest mean stddev;
var estimate;
run;

 Question 2

I am also trying to do the logistic regression with a glm to get a RR estimate so in place of logistic regression command I have written:

/*Run log-binomial regression to look at RR*/
proc genmod data=subset descending;
model event = treatment/dist=binomial link=log;
run;

however it seems to only run 1 sample simulation not the 10,100, or 1000 specified. AN ideas why it might be doing this?

Thanks heaps,

Ben

[FreelanceReinh]

You're welcome.

Re question 1: You can use a WHERE statement in the PROC FREQ step in order to restrict the p-values to those pertaining to variable treatment. Similarly, a CLASS statement in the PROC MEANS step would separate the estimates for intercept and treatment.

Example:

/* Create sample data for demonstration */

data subset2;
call streaminit(3141592);
do sampleno=1 to 1000;
  do _n_=1 to 100;
    treatment=_n_<=50;
    event=rand('bern',logistic(-1.789+1.234*treatment));
    output;
  end;
end;
run;

/* ... ODS and OPTIONS statements, PROC LOGISTIC and PROC FORMAT step as in your code ... */

proc freq data=logest;
where variable='treatment';
format ProbChiSq signif.;
tables ProbChiSq;
run;

proc means data=logest mean stddev;
class variable;
var estimate;
run;

Re question 2: Insert the "by sampleno;" statement into the PROC GENMOD step as you did with PROC PHREG and PROC LOGISTIC. (Also, note that the ParameterEstimates ODS table of PROC GENMOD has a different structure than that of PROC LOGISTIC. In particular, the variable containing 'Intercept', 'treatment', etc. is named Parameter, not Variable.)

[Ben12345]

Once again thankyou so much! Your a legend!