%macro icc9(data=, varlist=,  subject=, where=,   maxdec=4, exp=F,
            printsd=f,
	    method=ml, empcal=T, byvar=, notes=nonotes, outdat=, noprint=F, model=);
      /* to get reliability coefficients */
   
 /*********************
  by Ellen Hertzmark, Edmond Kabagambe, Eric Tchetgen, and Peter Gaccione
  **********************/
   
   
   
%local convstat i numvar  _fdl _nt ;
   
%let varlist=%upcase(&varlist);
%let subject=%upcase(&subject);
%let noprint=%upcase(&noprint);
%let empcal=%upcase(&empcal);
%let exp=%upcase(&exp);
%let printsd=%upcase(&printsd);
   
%let _fdl=%sysfunc(getoption(formdlim));
%let _nt=%sysfunc(getoption(notes));
   
   options &notes  nosyntaxcheck formdlim='-';
   
   
proc contents data=&data noprint  out=_d1_;  run;
data _d1_;  set _d1_ ;
   length varname $25  _lbl_ $40 ;
   varname=upcase(name);
   _lbl_=label;
   keep varname _lbl_;
   label _lbl_='label';
   run;
proc sort;  by varname;  run;
   
data _cucxxx_;  set &data;
%if %length(&where) ne 0 %then %do;  where &where;  %end;
   mervarbl=1;
   run;
proc sort data=_cucxxx_;  by &byvar mervarbl;  run;
   
%let numvar=%numargs(&varlist);
   
%do i=1 %to &numvar;
   %let thisvar = %scan(&varlist, &i, %str( ));
 /* data set with no missing values for thisvar */
data _cucuxx_;  set _cucxxx_;  where &thisvar gt .Z;  run;
   
proc freq noprint data=_cucuxx_;  tables &subject / out=_tabx_;
   by &byvar mervarbl;
   run;
data _tabx_;  set _tabx_;  where count gt 1;  run;
proc means noprint data=_tabx_;  var  count;
   output  out=_tab3_  n(count)=numw2 sum(count)=obsw2;
   by &byvar mervarbl;
   run;
data _tab3_;  set _tab3_;
   avgrep=round(obsw2/numw2, .1);
   label numw2='# subjects with >1 obs'
      obsw2='# obs these subjects have'
      avgrep='avg # meas/subj w >1 meas'
      ;
   
   call symput ('_numw2', trim(left(numw2)));
   call symput ('_obsw2', trim(left(obsw2)));
   call symput ('_avgrep', trim(left(avgrep)));
   run;
 /*
  proc print label data=_tab3_;  var numw2 obsw2;  run;
  */
   
   ods output
      convergencestatus=_cvstat_
      asycov=_asycov_
      Dimensions=_dims_
      NObs=_nobs_
      SolutionF=_solnf_
      CovParms=_cp_;
   ods listing close;
proc sort data=_cucuxx_;  by &subject;  run;
proc mixed data=_cucuxx_  method=&method asycov ; 
   model &thisvar= &model  / s  ddfm=bw ;
   random intercept / subject=&subject;
   by &byvar mervarbl;
   run;
   
data _cvstat_;  set _cvstat_;
   call symput ('convstat', status);
   run;
%if &convstat ne 0 %then %goto _bad;
   
data _yyy1_;  set _cp_;
   if upcase(covparm) eq "&subject" or upcase(covparm) eq 'CS' or upcase(covparm) eq 'INTERCEPT';
   sigb=estimate;  
   keep &byvar sigb mervarbl;
   run;
   
data _yyy2_;  set _cp_;  if upcase(covparm) eq 'RESIDUAL';
   sigw=estimate;
   keep &byvar sigw  mervarbl;
   run;
   
data _xxx1a_;  set _asycov_;
   covparm=upcase(covparm);
   if covparm eq "&subject" or covparm eq 'CS'  or covparm eq 'INTERCEPT';
   varsigb=covp1;  covbw=covp2;
   keep &byvar mervarbl varsigb covbw;
   run;
data _xxx1b_;  set _asycov_;
   covparm=upcase(covparm);
   if covparm eq 'RESIDUAL';
   varsigw=covp2;
   keep &byvar mervarbl varsigw;
   run;
   
data _xxx2a_;  set _nobs_;  if upcase(label) eq 'NUMBER OF OBSERVATIONS USED';
   nmeas=nobsused;
   keep &byvar nmeas mervarbl;
   run;
   
data _xxx2b_;  set _dims_;  if Descr eq 'Subjects';
   nsubj=value;
   keep &byvar nsubj mervarbl;
   run;
   
   /*V 6.12  name change from PARM to _EFFECT_*/
data _zzz_;  set _solnf_;  where upcase(effect) eq 'INTERCEPT';
   xbar=estimate;  sexbar=StdErr;  varxbar=StdErr*StdErr;
   keep &byvar xbar sexbar varxbar mervarbl;
   run;
   
data _aaa_;
   length varname $25;
   merge _yyy1_ _yyy2_ _xxx1a_ _xxx1b_ _xxx2a_ _xxx2b_ _zzz_  _tab3_ ;
   by &byvar mervarbl;
   varname = "&thisvar";
   totvar=sigb+sigw;
   rel=sigb/totvar;
   vlord=&i;
   if sigw gt 0 then sdw=sqrt(sigw);
   cvw=sdw/xbar;
   if . lt cvw lt 0 then do;  cvw=abs(cvw);  sign=-1;  end;
   else sign = 1;
   if totvar gt 0 then sdc=sqrt(totvar);
   
   
   /***************  peters code for variance of coeff of var  ******/
   varcov=varsigw/(4*xbar*xbar*sigw) + (sigw/(xbar**4))*varxbar;
   /*****************************************************************/
   
   /**  Make log transform for cv so that the cis are between 0 and 1 ****/
   lncvw=log(cvw);
   varlncvw=(1/cvw**2)*varcov;
   
   lnlclcvw=lncvw-1.96*sqrt(varlncvw);
   lnuclcvw=lncvw+1.96*sqrt(varlncvw);
   if sign eq -1 then do;  lx=lnuclcvw;  ux=lnlclcvw;
      lnlclcvw=lx;  lnuclcvw=ux;
      end;
   lclcvw=sign*exp(lnlclcvw);
   uclcvw=sign*exp(lnuclcvw);
   cvw=sign*cvw;
%makeestci(est=cvw, lcl=lclcvw, ucl=uclcvw, nameestci=estcicvw,
	   labelestci=%str(Estimated coefficient of within-subject variance (95% CI) ) );
   
   /**************peters code ****************/;
   varri=(sigb**2*varsigw+sigw**2*varsigb-2*sigb*sigw*covbw)/sigb**4 ;
   /*******************************************/
   
   logitrel=log(rel/(1-rel));
   vlogitre=((rel**2)*varri)/((1-rel)**2);
   selogitr=sqrt(vlogitre);
   lnlclrel=logitrel-1.96*selogitr;
   lnuclrel=logitrel+1.96*selogitr;
   lclrel=exp(lnlclrel)/(1+exp(lnlclrel));
   uclrel=exp(lnuclrel)/(1+exp(lnuclrel));
%makeestci(est=rel, lcl=lclrel, ucl=uclrel, nameestci=estcirel,
	   labelestci=%str(Estimated coefficient of reliability, ICC (95% CI) ) );
   
   drop sign lx ux;
   run;
   
%goto _good;
%_bad:
data _aaa_;
   length varname $25 ;
   rel=.;  avgrep=.;  cvw=.;  xbar=.;  sexbar=.;  varxbar=.;
   varname="&thisvar";  vlord=&i;
   lclcvw=.;  uclcvw=.;  lclrel=.;  uclrel=.;  uclrelo=.; lclrelo=.;
   lclcvwo=.; uclcvwo=.;
   nsubj=.;  nmeas=.;  numw2=.;
   sdc=.;  totvar=.;
   run;
   
%_good:
%if &i eq 1 %then %do;
   data _outpt_;  set _aaa_;  format _numeric_ 10.&maxdec;  run;
%end;
%else %do;
   data _outpt_;  set _outpt_ _aaa_;  run;
%end;
   
%end;  /* end of loop on vbls */



ods listing;


proc sort data=_outpt_;  by varname;  run;
data _outpt_;  merge _d1_ _outpt_ (in=inn);  by varname;  if inn;
run;
data _outpt_;  set _outpt_;
label nsubj='# subjects'  nmeas='# observations'  numw2='# subj with >1 meas'
  lclrel='lower 95% conf bound for reliability'
  uclrel='upper 95% conf bound for reliability'
  rel='reliability, ICC'
  lclcvw='lower 95% conf bound for w-in subj coeff of var'
  uclcvw='upper 95% conf bound for w-in subj coeff of var'
  cvw='within-subject coefficient of variation'
  vlord='vbl #'
  nsubj='# subjects'  nmeas='# observations'  xbar='Mean'
  sexbar='Standard Error'  sdc='Standard Deviation'
  vlord='order of vbl  '  sdw='Within-subj std dev'
  pctse='Std Err as percent of mean'  pctsd='Std Dev as percent of mean'
  emean='Mean on orig scale'  elcl='Lower 95% CI for mean'
  eucl='Upper 95% CI for mean'  elser='Lower end of SE range'
  euser='Upper end of SE range'
;
%if "&exp" eq "T" %then %do;
  emean=exp(xbar);  pctse=exp(sexbar);  pctsd=exp(sdc);
pctse=round(100*pctse, .1);  pctsd=round(100*pctsd, .1);
  elcl=exp(xbar-1.96*sexbar);  eucl=exp(xbar+1.96*sexbar);
  elser=exp(xbar-sexbar);  euser=exp(xbar+sexbar);
%makeestci(est=emean, lcl=elcl, ucl=eucl, nameestci=estciemean,
labelestci=%str(Exponentiated mean (95% SD range) ) );
%makeestci(est=emean, lcl=lser, ucl=user, nameestci=estciemse,
labelestci=%str(Exponentiated mean (95% CI) ) );
%end;
run;
proc sort  data=_outpt_;  by  vlord varname _lbl_ &byvar ;  run;

%if "&noprint" eq "F" %then %do;
%if %length(&where) ne 0 %then %do;
  title4 "Analysis run on subset WHERE &where..";
%end;
%if %length(&model) ne 0 %then %do;
  title5 "Adjusting for:  &model ";
%end;
title6 "Number of subjects with more than 1 measurement:  &_numw2";
title7 "These subjects have &_obsw2 measurements (&_avgrep on average)";
proc print label data=_outpt_;
var &byvar  nsubj nmeas estcirel estcicvw;
 by vlord varname _lbl_ ;
format nsubj nmeas  9.0;
run;

%if &printsd eq T %then %do;
%if &exp ne T %then %do;
proc print  label data=_outpt_;
format xbar sexbar sdc sdw 8.&maxdec ;
format nsubj nmeas numw2 9.0 ;
var &byvar  nsubj nmeas   xbar sexbar sdc sdw;
by vlord varname _lbl_;
run;
%end;

%if &exp eq T %then %do;
  proc print label data=_outpt_;
  format pctsd pctse 5.1 ;
  var &byvar nsubj nmeas estciemean estciemse pctsd pctse;
  by vloard varname _lbl_;
  run;
%end;
%end;  /* end of printsd eq t */
%end;  /* end of noprint eq f */

%if "&outdat" ne "" %then %do;
  data &outdat;  set _outpt_;  run;
%end;

proc datasets nolist;
delete _d1_ _cucuxx_  _xxx1a_ _xxx1b_ _xxx2a_ _xxx2b_ _cucxxx_
  _yyy2_ _solnf_ _cp_ _dims_  _asycov_ _cl_
_outpt_ _tab3_    _nobs_
_yyy1_ _aaa_
;

run;  quit;

options notes formdlim="&_fdl"  obs=max syntaxcheck ;
quit;
%mend ;


%*                                                                            ;
%* This macro was developed by Carrie Wager,Programmer,Channing Laboratory 1990;
%* Modified AMcD 1993, e hertzmark 1994 and L Chen 1996;
%*                                                                            ;

%macro numargs(arg, delimit);
   %if %quote(&arg)= %then %do;
        0
   %end;
   %else %do;
     %let n=1;
     %do %until (%qscan(%quote(&arg), %eval(&n), %str( ))=%str()); 
        %let n=%eval(&n+1);
        %end;
	%eval(&n-1)
   %end;
   %mend numargs;
%macro makeestci(est=, lcl=, ucl=, labelestci=%str(Estimate (95% CI)),
nameestci=estci);
length &nameestci $25;
length _vvv1 _vvv3 _vvv5 $6  _vvv2 _vvv6 $2 ;
  _vvv1=put(round(&est, .01), 5.2);  _vvv2=' (';
  _vvv3=put(round(&lcl, .01), 5.2);  _vvv4=', ';
  _vvv5=put(round(&ucl, .01), 5.2);  _vvv6=')';
   &nameestci=compress(left(_vvv1)) || _vvv2 || compress(left(_vvv3)) || _vvv4 || compress(left(_vvv5)) || _vvv6 ;
label &nameestci="&labelestci";
%mend;