/* macros for ML estimation of mean and variance observations */
/*   with below detection limit values */

/* These are research quality programs made available for public use */
/* No warranty or guarantee of any kind is made about correctness or */
/*   appropriateness */

/* written for version 9 of SAS */

/* Philip Dixon,  21 Jan 2004 */


%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, outdata);

/* Estimate mean and variance of censored log normal values */
/*   based on ml estimates of parameters of log transf. obs */
/* set up to process many groups of observations in sequence */
/*   using SAS BY group processing */

/* Input: dataset: name of SAS data set containing the data */
/*        chemid:  name of variable defining BY groups */
/*        valid:   name of variable containing log transformed obs */
/*           N.B. input must be log base e (i.e. natural log) transformed */
/*           prior to calling this procedure */
/*        censid:  name of variable indicating whether obs is censored */
/*        censval: value of &censid that indicates a censored obs */
/*        outdata: name of SAS data set to contain the output */
/*  a typical data set called WELLDATA might contain: */
/*  chemical logY belowDL */
/*    Ar      -3.1 0 */
/*    Ar      -2.0 1 */
/*    Ar       0.3 0 */
/*    Mo       1.2 0 */
/*    Mo       0.5 1 */
/* The second and fifth obs are censored.  The DL's are exp(-2.0) and exp(0.5) */
/* The appropriate call is then:
%lr(welldata, chemical, logY, belowDL, 1, wellmean);  */

/* Output: lmean, lsd, lse: estimates of mean, sd, and se for log response. */
/*         cv: coefficient of variation of responses  (not log response) */

/*         mean1, mean2, mean3: estimates of mean response after backtransformation */
/*         sd1, sd2, sd3: estimates of sd of response after backtransf. */
/*         se1,      se3: estimates of se of mean response, after backtransf */

/*  choice of estimator not obvious with censored data, so I calculate 3 */
/*    mean1, sd1, se1: naive backtransformation, assuming log normal */
/*    mean2, sd2: backtransf, adjusting for uncertainty in mean */
/*    mean3, sd3, se3: Finney's backtransformation, adjusting for uncert. in mean and variance */

/* I don't compute se for estimator 2 since I don't know how to do it. */
/* The use of effective N in Finney's method is a practical solution without */
/*   theoretical justification.  Sometime I'll write up and publish the details */
 
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 = intercept;
    sd =  _scale_;
    end;
  if (_type_ = "COV") and (_name_ = "Intercept") then do;
    se = sqrt(intercept);
    output;
    end;
    
  keep &chemid mean sd se;

data &outdata;
  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));
  
  se1 =  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));
  drop mean sd se;
   
%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;