*! version 2.0.2  6mar1998  TJS & NJC  STB-43 sg84
program define concord
 version 5.0

* Setup
 local if "opt"
 local in "opt"
 local weight "fweight opt"
 #delimit ;
 local options "Summary Graph(str) Level(real 95) BY(str) Pen(str)
   XScale(str) YScale(str) Connect(str) Symbol(str) XLAbel(string)
   YLAbel(string) L1title(str) T1title(str) T2title(str) SAving(str)
   SOrt *" ;
 #delimit cr
 local varlist "req ex min(2) max(2)"
 parse "`*'"
 parse "`varlist'", parse(" ,")

* Set temporary names
 tempvar d db dll dul m
 tempname dv
 tempvar touse byg rmaxis
 tempname zv k xb yb sx2 sy2 r rp sxy p u sep z zp
 tempname t set ll ul llt ult zt ztp sd1 sd2 sl

* Set up wgt and what to use
 if "`weight'" != "" { local wgt "[`weight'`exp']" }
 mark `touse' `if' `in' `wgt'
 markout `touse' `varlist'

* Generate CI z-value and label from Level()
 if `level' < 1 { local level = `level' * 100 }
 scalar `zv' = -1 * invnorm((1 - `level' / 100) / 2)
 local rl = `level'
 local level : di %7.0g `level'
 local level = ltrim("`level'")

* Generate BY groups
 if "`by'" != "" {
   confirm var `by'
   sort `by'
   qui by `by': gen byte `byg' = _n==1
   qui replace `byg' = sum(`byg')
   local byn = `byg'[_N]
 }
 else {
   qui gen byte `byg' = 1
   local byn = 1
 }

* Generate `by' labels -- if required
 if ("`by'" != "") {        
   tempvar kk km bylabel
   capture decode `by', gen(`bylabel')
   if _rc != 0 { 
     local type : type `by'
     qui gen `type' `bylabel' = `by'
   }
   qui by `by': gen `kk' = _n 
   qui egen `km' = max(_n), by(`by')
   qui replace `km' = . if `kk' != 1
 }

 if "`graph'" != "" {local graph = lower(substr("`graph'",1,1))}

* Do calculations
 local j = 1
 while `j' <= `byn' {  /* start of loop for each `by' group */
   di 
   if "`by'" != "" {
     sort `km'
     di in bl "--> `by' == " `bylabel'[`j']
     local byl : di "--> `by' == " `bylabel'[`j']
    }

* LOA (Bland & Altman) calculations
   qui gen `d' = `2' - `1'
   qui summ `d' if `touse' & `byg' == `j' `wgt'
   qui gen `db' = _result(3)
   scalar  `dv' = _result(4)
   scalar  `zv' = -1 * invnorm((1 - `level' / 100) / 2)
   qui gen `dll' = `db' - `zv' * sqrt(`dv')
   qui gen `dul' = `db' + `zv' * sqrt(`dv')
   qui gen `m' = (`1' + `2') / 2

* Concordance calculations
   qui summ `1' if `touse' & `byg' == `j' `wgt'
   scalar `k' = _result(1)
   scalar `yb' = _result(3)
   scalar `sy2' = _result(4) * (`k' - 1) / `k'
   scalar `sd1' = sqrt(_result(4))

   qui summ `2' if `touse' & `byg' == `j' `wgt'
   scalar `xb' = _result(3)
   scalar `sx2' = _result(4) * (`k' - 1) / `k'
   scalar `sd2' = sqrt(_result(4))

   qui corr `1' `2' if `touse' & `byg' == `j' `wgt'
   scalar `r' = _result(4)
   scalar `sl' = sign(`r') * `sd1' / `sd2' 
   #delimit ;
   scalar `rp' = min(tprob(_result(1) - 2, _result(4) *
     sqrt(_result(1) -  2) / sqrt(1 - _result(4)^2)) ,1) ;
   #delimit cr
   scalar `sxy' = `r' * sqrt(`sx2' * `sy2')
   scalar `p' = 2 * `sxy' / (`sx2' + `sy2' + (`yb' - `xb')^2)
   scalar `u' = (`yb' - `xb') / (`sx2' * `sy2')^.25

* --- variance, test, and CI for asymptotic normal approximation
   #delimit ;
   scalar `sep' = sqrt(((1 - ((`r')^2)) * (`p')^2 * (1 - 
     ((`p')^2)) / (`r')^2 + (4 * (`p')^3 * (1 - `p') * (`u')^2
     / `r') - 2 * (`p')^4 * (`u')^4 / (`r')^2 ) / (`k' - 2));
   #delimit cr
   scalar `z' = `p' / `sep'
   scalar `zp' = 2 * (1 - normprob(abs(`z')))
   scalar `ll' = `p' - `zv' * `sep'
   scalar `ul' = `p' + `zv' * `sep'

* --- statistic, variance, test, and CI for inverse hyperbolic
*      tangent transform to improve asymptotic normality
   scalar `t' = ln((1 + `p') / (1 - `p')) / 2
   scalar `set' = `sep' / (1 - ((`p')^2))
   scalar `zt' = `t' / `set'
   scalar `ztp' = 2 * (1 - normprob(abs(`zt')))
   scalar `llt' = `t' - `zv' * `set'
   scalar `ult' = `t' + `zv' * `set'
   scalar `llt' = (exp(2 * `llt') - 1) / (exp(2 * `llt') + 1)
   scalar `ult' = (exp(2 * `ult') - 1) / (exp(2 * `ult') + 1)

* Print output

   di in gr "Concordance correlation coefficient (Lin, 1989)" _n
   di in gr " rho_c   se(rho_c)   obs    [" _c
   if index("`level'",".") {
     di in gr %6.1f `level' "% CI  ]  p-value     CI type"
   }
   else { di in gr "   `level'% CI   ]  p-value     CI type" }
   di in gr _dup(63) "-"
   di in ye %6.3f `p' %10.3f `sep' %8.0f `k' %10.3f `ll' _c
   di in ye %7.3f `ul' %9.3f `zp' in gr "   asymptotic"
   di in ye _dup(24) " " %10.3f `llt' %7.3f `ult' %9.3f `ztp' _c
   di in gr "  z-transform"
  
   di _n in gr "Pearson's r =" in ye %7.3f `r' _c
   di    in gr "  Pr(r = 0) =" in ye %6.3f `rp' _c
   di    in gr "  C_b = rho_c/r =" in ye %7.3f `p' / `r'
   di    in gr "reduced major axis:   slope = " in ye %9.3f `sl' _c
   di    in gr "   intercept = " in ye %9.3f `yb'-`xb'*`sl'  

   di _n in gr "        difference           " _c
   if index("`level'",".") {
     di in gr %6.1f `level' "% limits of agreement"
   }
   else { di in gr "   `level'% limits of agreement" }
   di in gr "   average    std. dev.          (Bland & Altman, 1986)"                 
   di in gr _dup(63) "-"
   di in ye %10.3f `db' %12.3f sqrt(`dv') _c
   di in ye "         " %11.3f `dll' %11.3f `dul'
   
   if "`summary'" != "" {
     summ `1' `2' if `touse' & `byg' == `j' `wgt'
   }
* -----------------------------------------------------
* Graph options for loa plot
   if `j' == 1 & "`graph'" == "l" {

     if "`connect'" == "" { local connect "co(lll.)" }
     else {
       local lll = length("`connect'")
       if `lll' == 1 { local connect "co(`connect'l`connect'.)" }
       else if `lll' == 2 { 
         local connect = "co(`connect'" + substr("`connect'",1,1) + ".)"
       }
       else if `lll' == 3 { local connect "co(`connect'.)" }
       else { local connect = "co(" + substr("`connect'",1,4) +")" }
     }


     if "`symbol'" == "" { local symbol "sy(...o)" }
     else {
       local lll = length("`symbol'")
       if `lll' == 1 { local symbol "sy(`symbol'.`symbol'o)" }
       else if `lll' == 2 { 
         local symbol = "sy(`symbol'" + substr("`symbol'",1,1) + "o)"
       }
       else if `lll' == 3 { local symbol "sy(`symbol'o)" }
       else { local symbol = "sy(" + substr("`symbol'",1,4) +")" }
     }


     if "`pen'" == "" { local pen "pe(3532)" }
     else {
       local lll = length("`pen'")
       if `lll' == 1 { local pen "pe(`pen'5`pen'2)" }
       else if `lll' == 2 { 
         local pen = "pe(`pen'" + substr("`pen'",1,1) + "2)"
       }
       else if `lll' == 3 { local pen "pe(`pen'2)" }
       else { local pen = "pe(" + substr("`pen'",1,4) + ")" }
     }

     if "`l1title'" == "" {
       local l1title : variable label `1'
       if "`l1title'" == "" { local l1title "`1'" }
     }
     local l1title "l1(`l1title')"


     if "`t1title'" == "" {
       local t1title = "`level'% limits of agreement"
     }
     else if "`t1title'" == "." { /* "." means blank it out */
       local t1title ""  ""
     }
     local t1title "t1(`t1title')"


     if "`xlabel'" != "" { local xlabel "xlabel(`xlabel')" }
     else if index("`options'","xla") == 0 {local xlabel "xla" }

     if "`ylabel'" != "" { local ylabel "ylabel(`ylabel')" }
     else if index("`options'","yla") == 0 {local ylabel "yla" }


     if "`xscale'" == "" { local xscale "" }
     else { local xscale "xsc(`xscale')" }

     if "`yscale'" == "" { local yscale "" }
     else { local yscale "ysc(`yscale')" }


     if "`saving'" != "" { local saving "saving(`saving')" }
    }

   if "`graph'" == "l" {
     if "`t2title'" == "" {
       if "`byl'" != "" { local t2t "t2(`byl')" }
     }
     else if "`t2title'" == "." { /* "." means blank it out */
       local t2t ""
     }
     else local t2t "t2(`t2title')"

     gr `dll' `db' `dul' `d' `m' if `touse' & `byg' == `j' `wgt', /* 
        */ `pen' `connect' `symbol' `xlabel' `ylabel' yli(0)      /* 
        */ `xscale' `yscale' l2("difference of `2' and `1'")      /*
        */ b2("mean of `2' and `1'") `t1title' `t2t' `saving'     /*
        */ `options'
     more

* note: normal prob plot logic pilfered from qnorm
     tempvar Z Zv Psubi
     qui {
       sort `d' 
       gen `Psubi' = sum(`touse')
       replace `Psubi' = cond(`touse'==.,.,`Psubi'/(`Psubi'[_N] + 1))
       sum `d' if `touse' == 1, detail
       gen `Zv' = (`d' - _result(3)) / sqrt(_result(4))
       gen `Z'  = invnorm(`Psubi')
       label var `Z' "standard normal deviate"
     }
     graph `Z' `Zv' `d' if `touse' & `byg' == `j' `wgt', `t2t'    /*
        */ t1("Normal plot for differences") xla c(.l) gap(3)     /*
        */ s(oi) yla(-2,-1,0,1,2) b2("difference of `2' and `1'")

     if `byn' > 1 { more }
   } 
* ----------------------------------------------------  
* Graph options for concordance plot
   if `j' == 1 & "`graph'" == "c" {

     if "`connect'" == "" { local connect "co(.l.)" }
     else {
       local lll = length("`connect'")
       if `lll' == 1 { local connect "co(`connect'l.)" }
       else if `lll' == 2 { local connect "co(`connect'.)" }
       else { 
         local connect = "co(" + substr("`connect'",1,2) + ".)"
       }
     }


     if "`symbol'" == "" { local symbol "sy(o.i)" }
     else {
       local lll = length("`symbol'")
       if `lll' == 1 { local symbol "sy(`symbol'.i)" }
       else if `lll' == 2 { local symbol "sy(`symbol'i)" }
       else {
         local symbol = "sy(" + substr("`symbol'",1,2) + "i)"
       }
     }


     if "`pen'" == "" { local pen "pe(221)" }
     else {
       local lll = length("`pen'")
       if `lll' == 1 { local pen "pe(`pen'`pen'1)" }
       else if `lll' == 2 { local pen "pe(`pen'1)" }
       else { local pen = "pe(" + substr("`pen'",1,2) + "1)" }
     }


     if "`l1title'" == "" {
       local l1title : variable label `1'
       if "`l1title'" == "" { local l1title "`1'" }
     }
     local l1title "l1(`l1title')"

     if "`t1title'" == "" {
       local t1title = "Note: Data must" +                /*
     */  " overlay dashed line for perfect concordance"
     }
     else if "`t1title'" == "." { /* "." means blank it out */
       local t1title ""  ""
     }
     local t1title "t1(`t1title')"


     if "`t2title'" == "." { local t2title "" }
     else if "`t2title'" != "" { local t2title "t2(`t2title')"}


     if "`xlabel'" != "" { local xlabel "xlabel(`xlabel')" }
     else if index("`options'","xla") == 0 {local xlabel "xla" }

     if "`ylabel'" != "" { local ylabel "ylabel(`ylabel')" }
     else if index("`options'","yla") == 0 {local ylabel "yla" }


     if "`sort'" == "" { local sort "sort" }

     if "`saving'" != "" { local saving ", saving(`saving')" }
   }

* Graph concordance plot
   if "`graph'" == "c" {
     qui gen `rmaxis' = `sl' * (`2' - `xb') + `yb'
     gph open `saving'
       #delimit ;
       graph `1' `rmaxis' `2' `2' if `touse' & `byg' == `j'
         `wgt', `connect' `symbol' `t1title' `l1title'
         `t2title' `xlabel' `ylabel' `sort' `pen' `options';
       #delimit cr
     tempname ysca yloc xsca xloc xmin xstep
     scalar `ysca' = _result(5)
     scalar `yloc' = _result(6)
     scalar `xsca' = _result(7)
     scalar `xloc' = _result(8)
     scalar `xmin' = _result(3)
     scalar `xstep' = (_result(4) - `xmin') / 50
     local r1 = _result(2) * `ysca' + `yloc'
     local c1 = `xmin'     * `xsca' + `xloc'
     gph pen 1
     gph text `r1' `c1' 0 -1 `byl'
     local ig = 0
     while `ig' < 49 {
       local r1 = (`xmin' +  `ig'      * `xstep') * `ysca' + `yloc'
       local c1 = (`xmin' +  `ig'      * `xstep') * `xsca' + `xloc'
       local r2 = (`xmin' + (`ig' + 1) * `xstep') * `ysca' + `yloc'
       local c2 = (`xmin' + (`ig' + 1) * `xstep') * `xsca' + `xloc'
       gph line `r1' `c1' `r2' `c2'
       local ig = `ig' + 2
     }
     gph close
     if `byn' > 1 { more }
     drop `rmaxis'
   }
   if `byn' > 1 { capture drop `d' `db' `dll' `dul' `m' }
   local j = `j' + 1
 } /* end of loop for each `by' group */

* Save key values
 if `byn' == 1 {
   global S_1 = `k'
   global S_2 = `p'
   global S_3 = `sep'
   global S_4 = `ll'
   global S_5 = `ul'
   global S_6 = `llt'
   global S_7 = `ult'
   global S_8 = `p' / `r'
   global S_9 = `db'
   global S_10 = sqrt(`dv')
   global S_11 = `dll'
   global S_12 = `dul'
 }
end