* version 1.0.5 23set2004  MER  version 8.0
capture program drop kapci
program define kapci, rclass byable(recall)
	version 8.0


	if "`1'" == "" {
		di _n in gr "  Syntax is:" _n
		di    in wh "     kapci " in gr "[varlist] [if] [in] , [ " 		_c
		di    in wh "est" in gr "im(" in wh "an bc p n bsall" in gr ") "
		di    in wh _col(22) "w" in gr "gt" 					_c
                di    in gr "(" in wh "w w2 any_wgt" in gr ") " 			_c
                di    in wh "r" in gr "eps(" in wh "#" in gr ") " 			_c
                di    in wh "si" in gr "ze(" in wh "#" in gr ") "		
		di    in wh _col(22) "se" in gr "ed(" in wh "#" in gr ") "		_c
		di    in wh "ev" in gr "ery(" in wh "#" in gr ") "			_c		
		di    in wh "le" in gr "vel(" in wh "#" in gr ") " 			_c
		di    in wh "t" in gr "ab " in wh "w" in gr "ide"			
  		di    in wh _col(22) "sa" in gr "ving(" in wh "filename" in gr ") " 	_c
  		di    in wh "replace nom" in gr "sg ]"
	exit
	}


* Setup

	version 8
	syntax varlist [if] [in] [ , Reps(numlist) SIze(numlist) SEed(numlist)  	///
		EVery(numlist) Wgt(str) ESTim(str) Level(integer $S_level) 		///
	        SAving(str) REPLACE Tab WIDE NOMsg ]


	if "`options'" ~= "" { 
		local options = ", `options'" 
	}
	tokenize "`varlist'"
	marksample touse

	preserve

    	if "`if'"!="" {
       		keep `if'
    	} 
    	if ("`in'"!="") {
       		keep `in'
    	}


* Level value, etc.

	global zsc = invnorm(1-(1-`level'/100)/2)
	global N = _N


* Checking if estim=an is compatible with data

		qui inspect `1'
		local numlev `r(N_unique)'
		local nummeas : word count `varlist'	
		if `numlev' < 3  & `nummeas' < 3 { 
			local bs 0 
		}
			else local bs 1

	if "`estim'" == "an" & `bs' == 1 { 
		di in gr "  "
		di in gr "Note: Option " in wh "estim(an) " in gr "is only suitable for 2x2 data."
		di in wh "      bs " in gr " will be used."
	}

* Defaults

	if ("`estim'" == "") & (`bs' == 0) {
		local estim "an" 
	}
	if ("`estim'" == "") & (`bs' == 1) { 
		local estim "bc" 
	}	 
	

*************************************
* Displaying table if it is requested
*************************************

	if "`wide'" ~= "" {
		local wide_str ", wrap"
	}
	if "`wide'" == "" {
		local wide_str ""
	}	

	if ("`tab'" ~= "") & (`nummeas' < 3) {
		tab `varlist' if `touse' `wide_str'
	}
	if ("`tab'" ~= "") & (`nummeas' >= 3) {
		tab2 `varlist'if `touse' `wide_str'
	}


****************************************************
* Calculating analytical CI for kappa when estim=an
****************************************************

if (`bs' == 0) & ("`estim'" == "an") {		/*  Start of no bs situation */


* First ... extracting effective sample size used (noting byable!)

	qui summ `varlist' if `touse'
	local N = `r(N)'

	
* Call kappa and get midpoints

	qui kap `varlist' if `touse' `options'


* Saving scalars from kap as locals for return list

	local prop_e = r(prop_e)
	local prop_o = r(prop_o)


* Working macro ...

	local k = r(kappa)


* Extract table data 

	tempname tab2x2 a b c d agrN
	qui tab2 `varlist' , matcell(`tab2x2')

	scalar `a'=`tab2x2'[1,1]
	scalar `b'=`tab2x2'[1,2]
	scalar `c'=`tab2x2'[2,1]
	scalar `d'=`tab2x2'[2,2]
	scalar `agrN'=`a'+`b'+`c'+`d'


* Locals - marginals

	local p1 = (`a'+`b')/`agrN' 
	local p2 = (`a'+`c')/`agrN' 


* Quantity Q based on Fleiss, (1981), equations 13.15 - 13.18


 	local Q = ( ( ( (-1 + `k' ) * ( (-2*`k' *`p1' ) +          	/*    
*/	  	    ((`k' ^2)*`p1' ) + (4*`k' *(`p1'^2)) -        	/*
*/	  	    (2*(`k'^2)*(`p1'^2)) - (2*`k' *`p2' ) +     	/*
*/	  	    ((`k'^2)*`p2' ) - (4*`p1' *`p2' ) +           	/*
*/	  	    (8*`k' *`p1' *`p2' ) - (6*(`k'^2)*`p1' *`p2' ) +	/*
*/	  	    (4*(`p1'^2)*`p2' ) - (12*`k' *(`p1'^2)*`p2' ) +  	/*
*/	  	    (8*(`k'^2)*(`p1'^2)*`p2' ) + (4*`k' *(`p2'^2)) - 	/*
*/	  	    (2*(`k'^2)*(`p2'^2)) + (4*`p1' *(`p2'^2)) -     	/*
*/	  	    (12*`k' *`p1' *(`p2'^2)) +                       	/*
*/	  	    (8*(`k'^2)*`p1' *(`p2'^2)) -                   	/*
*/	  	    (4*(`p1'^2)*(`p2'^2)) +                        	/*
*/	  	    (12*`k' *(`p1'^2)*(`p2'^2)) -                  	/*
*/	  	    (8*(`k'^2)*(`p1'^2)*(`p2'^2)))                   	/*
*/	  	    ) /  ( (`p1'  + `p2'  - (2*`p1' *`p2' ) )^2 ) ) )


* Standard error, given kappa estimate  k_hat=`k' 

	local sek = (sqrt(`Q')) / (sqrt(`agrN') )


* CI

	local klow = `k'-($zsc*`sek')
	local kup = `k'+($zsc*`sek')

	if `kup' >= 1 {
		local kup = 1
	}
	if `klow' < -1 {
		local klow = -1
	}


* Display 

   	local type "analytical "
	local typeab "A"

	di _n in gr _col(42) "N=" `N' 
	di    in gr "{hline 48}" 
   	di    in gr " Kappa (" %2.0f `level' "% CI) = " in ye  %5.3f `k'	_c
   	di    in gr _col(24) " ("  in ye %5.3f `klow' in gr " - "  in ye %5.3f `kup' in gr ")" _c
   	di    in gr _col(44) "(" "`typeab'" ")"
        di    in gr "{hline 48}"  	
        di    in gr _col(2) "`typeab'" " = " "`type'"	


* Return list

	return scalar   N = `agrN'
	return scalar   z = $zsc
	return scalar   se = `sek'
	return scalar   prop_o = `prop_o'	
	return scalar   prop_e = `prop_e'	
	return scalar   ub_an = `kup'
	return scalar   lb_an = `klow'
	return scalar	kappa = `k'

		
} 	/* End of no bs situation */



*****************************************************
* Calculating analytical CI for kappa when estim!=an
*****************************************************

if (`bs' == 1)  | ((`bs' == 0) & ("`estim'" ~= "an")) {		/*  Start of bs situation */
 
      if "`estim'" ~= "an" {
	if "`estim'" ~= "" {
	  if "`estim'" ~= "bc" {
	    if "`estim'" ~= "n" {
	      if "`estim'" ~= "p" {
                if "`estim'" ~= "bsall" {
		     local estim "bc"
		     di in bl " "
		     di in bl "Note: Unknown bs CI type specified." 
		     di in wh "      bc " in gr "will be used." 
	        }
	      }
	    }
	  }
	}
      }
      if "`estim'" == "an" { 
      	  local estim "bc" 
      }
      

* Preparing ...

  tempfile mainfile 		
  qui save `mainfile', replace
  qui use `mainfile', clear
  
  tempfile tmpsave0 tmpsave1
  local byindex = _byindex()

  *----------------------------------------------------------------------------
  if ("`saving'" ~= "") & (_by()==1) { 
  	local sa_str "saving(`tmpsave1')" 
  }
  if ("`saving'" ~= "") & (_by()==0) { 
  	local sa_str "saving(`tmpsave0')" 
  }
  if "`saving'" == ""  & (_by()==0) { 
  	local sa_str "" 
  }
  *----------------------------------------------------------------------------
  if "`wgt'" ~= "" { 
  	local wgt_str "wgt(`wgt')"
   }
  local n : word count `varlist'
  if `n' > 2 & "`wgt'" ~= "" { 
	di _n in gr "Note: wgt() not allowed if varlist > 2. Option ignored." 
	local wgt_str ""
  }
  if "`wgt'" == "" { 
  	local wgt_str "" 
  }
  *----------------------------------------------------------------------------
  if "`reps'" ~= "" { 
  	local reps_str "reps(`reps')" 
  }
  if "`reps'" == "" { 
  	local reps_str "reps(5)" 
  	local reps 5
  	di _n in gr "Note: default number of bootstrap replications "		_c
	di    in gr "has been  						
  	di    in gr "      set to " in wh "5 " in gr "for syntax testing only."	_c
  	di    in wh "reps() " in gr "needs to "
  	di    in gr "      be increased when analysing real data." 		_n					
  }
  *----------------------------------------------------------------------------
  if "`seed'" ~= "" { 
  	set seed `seed' 
  	local seed_str "seed(`seed')"
  }
  *----------------------------------------------------------------------------
  if "`size'" ~= "" { 
  	if `size' < 5 { 
  		di in gr "Note: size() set to N"
  		local size $N
  	}
  	local size_str "size(`size')"
  }
  if "`size'" == "" { 
  	local size_str "" 
  	local size $N
  }
  *----------------------------------------------------------------------------
  if "`every'" ~= "" { 
  	local every_str "every(`every')" 
  }
  if "`every'" == "" { 
  	local every_str ""
  }
  *----------------------------------------------------------------------------
   if _by()==0 {
   	local first "kap `varlist' , `wgt_str'"
   }
   if _by()==1 {
   	local first "kap `varlist' if `touse' , `wgt_str'"
   }


* Calling bs
	
   if `reps' > 100 & "`nomsg'" == "" { 
  	di _n in gr "This may take quite a long time. Please wait ..." 
   }

   qui bs " `first' " r(kappa), `reps_str' `sa_str' level(`level') `size_str' `every_str' nowarn

   if ("`saving'" ~= "") & (_by()==0) {
  	qui use `tmpsave0'
	qui label data "kapci: varlist is `varlist'"
	qui rename _bs_1 _kapci_bs
	label var _kapci_bs "Options: `wgt_str' `reps_str' `seed_str' `size_str' `every_str'" 
	qui save `saving', `replace'
	restore
   }

      if ("`saving'" ~= "") & (_by()==1) {
  	qui use `tmpsave1'
	qui label data "kapci: varlist is `varlist'; byvars is `_byvars'; by-group is (`byindex')"
	qui rename _bs_1 _kapci_bs__`byindex'
	label var _kapci_bs__`byindex' "Options: `wgt_str' `reps_str' `seed_str' `size_str' `every_str'" 
	qui save `saving'`byindex', `replace'
	restore
   }


* Extracting sample size used (noting byable !)

   qui summ `varlist' if `touse'
   local N = `r(N)'


* Matrix extraction

   matrix tmp_mtx = e(b)
   local k = tmp_mtx[1,1]

   matrix tmp_mtx = e(ci_bc)
   local klow_bc = tmp_mtx[1,1]
   local kup_bc = tmp_mtx[2,1]   

   matrix tmp_mtx = e(ci_percentile)
   local klow_p = tmp_mtx[1,1]
   local kup_p = tmp_mtx[2,1]    

   matrix tmp_mtx = e(ci_normal)
   local klow_n = tmp_mtx[1,1]
   local kup_n = tmp_mtx[2,1]  

   matrix tmp_mtx = e(reps)
   local numreps = tmp_mtx[1,1]   

   matrix tmp_mtx = e(bias)
   local bias = tmp_mtx[1,1]

   matrix tmp_mtx = e(se)
   local se = tmp_mtx[1,1]
   matrix drop tmp_mtx

   
* Display 

   local dotdot "{hline 48}" 
   local col1 "_col(34)"
   local col2 "_col(43)"

   if "`estim'" ~= "bsall" {

	if "`estim'" == "bc" {
		local klow = `klow_bc'
		local kup  = `kup_bc'
	}
	
	if "`estim'" == "n" {	
		local klow = `klow_n'
		local kup  = `kup_n' 
	}

	if "`estim'" == "p" {	
		local klow = `klow_p'
		local kup  = `kup_p' 
	}

	if (`kup' >= 1) & (`kup' ~= .) {
		local kup = 1
	}
	if (`klow' < -1) & (`klow' ~= .) {
		local klow = -1
	}

   	if "`estim'" == "bc" { 
   		local type "bias corrected " 
   		local typeab "BC"
   		}

     	if "`estim'" == "n" { 
   		local type "normal " 
   		local typeab "N"
   		}

   	if "`estim'" == "p"  { 
   		local type "percentile " 
   		local typeab "P"
   		}


	di _n in gr _col(34) "B=" `numreps' _col(42) "N=" `N' 
	di    in gr "{hline 48}" 
   	di    in gr " Kappa (" %2.0f `level' "% CI) = " in ye  %5.3f `k'	_c
   	di    in gr _col(24) " ("  in ye %5.3f `klow' in gr " - "  in ye %5.3f `kup' in gr ")" _c
   	di    in gr _col(44) "(" "`typeab'" ")"
        di    in gr "{hline 48}"  	
        di    in gr _col(2) "`typeab'" " = " "`type'"

   }



   if "`estim'" == "bsall" {

	if `kup_n' >= 1 {
		local kup_n = 1
	}
	if `klow_n' < -1 {
		local klow_n = -1
	}

	local type1 "bias corrected" 
	local typeab1 "BC"
   	local type2 "percentile" 
	local typeab2 "P"
 	local type3 "normal"
 	local typeab3 "N" 


	di _n in gr _col(34) "B=" `numreps' _col(42) "N=" `N' 
	di    in gr "{hline 48}" 
   	di    in gr " Kappa (" %2.0f `level' "% CI) = " in ye  %5.3f `k'	_c
   	di    in gr _col(24) " ("  in ye %5.3f `klow_bc' in gr " - " in ye %5.3f `kup_bc' in gr ")"  _c
   	di    in gr _col(44) "(" "`typeab1'" ")"
   	di    in gr _col(24) " ("  in ye %5.3f `klow_p' in gr " - " in ye %5.3f `kup_p' in gr ")"    _c
   	di    in gr _col(44) "(" "`typeab2'" ")"  
   	di    in gr _col(24) " ("  in ye %5.3f `klow_n' in gr " - " in ye %5.3f `kup_n' in gr ")"    _c
   	di    in gr _col(44) "(" "`typeab3'" ")" 
        di    in gr "{hline 48}"  	
        di    in gr _col(2) "`typeab1'" " = " "`type1'" ", "					     _c
        di    in gr "`typeab2'" " = " "`type2'" ", "						     _c
        di    in gr "`typeab3'" " = " "`type3'"

   }   

* Return list

	return scalar   N_bs = `size'
	return scalar   N = $N
	return scalar   z = $zsc
	return scalar   reps = `numreps' 
	return scalar   bias = `bias' 
	return scalar   se = `se'   
	return scalar   lb_n = `klow_n' 
	return scalar   ub_n = `kup_n' 
	return scalar   lb_p = `klow_p' 
	return scalar   ub_p = `kup_p' 
	return scalar   lb_bc = `klow_bc'
	return scalar   ub_bc = `kup_bc'
	return scalar   kappa = `k'   

}	/*  End of no bs situation */


* Cleaning up

   macro drop zsc 
   macro drop N

end