**********************************************************************
*** LREGMC1 ***
*** last change: 8-14-96 ***
**********************************************************************
* SAS Macro for performing logistic regression when the *
* binary outcome is measured with imperfect sensitivity *
* sensitivity and specificity *
**********************************************************************
Written by
Larry Magder
Dept of Epidemiology and Preventive Medicine
University of Maryland at Baltimore
********************************************************************
;
%MACRO lregmc1
( DATA =_last_,
YVAR =y,
XVAR =x,
SENS = sensvar,
SPEC = specvar,
BETA =0,
ITER =30,
CRIT =0.00001,
MISSDEL=yes,
PRNITER=no);

%*** saving data set name for printed output;
%LET PRINDAT=&DATA;
%* Deleting observations with missing values, if requested *;
%IF %UPCASE(&MISSDEL)=YES
& %UPCASE(&DATA)^=MISSFREE
& %UPCASE(&DATA)^=WORK.MISSFREE
& %UPCASE(&DATA)^=SASWORK.MISSFREE %THEN %DO;
DATA missfree;
SET &DATA;
IF NMISS(OF &YVAR &XVAR &SENS &SPEC ) ^= 0 THEN Delete;
RUN;
%LET DATA=missfree;

%END;

%** Starting IML with size- and output specifications ***;
PROC IML ;
FILE print;
RESET noname nolog linesize=80 nocenter;

%** Handing the macro parameters over to IML ***;
USE &DATA;
SETIN &DATA NOBS nobs;

beta = { &BETA };
labely = { &YVAR };
labelx = { &XVAR }; labelx = labelx`;
labelsns = { &SENS };
labelspc = { &SPEC };
crit=1;
;
%************************ ESTB ROUTINE ********************************;
%*** An updated estimate for beta is calculated **;
%**********************************************************************;
START estb;
score=J(p,1,0);
inf=J(p,p,0);

%*---------------- Start of beta estimation loop ---------------------;
DO i=1 TO n;
%*** Reading the observation for person i ***;
READ VAR labely INTO yvar POINT i;
READ VAR labelx INTO xvar POINT i;
READ VAR labelsns INTO sens POINT i;
READ VAR labelspc INTO spec POINT i;
%*** Calculate phati;
ui=EXP(xvar*beta); phati=ui/(1+ui);
%*** Calculate yhati;
if yvar=1 then do;
yhati=sens*phati/(sens*phati+(1-spec)*(1-phati));
end;
if yvar=0 then do;
yhati=(1-sens)*phati/((1-sens)*phati+spec*(1-phati));
end;
%*** add to score and inf;
score=score+(yhati-phati)#xvar`;
inf=inf+phati#(1-phati)#xvar`*xvar;
END;
%*---------------- End of beta estimation loop ----------------------;
%*** Calculate beta;
beta=beta+INV(inf)*score;
%*** Print actual iterate for beta;
cname={"Estimate"};
%IF %UPCASE(&PRNITER)="YES" %THEN %DO;
PRINT labelx beta [COLNAME=cname];
%END;
FREE i score inf yhati xhati xvar yvar ui;
FINISH;

%*====================== MAIN PROGRAM ================================;
%************************** INIT ROUTINE *******************************
%* Initial printed Output is created. ****
** The individual observation numbers are extracted. ****
** The initial estimate for beta is calculated. ****
%**********************************************************************;
START;


PRINT / "Logistic Regression when outcome is measure with error ",
" and sensitivity and specificity are known",
"========================================================",,
"Data File: &PRINDAT",;

%*** Printing the variables to be analyzed; ***;
rname1 = {"Outcome variable:"};
rname2 = {"Predictors:"};
PRINT labely [ROWNAME=rname1],{&XVAR} [ROWNAME=RNAME2],;
rname1 = {"Sensitivity variable:"};
PRINT labelsns [ROWNAME=rname1], ;
rname1 = {"Specificity variable:"};
PRINT labelspc [ROWNAME=rname1], ;

%*--------------- Input first record ----------------------------;
p=NROW(labelx);

READ VAR labely INTO yvar;
READ VAR labelx INTO xvar;
READ VAR labelsns INTO sens ;
READ VAR labelspc INTO spec ;

IF NCOL(beta)=1 THEN DO;
xty=xvar`*yvar; xtx=xvar`*xvar;
END;
war=J(1,p,1); war=war#(war=xvar);
n=1;
vsum=yvar||xvar||sens||spec;
%*--------------- Start of input loop --------------------------;
DO j=2 TO nobs;

READ VAR labely INTO yvar POINT j;
READ VAR labelx INTO xvar POINT j;
READ VAR labelsns INTO sens POINT j;
READ VAR labelspc INTO spec POINT j;

war=war#(war=xvar);
vsum=vsum+(yvar||xvar||sens||spec);
n=n+1;

%*** Preparing for Initial estimate for beta (Least-Squares) ***;
IF NCOL(beta)=1 THEN DO;
xty=xty+xvar`*yvar;
xtx=xtx+xvar`*xvar;
END;
END;
%*------------------- End of input loop -------------------------------;

%*** Calculating observation means, maximum and minimum individual ***
%**** observation numbers and initial estimate for beta ***;
means=vsum/n;
IF NCOL(beta)=1 THEN beta=SOLVE(xtx,xty);
ELSE beta=SHAPE(beta,p,1,0);


%*** printing observation numbers in dataset ***;
rname1 = {"Number of non-missing obs used:"};
PRINT nobs [ROWNAME=rname1 FORMAT=8.0],;

cname = { &YVAR &XVAR &SENS &SPEC };

%*** Printing description of the dataset;
PRINT "Averages of Outcome variable and Predictors ",,
means [COLNAME=cname],,;
IF ALL(^war) THEN PRINT "*** WARNING: No intercept term in the model!";

PRINT "Initial estimate of regression coefficients:",,
labelx beta;
FREE means j vsum xtx xty xvar yvar;

FINISH;
RUN;
%******************* ESTIMATION routine ****************************
*** iteratively estimating beta, using ESTB ***
*********************************************************************;
START;
%*--------------------- Main iteration -------------------------------;
DO iter=1 TO &ITER WHILE(crit>&CRIT);
rname1={"===> Iteration:"};
%IF %UPCASE(&PRNITER)="YES" %THEN %DO;
PRINT iter [ROWNAME=rname1 FORMAT=3.0];
%END;
save=beta;
RUN estb;
crit=MAX(ABS(1-save/beta));
rname1={"Criterion of Convergence:"};
%IF %UPCASE(&PRNITER)="YES" %THEN %DO;
PRINT crit [ROWNAME=rname1 FORMAT=16.14],;
%END;
END;
%*------------------- End of iterations ------------------------------;

iter=iter-1;
IF iter>=&ITER & crit>&CRIT THEN DO;
PRINT ,," " {"WARNING: No convergence after"}
{&ITER} [FORMAT=3.0] {"iterations."},
"Criterion of convergence: "
crit [FORMAT=16.14] " > &CRIT",;
END;
ELSE DO;
PRINT ,," " {"NOTE: Convergence after"} iter [FORMAT=3.0]
{"iteration(s)."},,;
critw=0;
END;

FREE save ;
crit=1;
FINISH; RUN;
%*********************************************************************
*** Routine to estimate variances and log-likelihood ***
*********************************************************************;
START;
inf=J(p,p,0);
loglike=0;
DO i=1 TO n;
%*** Reading the observation for person i ***;
READ VAR labely INTO yvar POINT i;
READ VAR labelx INTO xvar POINT i;
READ VAR labelsns INTO sens POINT i;
READ VAR labelspc INTO spec POINT i;
%*** Calculate phati;
ui=EXP(xvar*beta); phati=ui/(1+ui);
%*** Calculate yhati;
if yvar=1 then do;
yhati=sens*phati/(sens*phati+(1-spec)*(1-phati));
end;
if yvar=0 then do;
yhati=(1-sens)*phati/((1-sens)*phati+spec*(1-phati));
end;
%*** add to inf;
inf=inf+(phati#(1-phati)-yhati#(1-yhati))#xvar`*xvar;
%*** add to loglike;
loglike=loglike+yvar#log(sens#phati+(1-spec)#(1-phati));
loglike=loglike+(1-yvar)#log((1-sens)#phati+spec#(1-phati));
END;
asympvar=INV(inf);
FREE i inf yhati xhati xvar yvar;
FINISH;
RUN;
%******************* OUTPUT routine ********************************
*** Printing results and ***
*** optionally writing output data set ***
*********************************************************************;
START;
sebeta=SQRT(VECDIAG(asympvar));
zbeta=beta/sebeta;
p_value=2#(1-probnorm(abs(zbeta)));


PRINT "Variance/covariance of estimates : ",,
asympvar [ROWNAME=labelx COLNAME=labelx],,;
or=exp(beta);
or_l=exp(beta-1.96#sebeta);
or_u=exp(beta+1.96#sebeta);

cname1={"Estimate"};
cname2={"SE"};
cname3={"z"};
cname4={"pvalue"};
cname5={"lower"};
cname6={"upper"};
cname7={"OR"};
PRINT
"Estimate with s.e., z-score, p-value, and confidence limits:",,
labelx beta [COLNAME=cname1]
sebeta [COLNAME=cname2 FORMAT=8.4]
zbeta [COLNAME=cname3 FORMAT=8.4]
p_value [COLNAME=cname4 FORMAT=8.4]
or [COLNAME=cname7 FORMAT=8.2]
or_l [COLNAME=cname5 FORMAT=8.2]
or_u [COLNAME=cname6 FORMAT=8.2],;
;
PRINT "Log-likelihood=" loglike;
IF ALL(^war) THEN PRINT "*** WARNING: No intercept term in the model!";



FINISH; RUN;
***_________________end of the program _____________________________***;
QUIT;
%MEND;