*! version 1.2.2 01feb2010
program define xfracpred
version 8
if "`e(fp_cmd)'"!="fracpoly" error 301
if "`e(cmd)'"=="stpm" | "`e(cmd2)'"=="stpm" {
local xb [xb]
local cmd stpm
}
else if "`e(cmd2)'"!="" local cmd "`e(cmd2)'"
else local cmd "`e(cmd)'"
if "`cmd'"=="mlogit" | "`cmd'"=="mprobit" {
di as err "`cmd' not yet supported by fracpred"
exit 198
}
local dist `e(fp_dist)'
local y "`e(fp_depv)'"
syntax newvarname [, For(str) Stdp Dresid All ]
if "`dresid'"!="" & ("`cmd'"=="clogit" | "`cmd'"=="stcrreg") {
di as err "dresid option of fracpred not supported for `cmd'"
exit 198
}
// Fix to cope with -mim-
if "`e(prefix2)'"=="mim" {
global S_esample _mim_e
}
else {
tempvar esample
gen byte `esample' = e(sample)
global S_esample `esample'
}
local all = cond("`all'"=="", "if $S_esample", "")
quietly if "`for'"=="" {
tempvar yhat
if `dist'==4 {
predict double `yhat', xb
if "`all'"!="" {
replace `yhat'=. if $S_esample==0
}
}
else _predict `yhat' `all', xb
// predicted values and SE
if "`dresid'"=="" {
if "`stdp'"=="" {
gen `typlist' `varlist'=`yhat'
lab var `varlist' "predicted `e(fp_depv)'"
}
else {
_predict `typlist' `varlist' `all', stdp
lab var `varlist' "SE of predicted `e(fp_depv)'"
}
}
// deviance residuals
else {
if "`e(fp_wgt)'"!="" & "`wts'"!="nowts" {
/* for weighted residual calc */
local exp "`e(fp_wexp)'"
}
gen `typlist' `varlist'=.
residual "`cmd'" `varlist' `yhat' "`exp'"
lab var `varlist' "deviance residuals for `e(fp_depv)'"
}
exit
}
// Deal with constant
tempname V
matrix `V' = e(V)
if ("`cmd'" == "stpm") matrix `V' = `V'["xb:", "xb:"]
fvexpand `e(fp_fvl)'
local nfvl: word count `r(varlist)'
local ncons = `nfvl' + 1
if (missing(`V'[`ncons',`ncons']) | "`cmd'"=="ologit" | "`cmd'"=="oprobit") local hascons
else local hascons constant
if ("`xb'" != "") local eq eq(xb)
else local eq
// Identify `xfor' as index of "for" variable in list of predictors
fvunab for : `for', max(1)
local i 1
local xfor 0
while `i' <= e(fp_nx) {
if "`for'" == "`e(fp_x`i')'" {
local xfor `i'
continue, break
}
local ++i
}
if `xfor'==0 {
noi di as err "`for' was not an xvar"
exit 198
}
if "`e(fp_k`xfor')'" == "." {
noi di as err "`for' not in selected model"
exit 198
}
// Identify transformed variable(s) representing `for', if present
if ("`e(fp_n`xfor')'" != "") local for `e(fp_n`xfor')'
// Compute (partial) prediction and SE
if "`stdp'" != "" {
xpredict `varlist' `all', with(`for') `hascons' `eq' stdp
lab var `varlist' "se(`for')"
}
else {
xpredict `varlist' `all', with(`for') `hascons' `eq'
lab var `varlist' "xb(`for')"
}
end
program define residual
/*
Compute deviance residuals.
Approach is to set output variable `r' to missing
if these residuals are undefined. This won't hurt fracplot.
*/
args cmd r eta exp
local dist `e(fp_dist)'
local y `e(fp_depv)'
if `dist'==0 { // Normal
replace `r' = `y'-`eta'
}
else if `dist'==1 { // works only for logit/probit; for other commands, set residuals to missing
cap drop `r'
cap predict `r', deviance
if _rc gen `r' = .
}
else if `dist'==2 { // Poisson
tempvar mu
gen `mu' = exp(`eta')
replace `r' = sign(`y'-`mu') /*
*/ *sqrt(cond(`y'==0, 2*`mu', /*
*/ 2*(`y'*log(`y'/`mu')-(`y'-`mu'))))
}
else if `dist'==3 { // refit Cox model on linear predictor only
local wgt [`e(fp_wgt)' `e(fp_wexp)']
local opts `e(fp_opts)'
tempvar touse
gen byte `touse' = $S_esample
tempname ests
cap drop `r'
_estimates hold `ests'
capture quietly cox `y' `eta' `wgt' if `touse', /*
*/ `opts' mgale(`r')
local rc = _rc
_estimates unhold `ests'
if `rc' exit `rc'
drop `touse'
// deviance residuals
local s=index("`opts'","dead(")
local dead=substr("`opts'",`s'+5,.)
local s=index("`dead'",")")
local dead=substr("`dead'",1,`s'-1)
replace `r'=sign(`r')*sqrt(-2*cond(`dead'==0, /*
*/ `r', `r'+log(1-`r') ) )
}
else if `dist'==4 { // glm
cap drop `r'
predict `r', deviance
}
else if `dist'==7 { // stcox, streg, stpm
if "`cmd'"=="stcox" {
local opts `e(fp_opts)'
tempvar touse mgale
gen byte `touse' = $S_esample
tempname ests
_estimates hold `ests'
cap drop `r'
capture quietly stcox `eta' if `touse', `opts' mgale(`mgale')
local rc = _rc
predict `r', deviance
_estimates unhold `ests'
if `rc' exit `rc'
drop `touse'
}
else if "`cmd'"=="streg" {
cap drop `r'
predict `r', deviance
}
}
/*
Note: other `dist' values neither accommodated nor trapped.
Deviance residuals will be missing by default.
*/
replace `r'=. if $S_esample==0
if `"`exp'"'!="" {
parse `"`exp'"', parse("=")
/*
Weighted residuals using weights standardized Stata-fashion
*/
sum `2' if $S_esample, meanonly
replace `r' = `r'*sqrt(`2'/r(mean))
}
end
program define xpredict
* Based on stand-alone xpredict version 1.0.8 PR 04apr2012
version 11.0
syntax newvarname [if] [in], WITH(varlist fv) ///
[ A(numlist) CONStant DOUble EQ(string) MEANzero mi Stdp xb ]
// Option xb is default so is redundant, included for backwards compatibility
if "`xb'" != "" {
local xb
}
if "`mi'"!="" {
confirm var _mj
}
fvexpand `with'
local with `r(varlist)'
if "`meanzero'" != "" & "`stdp'" != "" {
// Mark complete cases among `with' for later use
tempvar touse
mark `touse'
markout `touse' `with'
local origwith `with'
}
tempname tmp b V b2 V2
matrix `b'=e(b)
matrix `V'=e(V)
if "`constant'"!="" {
local with `with' _cons
}
if missing(`"`eq'"') {
local eq : coleq `b'
local eq : word 1 of `eq'
if `"`eq'"' == "-" local eq
}
if !missing(`"`eq'"') {
* include equation name with variable(s)
tokenize `with'
local with
while "`1'"!="" {
local with `with' `eq':`1'
mac shift
}
}
matselrc `b' `b2', row(1) col(`with')
matselrc `V' `V2', row(`with') col(`with')
local nc=colsof(`b2')
if "`a'"!="" {
local na: word count `a'
if `na'!=`nc' {
di in red "wrong number of elements in a(), should be `nc'"
exit 198
}
* construct matrix from `a' by hand
tempname A
matrix `A'=J(1,`nc',0)
tokenize `a'
local i 1
while `i'<=`nc' {
matrix `A'[1,`i']=``i''
local i=`i'+1
}
local cn: colnames `b2'
* Linear combination of elements in b2 with a.
* Multiply A by b2 elementwise and overwrite b2.
qui matewm `A' `b2' `b2'
matrix colnames `b2'=`cn'
* Form outer product of A with itself then multiply by V2 elementwise,
* overwriting V2.
matrix `A'=`A''*`A'
qui matewm `A' `V2' `V2'
matrix colnames `V2'=`cn'
matrix rownames `V2'=`cn'
}
_estimates hold `tmp'
capture ereturn post `b2' `V2'
local rc=_rc
if `rc' {
_estimates unhold `tmp'
error `rc'
}
qui if "`mi'"!="" {
* Partial prediction for all obs, then average over imputations
sum _mj, meanonly
local m=r(max)
tempvar f
if `"`if'"'=="" {
local If if _mj>0
}
else local If `if' & _mj>0
capture predict `double' `f' `If' `in', xb
local rc=_rc
_estimates unhold `tmp'
if `rc' {
error `rc'
}
sort _mi _mj
by _mi: gen `double' `varlist'=sum(`f')/`m' if _mj>0
by _mi: replace `varlist'=`varlist'[_N] if _n<_N & _mj>0
}
else {
tempvar newvar
predict `double' `newvar' `if' `in', `stdp' `options'
if "`meanzero'" != "" {
if "`stdp'" == "" {
qui sum `newvar' `if' `in'
qui gen `double' `varlist' = `newvar' - r(mean)
}
else {
matrix `V' = e(V)
mata: zcalc("`varlist'", "`V'", "`origwith'", "`touse'")
}
}
else rename `newvar' `varlist'
_estimates unhold `tmp'
}
end
* NJC 1.1.0 20 Apr 2000 (STB-56: dm79)
program def matselrc
* NJC 1.0.0 14 Oct 1999
version 6.0
gettoken m1 0 : 0, parse(" ,")
gettoken m2 0 : 0, parse(" ,")
if "`m1'" == "," | "`m2'" == "," | "`m1'" == "" | "`m2'" == "" {
di in r "must name two matrices"
exit 198
}
syntax , [ Row(str) Col(str) Names ]
if "`row'`col'" == "" {
di in r "nothing to do"
exit 198
}
tempname A B
mat `A' = `m1' /* this will fail if `matname' not a matrix */
local cols = colsof(`A')
local rows = rowsof(`A')
if "`col'" != "" {
if "`names'" != "" { local colnum 1 }
else {
capture numlist "`col'", int r(>0 <=`cols')
if _rc == 0 { local col "`r(numlist)'" }
else if _rc != 121 {
local rc = _rc
error `rc'
}
local colnum = _rc == 0
}
/* colnum = 1 for numbers, 0 for names */
tokenize `col'
local ncols : word count `col'
if `colnum' {
mat `B' = `A'[1..., `1']
local j = 2
while `j' <= `ncols' {
mat `B' = `B' , `A'[1..., ``j'']
local j = `j' + 1
}
}
else {
mat `B' = `A'[1..., "`1'"]
local j = 2
while `j' <= `ncols' {
mat `B' = `B' , `A'[1..., "``j''"]
local j = `j' + 1
}
}
mat `A' = `B'
local cols = colsof(`A')
}
if "`row'" != "" {
if "`names'" != "" { local rownum 0 }
else {
capture numlist "`row'", int r(>0 <=`rows')
if _rc == 0 { local row "`r(numlist)'" }
else if _rc != 121 {
local rc = _rc
error `rc'
}
local rownum = _rc == 0
}
/* rownum = 1 for numbers, 0 for names */
tokenize `row'
local nrows : word count `row'
if `rownum' {
mat `B' = `A'[`1', 1...]
local j = 2
while `j' <= `nrows' {
mat `B' = `B' \ `A'[``j'', 1...]
local j = `j' + 1
}
}
else {
mat `B' = `A'["`1'", 1...]
local j = 2
while `j' <= `nrows' {
mat `B' = `B' \ `A'["``j''", 1...]
local j = `j' + 1
}
}
mat `A' = `B'
}
mat `m2' = `A'
end
program define matewm
* 1.0.1 NJC 15 June 1999 STB-50 dm59
* 1.0.0 NJC 21 July 1998
version 5.0
parse "`*'", parse(" ,")
if "`3'" == "" | "`3'" == "," {
di in r "invalid syntax"
exit 198
}
matchk `1'
local A "`1'"
matchk `2'
local B "`2'"
matcfa `A' `B'
local nr = rowsof(matrix(`A'))
local nc = colsof(matrix(`A'))
local C "`3'"
mac shift 3
local options "Format(str)"
parse "`*'"
tempname D
mat `D' = J(`nr',`nc',1)
local i 1
while `i' <= `nr' {
local j 1
while `j' <= `nc' {
mat `D'[`i',`j'] = `A'[`i',`j'] * `B'[`i',`j']
local j = `j' + 1
}
local i = `i' + 1
}
if "`format'" == "" { local format "%9.3f" }
mat `C' = `D' /* allows overwriting of either `A' or `B' */
mat li `C', format(`format')
end
* 1.0.0 NJC 19 July 1998 STB-50 dm69
program def matcfa
* matrices conformable for addition?
* matcfa `1' `2'
version 5.0
if "`1'" == "" | "`2'" == "" | "`3'" != "" {
di in r "invalid syntax"
exit 198
}
tempname C
mat `C' = `1' + `2'
end
* 1.0.0 NJC 5 July 1998 STB-50 dm69
program def matchk
* matrix?
* matchk `1'
version 5.0
if "`1'" == "" | "`2'" != "" {
di in r "invalid syntax"
exit 198
}
tempname C
mat `C' = `1'
end
mata:
void zcalc(string scalar newvarname, string scalar cov_beta, string scalar xvarlist,
string scalar tousevar)
{
real colvector xbar
// Form views of data in Mata
V = st_matrix(cov_beta)
xvars = tokens(xvarlist)
st_view(X=., ., xvars, tousevar)
n = rows(X)
k = cols(X)
xbar = J(k, 1, 0)
for(j=1; j<=k; j++) {
xbar[j, 1] = sum(X[., j])
}
xbar = xbar :/ n
st_view(result=., ., st_addvar("float", newvarname), tousevar)
for(i=1; i<=n; i++) {
x1 = X[i, .]' - xbar
result[i] = sqrt(x1' * V * x1)
}
}
end