/* =================================================================================== Macro BLINPLUS - measurement error correction for SAS version 8. *** ALL DATASETS MUST BE SAS DATASETS - Put more info on macro here *** Last modified on 9/29/04 by strui, fix phreg problem ===================================================================================*/ %macro blinplus8( type = LOGISTIC, /* type can be LOGISTIC, REG, or PHREG */ valid =, /* Dataset with Validation Study data */ main_est =, /* Dataset containing est.s from main study regress. */ /* on the var.s given in err_var (below). */ weights =, /* Dataset with weights */ err_var =, /* List of the variables measured with error */ true_var =, /* List of the variables measured without error */ depend =, /* Name of the dependent variable */ change =, /* Percent Change Scale: T/F (need type=REG for T) */ internal = FALSE, /* If validation study is Internal, do you wish to */ /* combine estimates? T/F (if T, need to provide */ /* specify val_est argument */ val_est = /* Dataset containing Internal validation Est.s */ ); /* __________________________________________________________________________________ */ /* Check some arguments */; %if %numargs(&err_var) ^= %numargs(&true_var) %then %goto out1; %let type = %upcase(&type); %if &type = REG %then %let type = LS; %if &type = PHREG %then %let type = COX; %if &type ^= LOGISTIC and &type ^= LS and &type ^= COX %then %goto out2; %let change = %upcase(&change); %if &change = TRUE %then %let change = T; %if &type ^= LS & &change = T %then %do; %put; %put *** WARNING FROM MACRO BLINPLUS:; %put *** Cannot set CHANGE=T unless TYPE=REG (or LS).; %put *** Continuing as if CHANGE = F; %put; %let change= %str(); %end; %let internal = %upcase(&internal); %if &internal = TRUE %then %let internal = T; %if &internal = T and &val_est = %str() %then %goto out3; proc iml workspace=999; reset center noname; file print; *** Read in variable names and their "weights"; use &weights; read all var _char_ into NAM; read all var _num_ into WT; *** Read in validation data; use &valid; read all var { &err_var } into X; read all var { &true_var } into Y; close &weights &valid; n=nrow(X); p=ncol(X); *** Find the location of the variables that are measured with error; lev = loc({&true_var}^={&err_var}); /*vector of locations of the error variables*/ *** If robust=T, then read parameter estimates (and their cov matrix) from the output *** of %robust(), else read from data set of main study regression estimates; %let m_est=&main_est; use &m_est(where=(_name_^="Intercept")); read all var {&err_var} into BG; B = t(BG[1,]); /* px1 vector of main study regression param est.s */ VB = BG[2:p+1,]; /* VB is the Cov matrix of B; dim pxp */ seb= sqrt(vecdiag(VB)); /* seb[i] is the std error of B[i]; seb has dim px1 */ close; free BG; *** If internal option specified then read in parameter estimates and cov obtained *** from the primary regression on the true covariates in the internal validation data; %if &internal = T %then %do; use &val_est(where=(_name_^="Intercept")); read all var { &true_var } into Valid; BV = t(Valid[1,]); /* Est of regr. coeffs: from intern. validation study*/ VBV = Valid[2:p+1,]; /* Covariance matrix of vector BV; dim: pxp */ sebv = sqrt(vecdiag(VBV)); /* sebv[i] is the std. error of BV[i]; dim: px1 */ free Valid; close &val_est; %end; /* ____________________________________________ */ ***add column of ones to the matrix of covariates (for model intercept term); X=J(n,1,1) || X; *** Weighted Multivariate Regression (validation study regression); F=inv(t(X)*X); /* used to shorten following expressions */ GWI=F*t(X)*Y; /* Est. of parameter matrix, incl. intercept; dim: (p+1)xp */ GEV=t(GWI[,lev]); /* Intrcpt and coeff.s of error var.s in the validation regr.*/ ERR=(Y - X*GWI); /* Error (residuals), equals: Y - Yhat; dim nxp */ *** calc. regular covariance matrix for validation regression coeffiecients; S=(t(ERR)*ERR)/(n-p-1); /* S is the MSE; dim pxp */ VEV = S[lev,lev]; /* Variance matrix of the response error variables */ %if &type=COX %then %do; VG=S@(F[1:p,1:p]); /* VG is est. of Var(G); dim: p^2*p^2. It's a pxp matrix */ IGT=inv(t(GWI[1:(p),])); /* dim: pxp; Inverse G Transpose; G = GWI[1:p,] */ /* G is the matrix of param est.s (minus intercept) */ %end; %else %do; VG=S@(F[2:p+1,2:p+1]); /* VG is est. of Var(G); dim: p^2*p^2. It's a pxp matrix */ IGT=inv(t(GWI[2:(p+1),])); /* dim: pxp; Inverse G Transpose; G = GWI[2:(p+1),] */ /* G is the matrix of param est.s (minus intercept) */ %end; free S; /* of pxp matrices; (i,j)th matrix is cov(G[,i],G[,j]).*/ /* Cov[G[k,i],G[l,j]]=VG[(i-1)*p+k,(j-1)*p+l] */ free X Y F ERR; /* ____________________________________________ */ *** Create matrix (k) to rearrange elements of (IGT@t(IGT)) into the order needed; k = design(t( (1:p)@J(1,p,1) + repeat((0:p-1)#p,1,p) )); m_k = k*(IGT@t(IGT)); *** Calculate corrected betas and covariance matrix; BRC=t(IGT)*B; /* Corrected primary regression estimates; dim: px1 */ VBRC=t(IGT)*VB*IGT + (I(p)@t(B))*m_k*VG*t(m_k)*t(I(p)@t(B)); /* variance est. for BRC; dim: pxp */ sebrc=sqrt(vecdiag(VBRC)); /* Std err.s based on VBRC */ free GWI k m_k VG; *** Combine internal validation study estimates with corrected main study estimates *** (weighted average). Weighting done with inverse covariance of est. (properly scaled); %if &internal = T %then %do; EST1 = BRC; VEST1= VBRC; VE = BLOCK(VEST1,VBV); /* variance of vector of both estimates (BRC,BV) */ XC = repeat(I(p),2,1); /* X matrix for combining estimates */ M = inv(t(XC)*inv(VE)*XC)*t(XC)*inv(VE); BRCI = M*(EST1//BV); /* Combined estimate */ VBRCI = M*VE*t(M); /* Variance of BRCI */ sebrci = sqrt(vecdiag(VBRCI)); /* se's of components of BRCI */ free EST1 VEST1 VE XC M VBRCI; %end; ***-----------------------------------------------------; *** BEGIN FORMATTING AND OUTPUT ***; ***-----------------------------------------------------; %if &type=LOGISTIC or &type=COX %then %do; %let col_nam=" WT B SE OR 95% CI p"; %let bfunc=exp; %end; %else %do; /* if &type=LS */ %let bfunc=; %if &change ^= T %then %let col_nam=" WT B SE WT*B 95% CI p"; %else %let col_nam=" WT B SE % Change 95% CI p"// " in WT*B"; %end; %if &type = LOGISTIC %then %let head1 = Logistic Regression:; %else %if &type = COX %then %let head1 = Cox Regression:; %else %if &type = LS %then %let head1 = Linear Regression:; print "MEASUREMENT ERROR CORRECTION OF REGRESSION ESTIMATES VIA THE METHOD OF REGRESSION CALIBRATION", " References: Rosner BA, Spiegelman D, Willett WC. Amer J Epi 1990. 132:734-745. ", " Spiegelman D, Kipnis V, Carroll R. Stat in Med 2001. 20:139-160. ", "Programmers: Donna Spiegelman, Sc.D., Harvard School of Public Health, Boston MA", " Roger Logan, Ph.D., Harvard School of Public Health, Boston MA ", " Doug Grove, M.S., Aidan McDermott, M.A., Carrie Wager, B.S., ", " of The Channing Laboratory, Boston MA. ",; print "TYPE OF REGRESSION: &type"; %if &depend^= %then %do; print "Dependent Variable: &depend"; %end; *print "OPTIONS USED:"; options nocenter; print ,, "Validation Study Regression Coeffiecients:",, %if &type ^=COX %then %do; GEV[rowname=(NAM[lev]) colname=({intercpt}//NAM) format=best8.]; %end; %else %do; GEV[rowname=(NAM[lev]) colname=(NAM) format=best8.]; %end; print, "Variance matrix of error-variables:", (VEV)[rowname=(NAM[lev]) colname=(NAM[lev]) format=best8.],; titlu="Main study regression coeffiecients: Uncorrected"; %pr_table(est=b, se_est=seb, tab_tit=titlu, row_lab=NAM); titlc="Main study regression coefficients: Corrected"; %pr_table(est=BRC, se_est=sebrc, tab_tit=titlc, row_lab=NAM); %if &internal = T %then %do; titl_iv="Internal validation regression"; %pr_table(est=bv, se_est=sebv, tab_tit=titl_iv, row_lab = NAM); titl_ip="Pooled estimate: Internal"; %pr_table(est=brci, se_est=sebrci, tab_tit=titl_ip, row_lab = NAM); %end; %goto out; *** Error Messages ***; %out1: %err_msg(1); %goto out; %out2: %err_msg(2); %goto out; %out3: %err_msg(3); %goto out; %out4: %err_msg(4); %goto out; %out5: %err_msg(5); %goto out; %out6: %err_msg(6); %goto out; %out:; quit; run; %mend blinplus8; %macro err_msg(err_num); %put; %put; %put ********************** ERROR IN SAS MACRO BLINPLUS: **********************; %put ****** ******; %if &err_num = 1 %then %do; %put ****** The number of variables given in the ERR_VAR argument must ******; %put ****** be equal to the number given in the TRUE_VAR argument. ******; %end; %else %if &err_num = 2 %then %do; %put ****** Invalid option chosen for argument TYPE ******; %put ****** Possible values are: LS (or REG), COX (or PHREG), LOGISTIC ******; %end; %else %if &err_num = 3 %then %do; %put ****** Invalid argument combination: INTERNAL and VAL_EST ******; %put ****** When Internal=TRUE, a value for VAL_EST must be provided ******; %end; %else %if &err_num = 4 %then %do; %put ****** Invalid argument specification: When HETER_ME=TRUE the ******; %put ****** arg.s HET_REG, HET_VAR and HET_FUNC must also be provided ******; %end; %else %if &err_num = 5 %then %do; %put ****** Invalid argument specification. When ROBUST=TRUE ******; %put ****** the arguments PREDS and PRED must also be provided. ******; %end; %else %if &err_num = 6 %then %do; %put ****** Invalid argument combination: HETER_ME=T, HETER_CV=T ******; %put ****** Cannot run with both these options specified ******; %end; %put ****** ******; %put ******************************************************************************; %put; %put; %mend err_msg; %macro pr_table( est=, se_est=, tab_tit=, row_lab=, wt=wt, func = &bfunc, chng = &change, col_lab=&col_nam ); oddr= &func(&wt#&est); ub = &func(&wt#(&est + 1.96*&se_est)); lb = &func(&wt#(&est - 1.96*&se_est)); %if &chng=T %then %do; oddr= 100#(exp(oddr) - 1); ub = 100#(exp(ub) - 1); lb = 100#(exp(lb) - 1); %end; pval= 1 - probchi((&est/&se_est)#(&est/&se_est),1); obj1 = char(&wt,10,2)||char(&est||&se_est||oddr||lb,9,5); obj2 = left(concat("- ",right(char(ub,8,5))))||char(pval,9,5); out_obj= obj1||obj2; divider=" ________ _______________________________________________________________"; print , &tab_tit, divider,&col_lab, divider, out_obj[rowname=&row_lab]; print divider,,; %mend pr_table; %macro numargs(arg); %let n=1; %do %until (%scan(&arg,%eval(&n),%str( ))=%str()); %let n=%eval(&n+1); %end; %eval(&n-1) %mend numargs; %macro robust( cov = , /* Variance Covariance dataset */ preds = , /* Predicted Variable dataset */ xvar = , /* x variable names */ pred = , /* Predicted Variable name */ dep = , /* Dependent Variable name */ out = , /* Output dataset */ ); data _temp_; set &cov(where=(_name_="ESTIMATE")); data cov; set &cov(where=(_name_^="ESTIMATE")); run; proc iml; use cov; read all var{intercep &xvar } into cov; use &preds; read all var{ &xvar } into x; x = j(nrow(x),1,1) || x; read all var { &pred } into p; read all var { &dep } into d; rsq = (d - p)#(d - p); gg = j(ncol(x),ncol(x),0); do i = 1 to nrow(x); tmp = rsq[i] # t(x[i,])*(x[i,]); gg = gg + tmp; end; rob = cov*gg*cov; create &out from rob[colname={INTERCEP &xvar }]; append from rob; quit; run; data &out; merge cov &out; data &out; set _temp_ &out; run; %mend;