* version 4.3, DEC2013
* Artur Pokropek & Maciej Jakubowski
* regression with PIAAC data
cap program drop piaacreg
program define piaacreg, rclass
syntax [varlist(default=none)] [if] [in] ///
, save(string) [ weight(varlist max=1) rep(varlist) ///
countryid(varname) vemethodn(varname numeric) pvdep(string) pvindep1(string) pvindep2(string) pvindep3(string) ///
over(string) fast round(integer 2) cmd(string) cmdops(string) r2(string) cons or]
version 10.0
tempvar probka cntid jkvar proflevel1
tempname b se over_values
local oldover="`over'"
local longlist=0
local overpv=0
if "`over'"=="pvlit" | "`over'"=="pvnum" | "`over'"=="PVLIT" | "`over'"=="PVNUM" {
qui recode `over'1 (.=.) (0/175.9999=0) (176/225.9999=1) (226/275.9999=2) (276/325.9999=3) (326/375.9999=4) (376/999=5), gen(`proflevel1')
local overpv=1
local over="`proflevel1'"
}
else if "`over'"=="pvpsl" | "`over'"=="PVPLS" {
qui recode pvpsl1 (.=.) (0/240.9999=0) (241/290.9999=1) (291/340.9999=2) (341/999=3), gen(`proflevel1')
local overpv=1
local over="`proflevel1'"
}
if "`cmd'"=="" local cmd="reg"
if "`r2'"=="" {
if "`cmd'"=="reg" local r2="r2"
else if "`cmd'"=="logit" | "`cmd'"=="logistic" | "`cmd'"=="poisson" | "`cmd'"=="ologit" local r2="r2_p"
}
if "`pvdep'"=="" {
tokenize `varlist'
local outcome `1'
macro shift
local varlist `*'
local pv=0
}
else {
local outcome `pvdep'1
local pv=1
}
if "`cmd'"=="logit" | "`cmd'"=="logistic" {
sum `outcome', meanonly
if r(min)!=0 | r(mean)<=0 {
di as error "outcome has to vary and contain 0 and positive values only"
}
}
if "`cmd'"=="poisson" {
sum `outcome', meanonly
if r(mean)<=0 {
di as error "outcome has to vary and contain nonnegative values only"
}
}
if "`pvindep1'"!="" | "`pvindep2'"!="" | "`pvindep3'"!="" | `overpv'!=0 local pv=1
if "`countryid'"!="" {
capture confirm string variable `countryid'
if !_rc {
qui gen `cntid'=`countryid'
qui levelsof `cntid' `if' `in', local(cnt)
}
else {
qui tostring(`countryid'), gen(`cntid')
qui levelsof `cntid' `if' `in', local(cnt)
}
}
else {
cap confirm variable cntryid
if !_rc {
qui tostring(cntryid), gen(`cntid')
qui local cnt="36 40 124 203 208 233 246 250 276 372 380 392 410 528 578 616 703 724 752 840 SUB 56 826 OECD 196 643"
local longlist=1
}
else {
cap confirm variable CNTRYID
if !_rc {
qui tostring(CNTRYID), gen(`cntid')
qui local cnt="36 40 124 203 208 233 246 250 276 372 380 392 410 528 578 616 703 724 752 840 SUB 56 826 OECD 196 643"
local longlist=1
}
else {
di as error "Your dataset does not contain variable cntryid or CNTRYID. Please specify option countryid()"
exit 100
}
}
}
if "`weight'"=="" {
cap confirm variable spfwt0
if !_rc local weight="spfwt0"
else {
cap confirm variable SPFWT0
if !_rc local weight="SPFWT0"
else {
di as error "Your dataset does not contain main survey weight spfwt0 or SPFWT0. Please specify weight() option"
exit 100
}
}
}
local ok_small=0
local ok_big=0
if "`rep'"=="" {
forvalues i=1(1)80 {
cap confirm variable spfwt`i'
if !_rc local ok_small=`ok_small'+1
cap confirm variable SPFWT`i'
if !_rc local ok_big=`ok_big'+1
}
if `ok_small'==80 local rep="spfwt1-spfwt80"
else if `ok_big'==80 local rep="SPFWT1-SPFWT80"
else {
di as error "Your dataset does not contain all replicate weights spfwt1-spfwt80 or SPFWT1-SPFWT80. Please specify rep() option"
exit 100
}
}
cap confirm variable vemethodn
if !_rc {
qui gen `jkvar'=vemethodn
}
else {
cap confirm variable VEMETHODN
if !_rc {
qui gen `jkvar'=VEMETHODN
}
else {
if "`vemethodn'"=="" {
di as error "Your dataset does not contain variable vemethodn or VEMETHODN. Please specify option vemethodn()"
}
else gen `jkvar'=`vemethodn'
}
}
if "`if'"=="" local if=" if "
else local if=`"`if' "'+" & "
tempvar test
qui gen `test'=1
foreach var of local varlist {
capture confirm numeric variable `var'
if !_rc qui replace `test'=`test'*`var'
else di as error "`var' is not a numeric variable"
}
if "`pvdep'"=="" {
capture confirm numeric variable `outcome'
if !_rc qui replace `test'=`test'*`outcome'
else di as error "`outcome' is not a numeric variable"
}
else {
qui replace `test'=`test'*`pvdep'1
local outcome=""
}
if "`pvdep'"!="" qui replace `test'=`test'*`pvdep'1
if "`pvindep1'"!="" qui replace `test'=`test'*`pvindep1'1
if "`pvindep2'"!="" qui replace `test'=`test'*`pvindep2'1
if "`pvindep3'"!="" qui replace `test'=`test'*`pvindep3'1
if `overpv'==1 qui replace `test'=`test'*`oldover'1
local n_vars: word count `pvindep1' `pvindep2' `pvindep3' `varlist' `cons'
local n_cnt: word count `cnt'
if `longlist'==1 local n_cnt=`n_cnt'-2
if "`over'"!="" {
tab `over', nofreq matrow(over_values)
local n_over=r(r)
}
else {
local over="1"
matrix over_values=(1,1)
local n_over=1
}
tempname tempfile
file open `tempfile' using "`save'.html", write replace
file write `tempfile' `"
"' " | "
* TABLE HEADER - OVER CATEGORIES - if over() specified
if `n_over'>1 {
local col_span=2*`n_vars'+1
file write `tempfile' `" Over categories: | "'
forvalues i=1(1)`n_over' {
local numerek=over_values[`i',1]
if `overpv'!=1 {
local labelka: label (`over') `numerek'
file write `tempfile' `"`labelka' | "'
}
else {
if `numerek'==0 {
file write `tempfile' `"Below Level 1 | "'
}
else {
file write `tempfile' `"Level `numerek' | "'
}
}
}
}
* SECOND ROW - predictors
file write `tempfile' `" Variable | "'
forvalues i=1(1)`n_over' {
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
file write `tempfile' `"`var' | "'
}
file write `tempfile' `"`r2' | "'
}
* third row
file write `tempfile' "" "Country" " | "
forvalues i=1(1)`n_over' {
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
file write `tempfile' `""' "Coef." " | " `""' "S.E." " | "
local OECD_`var'_`i'=0
local OECD_SE_`var'_`i'=0
}
file write `tempfile' `""' "" " | "
local OECD_r2_`i'=0
}
local decimal=0.1^`round'
forvalues i=1(1)`n_over' {
local noc_`i'=0
* creating output matrices
mat b_over`i'=J(`n_cnt',`n_vars',.)
mat se_over`i'=J(`n_cnt',`n_vars',.)
mat r2_over`i'=J(`n_cnt',1,.)
}
local OECD_saved=0
local pvlist="`outcome'"
if "`pvindep1'"!="" local pvlist="`pvlist' `pvindep1'*"
if "`pvindep2'"!="" local pvlist="`pvlist' `pvindep2'*"
if "`pvindep3'"!="" local pvlist="`pvlist' `pvindep3'*"
if "`pvdep'"!="" local pvlist="`pvlist' `pvdep'*"
if `overpv'==1 local pvlist="`pvlist' `oldover'*"
local country_list=""
local cnt_i=0
foreach l of local cnt {
if "`l'"=="OECD" | "`l'"=="SUB" {
if "`l'"=="SUB" file write `tempfile' `"
Sub-national entities | "'
else {
local OECD_saved=1
file write `tempfile' `"
OECD Average | "'
forvalues i=1(1)`n_over' {
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
local mean=string(round(`OECD_`var'_`i''/`noc_`i'',`decimal'),"%20.`round'f")
local se=string(round(sqrt(`OECD_SE_`var'_`i'')/`noc_`i'',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | `se' | "'
}
local mean=string(round(`OECD_r2_`i''/`noc_`i'',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | "'
}
file write `tempfile' `"
Partners | "'
}
}
else {
_cnt `l'
local name=r(name)
local country_list="`country_list'`l' "
local cnt_i=`cnt_i'+1
di ""
di as result "`name'" _continue
file write `tempfile' "
" "`name'" " | "
sum `jkvar' if `cntid'=="`l'" , meanonly
if `r(N)'>0 {
local jk=1
sum `jkvar' if `cntid'=="`l'" , meanonly
cap local jk=`r(mean)'
if `jk'<1 | `jk'>2 {
di as error "Please specify values 1 or 2 in `jkvar'"
exit 100
}
forvalues i=1(1)`n_over' {
if 1!=`n_over' di as text " `i'" _c
sum `test' if `cntid'=="`l'" & `over'==over_values[`i',1], meanonly
if r(N)>2 {
return clear
cap drop `probka'
if `overpv'==0 qui gen `probka'=1 `if' `cntid'=="`l'" & `over'==over_values[`i',1] `in'
else qui gen `probka'=1 `if' `cntid'=="`l'" `in'
local temp=over_values[`i',1]
preserve
qui keep if `probka'==1
keep `pvlist' `varlist' `probka' `weight' `rep'
if `pv'==1 | `overpv'==1 {
if "`fast'"=="" _piaacregpv `outcome' `varlist', overval(`temp') overpv(`oldover') if(`probka') pvdep(`pvdep') pvindep1(`pvindep1') pvindep2(`pvindep2') pvindep3(`pvindep3') weight(`weight') rep(`rep') jk(`jk') cmd(`cmd') cmdops(`cmdops') r2(`r2') `cons' `or'
else _piaacregfastpv `outcome' `varlist', overval(`temp') overpv(`oldover') if(`probka') pvdep(`pvdep') pvindep1(`pvindep1') pvindep2(`pvindep2') pvindep3(`pvindep3') weight(`weight') cmd(`cmd') cmdops(`cmdops') r2(`r2') `cons' `or'
}
else {
if "`fast'"=="" _piaacreg `outcome' `varlist', if(`probka') weight(`weight') rep(`rep') jk(`jk') cmd(`cmd') cmdops(`cmdops') r2(`r2') `cons' `or'
else _piaacregfast `outcome' `varlist', if(`probka') weight(`weight') cmd(`cmd') cmdops(`cmdops') r2(`r2') `cons' `or'
}
matrix b=r(coef)
matrix se=r(se)
local noc_`i'=`noc_`i''+1
mat temp=b'
mat temp= temp[1,1..`n_vars']
mat b_over`i'[`cnt_i',1]=temp
mat temp=se'
mat temp= temp[1,1..`n_vars']
mat se_over`i'[`cnt_i',1]=temp
local row=1
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
local b=b[`row',1]
local se=se[`row',1]
local row=`row'+1
if `b'==0 & `se'==0 file write `tempfile' `" 0 | "' `" 0 | "'
else {
local mean=string(round(`b',`decimal'),"%20.`round'f")
local se=string(round(`se',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | "' `"`se' | "'
local OECD_`var'_`i'=`OECD_`var'_`i''+`mean'
local OECD_SE_`var'_`i'=`OECD_SE_`var'_`i''+`se'^2
}
}
local b=b[`row',1]
mat r2_over`i'[`cnt_i',1]=`b'
local mean=string(round(`b',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | "'
local OECD_r2_`i'=`OECD_r2_`i''+`mean'
restore
}
else {
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
file write `tempfile' `" . | "' `" . | "'
}
file write `tempfile' `" . | "'
}
}
}
}
}
if `OECD_saved'==0 {
file write `tempfile' `"
Average | "'
forvalues i=1(1)`n_over' {
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
local mean=string(round(`OECD_`var'_`i''/`noc_`i'',`decimal'),"%20.`round'f")
local se=string(round(sqrt(`OECD_SE_`var'_`i'')/`noc_`i'',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | `se' | "'
}
local mean=string(round(`OECD_r2_`i''/`noc_`i'',`decimal'),"%20.`round'f")
file write `tempfile' `" `mean' | "'
}
}
file write `tempfile' _n "
"
file close `tempfile'
di ""
di "Results saved in the `save'.html file"
forvalues i=1(1)`n_over' {
matrix colnames b_over`i' = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons'
matrix colnames se_over`i' = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons'
matrix rownames b_over`i' = `country_list'
matrix rownames se_over`i' = `country_list'
matrix rownames r2_over`i' = `country_list'
matrix colnames r2_over`i' = `r2'
}
if `n_over'==1 {
return matrix b=b_over1
return matrix se=se_over1
return matrix r2=r2_over1
}
else {
forvalues i=1(1)`n_over' {
return matrix b_over`i'=b_over`i'
return matrix se_over`i'=se_over`i'
return matrix r2_over`i'=r2_over`i'
}
}
end
***********
cap program drop _cnt
program define _cnt, rclass
args l
local name="`l'"
if "`l'"=="36" local name="Australia"
if "`l'"=="40" local name="Austria"
if "`l'"=="124" local name="Canada"
if "`l'"=="203" local name="Czech Republic"
if "`l'"=="276" local name="Germany"
if "`l'"=="208" local name="Denmark"
if "`l'"=="724" local name="Spain"
if "`l'"=="233" local name="Estonia"
if "`l'"=="246" local name="Finland"
if "`l'"=="250" local name="France"
if "`l'"=="826" local name="England/N. Ireland (UK)"
if "`l'"=="372" local name="Ireland"
if "`l'"=="380" local name="Italy"
if "`l'"=="392" local name="Japan"
if "`l'"=="410" local name="Korea"
if "`l'"=="528" local name="Netherlands"
if "`l'"=="578" local name="Norway"
if "`l'"=="616" local name="Poland"
if "`l'"=="643" local name="Russian Federation"
if "`l'"=="703" local name="Slovak Republic"
if "`l'"=="752" local name="Sweden"
if "`l'"=="840" local name="United States"
if "`l'"=="56" local name="Flanders (Belgium)"
if "`l'"=="196" local name="Cyprus"
return local name "`name'"
end
****************************************************
cap program drop _piaacreg
program define _piaacreg, rclass
syntax [varlist], weight(varlist max=1) rep(varlist) jk(integer) ///
cmd(string) if(varname numeric) [cmdops(string) r2(string) cons or]
tokenize `varlist'
local outcome `1'
macro shift
local varlist `*'
* SETUP
local num_rep : word count `rep'
local n_vars: word count `varlist'
if "`cons'"!="" {
local n_vars=`n_vars'+1
local stala="_cons"
}
matrix coef = J(`n_vars'+1,1,.)
matrix rownames coef =`varlist' `cons' `r2'
matrix se= J(`n_vars',1,.)
matrix rownames se =`varlist' `cons'
qui `cmd' `outcome' `varlist' [pw=`weight'] if `if'==1, `cmdops'
foreach var in `varlist' `stala' {
local fulstat_`var'=_b[`var']
local sumofsq_`var'=0
}
local fulstat_r2=e(`r2')
foreach replication of local rep {
qui `cmd' `outcome' `varlist' [pw=`replication'] if `if'==1, `cmdops'
foreach var in `varlist' `stala' {
local sq_`var' = (_b[`var'] - `fulstat_`var'')^2
local sumofsq_`var'=`sumofsq_`var''+ `sq_`var''
}
}
foreach var in `varlist' `stala' {
if `jk'==1 local sdr_`var'=sqrt(((`num_rep'-1)/`num_rep')*`sumofsq_`var'')
else local sdr_`var'=sqrt(1*`sumofsq_`var'')
}
local row=1
foreach var in `varlist' `stala' {
if "`or'"=="" {
matrix coef[`row',1] = `fulstat_`var''
matrix se[`row',1] = `sdr_`var''
}
else {
matrix coef[`row',1] = exp(`fulstat_`var'')
matrix se[`row',1] = exp(`fulstat_`var'')*`sdr_`var''
}
local row=`row'+1
}
matrix coef[`row',1] = `fulstat_r2'
matrix colnames coef = b
matrix colnames se = se
return matrix coef =coef
return matrix se =se
end
*************************************
cap program drop _piaacregfast
program define _piaacregfast, rclass
syntax [varlist], weight(varlist max=1) ///
cmd(string) if(varname numeric) [cmdops(string) r2(string) cons or]
tokenize `varlist'
local outcome `1'
macro shift
local varlist `*'
* SETUP
local n_vars: word count `varlist'
if "`cons'"!="" {
local n_vars=`n_vars'+1
local stala="_cons"
}
matrix coef = J(`n_vars'+1,1,.)
matrix rownames coef =`varlist' `cons' `r2'
matrix se= J(`n_vars',1,.)
matrix rownames se =`varlist' `cons'
qui `cmd' `outcome' `varlist' [pw=`weight'] if `if'==1, `cmdops'
local row=1
foreach var in `varlist' `stala' {
if "`or'"=="" {
matrix coef[`row',1] = _b[`var']
matrix se[`row',1] = _se[`var']
}
else {
matrix coef[`row',1] = exp(_b[`var'])
matrix se[`row',1] = exp(_b[`var'])*_se[`var']
}
local row=`row'+1
}
matrix coef[`row',1] = e(`r2')
matrix colnames coef = b
matrix colnames se = se
return matrix coef =coef
return matrix se =se
end
**************
cap program drop _piaacregpv
program define _piaacregpv, rclass
syntax [varlist(numeric default=none)], [overpv(string) overval(string) pvdep(string) pvindep1(string) pvindep2(string) pvindep3(string)] ///
weight(varlist max=1) rep(varlist) jk(integer) cmd(string) if(varname numeric) [cmdops(string) r2(string) cons or]
if "`overpv'"=="pvlit" | "`overpv'"=="pvnum" | "`overpv'"=="PVLIT" | "`overpv'"=="PVNUM" {
tempvar proflevel1 proflevel2 proflevel3 proflevel4 proflevel5 proflevel6 proflevel7 proflevel8 proflevel9 proflevel10
forvalues i=1(1)10 {
qui recode `overpv'`i' (.=.) (0/175.9999=0) (176/225.9999=1) (226/275.9999=2) (276/325.9999=3) (326/375.9999=4) (376/999=5), gen(`proflevel`i'')
}
local overpv=1
}
else if "`overpv'"=="pvpsl" | "`overpv'"=="PVPSL" {
tempvar proflevel1 proflevel2 proflevel3 proflevel4 proflevel5 proflevel6 proflevel7 proflevel8 proflevel9 proflevel10
forvalues i=1(1)10 {
qui recode pvpsl`i' (.=.) (0/240.9999=0) (241/290.9999=1) (291/340.9999=2) (341/999=3), gen(`proflevel`i'')
}
local overpv=1
}
else local overpv=0
tempvar probkaif
if "`pvdep'"=="" {
tokenize `varlist'
local outcome `1'
macro shift
local varlist `*'
local pv=0
}
else local pv=1
local num_rep : word count `rep'
local n_vars: word count `varlist' `pvindep1' `pvindep2' `pvindep3'
if "`cons'"!="" {
local n_vars=`n_vars'+1
local stala="_cons"
}
matrix coef = J(`n_vars'+1,10,.)
matrix rownames coef =`varlist' `pvindep1' `pvindep2' `pvindep3' `cons' `r2'
matrix se= J(`n_vars',10,.)
matrix rownames se =`varlist' `pvindep1' `pvindep2' `pvindep3' `cons'
forvalues ipv=1(1)10 {
di as text "." _c
if "`pvindep1'"!="" local pvindep="`pvindep1'`ipv'"
if "`pvindep2'"!="" local pvindep="`pvindep' `pvindep2'`ipv'"
if "`pvindep3'"!="" local pvindep="`pvindep' `pvindep3'`ipv'"
if `pv'!=0 local outcome="`pvdep'`ipv'"
cap drop `probkaif'
if `overpv'==1 qui gen `probkaif'=1 if `proflevel`ipv''==`overval'
else local probkaif=1
qui `cmd' `outcome' `pvindep' `varlist' [pw=`weight'] if `probkaif'==1, `cmdops'
foreach var in `pvindep' `varlist' `stala' {
local fulstat_`var'=_b[`var']
local sumofsq_`var'=0
}
local fulstat_r2=e(`r2')
if `fulstat_r2'==. local fulstat_r2=0
foreach replication of local rep {
qui `cmd' `outcome' `pvindep' `varlist' [pw=`replication'] if `probkaif'==1, `cmdops'
foreach var in `pvindep' `varlist' `stala' {
local sq_`var' = (_b[`var'] - `fulstat_`var'')^2
local sumofsq_`var'=`sumofsq_`var''+ `sq_`var''
}
}
foreach var in `pvindep' `varlist' `stala' {
if `jk'==1 local rVAR_`var'=(((`num_rep'-1)/`num_rep')*`sumofsq_`var'')
else local rVAR_`var'=(1*`sumofsq_`var'')
}
local row=1
foreach var in `pvindep' `varlist' `stala' {
matrix coef[`row',`ipv'] = `fulstat_`var''
matrix se[`row',`ipv'] = `rVAR_`var''
local row=`row'+1
}
matrix coef[`row',`ipv'] = `fulstat_r2'
}
mat coef=coef'
mat se=se'
mata: pvmat()
if "`or'"!="" {
local row=1
foreach var in `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' {
matrix newcoef[`row',1] = exp(newcoef[`row',1])
matrix newse[`row',1] = newcoef[`row',1]*newse[`row',1]
local row=`row'+1
}
}
matrix rownames newcoef = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' `r2'
matrix colnames newcoef = b
matrix rownames newse = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons'
matrix colnames newse = se
return matrix coef =newcoef
return matrix se =newse
end
**************
cap program drop _piaacregfastpv
program define _piaacregfastpv, rclass
syntax [varlist(numeric default=none)], [overpv(string) overval(string) pvdep(string) pvindep1(string) pvindep2(string) pvindep3(string)] ///
weight(varlist max=1) cmd(string) if(varname numeric) [cmdops(string) r2(string) cons or]
if "`overpv'"=="pvlit" | "`overpv'"=="pvnum" | "`overpv'"=="PVLIT" | "`overpv'"=="PVNUM" {
tempvar proflevel1 proflevel2 proflevel3 proflevel4 proflevel5 proflevel6 proflevel7 proflevel8 proflevel9 proflevel10
forvalues i=1(1)10 {
qui recode `overpv'`i' (.=.) (0/175.9999=0) (176/225.9999=1) (226/275.9999=2) (276/325.9999=3) (326/375.9999=4) (376/999=5), gen(`proflevel`i'')
}
local overpv=1
}
else if "`overpv'"=="pvpsl" | "`overpv'"=="PVPSL" {
tempvar proflevel1 proflevel2 proflevel3 proflevel4 proflevel5 proflevel6 proflevel7 proflevel8 proflevel9 proflevel10
forvalues i=1(1)10 {
qui recode pvpsl`i' (.=.) (0/240.9999=0) (241/290.9999=1) (291/340.9999=2) (341/999=3), gen(`proflevel`i'')
}
local overpv=1
}
else local overpv=0
tempvar probkaif
if "`pvdep'"=="" {
tokenize `varlist'
local outcome `1'
macro shift
local varlist `*'
local pv=0
}
else local pv=1
local n_vars: word count `varlist' `pvindep1' `pvindep2' `pvindep3'
if "`cons'"!="" {
local n_vars=`n_vars'+1
local stala="_cons"
}
matrix coef = J(`n_vars'+1,10,.)
matrix rownames coef =`varlist' `pvindep1' `pvindep2' `pvindep3' `cons' `r2'
matrix se= J(`n_vars',10,.)
matrix rownames se =`varlist' `pvindep1' `pvindep2' `pvindep3' `cons'
forvalues ipv=1(1)10 {
di as text "." _c
if "`pvindep1'"!="" local pvindep="`pvindep1'`ipv'"
if "`pvindep2'"!="" local pvindep="`pvindep' `pvindep2'`ipv'"
if "`pvindep3'"!="" local pvindep="`pvindep' `pvindep3'`ipv'"
if `pv'!=0 local outcome="`pvdep'`ipv'"
cap drop `probkaif'
if `overpv'==1 qui gen `probkaif'=1 if `proflevel`ipv''==`overval'
else local probkaif=1
qui `cmd' `outcome' `pvindep' `varlist' [pw=`weight'] if `probkaif'==1, `cmdops'
local row=1
foreach var in `pvindep' `varlist' `stala' {
if "`or'"=="" {
matrix coef[`row',`ipv'] = _b[`var']
matrix se[`row',`ipv'] = _se[`var']
}
else {
matrix coef[`row',`ipv'] = exp(_b[`var'])
matrix se[`row',`ipv'] = exp(_b[`var'])*_se[`var']
}
local row=`row'+1
}
matrix coef[`row',`ipv'] = e(`r2')
}
mat coef=coef'
mat se=se'
mata: pvmat()
matrix rownames newcoef = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons' `r2'
matrix colnames newcoef = b
matrix rownames newse = `pvindep1' `pvindep2' `pvindep3' `varlist' `cons'
matrix colnames newse = se
return matrix coef =newcoef
return matrix se =newse
end
***************
version 10
clear mata
mata:
void pvmat()
{
pvs_org=st_matrix("coef")
pvs=pvs_org[|1,1\10,cols(pvs_org)-1|]
repVAR=mean(st_matrix("se"))
pvdif=pvs:-mean(pvs)
sumofsqpv=pvdif'pvdif
VAR=1.1:*sumofsqpv:/9
pvVAR= diagonal(VAR)'
totVAR=repVAR+pvVAR
SDR=totVAR:^0.5
st_matrix("newse", SDR')
st_matrix("newcoef", mean(pvs_org)')
}
end