*! Date        : 29 Nov 2006
*! Version     : 1.03
*! Authors     : Adrian Mander
*! Email       : adrian.mander@mrc-hnr.cam.ac.uk
*! Description : Finding correlations and putting them into matrices

/*
version 1.03: Fixes bug due to missing data. Needed the novarlist option in marksample.
*/

pr matpwcorr
version 9.1
syntax [varlist] [if] [in] [, Gen]
marksample touse, novarlist

unab myvar : `varlist'
local nvar: list sizeof varlist

if "`gen'"~="" {
  qui gen var1=""
  qui gen var2=""
  qui gen corr=.
  qui gen pv=.
}

mat corr=I(`nvar')
mat pv=I(`nvar')-I(`nvar')
mat colnames corr= `varlist'
mat rownames corr= `varlist'
mat colnames pv= `varlist'
mat rownames pv= `varlist'

local line 1
local i 1
foreach v1 of local varlist {
  local j 1
  foreach v2 of local varlist {
     if `j'>`i' {
        qui corr `v1' `v2' if `touse'
        local t = abs(`r(rho)')* sqrt( (`r(N)'-2)/(1-`r(rho)'^2) )
        local p = 2*ttail(`r(N)'-2, `t')
        mat corr[`i',`j'] = `r(rho)'
        mat corr[`j',`i'] = `r(rho)'
        mat pv[`i',`j'] = `p'
        mat pv[`j',`i'] = `p'
        if "`gen'"~="" {
          if `line'>_N qui set obs `line'
          qui replace var1 = "`v1'" in `line'
          qui replace var2 = "`v2'" in `line'
          qui replace corr = `r(rho)' in `line'
          qui replace pv= `p'   in `line'
        }
        local `line++'
     }
     local `j++'     
  }
local `i++'
}

local ntests "`--line'"

if "`gen'"~="" {
  qui gen pv_bonf = pv * `ntests'
  qui replace pv_bonf=1 if pv_bonf>1
  qui gen pv_sidak = 1 - (1-pv)^`ntests'
  qui replace pv_sidak=1 if pv_sidak>1
  sort pv
  list var1 var2 corr pv pv_bonf pv_sidak if pv~=.
  di "A significant p-value (pv) under a Bonferroni correction would be 0.05/`ntests' tests = " 0.05/`ntests' 
  /*
  The Dunn-Sidak correction  is   ac = 1 - (1-ae)^(1/j)  
  The Bonferroni correction is    ac = ae/j
  */
  local sidak = 1 - (1-0.05)^(1/`line')
  di "A significant p-value (pv) under a Dunn-Sidak correction would be `sidak'"
}

mat pv_bonf = pv * `ntests'
mat pv_sidak = pv
forv i=1/`nvar' {
  forv j=1/`nvar' {
    mat pv_sidak[`i',`j'] = 1-(1-pv[`i',`j'])^`ntests'
    if pv_bonf[`i',`j']>1 mat pv_bonf[`i',`j']=1
    if pv_sidak[`i',`j']>1 mat pv_sidak[`i',`j']=1
  }
}

end