* Macro for the unified approach as presented in the JBS special issue on * agreement assessment titled 'Unified approach for Assessing Agreement'. * This macro requires SAS IML. No plotting will be produced. For binary * plotting, please use the macro specified in our JASA paper. Here, * target values are considered random only. * Written by Wenting Wu on 5/30/07, updated/validated by Lawrence Lin on * 6/14/07. updated/validated by Lawrence Lin on 5/28/08. * Data structure must be one observation per subject ID. The program will compute # of subjects. * Method/assay/rater and replicate must be arranged as m1_1-m1_r m2_1-m2_r, ..., mk_1-mk_r, etc., * Data set must be balanced. If not, delete the unbalanced data. * One might get error messages if the data are near perfect resulting * into near zero stadard errors; * Parameters are:; ********************************************************************* * dataset: specify the dataset name * var: dependent vaiable names to be evaluated. eg., m1_1-m1_r m2_1-m2_r, ..., mk_1-mk_r, etc. * k: the number of methods/raters/instruments/assay, etc * m: the number of replications for each k * CCC_a_: CCC allowance * CCC_a_intra: intra CCC allowance * CCC_a_inter: inter CCC allowance * CCC_a_total: total CCC allowance * CP_a: coverage probability (CP) allowance * TDI_a_: TDI allowance when m=1 * TDI_a_intra: intra TDI allowance * TDI_a_inter: inter TDI allowance * TDI_a_total: total TDI allowance * tran = 1 transformation used for inference for continuous data = 0 No transformation used for inference for categorical data * Error = "const" for constant error structure. Here, TDI is the absolute difference. = "prop" for proportional error structure. Here, TDI is the % change. Log transformation to data will be applied. * The above TDI_a_ values must be in % when Error="prop" is specified. * Dec: significant digits after the decimal point printed for TDI when the "const" error is specified (eg. dec=2). * Alpha: (1-alpha)% one-tailed confidence limit **********************************************************************; options mprint mlogic; %macro UnifiedAgreement(dataset=, var=, k=, m=, CCC_a_=, CCC_a_intra=, CCC_a_inter=, CCC_a_total=, CP_a=., TDI_a_=., TDI_a_intra=., TDI_a_inter=., TDI_a_total=., tran=, error="prop", dec=2, alpha=); **** Set ods escape character to insert simple formatting text into ODS output ****; ods escapechar='^'; %***** Create symbol for turning on subscript. ******; %LET SUBS=%str(^\sub); data a; set &dataset; id=_n_; conf=(1-&alpha)*100; call symputx('conf', conf); run; * Converts data with one observation per row containing method values of 1 to &k, ID values of 1 to n, and rep values of 1 to &m; data one; set a; array x{&k, &m} &var; do method=1 to &k; do rep=1 to &m; y=x(method, rep); output; end; end; run; data one; set one; * For proportional error; if &tran=1 and &error="prop" then y=log(y); * Construct k-1 dummy variables for k methods; %Do i=1 %to &k-1; k1=%eval(&k-1); call symputx('k1', k1); x&i=0; if method=&i then x&i=1; %end; * Create method by rep index; ind=(method-1)*&m+rep; *proc print data=one; run; data y ( keep = y ); set one; run; data x ( keep = x1-x&k1); set one; run; data id ( keep = id ); set one; run; data rating ( keep = ind ); set one; run; run; %if &m>1 %then %do; proc sort data=one; by id method rep; run; proc means noprint; var y; by id method; output out=outb stddev=stddev; run; data outb; set outb; cellvariance=stddev**2; run; proc means noprint; var cellvariance; output out=outc mean=meancellvar; run; data initialerror; set outc; keep meancellvar; run; data method (keep = method); set one; run; data rep (keep = rep); set one; run; data cellvar (keep = meancellvar); set initialerror; run; %end; proc iml worksize=5000; reset nolog noprint; USE Y; READ ALL INTO Y; USE ID; READ ALL INTO ID; USE X; READ ALL INTO X; nx = ncol(x); /* nx = # of marginal covariates */ initial = y || x; create best from initial; append from initial; close best; %macro colname(name,number); %do r=2 %to &number; &name&r %end; %mend colname; quit; /* uses PROC REG to get initial estimates for BETA */ proc reg data=best outest=esti noprint; model col1 = %colname(col,&k); run; data par(drop=_type_ _model_ _depvar_ _rmse_ col1); set esti; if (_type_ ne 'PARMS' ) then delete; run; %if &m=1 %then %do; proc iml worksize=5000; reset nolog noprint; USE PAR; READ ALL INTO BETA; beta = beta`; nbeta = nrow(beta); USE ID; READ ALL INTO ID; n = max(id); * total number of subjects; call symputx('obs',n); USE Y; READ ALL INTO Y; USE X; READ ALL INTO X; USE RATING; READ ALL INTO RATING; k = &k; m = &m; nc = k*m; nsigma = k + 3; sigmaalpha=j(1,1,0); sigmaerror=j(1,1,0); sigmabeta=j(1,1,0); U_i = j(1,1,0); sigma=j(nsigma,1,0); Nst = j(nc,nc,0); Do i=1 to n; * THIS LOOP CALCULATES Nst; free times T_i R_i; times = loc(id=i); T_i = ncol(times); R_i = RATING[times,]; Do p=1 to T_i; Do q=1 to T_i; s = R_i[p,1]; t = R_i[q,1]; Nst[s,t] = Nst[s,t] + 1; end; end; end; V = j(nc,nc,0); mu_i=j(nc,1,0); Do i=1 to n; * THIS LOOP CALCULATES V; free times T_i Y_i R_i X_i mu_i; times = loc(id=i); T_i = ncol(times); Y_i = Y[times,]; R_i = RATING[times,]; X_i = X[times,]; X_i = j(T_i,1,1) || X_i; mu_i = X_i*beta; Do p=1 to T_i; Do q=1 to T_i; s = R_i[p,1]; t = R_i[q,1]; V[s,t] = V[s,t] + (Y_i[p,1] - mu_i[p,1])*(Y_i[q,1] - mu_i[q,1]) / Nst[s,t]; end; end; end; * This ends the loop over each person for V; *print, V mu_i k; muk = j(k,1,0); *mean of method; Do s = 1 to k; muk[s,1] = mu_i[s,1]; end; *initial value for sigma; Do s=1 to k; sigma[s,1]=muk[s,1]; End; Do s = 1 to k-1; Do t = s+1 to k; sigmabeta= sigmabeta + (muk[s,1]-muk[t,1])**2/(k*(k-1)); end; end; Do s = 1 to k-1; Do t = s+1 to k; sigmaalpha = sigmaalpha + 2*V[s, t]/(k*(k-1)); end; end; A=j(1,1,0); Do s = 1 to k; A = A + V[s,s]/k; end; sigmaerror = A - sigmaalpha; sigma[k+1,1]=sigmabeta; sigma[k+2,1]=sigmaalpha; sigma[k+3,1]=sigmaerror; *print, sigma muk; crit1=1; Do it=1 to 35 while (crit1 > 0.00001); U2 = j( nsigma, 1, 0); FHF = j( nsigma, nsigma, 0); u2sq = fhf; Do i=1 to n; free times T_i Y_i R_i X_i mu_i YY_i H_i V_i; times = loc(id=i); T_i = ncol(times); Y_i = Y[times,]; R_i = RATING[times,]; X_i = X[times,]; X_i = j(T_i,1,1) || X_i; mu_i = X_i*beta; YY_i = j(nsigma,1,0); Do s = 1 to k; * YY_i[1,1] to YY_i[k,1]; YY_i[s,1] = YY_i[s,1] + Y_i[s,1];*the 1st and 2nd term in y; end; Do s = 1 to k-1; *YY_i[k+1,1] for sigmabeta+sigmaerror; do t = 2 to k; YY_i[k+1,1]= YY_i[k+1,1] + (Y_i[s,1]-Y_i[t,1])**2/(k*(k-1)); end; end; Do s = 1 to k; *YY_i[k+3, 1] for sigmaerror+sigmaalpha; YY_i[k+2,1] = YY_i[k+2,1] + (Y_i[s,1]- muk[s,1])**2/(k); end; Do s = 1 to k-1; *YY_i[k+2,1] for sigmaalpha; do t = s+1 to k; YY_i[k+3,1] = YY_i[k+3,1] + 2*((Y_i[s,1] - muk[s,1])*(Y_i[t,1] - muk[t,1]))/(k*(k-1)); end; end; delta_i = j(nsigma,1,0); H_i = j(nsigma,nsigma,0); F_i = j(nsigma,nsigma,0); *delta matrix; Do s=1 to k; delta_i[s,1] = muk[s,1]; end; delta_i[k+1,1] = sigmabeta+sigmaerror; delta_i[k+2,1] = sigmaalpha+sigmaerror; delta_i[k+3,1] = sigmaalpha; Do s=1 to k; H_i[s,s] = sigmaalpha+sigmaerror; end; H_i[k+1,k+1] = (4*sigmaerror**2+8*sigmabeta*sigmaerror)/(k*(k-1)); H_i[k+2,k+2] = (2*k*sigmaalpha**2+2*sigmaerror**2+4*sigmaalpha*sigmaerror)/k; H_i[k+3,k+3] = (2*k*(k-1)*sigmaalpha**2+2*k*sigmaerror*sigmaalpha+2*sigmaerror**2)/(k*(k-1)); *F; Do s = 1 to k; F_i[s,s]=1; end; F_i[k+1,k+1] = 1; F_i[k+1,k+3] = 1; F_i[k+2,k+2] = 1; F_i[k+2,k+3] = 1; F_i[k+3,k+2] = 1; *print, H_i F_i delta_i YY_i ; u2_i = j(nsigma,1,0); fhf_i = j(nsigma,1,0); u2_i = F_i`*inv(H_i)*(YY_i-delta_i); fhf_i = F_i`*inv(H_i)*F_i; *print, u2_i fhf_i; U2 = U2 + u2_i; FHF = FHF + fhf_i; u2sq = u2sq + u2_i*u2_i`; *print, U2 FHF u2sq; end; DELTA2 = solve( FHF, U2 ); sigma = sigma + DELTA2; crit1 = max(ABS(DELTA2)); end; bread = j(nsigma,nsigma,0); meat = bread; bread = FHF; meat = u2sq; var = inv(bread) * meat * inv(bread)`; var_sigmabeta = var[k+1,k+1]; var_sigmaalpha = var[k+2,k+2]; var_sigmaerror = var[k+3,k+3]; cov_betaalpha = var[k+1,k+2]; cov_betaerror = var[k+1,k+3]; cov_alphaerror = var[k+2,k+3]; *indexes: precision, accuracy, ccc, TDI, CP and their transformations; rho = j(1,1,0); rho = sigmaalpha/(sigmaalpha + sigmaerror); rho_z = j(1,1,0); rho_z = log((1+rho)/(1-rho))/2; *z-transformation; accuracy = j(1,1,0); accuracy = (sigmaalpha + sigmaerror)/(sigmaalpha + sigmabeta + sigmaerror); accuracy_l = j(1,1,0); accuracy_l = log(accuracy/(1-accuracy)); *logit-transformation; ccc = j(1,1,0); ccc = sigmaalpha/(sigmaalpha + sigmabeta + sigmaerror); ccc_z = j(1,1,0); ccc_z = log((1+ccc)/(1-ccc))/2; *z-transformation; * Not for categorical data; %if &tran=1 %then %do; esquare = (2*sigmaerror+2*sigmabeta); * MSD; TDI = j(1,1,0); TDI = probit(1-(1-&CP_a)/2)*sqrt(esquare); TDI_w = j(1,1,0); TDI_w = log(esquare); *ln(e^2); CP = j(1,1,0); kp0_=&TDI_a_; if &error="prop" then kp0_=log(&TDI_a_/100+1); CP = 1-2*(1-probnorm(kp0_/sqrt(esquare))); cp_l = j(1,1,0); CP_l = log(cp/(1-cp)); *logit-transformation; %end; *s.e. for each index; var_rho = j(1,1,0); var_rho = ((1-rho)**2*var_sigmaalpha+(rho)**2*var_sigmaerror-2*(1-rho)*rho*cov_alphaerror)/ (sigmaalpha + sigmaerror)**2; se_rho = sqrt(var_rho); rho_notran_lower=rho-probit(1-&alpha)*se_rho; *lower bound for rho without z-tranform; se_rho_z=se_rho/(1-rho**2); *z-transformation; rho_z_low=rho_z-probit(1-&alpha)*se_rho_z; *low bound for z; rho_tran_lower=(exp(2*rho_z_low)-1)/(exp(2*rho_z_low)+1); *lower bound for rho based on z-tranform; var_ccc = j(1,1,0); var_ccc = ((1-ccc)**2*var_sigmaalpha+ccc**2*(var_sigmabeta+var_sigmaerror+2*cov_betaerror) -(2*(1-ccc)*ccc*(cov_betaalpha+cov_alphaerror)))/(sigmaalpha + sigmabeta + sigmaerror)**2; se_ccc = sqrt(var_ccc); ccc_notran_lower=ccc-probit(1-&alpha)*se_ccc; *lower and upper confidence interval without transformation; se_ccc_z = se_ccc/(1-(ccc)**2); *z-transformation; ccc_z_low=ccc_z-probit(1-&alpha)*se_ccc_z; ccc_tran_lower=(exp(2*ccc_z_low)-1)/(exp(2*ccc_z_low)+1); *lower and upper confidence interval based on z-transformation; var_accuracy = j(1,1,0); var_accuracy = ((1-accuracy)**2*(var_sigmaalpha+var_sigmaerror+2*cov_alphaerror)+ (accuracy)**2*var_sigmabeta-2*(1-accuracy)*accuracy*(cov_betaalpha+cov_betaerror)) /(sigmaalpha + sigmabeta + sigmaerror)**2; se_accuracy = sqrt(var_accuracy); accuracy_notran_lower=accuracy-probit(1-&alpha)*se_accuracy; *lower and upper confidence interval without transformation; se_accuracy_l = se_accuracy/(accuracy*(1-accuracy)); accuracy_l_low = accuracy_l-probit(1-&alpha)*se_accuracy_l; *low bound for l; accuracy_tran_lower = exp(accuracy_l_low)/(1+exp(accuracy_l_low)); var_esqua = j(1,1,0); var_esqua = 4*(var_sigmabeta+var_sigmaerror+2*cov_betaerror); se_esqua = sqrt(var_esqua); %if &tran=1 %then %do; se_TDI_w= sqrt(var_esqua)/esquare; *log-transformation, for TDI only consider tranformation; TDI_w_up=TDI_w+probit(1-&alpha)*se_TDI_w; *up bound for w Problems?; TDI_e_up=sqrt(exp(TDI_w_up)); *up bound for e; TDI_tran_upper=probit(1-(1-&CP_a)/2)*TDI_e_up; *up bound for TDI; %end; esquare = j(1,1,0); esquare = (2*sigmaerror+2*sigmabeta); var_esquare=j(1,1,0); var_esquare = 4*(var_sigmabeta+var_sigmaerror+2*cov_betaerror); %if &tran=1 %then %do; var_cp=j(1,1,0); var_cp=exp(-(kp0_)**2/(esquare))*(1+(kp0_)**2/(esquare))**2*var_esquare/(8*3.1415926*(kp0_)**2*esquare); se_cp=j(1,1,0); se_cp=sqrt(var_cp); cp_notran_lower=cp-probit(1-&alpha)*se_cp; se_cp_l = se_cp/(cp*(1-cp)); *se based on logit-transformation; cp_l_low = cp_l-probit(1-&alpha)*se_cp_l; cp_tran_lower = exp(cp_l_low)/(1+exp(cp_l_low)); RBS=2*sigmabeta/esquare; %end; %if &tran=0 %then %do; estimates = ccc || rho || accuracy || ccc_notran_lower||rho_notran_lower|| accuracy_notran_lower; create output from estimates; append from estimates; %end; %if &tran=1 %then %do; if &error="prop" then do; TDI=100*(exp(TDI)-1); TDI_tran_upper=100*(exp(TDI_tran_upper)-1); end; estimates = ccc || rho || accuracy || TDI || cp || ccc_tran_lower||rho_tran_lower|| accuracy_tran_lower|| TDI_tran_upper || cp_tran_lower || RBS; create output from estimates; append from estimates; %end; *proc print data=output; run; data r1; set output; S='Estimate '; CCC=Col1; Precision=Col2; Accuracy=Col3; TDI=Col4; CP=Col5; RBS=col11; output; S="&conf.% Conf. Limit"; CCC=Col6; Precision=Col7; Accuracy=Col8; TDI=Col9; CP=Col10; RBS=.; output; S='Allowance'; CCC=&CCC_a_; Precision=.; Accuracy=.; TDI=&TDI_a_; CP=&CP_a; RBS=.; output; keep s ccc Precision Accuracy TDI CP RBS; run; proc report data=r1 headline headskip split='|' nowindows; %if &tran=1 %then %do; column s CCC Precision Accuracy TDI CP RBS; %end; %else %do; column s CCC Precision Accuracy; %end; define s/ 'Statistics' format=$20. left; define CCC/ 'CCC' format=10.4 right; define Precision / 'Precision' format=10.4 right; define Accuracy / 'Accuracy' format=10.4 right; %if &tran=1 %then %do; %if &error="const" %then %do; define TDI / "TDI &subs &CP_a" format=10.&dec right; define CP / "CP &subs &TDI_a_" format=10.4 right; %end; %else %do; define TDI / "TDI% &subs &CP_a" format=10.2 right; define CP / "CP &subs &TDI_a_.%" format=10.4 right; %end; define RBS / "RBS*" format=10.2 right; compute after _page_; line "For k=&k, n=&obs, and m=&m"; line "*: The relative bias squared (RBS) must be less than 1 or 8 for CP_a of 0.9 or 0.8,"; line "respectively, in order for the approximated TDI and CP values to be valid."; endcomp; %end; %else %do; compute after _page_; line "For k=&k, n=&obs, and m=&m"; endcomp; %end; run; quit; ods rtf close; %end; %if &m>1 %then %do; proc iml worksize=5000; reset nolog noprint; USE PAR; READ ALL INTO BETA; beta = beta`; nbeta = nrow(beta); USE ID; READ ALL INTO ID; n = max(id); * total number of persons (clusters); call symputx('obs',n); USE Y; READ ALL INTO Y; USE X; READ ALL INTO X; USE RATING; READ ALL INTO RATING; USE METHOD; READ ALL INTO METHOD; USE CELLVAR; READ ALL INTO CELLVAR; USE rep; READ ALL INTO rep; k = &k; m = &m; nc = k*m; sigmaalpha=j(1,1,0); sigmaerror=j(1,1,0); sigmaerror=cellvar; sigmabeta=j(1,1,0); sigmagamma=j(1,1,0); U_i = j(1,1,0); nsigma = k + 4; sigma=j(nsigma,1,0); Nst = j(nc,nc,0); Do i=1 to n; * THIS LOOP CALCULATES Nst; free times T_i R_i; times = loc(id=i); T_i = ncol(times); R_i = RATING[times,]; Do j=1 to T_i; Do l=1 to T_i; s = R_i[j,1]; t = R_i[l,1]; Nst[s,t] = Nst[s,t] + 1; end; end; end; *print, Nst; Ytotal = j(nc,1,0); Do i=1 to n; free times Y_i T_i R_i; times = loc(id=i); T_i = ncol(times); R_i = RATING[times,]; Y_i = Y[times,]; Ytotal = Ytotal + Y_i; end; *print, Ytotal; Ybar = j(nc,1,0); Do s=1 to nc; Ybar[s,1] = Ytotal[s,1]/Nst[s,s]; end; *print, Ybar ; *the covariance structure for 1,2,... k*m columns* *each combination of (j,l) will be treated as a column in order to calculate the sigmas; VC = j(nc,nc,0); *variance of column; Do i = 1 to n; free times Y_i T_i R_i; times = loc(id=i); T_i = ncol(times); R_i = RATING[times,]; Y_i = Y[times,]; VC = VC + (Y_i - Ybar)*(Y_i - Ybar)`; end; Do s=1 to nc; Do t=1 to nc; VC[s,t] = VC[s,t]/ Nst[s,t]; end; end; * VM is the covariance structure for methods 1,2, ...K; * within each method, all obs. from different replications will be pooled together; * muk is the vector of methods mean; V = j(nc,nc,0); Do i=1 to n; free times T_i Y_i R_i; times = loc(id=i); T_i = ncol(times); Y_i = Y[times,]; R_i = rep[times,]; X_i = X[times,]; X_i = j(T_i,1,1) || X_i ; mu_i = X_i*beta; V = V + (Y_i - mu_i) * (Y_i - mu_i)` ; end; muk = j(k,1,0); *mean of method; VM = j(k,k,0); NM = j(k,k,0); *variance of method; Do s=1 to k; do l=1 to m; VM[s,s] = VM[s,s] + V[(s-1)*m+l,(s-1)*m+l]; NM[s,s] = NM[s,s] + Nst[(s-1)*m+l,(s-1)*m+l]; muk[s,1] = mu_i[(s-1)*m+1,1]; end; VM[s,s] = VM[s,s]/NM[s,s]; end; *initial value for sigma; do s=1 to k; sigma[s,1]=muk[s,1]; end; do s = 1 to k-1; do t = s+1 to k; sigmabeta= sigmabeta + (muk[s,1]-muk[t,1])**2/(k*(k-1)); end; end; do s = 1 to k-1; do t = s+1 to k; do l = 1 to m; do r = 1 to m; sigmaalpha = sigmaalpha + 2*VC[(s-1)*m+l, (t-1)*m+r]/(m*m*k*(k-1)); end; end; end; end; VC1 = j(1,1,0); do s = 1 to k; do l = 1 to m; VC1 = VC1 + VC[(s-1)*m+l,(s-1)*m+l]/(m*m*k); end; end; VC2 = j(1,1,0); do s = 1 to k; do l = 1 to m-1; do r = l+1 to m; VC2 = VC2 + 2*VC[(s-1)*m+l,(s-1)*m+r]/(m*m*k); end; end; end; sigmagamma = VC1 + VC2 - sigmaerror/m - sigmaalpha; sigma[k+1,1]=sigmabeta; sigma[k+2,1]=sigmaalpha; sigma[k+3,1]=sigmaerror; sigma[k+4,1]=sigmagamma; *print, sigmaalpha sigmabeta sigmaerror sigmagamma sigma muk; crit1=1; Do it=1 to 35 while (crit1 > 0.0000000001); U2 = j( nsigma, 1, 0); FHF = j( nsigma, nsigma, 0); u2sq = fhf; Do i=1 to n; free times T_i Y_i R_i X_i mu_i; times = loc(id=i); T_i = ncol(times); Y_i = Y[times,]; R_i = rep[times,]; YY_i = j(nsigma,1,0); Do s = 1 to k; * YY_i[1,1] to YY_i[k,1]; do l = 1 to m; YY_i[s,1] = YY_i[s,1] + Y_i[(s-1)*m+l,1];*the 1st and 2nd term in y; end; YY_i[s,1] = YY_i[s,1]/m; *check: l may not be equal to m for each i; end; Do s = 1 to k-1; *YY_i[k+1,1] for sigmabeta; do t = 2 to k; YY_i[k+1,1]= YY_i[k+1,1] + (YY_i[s,1]-YY_i[t,1])**2/(k*(k-1)); end; end; Do s = 1 to k-1; *YY_i[k+2,1] for sigmaalpha; do t = s+1 to k; do l = 1 to m; do r = 1 to m; YY_i[k+2,1] = YY_i[k+2,1] + 2*((Y_i[(s-1)*m+l,1] - Ybar[(s-1)*m+l,1])*(Y_i[(t-1)*m+r,1] - Ybar[(t-1)*m+r,1]))/(m*m*k*(k-1)); end; end; end; end; do s = 1 to k; *YY_i[k+3, 1] for sigmaerror; do l = 1 to m; YY_i[k+3,1] = YY_i[k+3,1] + (Y_i[(s-1)*m+l,1]- YY_i[s,1])**2/(k*(m-1)); end; end; *YY_i[k+4,1] for sigmagamma; A_i = j(1,1,0); do s = 1 to k; do l = 1 to m; A_i[1,1] = A_i[1,1] + (Y_i[(s-1)*m+l,1] - ybar[(s-1)*m+l,1])**2/(m*m*k); end; end; *print, A_i; B_i = j(1,1,0); do s = 1 to k; do l = 1 to m-1; do r = l+1 to m; B_i[1,1] = B_i[1,1] + 2*(Y_i[(s-1)*m+l,1] - ybar[(s-1)*m+l,1]) *(Y_i[(s-1)*m+r,1] - ybar[(s-1)*m+r,1])/(m*m*k); end; end; end; *print, B_i; YY_i[k+4,1] = YY_i[k+4,1] + A_i + B_i; * * end of YY_i[k+4,1]; delta_i = j(nsigma,1,0); H_i = j(nsigma,nsigma,0); F_i = j(nsigma,nsigma,0); *delta matrix; do s=1 to k; delta_i[s,1] = muk[s,1]; end; delta_i[k+1,1] = sigmabeta+sigmaerror/m+sigmagamma; delta_i[k+2,1] = sigmaalpha; delta_i[k+3,1] = sigmaerror; delta_i[k+4,1] = sigmagamma+sigmaalpha+sigmaerror/m; *end of delta; *H matrix; Do s=1 to k; H_i[s,s] = sigmagamma+sigmaalpha+sigmaerror/m; end; H_i[k+1, k+1] = (k*(k-1)*(3*k-2)/(m*m))*sigmaerror**2+k*(k-1)*(3*k-2)*sigmagamma**2 +(4*k*k*(k-1)/m)*sigmagamma*sigmaerror +(8*k*(k-1)/m)*sigmabeta*sigmaerror+8*k*(k-1)*sigmabeta*sigmagamma; H_i[k+2, k+2] = m**4*k*(k-1)*(2*k-3)*sigmaalpha**2+m**4*k*(k-1)*(2*k-3)/2*sigmagamma**2 +(k*(k-1)*m*m/2)*sigmaerror**2 + m**4*k*(k-1)*(2*k-3)*sigmaalpha*sigmagamma +((2+(m-1)**2+2*m*(m-1)*(k-2))*k*(k-1)*m**2/2)*sigmaalpha*sigmaerror +((2+(m-1)**2+2*m*(m-1)*(k-2))*k*(k-1)*m**2/2)*sigmagamma*sigmaerror; H_i[k+3, k+3] = (2*k/(m-1))*sigmaerror**2; H_i[k+4, k+4] = (2*m**2*k**2+k*m*(m-1)*(2*m-3+(k-1)*m*(m-1)/2)+2*(3*m-5)/m**2)*sigmaalpha**2 +(2*k*m**2+k*m*(m-1)*(2*m-3)+k*m**2*(m-1)*(k+2)/2-m**2*k**2)*sigmagamma**2 +(2*k*m+k*m*(m-1)/2)*sigmaerror**2 +(4*k*m**2+2*k*m*(m-1)*(2*m-3)+k*m*(m-1)*(m*k+m+2)-2*m**2*k**2)*sigmaalpha*sigmagamma +(4*k*m+k*m*(m-1)**2+m*k*(m-1)*(m*k+4)/2-m**2*k**2)*sigmagamma*sigmaerror +(4*k*m+k*m*(m-1)**2+m**2*k**2*(m-1)/2-m**2*k**2)*sigmaalpha*sigmaerror; *end of H; *F matrix; do s = 1 to k; F_i[s,s]=1; end; F_i[k+1,k+1] = 1; F_i[k+1,k+3] = 1/m; F_i[k+1,k+4] = 1; F_i[k+2,k+2] = 1; F_i[k+3,k+3] = 1; F_i[k+4,k+2] = 1; F_i[k+4,k+3] = 1/m; F_i[k+4,k+4] = 1; *end of F; *print, H_i F_i delta_i YY_i ; u2_i = j(k+4,1,0); fhf_i = j(k+4,1,0); u2_i = F_i`*inv(H_i)*(YY_i-delta_i); fhf_i = F_i`*inv(H_i)*F_i; U2 = U2 + u2_i; FHF = FHF + fhf_i; u2sq = u2sq + u2_i*u2_i`; *print, U2 FHF u2sq; end; DELTA2 = solve( FHF, U2 ); sigma = sigma + DELTA2; crit1 = max(ABS(DELTA2)); *print, H_i fhf u2sq sigma; end; crit = max(crit1); *print, it crit; bread = j(nsigma,nsigma,0); meat = bread; bread = FHF; meat = u2sq; var = inv(bread) * meat * inv(bread)`; var_sigmabeta = var[k+1,k+1]; var_sigmaalpha = var[k+2,k+2]; var_sigmaerror = var[k+3,k+3]; var_sigmagamma = var[k+4,k+4]; cov_betaalpha = var[k+1,k+2]; cov_betaerror = var[k+1,k+3]; cov_betagamma = var[k+1,k+4]; cov_alphaerror = var[k+2,k+3]; cov_alphagamma = var[k+2,k+4]; cov_errorgamma = var[k+3,k+4]; sigmaadjust=sigmabeta||sigmaalpha||sigmaerror||sigmagamma; create adjust from sigmaadjust; append from sigmaadjust; close adjust; varcov=var_sigmabeta||var_sigmaalpha||var_sigmaerror||var_sigmagamma||cov_betaalpha|| cov_betaerror||cov_betagamma||cov_alphaerror||cov_alphagamma||cov_errorgamma; create varcov from varcov; append from varcov; close varcov; data adjust2; *change non-positive estimate to zero; set adjust; if col2<0 then col2=0.00000001; if col3<0 then col3=0.00000001; if col4<0 then col4=0.00000001; run; data sigmabeta (keep = col1); set adjust2; run; data sigmaalpha (keep = col2); set adjust2; run; data sigmaerror (keep = col3); set adjust2; run; data sigmagamma (keep = col4); set adjust2; run; data var_sigmabeta (keep = col1); set varcov; run; data var_sigmaalpha (keep = col2); set varcov; run; data var_sigmaerror(keep = col3); set varcov; run; data var_sigmagamma (keep = col4); set varcov; run; data cov_betaalpha (keep = col5); set varcov; run; data cov_betaerror (keep = col6); set varcov; run; data cov_betagamma (keep = col7); set varcov; run; data cov_alphaerror (keep = col8); set varcov; run; data cov_alphagamma (keep = col9); set varcov; run; data cov_errorgamma(keep = col10); set varcov; run; proc iml worksize=5000; reset nolog noprint; USE sigmabeta; READ ALL INTO sigmabeta; USE sigmaalpha; READ ALL INTO sigmaalpha; USE sigmagamma; READ ALL INTO sigmagamma; USE sigmaerror; READ ALL INTO sigmaerror; USE var_sigmabeta; READ ALL INTO var_sigmabeta; USE var_sigmaalpha ; READ ALL INTO var_sigmaalpha ; USE var_sigmagamma; READ ALL INTO var_sigmagamma; USE var_sigmaerror; READ ALL INTO var_sigmaerror; USE cov_betaalpha; READ ALL INTO cov_betaalpha; USE cov_betaerror; READ ALL INTO cov_betaerror; USE cov_betagamma; READ ALL INTO cov_betagamma; USE cov_alphaerror; READ ALL INTO cov_alphaerror; USE cov_alphagamma; READ ALL INTO cov_alphagamma; USE cov_errorgamma; READ ALL INTO cov_errorgamma; k = &k; m = &m; *print, sigmabeta sigmaalpha sigmagamma sigmaerror var_sigmabeta var_sigmaalpha var_sigmagamma var_sigmaerror ; *print, cov_betaalpha cov_betaerror cov_betagamma cov_alphaerror cov_alphagamma cov_errorgamma; *all indexes; *intra; rho_intra = j(1,1,0); rho_intra = (sigmaalpha + sigmagamma)/(sigmaalpha + sigmagamma + sigmaerror); rho_intra_z = j(1,1,0); rho_intra_z = log((1+rho_intra)/(1-rho_intra))/2; ccc_intra = j(1,1,0); ccc_intra = (sigmaalpha + sigmagamma)/(sigmaalpha + sigmagamma + sigmaerror); ccc_intra_z = j(1,1,0); ccc_intra_z = log((1+ccc_intra)/(1-ccc_intra))/2; esquare_intra = j(1,1,0); esquare_intra = 2*sigmaerror; %if &tran=1 %then %do; TDI_intra = j(1,1,0); TDI_intra = probit(1-(1-&CP_a)/2)*sqrt(esquare_intra); TDI_intra_w = j(1,1,0); TDI_intra_w = log(2*sigmaerror); kp0_intra=&TDI_a_intra; if &error="prop" then kp0_intra=log(&TDI_a_intra/100+1); CP_intra = j(1,1,0); CP_intra = 1-2*(1-probnorm(kp0_intra/sqrt(esquare_intra))); CP_intra_l=j(1,1,0); CP_intra_l = log(CP_intra/(1-CP_intra)); %end; *inter; rho_inter = j(1,1,0); rho_inter = sigmaalpha/(sigmaalpha + sigmagamma + sigmaerror/m); rho_inter_z = j(1,1,0); rho_inter_z = log((1+rho_inter)/(1-rho_inter))/2; ccc_inter_1 = j(1,1,0); ccc_inter_1 = sigmaalpha/(sigmaalpha + sigmagamma + sigmaerror/m + sigmabeta); ccc_inter_1_z = j(1,1,0); ccc_inter_1_z=log((1+ccc_inter_1)/(1-ccc_inter_1))/2; ccc_inter_2 = j(1,1,0); ccc_inter_2 = sigmaalpha/(sigmaalpha + sigmagamma + sigmabeta); ccc_inter_2_z = j(1,1,0); ccc_inter_2_z=log((1+ccc_inter_2)/(1-ccc_inter_2))/2; accuracy_inter = j(1,1,0); accuracy_inter = (sigmaalpha + sigmagamma + sigmaerror/m)/(sigmaalpha + sigmagamma + sigmaerror/m + sigmabeta); accuracy_inter_l = j(1,1,0); accuracy_inter_l = log(accuracy_inter/(1-accuracy_inter)); *logit-transformation; esquare_inter = j(1,1,0); esquare_inter = 2*(sigmaerror/m+sigmagamma+sigmabeta); %if &tran=1 %then %do; TDI_inter = j(1,1,0); TDI_inter = probit(1-(1-&CP_a)/2)*sqrt(esquare_inter); TDI_inter_w = j(1,1,0); TDI_inter_w = log(esquare_inter); *ln(e^2); kp0_inter=&TDI_a_inter; if &error="prop" then kp0_inter=log(&TDI_a_inter/100+1); CP_inter = j(1,1,0); CP_inter = 1-2*(1-probnorm(kp0_inter/sqrt(esquare_inter))); CP_inter_l=j(1,1,0); CP_inter_l = log(CP_inter/(1-CP_inter)); %end; *total; rho_total = j(1,1,0); rho_total = sigmaalpha/(sigmaalpha + sigmagamma + sigmaerror); rho_total_z = j(1,1,0); rho_total_z = log((1+rho_total)/(1-rho_total))/2; ccc_total = j(1,1,0); ccc_total = sigmaalpha/(sigmaalpha + sigmagamma + sigmaerror + sigmabeta); ccc_total_z = j(1,1,0); ccc_total_z=log((1+ccc_total)/(1-ccc_total))/2; accuracy_total = j(1,1,0); accuracy_total = (sigmaalpha + sigmagamma + sigmaerror)/(sigmaalpha + sigmagamma + sigmaerror + sigmabeta); accuracy_total_l = j(1,1,0); accuracy_total_l = log(accuracy_total/(1-accuracy_total)); *logit-transformation; esquare_total = j(1,1,0); esquare_total = 2*(sigmaerror+sigmagamma+sigmabeta); %if &tran=1 %then %do; TDI_total = j(1,1,0); TDI_total = probit(1-(1-&CP_a)/2)*sqrt(esquare_total); TDI_total_w = j(1,1,0); TDI_total_w = log(esquare_total); *ln(e^2); kp0_total=&TDI_a_total; if &error="prop" then kp0_total=log(&TDI_a_total/100+1); CP_total = j(1,1,0); CP_total = 1-2*(1-probnorm(kp0_total/sqrt(esquare_total))); CP_total_l=j(1,1,0); CP_total_l = log(CP_total/(1-CP_total)); %end; *s.e. of intra; *-ccc and rhoc; rhoc0 = j(1,1,0); rhoc0 = ccc_intra; var_ccc_intra = ((1-rhoc0)**2*(var_sigmaalpha+var_sigmagamma+2*cov_alphagamma)+(rhoc0)**2*var_sigmaerror -2*(1-rhoc0)*rhoc0*(cov_alphaerror+cov_errorgamma))/(sigmaalpha + sigmagamma + sigmaerror)**2; se_ccc_intra = sqrt(var_ccc_intra); *before z-transformation; ccc_intra_notran_lower = ccc_intra-probit(1-&alpha)*se_ccc_intra; se_ccc_intra_z = se_ccc_intra/(1-(ccc_intra)**2); *z-transformation; ccc_intra_z_low=ccc_intra_z-probit(1-&alpha)*se_ccc_intra_z; ccc_intra_tran_lower = (exp(2*ccc_intra_z_low)-1)/(exp(2*ccc_intra_z_low)+1); se_rho_intra = se_ccc_intra; se_rho_intra_z = se_ccc_intra_z; rho_intra_z_low=ccc_intra_z_low; rho_intra_notran_lower = ccc_intra_notran_lower; rho_intra_tran_lower=ccc_intra_tran_lower; var_esqua_intra = j(1,1,0); var_esqua_intra = 4*(var_sigmaerror); se_esqua_intra = sqrt(var_esqua_intra); %if &tran=1 %then %do; *-TDI; se_TDI_intra_w= sqrt(var_esqua_intra)/(2*sigmaerror); *log-transformation; TDI_intra_w_up=TDI_intra_w+probit(1-&alpha)*se_TDI_intra_w; *up bound for w; TDI_intra_e_up=sqrt(exp(TDI_intra_w_up)); *up bound for e; TDI_intra_tran_upper=probit(1-(1-&CP_a)/2)*TDI_intra_e_up; *up bound for TDI; *-CP; var_notran_cpintra = j(1,1,0); var_notran_cpintra = (1/(exp((kp0_intra)**2/esquare_intra))*(1/kp0_intra+kp0_intra/esquare_intra)**2*var_esqua_intra)/(8*3.1415926*esquare_intra); se_notran_cpintra = j(1,1,0); se_notran_cpintra = sqrt (var_notran_cpintra); var_tran_cpintra = var_notran_cpintra/((cp_intra*(1-cp_intra))**2); se_tran_cpintra = se_notran_cpintra/(cp_intra*(1-cp_intra)); *se based on logit-transformation; cp_intra_notran_lower = cp_intra - probit(1-&alpha)*se_notran_cpintra; cp_intra_l_low = cp_intra_l-probit(1-&alpha)*se_tran_cpintra; cp_intra_tran_lower = exp(cp_intra_l_low)/(1+exp(cp_intra_l_low)); %end; *s.e. of inter; *-ccc_inter_1: with sigmaerror/m; rhoc1 = j(1,1,0); rhoc1 = ccc_inter_1; var_ccc_inter_notran = ((1-rhoc1)**2*var_sigmaalpha+(rhoc1)**2*(var_sigmabeta+var_sigmagamma+var_sigmaerror/(m**2) +2*cov_betagamma+2*cov_betaerror/m+2*cov_errorgamma/m)-(2*(1-rhoc1)*rhoc1*(cov_betaalpha+cov_alphagamma +cov_alphaerror/m)))/(sigmaalpha + sigmabeta + sigmagamma + sigmaerror/m)**2; se_ccc_inter_notran = sqrt(var_ccc_inter_notran); ccc_inter_notran_lower = ccc_inter_1 -probit(1-&alpha)*se_ccc_inter_notran; se_ccc_inter_1_z = se_ccc_inter_notran/(1-(ccc_inter_1)**2); *z-transformation; ccc_inter_z_low=ccc_inter_1_z-probit(1-&alpha)*se_ccc_inter_1_z; ccc_inter_tran_lower=(exp(2*ccc_inter_z_low)-1)/(exp(2*ccc_inter_z_low)+1); *-rho_inter; rhoc2b = j(1,1,0); rhoc2b = rho_inter; *precision; var_rho_inter = ((1-rhoc2b)**2*var_sigmaalpha+(rhoc2b)**2*(var_sigmagamma +var_sigmaerror/(m**2)+2*cov_errorgamma/m)-(2*(1-rhoc2b)*rhoc2b*(cov_alphagamma +cov_alphaerror/m)))/(sigmaalpha + sigmagamma + sigmaerror/m)**2; se_rho_inter = sqrt(var_rho_inter); rho_inter_notran_lower = rho_inter - probit(1-&alpha)* se_rho_inter; se_rho_inter_z=se_rho_inter/(1-rho_inter**2); *z-transformation; rho_inter_z_low=rho_inter_z-probit(1-&alpha)*se_rho_inter_z; *low bound for z; rho_inter_tran_lower=(exp(2*rho_inter_z_low)-1)/(exp(2*rho_inter_z_low)+1); *low bound for rho; *-accuracy_inter; var_accuracy_inter = j(1,1,0); var_accuracy_inter = ((1-accuracy_inter)**2*(var_sigmaalpha+var_sigmagamma+var_sigmaerror/(m**2) +2*cov_alphaerror/m+2*cov_alphagamma+2*cov_errorgamma/m)+(accuracy_inter)**2*var_sigmabeta -2*(1-accuracy_inter)*accuracy_inter*(cov_betaalpha+cov_betaerror/m+cov_betagamma)) /(sigmaalpha + sigmabeta + sigmaerror/m+sigmagamma)**2; se_accuracy_inter = sqrt(var_accuracy_inter); accuracy_inter_notran_lower = accuracy_inter - probit(1-&alpha)*se_accuracy_inter; se_accuracy_inter_l = se_accuracy_inter/(accuracy_inter*(1-accuracy_inter)); accuracy_inter_l_low = accuracy_inter_l-probit(1-&alpha)*se_accuracy_inter_l; *low bound for l; accuracy_inter_tran_lower = exp(accuracy_inter_l_low)/(1+exp(accuracy_inter_l_low)); *low bound for accuracy; var_esqua_inter = j(1,1,0); var_esqua_inter = 4*(var_sigmabeta+var_sigmagamma+var_sigmaerror/(m**2) +2*cov_betagamma+2*cov_betaerror/m+2*cov_errorgamma/m); se_esqua_inter = sqrt(var_esqua_inter); %if &tran=1 %then %do; *-TDI_inter; se_TDI_inter_w= sqrt(var_esqua_inter)/(2*sigmaerror/m+2*sigmabeta+2*sigmagamma); *log-transformation; TDI_inter_w_low=TDI_inter_w-probit(1-&alpha)*se_TDI_inter_w; *low bound for w; TDI_inter_w_up=TDI_inter_w+probit(1-&alpha)*se_TDI_inter_w; *up bound for w; TDI_inter_e_low=sqrt(exp(TDI_inter_w_low)); TDI_inter_e_up=sqrt(exp(TDI_inter_w_up)); TDI_inter_tran_lower=probit(1-(1-&CP_a)/2)*TDI_inter_e_low; TDI_inter_tran_upper=probit(1-(1-&CP_a)/2)*TDI_inter_e_up; *-CP_inter; var_notran_cpinter = j(1,1,0); var_notran_cpinter = (1/(exp((kp0_inter)**2/esquare_inter))*(1/kp0_inter+kp0_inter/esquare_inter)**2*var_esqua_inter)/(8*3.1415926*esquare_inter); se_notran_cpinter = j(1,1,0); se_notran_cpinter = sqrt (var_notran_cpinter); cp_inter_notran_lower = cp_inter- probit(1-&alpha)*se_notran_cpinter; var_tran_cpinter = var_notran_cpinter/((cp_inter*(1-cp_inter))**2); se_tran_cpinter = se_notran_cpinter/(cp_inter*(1-cp_inter)); cp_inter_l_low = cp_inter_l-probit(1-&alpha)*se_tran_cpinter; cp_inter_tran_lower = exp(cp_inter_l_low)/(1+exp(cp_inter_l_low)); RBS_inter= sigmabeta/(sigmaerror/m+sigmagamma); %end; *s.e. of total; *-ccc_total; rhoc3 = j(1,1,0); rhoc3 = ccc_total; var_ccc_total = ((1-rhoc3)**2*var_sigmaalpha+(rhoc3)**2*(var_sigmabeta+var_sigmagamma+var_sigmaerror) -2*(1-rhoc3)*rhoc3*(cov_betaalpha+cov_alphagamma+cov_alphaerror)+2*(rhoc3)**2*(cov_betaerror+cov_betagamma +cov_errorgamma))/(sigmaalpha + sigmabeta + sigmagamma + sigmaerror)**2; se_ccc_total_notran = sqrt(var_ccc_total); ccc_total_notran_lower = ccc_total - probit(1-&alpha)*se_ccc_total_notran; se_ccc_total_z = se_ccc_total_notran/(1-(ccc_total)**2); *z-transformation; ccc_total_z_low=ccc_total_z-probit(1-&alpha)*se_ccc_total_z; ccc_total_tran_lower=(exp(2*ccc_total_z_low)-1)/(exp(2*ccc_total_z_low)+1); *-rho_total; rhoc3b = j(1,1,0); rhoc3b = rho_total; var_rho_total = ((1-rhoc3b)**2*var_sigmaalpha+(rhoc3b)**2*(var_sigmagamma +var_sigmaerror+2*cov_errorgamma)-(2*(1-rhoc3b)*rhoc3b*(cov_alphagamma +cov_alphaerror)))/(sigmaalpha + sigmagamma + sigmaerror)**2; se_rho_total = sqrt(var_rho_total); rho_total_notran_lower = rho_total - probit(1-&alpha)*se_rho_total; se_rho_total_z=se_rho_total/(1-rho_total**2); *z-transformation; rho_total_z_low=rho_total_z-probit(1-&alpha)*se_rho_total_z; *low bound for z; rho_total_tran_lower=(exp(2*rho_total_z_low)-1)/(exp(2*rho_total_z_low)+1); *low bound for rho; *-accuracy_total; var_accuracy_ = j(1,1,0); var_accuracy_total = ((1-accuracy_total)**2*(var_sigmaalpha+var_sigmagamma+var_sigmaerror +2*cov_alphaerror+2*cov_alphagamma+2*cov_errorgamma)+(accuracy_total)**2*var_sigmabeta -2*(1-accuracy_total)*accuracy_total*(cov_betaalpha+cov_betaerror+cov_betagamma)) /(sigmaalpha + sigmabeta + sigmaerror+sigmagamma)**2; se_accuracy_total = sqrt(var_accuracy_total); accuracy_total_notran_lower = accuracy_total - probit(1-&alpha)*se_accuracy_total; se_accuracy_total_l = se_accuracy_total/(accuracy_total*(1-accuracy_total)); accuracy_total_l_low = accuracy_total_l-probit(1-&alpha)*se_accuracy_total_l; *low bound for l; accuracy_total_tran_lower = exp(accuracy_total_l_low)/(1+exp(accuracy_total_l_low)); *low bound for accuracy; var_esqua_total = j(1,1,0); var_esqua_total = 4*(var_sigmabeta+var_sigmagamma+var_sigmaerror +2*cov_betagamma+2*cov_betaerror+2*cov_errorgamma); se_esqua_total = sqrt(var_esqua_total); %if &tran=1 %then %do; *-TDI_total; se_TDI_total_w= sqrt(var_esqua_total)/(2*sigmaerror+2*sigmabeta+2*sigmagamma); *log-transformation; TDI_total_w_up=TDI_total_w+probit(1-&alpha)*se_TDI_total_w; *up bound for w; TDI_total_e_up=sqrt(exp(TDI_total_w_up)); TDI_total_tran_upper=probit(1-(1-&CP_a)/2)*TDI_total_e_up; *-CP_total; var_notran_cptotal = j(1,1,0); var_notran_cptotal = (1/(exp((kp0_total)**2/esquare_total))*(1/kp0_total+kp0_total/esquare_total)**2*var_esqua_total)/(8*3.1415926*esquare_total); se_notran_cptotal = j(1,1,0); se_notran_cptotal = sqrt (var_notran_cptotal); var_tran_cptotal = var_notran_cptotal/((cp_total*(1-cp_total))**2); se_tran_cptotal = se_notran_cptotal/(cp_total*(1-cp_total)); cp_total_notran_lower = cp_total - probit(1-&alpha)*se_notran_cptotal; cp_total_l_low = cp_total_l-probit(1-&alpha)*se_tran_cptotal; cp_total_tran_lower = exp(cp_total_l_low)/(1+exp(cp_total_l_low)); RBS_total= sigmabeta/(sigmaerror+sigmagamma); %end; %if &tran=0 %then %do; estimates = ccc_intra || ccc_inter_1 || ccc_total || rho_intra || rho_inter || rho_total || accuracy_inter || accuracy_total || ccc_intra_notran_lower || ccc_inter_notran_lower || ccc_total_notran_lower || rho_intra_notran_lower || rho_inter_notran_lower || rho_total_notran_lower || accuracy_inter_notran_lower || accuracy_total_notran_lower ; create output from estimates; append from estimates; %end; %if &tran=1 %then %do; if &error="prop" then do; TDI_intra=100*(exp(TDI_intra)-1); TDI_inter=100*(exp(TDI_inter)-1); TDI_total=100*(exp(TDI_total)-1); TDI_intra_tran_upper=100*(exp(TDI_intra_tran_upper)-1); TDI_inter_tran_upper=100*(exp(TDI_inter_tran_upper)-1); TDI_total_tran_upper=100*(exp(TDI_total_tran_upper)-1); end; estimates = ccc_intra || ccc_inter_1 || ccc_total || rho_intra || rho_inter || rho_total ||accuracy_inter || accuracy_total || TDI_intra || TDI_inter || TDI_total || CP_intra || CP_inter || CP_total || ccc_intra_tran_lower || ccc_inter_tran_lower || ccc_total_tran_lower || rho_intra_tran_lower || rho_inter_tran_lower || rho_total_tran_lower || accuracy_inter_tran_lower || accuracy_total_tran_lower || TDI_intra_tran_upper || TDI_inter_tran_upper || TDI_total_tran_upper || cp_intra_tran_lower || cp_inter_tran_lower || cp_total_tran_lower || RBS_inter||RBS_total; create output from estimates; append from estimates; %end; *proc print data=output; run; data r; set output; %if &tran=1 %then %do; t='Intra'; S='Estimate '; CCC=Col1; Precision=Col4; Accuracy=.; TDI=Col9; CP=Col12; RBS=.; output; S="&conf.% Conf. Limit"; CCC=Col15; Precision=Col18; Accuracy=.; TDI=Col23; CP=Col26; RBS=.; output; S='Allowance'; CCC=&CCC_a_intra; Precision=&CCC_a_intra; Accuracy=.; TDI=&TDI_a_intra; CP=&CP_a; RBS=.; output; t='Inter'; S='Estimate '; CCC=Col2; Precision=Col5; Accuracy=Col7; TDI=Col10; CP=Col13; RBS=col29; output; S="&conf.% Conf. Limit"; CCC=Col16; Precision=Col19; Accuracy=Col21; TDI=Col24; CP=Col27; RBS=.; output; S='Allowance'; CCC=&CCC_a_Inter; Precision=.; Accuracy=.; TDI=&TDI_a_Inter; CP=&CP_a; RBS=.; output; t='Total'; S='Estimate '; CCC=Col3; Precision=Col6; Accuracy=Col8; TDI=Col11; CP=Col14; RBS=col30; output; S="&conf.% Conf. Limit"; CCC=Col17; Precision=Col20; Accuracy=Col22; TDI=Col25; CP=Col28; RBS=.; output; S='Allowance'; CCC=&CCC_a_Total; Precision=.; Accuracy=.; TDI=&TDI_a_Total; CP=&CP_a; RBS=.; output; keep t s ccc Precision Accuracy TDI CP RBS; %end; %else %do; t='Intra'; S='Estimate '; CCC=Col1; Precision=Col4; Accuracy=.; output; S="&conf.% Conf. Limit"; CCC=Col9; Precision=Col12; Accuracy=.; output; S='Allowance'; CCC=&CCC_a_intra; Precision=&CCC_a_intra; Accuracy=.; output; t='Inter'; S='Estimate '; CCC=Col2; Precision=Col5; Accuracy=Col7; output; S="&conf.% Conf. Limit"; CCC=Col10; Precision=Col13; Accuracy=Col15; output; S='Allowance'; CCC=&CCC_a_Inter; Precision=.; Accuracy=.; output; t='Total'; S='Estimate '; CCC=Col3; Precision=Col6; Accuracy=Col8; output; S="&conf.% Conf. Limit"; CCC=Col11; Precision=Col14; Accuracy=Col16; output; S='Allowance'; CCC=&CCC_a_Total; Precision=.; Accuracy=.; output; keep t s ccc Precision Accuracy; %end; run; proc report data=r headline headskip split='|' nowindows; %if &tran=1 %then %do; column t s CCC Precision Accuracy TDI CP RBS; %end; %else %do; column t s CCC Precision Accuracy; %end; define t/ group order=data 'Type' format=$5. left; define s/ 'Statistics' format=$20. left; define CCC/ 'CCC' format=10.4 right; define Precision / 'Precision' format=10.4 right; define Accuracy / 'Accuracy' format=10.4 right; %if &tran=1 %then %do; %if &error="const" %then %do; define TDI / "TDI &subs CP_a" format=10.&dec right; define CP / "CP &subs TDI_a" format=10.4 right; %end; %else %do; define TDI / "TDI% &subs CP_a" format=10.2 right; define CP / "CP &subs TDI_a%" format=10.4 right; %end; define RBS / "RBS*" format=10.2 right; compute after _page_; line "For k=&k, n=&obs, and m=&m"; line "*: The relative bias squared (RBS) must be less than 1 or 8 for CP_a of 0.9 or 0.8,"; line "respectively, in order for the approximated TDI and CP values to be valid."; endcomp; %end; %else %do; compute after _page_; line "For k=&k, n=&obs, and m=&m"; endcomp; %end; run; ods rtf close; quit; %end; ods rtf close; %mend UnifiedAgreement; *************** The followings are examples ************************* * of running the macro using the two examples in the paper. Again, * data structure must have a dependent variable (Y), subject (sample) ID, * method/assay/rater and replicate ID with values of 1, 2, 3, ..., etc., * with one observation per row.; * Give ODS output detination then run the macro; /* Eample 1 (continuous data): ods rtf file='x:\xxxx\xxxx\EX1.rtf' style=styles.TablesRTF NOTOC_DATA; %UnifiedAgreement(dataset=xxx, var=m1_1 m1_2 m2_1 m2_2, k=2, m=2, CCC_a_intra=0.9943, CCC_a_inter=0.9775, CCC_a_total=0.9775, CP_a=0.9, tran=1, TDI_a_intra=75, TDI_a_inter=150, TDI_a_total=150, error="const", dec=1, alpha=0.025); Eample 2 (categorical data): ods rtf file='x:\xxxx\xxxx\EX2.rtf' style=styles.TablesRTF NOTOC_DATA; %UnifiedAgreement(dataset=xxxx, var=m1_1 m1_2 m2_1 m2_2, k=2, m=2, CCC_a_intra=0.75, CCC_a_inter=0.4375, CCC_a_total=0.4375, tran=0, alpha=0.025); */