* Multivariate Reference Interval Based on Global Smoothing;



%macro GlobalQSmooth(dataset,grid,gamma,p,outdata);

/*

Inputs:

 

DATASET = Data set to be analyzed

 

GRID    = Data set with grid points (X variable)

 

GAMMA   = Quantile of the conditional distribution function

 

P       = Degree of the polynomial approximation

 

OUTDATA = Output data set containing the polynomial (POLY variable),

          its first derivative (FIRSTDER variable) and second derivative

          (SECDER variable) evaluated at the grid points	

*/



data variable;

    set &dataset;

    array cov{*} cov1-cov&p;

    do i=1 to &p; cov{i}=x**i; end;

proc mixed data=variable;

    model y=cov1-cov&p/s;

    ods output SolutionF=param;

data equation;

    set param nobs=m;

    retain c0-c&p;

    if effect="Intercept" then c0=estimate;

    call symput("c0",compress(put(c0,best8.)));

    %do i=1 %to &p;

        if effect="cov&i" then c&i=estimate;

        call symput("c&i",compress(put(c&i,best8.)));

    %end;

    if _n_=m;

    run;

proc nlin data=variable nohalve maxiter=500 converge=0.0005;

    parms c0=&c0 %do i=1 %to &p; c&i=&&c&i %end;;

    model y=c0 %do i=1 %to &p; +c&i*cov&i %end;;

    der.c0=1;

    %do i=1 %to &p; der.c&i=cov&i; %end;

    resid=y-model.y;

    if resid>0 then _weight_=&gamma/resid;

    if resid<0 then _weight_=(&gamma-1)/resid;

    if resid=0 then _weight_=0;

    ods output ParameterEstimates=est;

data coef;

    set est nobs=m;

    retain d0-d&p;

    %do i=0 %to &p;

        if parameter="c&i" then d&i=estimate;

    %end;

    if _n_=m;

data &outdata;

    set &grid;

    if _n_=1 then set coef;

    array d{*} d1-d&p;

    poly=d0; firstder=d1;

    if &p>=2 then secder=2*d2; else secder=0;

    do i=1 to &p;

        poly=poly+d{i}*x**i;

        if i>=2 then firstder=firstder+i*d{i}*x**(i-1);

        if i>=3 then secder=secder+(i-1)*i*d{i}*x**(i-2);

    end;

    run;

%mend GlobalQSmooth;