FILENAME goptions URL "http://192.38.117.59/~linearpredictors/datafiles/goptions.sas";
%include goptions;

FILENAME url URL "http://192.38.117.59/~linearpredictors/datafiles/readFibrosis.sas";
%include url;

ods output Estimates=or23vs01;

data fibrosis;
set fibrosis;    

if stage<0 then delete;
if ha<0 then delete;
if ykl40<0 then delete;
if p3np<0 then delete;

log2ykl40=log2(ykl40);
log2p3np=log2(p3np);
log2ha=log2(ha);

/* Definition of standardized covariates */
/* PI/sqrt(3) is the SD of the Bernoulli distribution */
stlog2ykl40=log2(ykl40)/(CONSTANT('PI')/sqrt(3));
stlog2p3np=log2(p3np)/(CONSTANT('PI')/sqrt(3));
stlog2ha=log2(ha)/(CONSTANT('PI')/sqrt(3));

fibrosis3=(stage>2);
fibrosis23=(stage>1);
fibrosis123=(stage>0);
RUN;

title1 'SD of standardized covariates';
title2 'Note: Explanation of 1SD effect in program';
PROC TABULATE DATA=fibrosis; 
    VAR stlog2ha stlog2p3np stlog2ykl40;
    TABLE (stlog2ha stlog2p3np stlog2ykl40)*(N*(f=5.0) STDDEV*(f=8.3));
RUN;

ods output Estimates=propodds;

title1 'Estimates from table on the last page';
proc genmod data=fibrosis descending;
  model stage=log2ha log2p3np log2ykl40
              / dist=multinomial link=cumlogit;
estimate 'ha doubling OR' log2ha 1 / exp;
estimate 'p3np doubling OR' log2p3np 1 / exp;
estimate 'ykl40 doubling OR' log2ykl40 1 / exp;
estimate 'log2ha 1SD OR' log2ha 1.073 / exp;
estimate 'log2p3np 1SD OR' log2p3np 0.636 / exp;
estimate 'log2ykl40 1SD OR' log2ykl40 0.691 / exp;
run;

DATA propodds;
SET propodds;
    
type=substr(Label,1,3);
RUN;

title 'Estimates with confidence intervals';
PROC PRINT DATA=propodds; where type ne 'Exp';
    VAR Label  ProbChiSq;
RUN;

title1 'ORs with confidence intervals';
title2 'Note: Explanation of 1SD effect in program';
PROC PRINT DATA=propodds NOOBS; where type='Exp';
    VAR Label LBetaEstimate LBetaLowerCL LBetaUpperCL;
RUN;

title1 'and P-values';
title2;
PROC PRINT DATA=propodds NOOBS; where type ne 'Exp';
    VAR Label  ProbChiSq;
RUN;