%macro finney(n, t, outval);
/* compute Finney's factor, using Gilbert equ 13.4, p 165 */
/*  infinite sum continued far enough so that last term is < 1E-7 */
 
XXfinn =  1 + (&n-1)*&t/&n;
XXdelta =  (&n-1)**3 * &t*&t/(2*&n*&n*(&n+1));

XXi =  3;

do until(XXdelta < 0.0000001);
  XXfinn = XXfinn + XXdelta;
  XXdelta = XXdelta * (&n-1)*(&n-1)*&t/(XXi*&n*(&n+2*XXi-3));
  
  XXi =  XXi + 1;
  end;

&outval =  XXfinn;
%mend;
  
%macro lr(dataset, chemid, valid, censid, censval);

data XXnew;
  set &dataset;
    
  if &censid =  &censval then do;
    lower = .;
    upper = &valid;
    end;
   else do;
    lower = &valid;
    upper = &valid;
    end;
    
proc sort data = XXnew;
   by &chemid;
   
proc lifereg noprint outest = XXest covout;
   by &chemid;

   model (lower,upper) = /d = normal ;

data XXest2;
  set XXest;
  by &chemid;
  
  retain mean sd se;
  
  if first.&chemid then do;
    mean =  .;
    sd =  .;
    se = .;
    end;
    
  if _type_ =  "PARMS" then do;
    mean = intercep;
    sd =  _scale_;
    end;
  if (_type_ = "COV") and (_name_ = "INTERCPT") then do;
    se = sqrt(intercep);
    output;
    end;
    
  keep &chemid mean sd se;

data XXest3;
  set XXest2;
  
  lmean =  mean;
  lsd =  sd;
  lse =  se;
  
/* calculate Finney's adjustment using effective N =  (sd/se) ^ 2 */

  effN =  (lsd/lse)**2;
  t =  lsd*lsd/2;
  %finney(effN,t, psi);
    
  mean1 = exp(lmean+lsd**2/2);
  mean2 =  exp(lmean+lsd**2/2 - lse**2/2);
  mean3 = exp(lmean)*psi;
  
  gmean =  exp(lmean);
  
  cv =  100*sqrt(exp(lsd**2)-1);

  sd1 = sqrt(exp(2*lmean + 2*lsd**2) - exp(2*lmean + lsd**2));
  sd2 = sqrt(exp(2*lmean+2*lsd**2 - lse**2) - exp(2*lmean + lsd**2 - lse**2));
  
  t2 = lsd*lsd*2;
  %finney(effN,t2, psi2);

  t3 = lsd*lsd*(effN-2)/(effN-1);
  %finney(effN,t3,psi3);
    
  sd3 = sqrt(exp(2*lmean)*(psi2-psi3));
  
  se =  sqrt(exp(2*lmean + 2*lse**2)-exp(2*lmean +  lse**2));  
/*  se3 = sqrt(exp(2*lmean)*(psi** 2 - psi2));
     /* Gilbert's version of Bradu and Mundlak enbiased est of var of mean3 */

  t4 = effN*(lsd*lsd-lse*lse)/(2*(effN-1));
  %finney(effN,t4, psi4);
  t5 = effN*(lsd*lsd-lse*lse)/(effN-1);
  %finney(effN,t5, psi5);

  se3 = sqrt(exp(2*lmean)*(psi4*psi4-psi5));
   
%mend;
       

%macro logn(datain, dataout);
/* compute backtransformed estimates of mean, etc. */

/*  datain: name of input data set.  */
/*    Must contain variables: mean sd se */

/*  dataout: name of output data set.  */
/*    Contains all input variables plus: */
/*      lmean, lsd, lse: old mean, sd, and se */;
/*      mean: arithmetic mean.  estimated by exp(lmean+(lsd**2 -lse**2)/2) */
/*      gmean: geometric mean.  estimated by exp(lmean - lse**2/2) */
/*      cv.: coef. of var. estimated by 100*sqrt(exp(lsd**2)-1)) */
/*      sd: estimated by sqrt(exp(
/*        2 lmean + 2 lsd**2 - lse**2) - exp(2 lmean + lsd**2 - lse**2)))) */
/*      s.e: est.by sqrt(exp(2 lmean + 2 lse**2)-exp(2 lmean +  lse**2)))) */

data &dataout;
  set &datain;
  
  lmean =  mean;
  lsd =  sd;
  lse =  se;
  

  mean =  exp(lmean+lsd**2/2 - lse**2/2);
  gmean =  exp(lmean);
  
  cv =  100*sqrt(exp(lsd**2)-1);

  sd = sqrt(exp(2*lmean + 2*lsd**2 - lse**2) - exp(2*lmean + lsd**2 - lse**2));
  
  se =  sqrt(exp(2*lmean + 2*lse**2)-exp(2*lmean +  lse**2));
  

%mend;