|
|
|
||||
|
|
|||||
Figure 11.1, page 359
proc sort data = 'c:\alda\firstsex';
by pt;
run;
proc lifetest data = 'c:\alda\firstsex';
time time*censor(1);
by pt;
ods output ProductLimitEstimates = tsex;
run;
data fig11_1b;
set tsex (where = (survival~=.));
if time = 0 then time = 6;
lags = lag(survival);
by pt;
if ~ first.pt then do;
hazard = 1 - survival/lags;
end;
keep survival hazard time pt;
run;
goptions reset=all device=gif570;
symbol1 c=red i=join;
symbol2 c=blue i=join width=2;
axis1 order=(0 to .5 by 0.1) label=(a=90 'Estimated hazard probability');
axis2 order = (0 to 1 by .5) label=(a=90 'Estimated survival probability');
axis3 order = (6 to 12 by 1) minor = none label=('Grade');
proc gplot data=fig11_1b;
format time 2.0 survival hazard 3.1;
plot hazard*time = pt / vaxis=axis1 haxis=axis3;
plot survival*time = pt / vaxis=axis2 haxis=axis3 vref=0.5 lvref=21;
run;
quit;
Table 11.1, page 360. We will do create the table in two parts. The first part is based on the result from previous example. The data set tsex came from proc lifetest from previous example.
options nocenter nodate; data table11_1by_pt; set tsex (where = (survival~=.)); format left 4.0; myleft = lag(left); myfail = lag(failed); if time = 0 then time = 6; lags = lag(survival); by pt; if ~ first.pt then do; failed = failed - myfail; hazard = 1 - survival/lags; end; keep survival hazard time pt myleft failed; run;
proc print data = table11_1by_pt noobs; where hazard ~=.; format time 3.0; var time myleft failed hazard survival; run;
TIME myleft Failed hazard Survival
7 72 2 0.02778 0.9722 8 70 2 0.02857 0.9444 9 68 8 0.11765 0.8333 10 60 8 0.13333 0.7222 11 52 10 0.19231 0.5833 12 42 8 0.19048 0.4722 7 108 13 0.12037 0.8796 8 95 5 0.05263 0.8333 9 90 16 0.17778 0.6852 10 74 21 0.28378 0.4907 11 53 15 0.28302 0.3519 12 38 18 0.47368 0.1852
Now, we will create the part of Table 11.1 for overall survival regardless the parenting transitions.
proc lifetest data = alda.firstsex; time time*censor(1); ods output ProductLimitEstimates = tall; run; data table11_1all; set tall (where = (survival~=.)); format left 4.0; myleft = lag(left); myfail = lag(failed); if time = 0 then time = 6; lags = lag(survival); failed = failed - myfail; hazard = 1 - survival/lags; keep survival hazard time pt myleft failed; run; proc print data = table11_1all noobs; where hazard ~=.; format time 3.0; var time myleft failed hazard survival; run;
TIME myleft Failed hazard Survival
7 180 15 0.08333 0.9167 8 165 7 0.04242 0.8778 9 158 24 0.15190 0.7444 10 134 29 0.21642 0.5833 11 105 25 0.23810 0.4444 12 80 26 0.32500 0.3000
Figure 11.2, page 363
Manipulate the data set fig11_1b containing the hazard function created for Figure 11.1 in order to create the estimated odds and the estimated logit(hazard).
data fig11_2;
set fig11_1b;
odds = hazard/(1-hazard);
logith = log(odds);
run;
goptions reset=all ;
symbol1 c=red i=join;
symbol2 c=blue i=join width=2;
axis1 order=(0 to 1 by 0.2) minor = none label=(a=90 'Estimated hazard');
axis2 order=(0 to 1 by 0.2) minor = none label=(a=90 'Estimated odds');
axis3 order=(-4 to 0 by 1) minor = none label=(a=90 'Estimated logit(hazard)');
axis4 order =(7 to 12 by 1) minor = none label=('Grade');
proc gplot data=fig11_2;
format time 3.0;
plot hazard*time = pt / vaxis=axis1 haxis=axis4;
plot odds*time = pt / vaxis=axis2 haxis = axis4;
plot logith*time = pt / vaxis=axis3 haxis = axis4;
run;
quit;
Figure 11.3 on page 366.Panel A: logit hazard is horizontal with time.
proc logistic data = alda.firstsex_pp descending; model event = pt; run;
Analysis of Maximum Likelihood Estimates
Standard Wald Parameter DF Estimate Error Chi-Square Pr > ChiSq
Intercept 1 -2.1493 0.1714 157.2166 <.0001
PT 1 0.7131 0.2084 11.7049 0.0006
data fig11_3a;
set fig11_1b;
if pt = 0 then do;
lp0 = -2.14931;
logit_h0 = log(hazard/(1-hazard));
end;
if pt = 1 then do;
lp1 = -2.14931 + 0.71314;
logit_h1 = log(hazard/(1-hazard));
end;
run;
goptions reset = all device=gif570;
axis1 order = (6 to 12 by 1) minor = none;
axis2 order = (-4 to 0 by 1) minor = none label=none;
symbol1 c=blue v=dot i=none;
symbol2 c=red v=circle i= none;
symbol3 c=blue i=j value=none;
symbol4 c=red i=j l=4 value=none;
proc gplot data=fig11_3a;
plot (logit_h0 logit_h1 lp0 lp1)*time
/overlay legend vaxis=axis2 haxis=axis1;
run;
quit;
Panel B: logit hazard is linear with time.
proc logistic data = alda.firstsex_pp descending; model event = pt period; run;
Analysis of Maximum Likelihood Estimates
Standard Wald Parameter DF Estimate Error Chi-Square Pr > ChiSq
Intercept 1 -6.3048 0.6703 88.4689 <.0001
PT 1 0.8753 0.2169 16.2771 <.0001
PERIOD 1 0.4300 0.0641 45.0418 <.0001
data fig11_3b;
set fig11_1b;
if pt = 0 then do;
lp0 = -6.3048 +.43*time;
logit_h0 = log(hazard/(1-hazard));
end;
if pt = 1 then do;
lp1 = -6.3048 + .8753 +.43*time;
logit_h1 = log(hazard/(1-hazard));
end;
run;
proc gplot data=fig11_3b;
plot (logit_h0 logit_h1 lp0 lp1)*time
/overlay legend vaxis=axis2 haxis=axis1;
run;
quit;
Panel C: Logit hazard is completely general with time.
proc logistic data = alda.firstsex_pp descending;
model event = pt d7-d12;
output out = pred xbeta=g;
run;
proc sort data = pred;
by pt period;
run;
proc sort data = fig11_1b;
by pt time;
run;
data fig11_3c;
merge pred (rename = (period = time)) fig11_1b;
by pt time;
run;
data fig11_3c2;
set fig11_3c;
if pt = 0 then do;
lp0 = g;
logit_h0 = log(hazard/(1-hazard));
end;
if pt = 1 then do;
lp1 = g;
logit_h1 = log(hazard/(1-hazard));
end;
run;
proc gplot data=fig11_3c2; plot (logit_h0 logit_h1 lp0 lp1)*time /overlay legend vaxis=axis2 haxis=axis1; run; quit;
Figure 11.4, page 374
The model that we ran for Panel C of Figure 11.3 will produce the same plot for Panel A here. We will continue to use the data set called pred created in the previous example.
data fig11_4; set pred; if pt = 0 then do; lhz0 = g; odds0 = exp(g); hazard0 = odds0/(1+odds0); end; if pt = 1 then do; lhz1 = g; odds1 = exp(g); hazard1 = odds1/(1+odds1); end; rename period = time; run;
axis1 order = (6 to 12 by 1) minor = none label=('Grade');
axis2 order = (-4 to 1 by 1) minor = none label=(a=90 'Logit(hazard)');
axis3 order = (0.0 to 0.8 by .2) minor = none label=('Odds');
axis4 order = (0.0 to 0.5 by .1) minor = none label=('Hazard');
symbol1 c=blue i=j value=none;
symbol2 c=red i=j l=4 value=none;
legend label=none value=(height=1 font=swiss 'Pt=0' 'Pt=1' )
position=(top left inside) mode=share cborder=black;
proc gplot data=fig11_4;
plot (lhz0 lhz1)*time
/ legend = legend1 overlay vaxis=axis2 haxis=axis1;
plot (odds0 odds1)*time
/ legend = legend1 overlay vaxis=axis3 haxis=axis1;
plot (hazard0 hazard1)*time
/ legend = legend1 overlay vaxis=axis4 haxis=axis1;
run;
quit;
Table 11.3, page 386
* Model A ; proc logistic data="c:\alda\firstsex_pp" descending ; model event = d7 d8 d9 d10 d11 d12 / noint ; run;
Model Fit Statistics
Without With
Criterion Covariates Covariates
AIC 1139.534 663.955
SC 1139.534 692.226
-2 Log L 1139.534 651.955
Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 487.5786 6 <.0001
Score 421.4842 6 <.0001
Wald 277.1388 6 <.0001
Analysis of Maximum Likelihood Estimates
Standard Wald
Parameter DF Estimate Error Chi-Square Pr > ChiSq
D7 1 -2.3979 0.2697 79.0611 <.0001
D8 1 -3.1167 0.3862 65.1115 <.0001
D9 1 -1.7198 0.2217 60.2016 <.0001
D10 1 -1.2867 0.2098 37.6195 <.0001
D11 1 -1.1632 0.2291 25.7699 <.0001
D12 1 -0.7309 0.2387 9.3751 0.0022
* Model B ;
proc logistic data="c:\alda\firstsex_pp" descending ;
model event = d7 d8 d9 d10 d11 d12 pt / noint ;
pt: test pt;
run;
Model Fit Statistics
Without With
Criterion Covariates Covariates
AIC 1139.534 648.662
SC 1139.534 681.644
-2 Log L 1139.534 634.662
Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 504.8722 7 <.0001
Score 429.6738 7 <.0001
Wald 271.0944 7 <.0001
Analysis of Maximum Likelihood Estimates
Standard Wald
Parameter DF Estimate Error Chi-Square Pr > ChiSq
D7 1 -2.9943 0.3175 88.9379 <.0001
D8 1 -3.7001 0.4206 77.4055 <.0001
D9 1 -2.2811 0.2724 70.1309 <.0001
D10 1 -1.8226 0.2585 49.7268 <.0001
D11 1 -1.6542 0.2691 37.7872 <.0001
D12 1 -1.1791 0.2716 18.8484 <.0001
PT 1 0.8736 0.2174 16.1471 <.0001
Linear Hypotheses Testing Results
Wald
Label Chi-Square DF Pr > ChiSq
pt 16.1471 1 <.0001* Model C ; proc logistic data="c:\alda\firstsex_pp" descending ; model event = d7 d8 d9 d10 d11 d12 pas / noint ; pas: test pas; run;
Model Fit Statistics
Without With
Criterion Covariates Covariates
AIC 1139.534 651.169
SC 1139.534 684.151
-2 Log L 1139.534 637.169
Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 502.3654 7 <.0001
Score 428.8665 7 <.0001
Wald 272.9013 7 <.0001
Standard Wald
Parameter DF Estimate Error Chi-Square Pr > ChiSq
D7 1 -2.4646 0.2741 80.8416 <.0001
D8 1 -3.1591 0.3890 65.9498 <.0001
D9 1 -1.7297 0.2245 59.3638 <.0001
D10 1 -1.2851 0.2127 36.5178 <.0001
D11 1 -1.1360 0.2324 23.8998 <.0001
D12 1 -0.6421 0.2428 6.9921 0.0082
PAS 1 0.4428 0.1140 15.0984 0.0001
Linear Hypotheses Testing Results
Wald
Label Chi-Square DF Pr > ChiSq
pas 15.0984 1 0.0001* Model D ; proc logistic data="c:\alda\firstsex_pp" descending ; model event = d7 d8 d9 d10 d11 d12 pt pas / noint ; pt: test pt; pas: test pas; run;
Model Fit Statistics
Without With
Criterion Covariates Covariates
AIC 1139.534 645.147
SC 1139.534 682.841
-2 Log L 1139.534 629.147
Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 510.3870 8 <.0001
Score 432.4782 8 <.0001
Wald 269.8099 8 <.0001
Analysis of Maximum Likelihood Estimates
Standard Wald
Parameter DF Estimate Error Chi-Square Pr > ChiSq
D7 1 -2.8932 0.3206 81.4252 <.0001
D8 1 -3.5847 0.4231 71.7689 <.0001
D9 1 -2.1502 0.2775 60.0588 <.0001
D10 1 -1.6932 0.2647 40.9314 <.0001
D11 1 -1.5177 0.2757 30.2938 <.0001
D12 1 -1.0099 0.2811 12.9040 0.0003
PT 1 0.6605 0.2367 7.7855 0.0053
PAS 1 0.2964 0.1254 5.5872 0.0181
Linear Hypotheses Testing Results
Wald
Label Chi-Square DF Pr > ChiSq
pt 7.7855 1 0.0053
pas 5.5872 1 0.0181
Table 11.4, page 388Interpreting the results of fitting an initial discrete-time hazard model including the main effect of time. Expressing parameter estimates as fitted odds and fitted hazard probabilities. Note that SAS, by default, will produce the parameter estimate and fitted odds without any extra effort. However, you can use the steps below to capture the output using ODS and then manipulate it to create the table below like Table 14.1 which includes the Fitted Hazard.
* Model A, repeated but capturing parameter estimates in file parest ; ods output ParameterEstimates=parest ; proc logistic data="c:\alda\firstsex_pp" descending ; model event = d7 d8 d9 d10 d11 d12 / noint ; run; ods output close;
<Output suppressed to save space>
* Computing odds and hazard;
data table11_4;
set parest;
FittedOdds = exp(Estimate);
Hazard = 1 / (1 + exp(-Estimate)) ;
run;
proc print data=table11_4;
var Variable Estimate FittedOdds Hazard ;
run;
Fitted
Obs Variable Estimate Odds Hazard
1 D7 -2.3979 0.09091 0.08333
2 D8 -3.1167 0.04430 0.04242
3 D9 -1.7198 0.17910 0.15190
4 D10 -1.2867 0.27619 0.21642
5 D11 -1.1632 0.31250 0.23810
6 D12 -0.7309 0.48148 0.32500
Figure 11.6, page 393
proc logistic data = alda.firstsex_pp descending;
model event = pt d7-d12 /noint;
output out = pred (keep= period pt g) xbeta=g;
run;
proc sort data = pred nodup out=fig11_6a;
by pt period g;
run;
data fig11_6all;
if _n_ = 1 then do;
period = 6;
pt = 0;
output;
period = 6;
pt = 1;
output; end;
set fig11_6a;
output;
run;
proc sort data = fig11_6all;
by pt period;
run;
data fig11_6a1;
retain survivor0 survivor1 ;
set fig11_6all;
if pt = 0 and period = 6 then
survivor0 = 1 ;
if pt = 0 and period ~= 6 then do;
lhz0 = g;
odds0 = exp(g);
hazard0 = odds0/(1+odds0);
survivor0 = survivor0*(1-hazard0);
end;
if pt = 1 and period = 6 then do;
survivor1 = 1 ;
survivor0 = .; end;
if pt = 1 and period ~=6 then do;
lhz1 = g;
odds1 = exp(g);
hazard1 = odds1/(1+odds1);
survivor1 = survivor1*(1-hazard1);
end;
run;
goptions reset=all;
symbol1 c=blue i=j width=2;
symbol2 c=red i=join;
axis1 order=(-4 to 0 by 1) label=(a=90 'Logit(hazard)');
axis2 order=(0 to 0.5 by 0.1) label=(a=90 'Hazard');
axis3 order=(0 to 1 by 0.5) label=(a=90 'Fitted survival probability');
axis4 order=(6 to 12 by 1) label=('Grade') minor = none;
legend1 label=none value=(height=1 font=swiss 'Pt=0' 'Pt=1' )
position=(top left inside) mode=share cborder=black;
legend2 label=none value=(height=1 font=swiss 'Pt=0' 'Pt=1' )
position=(bottom left inside) mode=share cborder=black;
proc gplot data=fig11_6a1;
plot (lhz0 lhz1)*period / overlay vaxis=axis1 haxis=axis4 legend=legend1;
plot (hazard0 hazard1)*period / overlay vaxis=axis2 haxis=axis4 legend=legend1;
plot (survivor0 survivor1)*period / overlay vaxis=axis3 haxis=axis4 legend=legend2
vref=.5 lvref=21 href=9.9 11.8 lhref=21;
run;
quit;
Figure 11.7, page 395Presents the prototypical hazard and survivor functions for model D in Table 11.3.
Fitting the logistic model with main effects pt and pas and outputting the parameter estimates to the data set estimate.
proc logistic data=datasets.firstsex_pp descending out=estimate noprint; model event = d7-d12 PT PAS/noint; run;
Reconstructing fitted hazard and survivor functions.
data estimate (replace=yes);
set estimate;
array time[6] d7-d12;
do group = 0 to 1;
do antisocial = -1 to 1;
survivor = 1;
do period = 1 to 6;
grade = period + 6;
x = time[period] + group*PT +antisocial*PAS;
hazard = 1/(1 + exp(-x));
survivor = (1 - hazard)*survivor;
output;
end;
end;
end;
keep group antisocial grade x hazard survivor;
run;
Adding initial observations where grade=6 for each combination of pt and pas; generating a variable called ag for each combination of levels for pt and pas for graphing purposes and finally sorting the data.
data temp;
input survivor antisocial group grade;
cards;
1 1 0 6
1 0 0 6
1 -1 0 6
1 1 1 6
1 0 1 6
1 -1 1 6
;
run;
data estimate;
set temp estimate;
if antisocial=1 and group=0 then ag = 1;
if antisocial=0 and group=0 then ag = 2;
if antisocial=-1 and group=0 then ag = 3;
if antisocial=1 and group=1 then ag = 4;
if antisocial=0 and group=1 then ag = 5;
if antisocial=-1 and group=1 then ag = 6;
run;
proc sort data=estimate;
by group antisocial;
run;
goptions reset=all;
symbol1 c=blue i=j width=2;
symbol2 c=blue i=j width=2 l=4;
symbol3 c=blue i=j width=2 l=3;
symbol4 c=red i=j;
symbol5 c=red i=j l=4;
symbol6 c=red i=j l=3;
axis1 order=(0 to 0.5 by 0.1) label=(a=90 'Fitted hazard');
axis2 order=(0 0.5 1.0) label=(a=90 'Fitted survival probabilities');
legend1 label=none value=(height=1 font=swiss
'pas=1, pt=0' 'pas=0, pt=0' 'pas=-1, pt=0' 'pas=1, pt=1' 'pas=0, pt=1' 'pas=-1, pt=1')
position=(top left inside) mode=share cborder=black;
legend2 label=none value=(height=1 font=swiss
'pas=1, pt=0' 'pas=0, pt=0' 'pas=-1, pt=0' 'pas=1, pt=1' 'pas=0, pt=1' 'pas=-1, pt=1')
position=(bottom left inside) mode=share cborder=black;
proc gplot data=estimate;
plot hazard*grade=ag / overlay vaxis=axis1 legend=legend1;
plot survivor*grade=ag / overlay vaxis=axis2 legend=legend2 vref=0.5 lvref=21;
run;
quit;UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services
