You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
496 lines
12 KiB
Plaintext
496 lines
12 KiB
Plaintext
7 months ago
|
* 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
|
||
|
|
|