data calib;
    set elisa_cs;
    if plate='A' and series=1;
    array y(3) rep1 rep2 rep3;
    do j=1 to 3;
        signal=y(j);
        w=1/signal**1.3;
    output;
    end;
proc mixed data=calib method=reml;
    model signal=concentration/solution outpredm=predict;
    weight w;    
    ods output solutionf=parameter_estimates;
proc sort data=predict;
    by concentration;
axis1 minor=none label=(angle=90 'Signal') order=(0 to 4 by 1);
axis2 minor=none label=('Concentration') logbase=10 logstyle=expand;    
symbol1 value=none i=join color=black line=1 width=3;
symbol2 value=dot i=none color=black height=5;
proc gplot data=predict;
    plot (pred signal)*concentration/overlay frame haxis=axis2 vaxis=axis1;
    run;
    quit;



proc nlin data=calib method=marquardt outest=parameter_estimates;
    parameters top=3 bottom=0.2 c50=250 slope=1;
    model signal=top+(bottom-top)/(1+(concentration/c50)**slope);
    _weight_=w;
    output out=fitted_values predicted=pred;
proc sort data=fitted_values;
    by concentration;
axis1 minor=none label=(angle=90 'Signal') order=(0 to 4 by 1);
axis2 minor=none label=('Concentration') logbase=10 logstyle=expand;    
symbol1 value=none i=join color=black line=1 width=3;
symbol2 value=dot i=none color=black height=5;
proc gplot data=fitted_values;
    plot (pred signal)*concentration/overlay frame haxis=axis2 vaxis=axis1;
    run;
    quit;



proc nlmixed data=calib technique=nrridg;
    parms top=3 bottom=0.2 c50=250 slope=1 theta=1 var=0.0001;
    expect=top+(bottom-top)/(1+(concentration/c50)**slope);
    model signal~normal(expect,(expect**theta)*var);
    predict top+(bottom-top)/(1+(concentration/c50)**slope) out=fitted_values;
    run;



proc nlmixed data=calib technique=nrridg;
    parms top=3 bottom=0.2 c50=250 slope=1 theta=1 var=0.0001 g=1;
    bounds g>=1, g<=1, theta>=0, var>0; * Constraints on model parameters;
    expect=top+(bottom-top)/((1+(concentration/c50)**slope)**g);
    model signal~normal(expect,(expect**theta)*var);
    predict top+(bottom-top)/((1+(concentration/c50)**slope)**g) out=fitted_values;
    ods output ParameterEstimates=parm_est_repl;
    ods output CovMatParmEst=cov_parm;
* Data set containing parameter estimates;
data b;
    set parm_est_repl;
    where parameter in ('top','bottom','c50','slope');
    keep estimate;
* Data set containing the covariance matrix of the parameter estimates;
data covb;
    set cov_parm;
    where parameter in ('top','bottom','c50','slope');
    keep top bottom c50 slope;
proc sql noprint;
    select distinct estimate into: sigma_sq from parm_est_repl where parameter='var';
    select distinct estimate into: theta from parm_est_repl where parameter='theta';
    select min(concentration) into: minconc from calib where concentration>0;
    select max(concentration) into: maxconc from calib;
    run;
    quit;
%let m=3; * Number of replicates at each concentration level;
proc iml;
    * Import the data sets;
    use b; read all into b;
    use covb; read all into covb;
    * Initialize matrices;
    y=j(101,1,0);
    h=j(101,1,0);
    hy=j(101,1,0);
    hb=j(101,4,0);
    varx0=j(101,1,0);
    pp=j(101,1,0);
    top=b[1];
    bottom=b[2];
    c50=b[3];
    slope=b[4];
    * Calculate the precision profile;
    do i=1 to 101;
        h[i,1]=10**(log10(&minconc)+(i-1)*(log10(&maxconc)-log10(&minconc))/100);
        y[i,1]=top+(bottom-top)/(1+(h[i,1]/c50)**slope);
        hy[i,1]=h[i,1]*(top-bottom)/(slope*(bottom-y[i,1])*(y[i,1]-top));
        hb[i,1]=h[i,1]/(slope*(y[i,1]-top));
        hb[i,2]=h[i,1]/(slope*(bottom-y[i,1]));
        hb[i,3]=h[i,1]/c50;
        hb[i,4]=-h[i,1]*log((bottom-y[i,1])/(y[i,1]-top))/(slope**2);
        varx0[i,1]=((hy[i,1]**2)*&sigma_sq*(y[i,1]**(2*&theta))/&m)+hb[i,]*covb*hb[i,]`;
        pp[i,1]=100*sqrt(varx0[i,1])/h[i,1];
    end;
    create plot var{h hy y pp};
    append;
    quit;
axis1 minor=none label=(angle=90 'CV (%)') order=(0 to 30 by 10);
axis2 minor=none label=('Concentration') logbase=10 logstyle=expand;    
symbol1 value=none i=join color=black line=1 width=3;
proc gplot data=plot;
    plot pp*h/frame haxis=axis2 vaxis=axis1 vref=20 lvref=34 href=4.2 lhref=34;
    run;
    quit;


data calib;
    set elisa_cs;
    array y(3) rep1 rep2 rep3;
    do j=1 to 3;
        signal=y(j);
    output;
    end;
proc datasets nolist;
    delete parm_est;
%macro calib(series, plate);
    proc nlmixed data=calib technique=nrridg;
        parms top=3 bottom=0.2 c50=250 slope=1 theta=1 var=0.0001;
        expect=top+(bottom-top)/(1+(concentration/c50)**slope);
        model signal~normal(expect,(expect**(theta*2))*var);
        where series=&series and plate=&plate;
        ods output ParameterEstimates=parm_est_tmp;
    data parm_est_tmp;
        set parm_est_tmp;
        series=&series;
        plate=&plate;
    proc append data=parm_est_tmp base=parm_est;
        run;
%mend;
%calib(series=1, plate='A');
%calib(series=1, plate='B');
%calib(series=2, plate='A');
%calib(series=2, plate='B');
%calib(series=3, plate='A');
%calib(series=3, plate='B');
%calib(series=4, plate='A');
%calib(series=4, plate='B');
data valid;
    set elisa_vs;
    array y(3) rep1 rep2 rep3;
    do j=1 to 3;
        signal=y(j);
        output;
    end;
proc sort data=valid;
    by series plate;
proc transpose data=parm_est out=parm_est_t;
    var estimate;
    id parameter;
    idlabel parameter;
    by series plate;
data calc_conc;
    merge valid parm_est_t;
    by series plate;
    drop _name_ theta var;
data calc_conc;
    set calc_conc;
    if plate='A' then run=series;
    else run=series+4;
    calc_conc=c50*((((bottom-top)/(signal-top))-1)**(1/slope));
    run;


proc sql;
    create table linprof as
    select concentration, calc_conc, 't1' as type
           from calc_conc union
    select distinct concentration, concentration*0.7 as calc_conc, 't2' as type
           from calc_conc union
    select distinct concentration, concentration*1.3 as calc_conc, 't3' as type
           from calc_conc union
    select distinct concentration, concentration as calc_conc, 't4' as type
           from calc_conc;
    run;
    quit;
proc sort data=linprof;
    by type concentration;
axis1 label=("Nominal concentration") order=(0 to 600 by 200);
axis2 label=(angle=90 "Calculated concentration") order=(0 to 600 by 200);
symbol1 value=circle i=none color=black height=5;
symbol2 value=none i=join color=black line=34 width=3;
symbol3 value=none i=join color=black line=34 width=3;
symbol4 value=none i=join color=black line=1 width=3;
proc gplot data=linprof;
    plot calc_conc*concentration=type/haxis=axis1 vaxis=axis2 nolegend;
    run;
    quit; 



proc sort data=calc_conc;
    by concentration;
proc mixed data=calc_conc;
    class run;
    model calc_conc=/solution;
    random run;
    by concentration;
    ods output dimensions=dim;
    ods output solutionf=solf;
    ods output CovParms=var_comp;
proc sort data=calc_conc;
    by concentration run plate;
proc univariate data=calc_conc noprint;
    var calc_conc;
    by concentration run;
    output out=stat_by_conc_run n=n;
proc sql;
    create table stat_by_conc as
    select t1.concentration, cap_n, sum_n_sq/cap_n as n, n_run, mean, 
         sqrt(sigma_w_sq) as sigma_w, 
         100*sqrt(sigma_w_sq)/t1.concentration as cv_w, 
         sqrt(sigma_b_sq) as sigma_b, 
         100*sqrt(sigma_b_sq)/t1.concentration as cv_b, 
         sqrt(sigma_w_sq+sigma_w_sq) as sigma_t, 
         100*sqrt(sigma_w_sq+sigma_w_sq)/t1.concentration as cv_t
    from (select concentration, value as cap_n from dim where 
         descr='Observations Used') as t1,
         (select concentration, estimate as mean from solf) as t2,
         (select concentration, sum(n**2) as 
         sum_n_sq from stat_by_conc_run group by concentration) as t3,
         (select concentration, estimate as sigma_b_sq from var_comp 
         where covparm='run') as t4,
         (select concentration, estimate as sigma_w_sq from var_comp 
         where covparm='Residual') as t5,
         (select count(run) as n_run from (select distinct run from 
         stat_by_conc_run)) as t6 where t1.concentration=t2.concentration=
         t3.concentration=t4.concentration=t5.concentration;
    run;
    quit;
* Calculation of validation criteria;
data stat_by_conc;
    set stat_by_conc;
    format mean re cv_w cv_b cv_t te 4.1 low_tol_lim_rel upp_tol_lim_rel 5.1;
    re=100*(mean-concentration)/concentration;
    te=abs(re)+cv_t;
    r=(sigma_b/sigma_w)**2;
    b=sqrt((r+1)/(n*r+1));
    ip=sigma_b**2+sigma_w**2;
    ddl=((r+1)**2)/(((r+1/n)**2)/(n_run-1)+(1-1/n)/(cap_n));
    low_tol_lim_abs=mean-tinv(0.975,ddl)*sqrt(ip*(1+1/(cap_n*b**2)));
    upp_tol_lim_abs=mean+tinv(0.975,ddl)*sqrt(ip*(1+1/(cap_n*b**2)));
    low_tol_lim_rel=(low_tol_lim_abs-concentration)*100/concentration;
    upp_tol_lim_rel=(upp_tol_lim_abs-concentration)*100/concentration;
    label concentration='Nominal*concentration'
          mean='Calculated*concentration'
          re='Relative*error (%)' 
proc sort data=calc_conc;
    by concentration;
proc mixed data=calc_conc;
    class run;
    model calc_conc=/solution;
    random run;
    by concentration;
    ods output dimensions=dim;
    ods output solutionf=solf;
    ods output CovParms=var_comp;
proc sort data=calc_conc;
    by concentration run plate;
proc univariate data=calc_conc noprint;
    var calc_conc;
    by concentration run;
    output out=stat_by_conc_run n=n;
proc sql;
    create table stat_by_conc as
    select t1.concentration, cap_n, sum_n_sq/cap_n as n, n_run, mean, 
         sqrt(sigma_w_sq) as sigma_w, 
         100*sqrt(sigma_w_sq)/t1.concentration as cv_w, 
         sqrt(sigma_b_sq) as sigma_b, 
         100*sqrt(sigma_b_sq)/t1.concentration as cv_b, 
         sqrt(sigma_w_sq+sigma_w_sq) as sigma_t, 
         100*sqrt(sigma_w_sq+sigma_w_sq)/t1.concentration as cv_t
    from (select concentration, value as cap_n from dim where 
         descr='Observations Used') as t1,
         (select concentration, estimate as mean from solf) as t2,
         (select concentration, sum(n**2) as 
         sum_n_sq from stat_by_conc_run group by concentration) as t3,
         (select concentration, estimate as sigma_b_sq from var_comp 
         where covparm='run') as t4,
         (select concentration, estimate as sigma_w_sq from var_comp 
         where covparm='Residual') as t5,
         (select count(run) as n_run from (select distinct run from 
         stat_by_conc_run)) as t6 where t1.concentration=t2.concentration=
         t3.concentration=t4.concentration=t5.concentration;
    run;
    quit;
* Calculation of validation criteria;
data stat_by_conc;
    set stat_by_conc;
    format mean re cv_w cv_b cv_t te 4.1 low_tol_lim_rel upp_tol_lim_rel 5.1;
    re=100*(mean-concentration)/concentration;
    te=abs(re)+cv_t;
    r=(sigma_b/sigma_w)**2;
    b=sqrt((r+1)/(n*r+1));
    ip=sigma_b**2+sigma_w**2;
    ddl=((r+1)**2)/(((r+1/n)**2)/(n_run-1)+(1-1/n)/(cap_n));
    low_tol_lim_abs=mean-tinv(0.975,ddl)*sqrt(ip*(1+1/(cap_n*b**2)));
    upp_tol_lim_abs=mean+tinv(0.975,ddl)*sqrt(ip*(1+1/(cap_n*b**2)));
    low_tol_lim_rel=(low_tol_lim_abs-concentration)*100/concentration;
    upp_tol_lim_rel=(upp_tol_lim_abs-concentration)*100/concentration;
    label concentration='Nominal*concentration'
          mean='Calculated*concentration'
          re='Relative*error (%)' 
          cv_w='Within-run*CV (%)'
          cv_b='Run-to-run*CV (%)'
          cv_t='Intermediate*precision*CV (%)'
          te='Total*error(%)' 
          low_tol_lim_rel='Upper 95%*tolerance*limit'
          upp_tol_lim_rel='95% TI*upper*limit';
proc print data=stat_by_conc noobs split='*';
    format mean re cv_w cv_b cv_t te 4.1 low_tol_lim_rel upp_tol_lim_rel 5.1;
    var concentration mean re cv_w cv_b cv_t te;    
    run;
          cv_w='Within-run*CV (%)'
          cv_b='Run-to-run*CV (%)'
          cv_t='Intermediate*precision*CV (%)'
          te='Total*error(%)' 
          low_tol_lim_rel='Upper 95%*tolerance*limit'
          upp_tol_lim_rel='95% TI*upper*limit';
proc print data=stat_by_conc noobs split='*';
    format mean re cv_w cv_b cv_t te 4.1 low_tol_lim_rel upp_tol_lim_rel 5.1;
    var concentration mean re cv_w cv_b cv_t te;    
    run;





proc print data=stat_by_conc noobs split='*';    
    var concentration mean low_tol_lim_rel upp_tol_lim_rel;    
    run;
axis1 minor=none label=(angle=90 "Total error (%)") order=(0 to 70 by 10);
axis2 minor=none label=("Nominal concentration") logbase=10 logstyle=expand;
symbol1 value=dot i=join color=black line=1 width=3 height=3;
proc gplot data=stat_by_conc;
    format te 4.0;
    plot te*concentration/frame haxis=axis2 vaxis=axis1 vref=30 lvref=34;
    run;
    quit;
axis1 minor=none label=(angle=90 "Relative error (%)") order=(-120 to 100 by 20);
axis2 minor=none label=("Nominal concentration") logbase=10 logstyle=expand;
symbol1 value=dot i=join color=black line=1 width=3 height=3;
symbol2 value=none i=join color=black line=8 width=3;
symbol3 value=dot i=join color=black line=1 width=3 height=3;
proc gplot data=stat_by_conc;
    format low_tol_lim_rel re upp_tol_lim_rel 4.0;
    plot (low_tol_lim_rel re upp_tol_lim_rel)*concentration/overlay frame 
        haxis=axis2 vaxis=axis1 vref=-30,30 lvref=34;
    run;
    quit;


proc sql;
    create table te as
    select distinct 2 as type, calc_conc.concentration, 
        abs(calc_conc-calc_conc.concentration)*100/calc_conc.concentration as te
    from calc_conc, stat_by_conc
    where calc_conc.concentration=stat_by_conc.concentration
    union
    select 1 as type, concentration, te from stat_by_conc;
    run;
    quit;
axis1 minor=none label=(angle=90 "Total error (%)") order=(0 to 100 by 20);
axis2 minor=none label=("Nominal concentration") logbase=10 logstyle=expand;
symbol1 value=dot i=join color=black line=1 width=3;
symbol2 value=circle i=none color=black height=5;
proc gplot data=te;
    plot te*concentration=type/vaxis=axis1 haxis=axis2 vref=30 lvref=34 nolegend;
    run;
    quit;
data me;
    set stat_by_conc;
    array y(3) low_tol_lim_rel re upp_tol_lim_rel;
    do j=1 to 3;
        me=y(j); type=j;
    output;
    end;
proc sql;
    create table me as
    select distinct 4 as type, calc_conc.concentration, 
        (calc_conc-calc_conc.concentration)*100/calc_conc.concentration as me
    from calc_conc, stat_by_conc
    where calc_conc.concentration=stat_by_conc.concentration
    union
    select type, concentration, me from me;
    run;
    quit;
axis1 minor=none label=(angle=90 "Relative error (%)") order=(-120 to 100 by 20);
axis2 minor=none label=("Nominal concentration") logbase=10 logstyle=expand;
symbol1 value=dot i=join color=black line=1 width=3 height=5;
symbol2 value=none i=join color=black line=20 width=3;
symbol3 value=dot i=join color=black line=1 width=3 height=5;
symbol4 value=circle i=none color=black height=5;
proc gplot data=me;
    plot me*concentration=type/frame haxis=axis2 vaxis=axis1 nolegend vref=-30,30 lvref=34;
    run;
    quit;