📄 boxcox_6.ado
字号:
*! version 3.0.15 16sep2004
*! made boxcox_6 on 22may2000
program define boxcox_6, eclass
version 6, missing
if replay() {
if "`e(cmd)'"!="boxcox_6" { error 301 }
}
else {
/*
Global macros used, other than S_*:
S_B_exp, S_B_ix, S_B_ropt, S_B_rhs, S_B_lmax, S_B_rms,
S_B_gm.
*/
syntax varlist(numeric) [if] [in] [aw fw] [, /*
*/ DELta(real .01) Generate(string) GRaph GRVars /*
*/ Iterate(integer 0) Level(integer 0) noLOg /*
*/ LStart(real 1e20) MEAn MEDian /*
*/ Ropt(string) SAVing(string) Zero(real .001) ]
if "`generat'"!="" {
conf new var `generat'
}
if "`grvars'"!="" {
conf new var _ll _lx
}
tempvar x lx ll s xtsf
marksample touse
est clear
if `level'<10 | `level'>99 {
local level 95
local ci 0
}
else local ci=`level'/100
global S_B_exp `"[`weight'`exp']"'
global S_B_ropt "`ropt'"
tokenize `varlist'
local lhs "`1'"
mac shift
global S_B_rhs "`*'"
qui gen double `x' = `lhs' if `touse'
qui summ `x'
if r(N)<3 { error 2001 }
if r(min) == r(max) {
di in red "`lhs' has no variance"
exit 498
}
if (`zero' == 0) & (`ci' > 0) {
di in red "no confidence interval when zero = 0"
exit 198
}
if `iterate' <= 0 { local iterate 20 }
local noisy = ("`log'" == "")
local finalrc 0
/*
Calc geometric mean of `x' or its transform
*/
if "`mean'"!="" | "`median'"!="" {
if "`mean'"!="" & "`median'"!="" { error 198 }
if "`median'"!="" {
qui summ `x', d
local cent = r(p50)
}
else {
qui summ `x', meanonly
local cent = r(mean)
}
local sym 1
quietly {
gen `s' = sign(`x'-`cent')
replace `x' = abs(`x'-`cent')+1
gen double `xtsf' = ln(`x')
summ `xtsf', meanonly
global S_B_gm = exp(r(mean))
}
}
else {
local sym 0
capture assert `x'>0
if _rc {
di in red "`lhs' has zero or negative values"
exit 498
}
quietly {
gen `s' = 1
gen double `xtsf' = ln(`x')
summ `xtsf', meanonly
global S_B_gm = exp(r(mean))
}
}
/*
Calculate LL for untransformed data.
*/
_crcbcll `xtsf' `x' `s' 1 1
local lraw $S_B_lmax
_crcbcll `xtsf' `x' `s' 0 $S_B_gm
local lln $S_B_lmax
_crcbcll `xtsf' `x' `s' -1 $S_B_gm
local lneg1 $S_B_lmax
local k = cond(`lstart'<.999e20, `lstart', 1)
if `noisy' & `zero'!=0 {
di in gr "(note: iterations performed using zero =" /*
*/ in ye `zero' in gr ")" _n
}
qui gen `lx' = .
qui gen `ll' = .
global S_B_ix 0
local kl .
local kh .
local target 0
local deriv 1
global S_1 0
global S_5 0
_crcbcrt `x' `lx' `ll' `s' `xtsf' `k' `zero' `delta' /*
*/ `target' `iterate' `noisy' `deriv'
if `zero'==0 {
di in gr "Log likelihood(L=" $S_1 ") = " in ye $S_5
global S_4 $S_5
exit
}
local k $S_1 /* lambda */
local ltsf $S_5 /* maximized log likelihood */
if ($S_2 == -1) {
capture noisily error 430
local finalrc 430
}
local x2 0
if (`ci' > 0) & (`k' < .) & `finalrc'==0 {
/*
Find support interval for lambda.
First calculate LL for untransformed data.
*/
local x2 = 0.5*invnorm(0.5+0.5*`ci')^2
/*
Lower confidence limit (if it exists).
If not, it is taken as missing (minus infinity).
*/
local target = -`x2'
local deriv 0
global S_1 `k'
/*
Starting value for lower limit for lambda
is lambda-0.1. Pass max LL via S_5
and lambda via S_1.
*/
local kl = `k'-.5
if `noisy' {
di _n in gr /*
*/ "Iterations for lower confidence interval:"
}
_crcbcrt `x' `lx' `ll' `s' `xtsf' `kl' `zero' /*
*/ `delta' `target' `iterate' `noisy' `deriv'
local kl $S_1
if ($S_2 == -1) {
local kl .
di in red /*
*/ "lower limit: convergence not achieved"
}
local target = -`target'
global S_1 `k'
global S_5 `ltsf'
local kh = `k'+0.5
if `noisy' {
di _n in gr /*
*/ "Iterations for upper confidence interval:"
}
_crcbcrt `x' `lx' `ll' `s' `xtsf' `kh' `zero' /*
*/ `delta' `target' `iterate' `noisy' `deriv'
local kh $S_1
if ($S_2 == -1) {
local kh .
di in red /*
*/ "upper limit: convergence not achieved"
}
}
lab var `ll' "Log likelihood"
lab var `lx' "Lambda"
if "`graph'"!="" {
if `"`saving'"'!="" {
local saving `"saving(`"`saving'"')"'
}
if (`ltsf' < .) & (`x2' > 0) {
local yl = `ltsf'-`x2'
gr7 `ll' `lx', c(s) s(O) xla yla /*
*/ yline(`yl') sort `saving'
}
else {
gr7 `ll' `lx', c(s) s(O) xla yla sort `saving'
}
}
if `finalrc'==0 {
if "`generat'"!="" | "$S_B_rhs"!="" {
if "`generat'"=="" {
tempvar generat
local depname "depname(_boxcox)"
}
if abs(`k')>1e-10 {
qui gen `generat' = `s'*((`x')^`k'-1)/`k'
version 2.1, missing
local kr=round(`k',.0001)
version 6, missing
if "`mean'"!="" {
lab var `generat'/*
*/ "s*BC(|`lhs'-mn+1|,`kr')"
}
else if "`median'"!="" {
lab var `generat'/*
*/ "s*BC(|`lhs'-md+1|,`kr')"
}
else lab var `generat' /*
*/ "BC(`lhs',`kr')"
}
else {
qui gen `generat' = `s'*ln(`x')
if "`mean'"!="" {
lab var `generat' /*
*/ "s*ln(|`lhs'-mn+1|)"
}
else if "`median'"!="" {
lab var `generat' /*
*/ "s*ln(|`lhs'-md+1|)"
}
else lab var `generat' "ln(`lhs')"
}
if "$S_B_rhs"!="" {
qui reg `generat' $S_B_rhs $S_B_exp, /*
*/ $S_B_ropt `depname'
global S_E_reg 1 /* double save */
est scalar reg = 1
}
}
}
if "`grvars'"!="" {
rename `ll' _ll
rename `lx' _lx
}
/* double save results in e() and S_E_ */
est scalar L = `k'
global S_E_L `k'
if `ci' {
est scalar lb_L = `kl'
est scalar ub_L = `kh'
global S_E_L0 `kl'
global S_E_L1 `kh'
}
else {
est scalar lb_L = .
est scalar ub_L = .
global S_E_L0
global S_E_L1
}
if "`mean'`median'" != "" {
est local center "`mean'`median'"
if "`mean'" != "" { est scalar mean = `cent' }
else { est scalar median = `cent' }
global S_E_xf "`mean'`median'"
global S_E_ctr `cent'
}
est local ll_0 /* clear from regress */
est scalar ll = `ltsf'
est local depvar "`lhs'"
est scalar level = `level'
est scalar chi2_m1 = 2*(`ltsf'-`lneg1')
est scalar chi2_0 = 2*(`ltsf'-`lln')
est scalar chi2_1 = 2*(`ltsf'-`lraw')
est scalar rc = `finalrc'
global S_E_ll `ltsf'
global S_E_depv "`lhs'"
global S_E_lvl `level'
global S_E_llm1 = e(chi2_m1)
global S_E_ll0 = e(chi2_0)
global S_E_ll1 = e(chi2_1)
global S_E_rc `finalrc'
est local predict boxcox_p
est local cmd "boxcox_6"
global S_E_cmd "boxcox_6"
mac drop S_B_ropt S_B_exp S_B_ix S_B_rhs S_B_lmax /*
*/ S_B_rms S_B_gm
} /* end estimate */
/* double save also in S_# (in addition to S_E_ and e()) */
global S_1 $S_E_L
global S_2 $S_E_L0
global S_3 $S_E_L1
global S_4 $S_E_ll
global S_5 $S_E_llm1
global S_6 $S_E_ll0
global S_7 $S_E_ll1
if "`e(center)'" != "" {
di _n in gr "Transform: " /*
*/ "sign(`e(depvar)'-`e(center)')*((|`e(depvar)'-`e(center)'|)+1^L-1)/L, `e(center)' = " /*
*/ e(`e(center)')
}
else di _n in gr "Transform: (`e(depvar)'^L-1)/L"
di _n in gr /*
*/ _col(18) "L" /*
*/ _col(26) "[`e(level)'% Conf. Interval]" /*
*/ _col(51) "Log Likelihood" _n /*
*/ _col(13) _dup(52) "-"
di _col(13) %9.4f e(L) " " _c
if e(lb_L)<. {
di %9.4f e(lb_L) " " %9.4f e(ub_L) _c
}
else di in gr " (not calculated) " _c
di " " %10.0g e(ll) _n
di in gr _col(6) "Test: L == -1" _col(26) "chi2(1) = " /*
*/ in ye %8.2f e(chi2_m1) /*
*/ in gr _col(48) "Pr>chi2 =" /*
*/ in ye %8.4f chiprob(1,e(chi2_m1))
di in gr _col(13) "L == 0" _col(26) "chi2(1) = " /*
*/ in ye %8.2f e(chi2_0) /*
*/ in gr _col(48) "Pr>chi2 =" /*
*/ in ye %8.4f chiprob(1,e(chi2_0))
di in gr _col(13) "L == 1" _col(26) "chi2(1) = " /*
*/ in ye %8.2f e(chi2_1) /*
*/ in gr _col(48) "Pr>chi2 =" /*
*/ in ye %8.4f chiprob(1,e(chi2_1))
if "`e(reg)'"!="" {
di _n in gr "(type " in ye "regress" in gr /*
*/ " without arguments for regression estimates conditional on L)"
}
error `e(rc)'
end
program define _crcbcrt
/*
Right-hand side variables for regression or anova passed in $S_B_rhs.
Returns lambda in $S_1, #iterations in $S_2, max LL in $S_5.
*/
local k0 = $S_1 /* ML est of lambda, if already calculated */
local ltsf = $S_5 /* Max log likelihood, if already calculated */
local x "`1'"
local lx "`2'"
local ll "`3'"
local s "`4'"
local xtsf "`5'"
local k `6' /* Starting value for lambda */
local zero `7' /* def of zero, 0 -> calc. just once */
local delta `8'
local target `9' /* target LL */
local itmax `10' /* max allowed iterations */
local noisy `11' /* 1 -> noisy mode */
local deriv `12'
local iter 0
/*
Transformation is standardized Box-Cox (Atkinson p.87)
*/
_crcbcll `xtsf' `x' `s' `k' $S_B_gm
local lmax $S_B_lmax
local rms $S_B_rms
if `noisy' & `zero'!=0 {
di in gr _col(4) /*
*/ "Iteration Lambda Zero Variance LL" /*
*/ _n _col(4) _dup(60) "-"
}
if `zero' == 0 { local f0 0 }
else {
if `deriv' {
local e = `k'+`delta'
_crcbcll `xtsf' `x' `s' `e' $S_B_gm
local f0 = ($S_B_lmax-`lmax')/`delta'
}
else {
local f0 = `ltsf'-$S_B_lmax
if `k' < `k0' { local f0 = -`f0' }
local f0 = `f0'-`target'
}
/* wwg */
local f0m2 .
local f0m1 .
local k2 .
local k1 .
local reason
/* end wwg */
while (abs(`f0') > `zero') & (`iter' < `itmax') {
local f0m2 `f0m1'
local k2 `k1'
local f0m1 `f0'
local k1 `k'
if `noisy' {
di in ye _col(4) %7.0f `iter' /*
*/ _col(16) %9.4f `k' /*
*/ %13.5f `f0' /*
*/ " " %11.0g $S_B_rms %13.5f $S_B_lmax /*
*/ in gr " `reason'"
}
local reason
global S_B_ix = $S_B_ix+1
if $S_B_ix<=_N {
qui replace `lx' = `k' in $S_B_ix
qui replace `ll' = `lmax' in $S_B_ix
}
local iter = `iter'+1
local e = `k'+`delta'
_crcbcll `xtsf' `x' `s' `e' $S_B_gm
if `deriv' {
local f1 $S_B_lmax
local e = `e'+`delta'
_crcbcll `xtsf' `x' `s' `e' $S_B_gm
local f1 = ($S_B_lmax-`f1')/`delta'
}
else {
local f1 = `ltsf'-$S_B_lmax
if `k' < `k0' { local f1 = -`f1' }
local f1 = `f1'-`target'
}
local m = (`f1'-`f0')/`delta'
if `m' == 0 {
noisily di in red /*
*/ "Convergence problem, doubling the value of delta."
local delta = `delta'*2
}
else {
local k = `k'-`f0'/`m'
/* BEGIN INSERT */
if `f0m2'<. & `f0m1'<. & `f0m2'*`f0m1'<0 {
/* k1 k k2 or k2 k k1 */
if !(`k1'<`k' & `k'<`k2' | /*
*/ `k2'<`k' & `k'<`k1') {
local k=(`k1'+`k2')/2
local reason "(*)"
local hasmp "yes"
}
}
if abs(`k')>20 {
local k=uniform()*2 - 1
local reason "(**)"
local hasrnd "yes"
local iter 0 /* reset */
}
/* END INSERT */
_crcbcll `xtsf' `x' `s' `k' $S_B_gm
local lmax $S_B_lmax
local rms $S_B_rms
if `deriv' {
local e = `k'+`delta'
_crcbcll `xtsf' `x' `s' `e' $S_B_gm
local f0 = ($S_B_lmax-`lmax')/`delta'
}
else {
local f0 = `ltsf'-$S_B_lmax
if `k' < `k0' { local f0 = -`f0' }
local f0 = `f0'-`target'
}
}
if $S_B_lmax >= . { local iter `itmax' }
}
if abs(`f0') > `zero' {
/*
Convergence failure - iterations set to -1.
*/
global S_2 -1
}
else global S_2 `iter'
}
if `noisy' & `zero'!=0 {
di in ye _col(4) %7.0f `iter' /*
*/ _col(16) %9.4f `k' /*
*/ %13.5f `f0' " " %11.0g $S_B_rms %13.5f $S_B_lmax /*
*/ _n _col(4) in gr _dup(60) "-"
if "`hasmp'"=="yes" {
di in gr /*
*/ " * (lambda obtained as midpoint of previous two iterations)"
}
if "`hasrnd'"=="yes" {
di in gr /*
*/ " ** (lambda obtained by randomly choosing a new starting value)"
}
if ("`hasrnd'"=="yes" | "`hasmp'"=="yes") & $S_2!=-1 {
di in gr /*
*/ " (despite this, convergence was achieved)"
}
}
global S_1 `k'
global S_5 `lmax'
end
program define _crcbcll
/*
Either regress `xtsf' on $S_B_rhs or summarize `xtsf'.
RMS, in $S_B_rms, and partially maximized likelihood in $S_B_lmax.
*/
local xtsf "`1'" /* var to be (re)created */
local x "`2'"
local sign "`3'" /* sign of `x' */
local k `4' /* Lambda */
local gmean `5' /* geometric mean of `x' */
cap drop `xtsf'
quietly {
if abs(`k')>1e-10 {
gen double `xtsf'=/*
*/ `sign'*((`x')^`k'-1)/(`k'*(`gmean')^(`k'-1))
}
else {
gen double `xtsf'=`sign'*`gmean'*ln(`x')
}
if "$S_B_rhs" == "" {
summ `xtsf' $S_B_exp, detail
local rss = (r(N)-1)*r(Var)
global S_B_rms = `rss'/r(N)
if $S_B_rms == 0 { global S_B_rms = 1e-7 }
global S_B_lmax = -0.5*r(N)*ln($S_B_rms)
}
else {
reg `xtsf' $S_B_rhs $S_B_exp, $S_B_ropt
local rss = e(rss)
global S_B_rms = `rss'/e(N)
if $S_B_rms == 0 { global S_B_rms = 1e-7 }
global S_B_lmax = -0.5*e(N)*ln($S_B_rms)
}
}
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -