options ls=80 nodate;

*change name of path to read in fulldata.dat;
data fulldata;
 infile 'E\GeneralGamma\fulldata.dat';
 input pubid entry exit event p1 p2 p3 p4 per1234;
run;

proc sort;
 by per1234;
run;

* General gamma model with left truncation and right censoring.
   event = 0 for right censored observations, event = 1 for event times
   entry >= 0, time of entry for left truncated observations
   exit > 0, time of event or right censoring;

title 'Final Generalized Gamma Model';

proc nlmixed cov fd;

 parms beta=0  sigma=1;

 if (per1234 = 1) then lambda = 1/sigma;
 if (per1234 = 2) then lambda = sigma;
 if (per1234 = 3) then lambda = 1/sigma;
 if (per1234 = 4) then lambda = 1;

* Log likelihood for uncensored observations, no late entry: log[f(t)];
 if (entry = 0 and event = 1) then
    ll = log(     pdf('gamma',exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) )
                                + lambda*log(exit)/sigma  + log(lambda)  - log(sigma);

* Log likelihood for censored observations, no late entry: log[S(t)];
 if (entry = 0 and event = 0) then
    ll = log( 1 - cdf('gamma',exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );

* Log likelihood for uncensored observations, with late entry: log[f(t)/S(w)];
 if (entry > 0 and event = 1) then
    ll = log(     pdf('gamma', exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) )
                                + lambda*log(exit)/sigma  + log(lambda)  - log(sigma)
       - log( 1 - cdf('gamma',entry**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );

* Log likelihood for censored observations, with late entry: log[S(t)/S(w)];
 if (entry > 0 and event = 0) then
    ll = log( 1 - cdf('gamma', exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) )
       - log( 1 - cdf('gamma',entry**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );

 model exit ~ general(ll);

 estimate 'logsigma' log(sigma);

 by per1234;
run;