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 'Full Generalized Gamma Model';

proc nlmixed cov fd;

 parms beta=0  sigma=1  lambda=0.5;

* Log likelihood for uncensored observations, no late entry: log[f(t)];
 if (entry = 0 and event = 1 and lambda < 0) 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);
 if (entry = 0 and event = 1  and  lambda > 0) 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 and lambda < 0) then
    ll = log(     cdf('gamma',exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );
 if (entry = 0 and event = 0  and  lambda > 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 and lambda < 0) 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(     cdf('gamma',entry**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );
 if (entry > 0 and event = 1  and  lambda > 0) 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 and lambda < 0) then
    ll = log(     cdf('gamma', exit**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) )
       - log(     cdf('gamma',entry**(lambda/sigma), 1/(lambda*lambda), lambda*lambda*exp(beta*lambda/sigma)) );
 if (entry > 0 and event = 0  and  lambda > 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;