*! version 2.2.3  8sep2000  TJS & NJC  bug fixes & mods  (STB-58: sg84.3)
program define concord, rclass
*  version 2.2.0 16dec1999  TJS & NJC  version 6 changes (STB-54: sg84.2)
*  version 2.1.6 18jun1998  TJS & NJC  STB-45 sg84.1
*  version 2.0.2  6mar1998  TJS & NJC  STB-43 sg84
*
* syntax: concord vary varx [fw] [if] [in] [,
*         BY(varname) Summary LEvel(real 95) Graph(str)
*               noREF REG NPsaving(str) SND(str) graph_options]

* Syntax help

 if "`1'" == "" {
   di _n in gr "Syntax is:" _n
   di in wh "  concord" in gr " vary varx [weight] [" _c
   di in wh "if" in gr " exp] [" in wh "in" in gr " range]"
   di in gr "         [ " in wh ", by(" in gr "byvar" _c
   di in wh ") s" in gr "ummary" in wh " le" in gr "vel" _c
   di in wh "(" in gr "level" in wh ") g" in gr "raph" _c
   di in wh "(c" in gr "cc|" in wh "l" in gr "oa" in wh ")"
   di in wh "             noref reg np" in gr "saving" in wh "(" _c
   di in gr "filename[" in wh ", replace" in gr "]" in wh ")"
   di in wh "             snd(" in gr "sndvar[" in wh ", replace" _c
   di in gr "]" in wh ") " in gr "graph_options ]"
   exit
 }

* Setup

 version 6
 syntax varlist(numeric min=2 max=2) [fw] [if] [in]    /*
 */ [ ,  BY(varname) Summary LEvel(real 95) Graph(str) /*
 */      noREF REG NPsavin(str) SND(str) * ]
 marksample touse
 tokenize `varlist'

* Set temporary names

 tempvar d db dll dul m byg kk bylabel
 tempname dsd 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

 if "`weight'" != "" { local wgt "[`weight'`exp']" }

* 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

 sort `touse' `by'
 qui by `touse' `by': gen byte `byg' = _n == 1 if `touse'
 if "`by'" != "" { qui gen `kk' = _n if `byg' == 1 }
 qui replace `byg' = sum(`byg')
 local byn = `byg'[_N]

* Generate `by' labels -- if required

 if "`by'" != "" {
   capture decode `by', gen(`bylabel')
   if _rc != 0 {
     local type : type `by'
     qui gen `type' `bylabel' = `by'
   }
 }

 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 `kk'
     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 `byg' == `j' `wgt'
   qui gen  `db'  = r(mean)
   scalar   `dsd' = sqrt(r(Var))
   qui gen  `dll' = `db' - `zv' * `dsd'
   qui gen  `dul' = `db' + `zv' * `dsd'
   qui gen  `m'   = (`1' + `2') / 2
   if "`reg'" == "reg" {
     tempvar fit
     qui regr `d' `m' if `byg' == `j' `wgt'
     qui predict `fit'
   }

* Concordance calculations

   qui summ `1' if `byg' == `j' `wgt'
   scalar `k'   = r(N)
   scalar `yb'  = r(mean)
   scalar `sy2' = r(Var) * (`k' - 1) / `k'
   scalar `sd1' = sqrt(r(Var))

   qui summ `2' if `byg' == `j' `wgt'
   scalar `xb'  = r(mean)
   scalar `sx2' = r(Var) * (`k' - 1) / `k'
   scalar `sd2' = sqrt(r(Var))

   qui corr `1' `2' if `byg' == `j' `wgt'
   scalar `r'   = r(rho)
   scalar `sl'  = sign(`r') * `sd1' / `sd2'
   #delimit ;
   scalar `rp'  = min(tprob(r(N) - 2, r(rho) *
     sqrt(r(N) -  2) / sqrt(1 - r(rho)^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        CI type"
   }
   else di in gr "   `level'% CI   ]     P        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'

   local dl  = length("Difference (`2' - `1')")
   local sk1 = 32 - `dl'
   local sk  = min(int(`sk1' / 2), 3)
   local sk1 = `sk1' - `sk'
   di _n in gr _sk(`sk') "Difference (`2' - `1')" _sk(`sk1') _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 `dsd' _c 
   di in ye "            " %11.3f `dll' %11.3f `dul'

   if "`summary'" != "" { summ `1' `2' if `byg' == `j' `wgt' }

* setup local options for passing to graph routines

   if "`byl'"      != "" { local byls byl("`byl'") }
   if "`snd'"      != "" { local snds snd("`snd'") }
   if "`level'"    != "" { local levs level(`level') }
   if "`npsavin'" != "" { local npsavin "npsavin(`npsavin')" }

* loa graphs

   if "`graph'" == "l" {
     #delimit ;
     gphloa `1' `2' `dll' `db' `dul' `d' `m' `byg'
       `fit' `wgt', j(`j') byn(`byn') `byls' `reg' `ref'
       `snds' `levs' `npsavin' `options' ;
     #delimit cr
   }

* ccc graphs

   local sll=`sl'
   local xbl=`xb'
   local ybl=`yb'
   if "`graph'" == "c" {
     #delimit ;
     gphccc `1' `2' `byg' `wgt', j(`j') xb(`xbl')
        yb(`ybl') sl(`sll') byn(`byn') `byls' `options' ;
     #delimit cr
   }

   if `byn' > 1 { capture drop `d' `db' `dll' `dul' `m' }
   local j = `j' + 1
 } /* end of loop for each `by' group */


* save globals

 if `byn' == 1 {
   return scalar        N = `k'
   return scalar    rho_c = `p'
   return scalar se_rho_c = `sep'
   return scalar  asym_ll = `ll'
   return scalar  asym_ul = `ul'
   return scalar  z_tr_ll = `llt'
   return scalar  z_tr_ul = `ult'
   return scalar      C_b = `p' / `r'
   return scalar     diff = `db'
   return scalar  sd_diff = `dsd'
   return scalar   LOA_ll = `dll'
   return scalar   LOA_ul = `dul'

* double save globals

   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 = `dsd'
   global S_11 = `dll'
   global S_12 = `dul'
 }
end


program define gphloa
 version 6
*-----------------------------------------------------
* loa graphs
* ----------------------------------------------------
 syntax varlist(numeric min=2) [fw] [ , J(int 1) BYN(int 1)  /*
     */ BYL(str) REG noREF SND(str)Pen(str) Connect(str)     /*
     */ Symbol(str) LEvel(real 95) L2title(str) T1title(str) /*
     */ T2title(str) XLAbel(str) YLAbel(str) B2title(str)    /*
     */ SAving(str) NPsaving(str) * ]
 tokenize `varlist'

 tempvar dll db dul d m byg Z Zv Psubi touse2
 local dll `3'
 local db  `4'
 local dul `5'
 local d   `6'
 local m   `7'
 local byg `8'
 if "`reg'" == "reg" {
   tempvar fit
   local fit `9'
 }
 if "`weight'" != "" { local wgt "[`weight'`exp']" }

* Graph options for loa plot
 if "`reg'" == "reg" & length("`connect'") < 5 { local fc "l" }
 if "`connect'" == "" { local connect "co(lll.`fc')" }
 else {
   local lll = length("`connect'")
   if `lll' == 1 { local connect "co(`connect'l`connect'.`fc')" }
   if `lll' == 2 { 
     local connect = "co(`connect'" + substr("`connect'",1,1) + ".`fc')"
   }
   if `lll' == 3 { local connect "co(`connect'.`fc')" }
   if `lll' == 4 { 
     local connect = "co(" + substr("`connect'",1,4) + "`fc')"
   }
   if `lll' >= 5 { 
     if "`reg'" == "reg" { 
       local connect = "co(" + substr("`connect'",1,5) + ")"
     }
     else { local connect = "co(" + substr("`connect'",1,4) + ")" }
   }
 }

 local fc 
 if "`reg'" == "reg" & length("`symbol'") < 5 { local fc "i" }
 if "`symbol'" == "" { local symbol "sy(iiio`fc')" }
 else {
   local lll = length("`symbol'")
   if `lll' == 1 { local symbol "sy(`symbol'i`symbol'o`fc')" }
   else if `lll' == 2 {
     local symbol = "sy(`symbol'" + substr("`symbol'",1,1) + "o`fc')"
   }
   else if `lll' == 3 { local symbol "sy(`symbol'o`fc')" }
   else if `lll' == 4 { 
     local symbol = "sy(" + substr("`symbol'",1,4) + "`fc')"
   }
   else if `lll' >= 5 { 
     if "`reg'" == "reg" {
       local symbol = "sy(" + substr("`symbol'",1,5) + ")"
     }
     else { local symbol = "sy(" + substr("`symbol'",1,4) + ")" }
   }
 }

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

 if `"`t1title'"' == `""' { local t1title "`level'% Limits Of Agreement" }
 local t1title "t1(`"`t1title'"')"

 if `"`t2title'"' == `""' {
   if "`byl'" != "" { local t2title "`byl'" }
 }
 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" }

 local name2 : variable label `2'
 if `"`name2'"' == `""' { local name2 "`2'" }

 local name1 : variable label `1'
 if `"`name1'"' == `""' { local name1 "`1'" }

 local l2t `"`l2title'"'
 if `"`l2title'"' == `""' /*
 */ { local l2title `"Difference of `name2' and `name1'"' }
 local l2title "l2t(`"`l2title'"')"

 if `"`b2title'"' == ""  /*
 */ { local b2title `"Mean of `name2' and `name1'"' }
 local b2title "b2t(`"`b2title'"')"

 if "`ref'" != "noref" { local refline "yli(0)" }

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

* graph loa plot
 #delimit ;
 gr `dll' `db' `dul' `d' `fit' `m' if `byg' == `j'
    `wgt', `pen' `connect' `symbol' `xlabel' `ylabel' `refline'
    `xscale' `yscale' `l1title' `l2title' `b2title' `t1title'
    `t2title' `saving' `options' ;
 #delimit cr
 more

* normal prob plot
* note: logic pilfered from qnorm

 mark `touse2' if `byg' == `j'
 qui {
   gsort -`touse2' `d'
   gen `Psubi' = sum(`touse2')
   replace `Psubi' = cond(`touse2'==0,.,`Psubi'/(`Psubi'[_N] + 1))
   sum `d' if `touse2' == 1, detail
   gen `Zv' = (`d' - r(mean)) / sqrt(r(Var))
   gen `Z'  = invnorm(`Psubi')
   label var `Z' "Standard Normal Deviate"
 }
 local l2title ""
 if `"`l2t'"' != `""' { local l2title "l2t(`"`l2t'"')" }

 local pen2 = "pe(" + substr("`pen'",7,2)
 if length("`pen2'") < 6 { local pen2 = substr("`pen2'",1,4) + "4)" }

 local sym2 = "sy(" + substr("`symbol'",7,2)
 if length("`sym2'") < 6 { local sym2 = substr("`sym2'",1,4) + "o)" }

 local con2 = "co(" + substr("`connect'",7,2)
 if length("`con2'") < 6 { local con2 = substr("`con2'",1,4) + "l)" }

 #delimit ;
 graph `Z' `Zv' `d' if `byg' == `j' `wgt', `t2title' `pen2' `sym2'
       t1("Normal plot for differences") xla `con2' gap(3)
       yla(-2,-1,0,1,2) b2(`"Difference of `name2' and `name1'"')
       `l2title' `b1title' `npsavin' `options' ;
 #delimit cr

* save snd var?

 if "`snd'" != "" {
   local c = index("`snd'",",")
   if `c' != 0 {
     local snd = substr("`snd'",1,`c'-1) + substr("`snd'",`c'+1, .)
   }
   local snd1 : word 1 of `snd'
   if `byg' == 1 {
     local replace : word 2 of `snd'
     if "`replace'" == "replace" { capture drop `snd1' }
     capture confirm new var `snd1'
     if _rc == 0 { qui gen `snd1' = . }
     else {
       local rc = _rc
       di in re "`snd1' exists. Use 'replace' option: snd(snd_var, replace)." 
       exit `rc'
     }
   }
   local snd `snd1'
   qui replace `snd' = `Z' if `byg' == `j'
   label var `snd' "standard normal deviate"
 }

 if `byn' > 1 { more }
end


program define gphccc
 version 6
*-----------------------------------------------------
* ccc graph
* ----------------------------------------------------
 syntax varlist(numeric min=2) [fw] [ , J(int 1) XB(real 0) /*
     */ YB(real 0) SL(real 0) BYN(int 1) BYL(str) Pen(str)  /*
     */ Connect(str) Symbol(str) L1title(str) T1title(str)  /*
     */ T2title(str) XLAbel(string) YLAbel(string) SAving(str) * ]
 tokenize `varlist'

 tempvar byg rmaxis
 tempname ysca yloc xsca xloc xmin xstep
 local byg   `3'

 if "`weight'" != "" { local wgt "[`weight'`exp']" }

* Graph options for concordance plot

 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"
 }
 local t1title "t1(`"`t1title'"')"

 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 "`saving'" != "" { local saving ", saving(`saving')" }

* Graph concordance plot

 qui gen `rmaxis' = `sl' * (`2' - `xb') + `yb'

 gph open `saving'
 #delimit ;
 graph `1' `rmaxis' `2' `2' if `byg' == `j'
       `wgt', `connect' `symbol' `t1title' `l1title'
       `t2title' `xlabel' `ylabel' sort `pen' `options';
 #delimit cr
 scalar `ysca'  = r(ay)
 scalar `yloc'  = r(by)
 scalar `xsca'  = r(ax)
 scalar `xloc'  = r(bx)
 scalar `xmin'  = r(xmin)
 scalar `xstep' =(r(xmax) - `xmin') / 50
 local r1 = r(ymax) * `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 }
end