*! version 1.0.0 PR 13oct2021
program define metatef, rclass
version 11.0
local cmdline : copy local 0
mata: _parse_colon("hascolon", "rhscmd")
if !`hascolon' error 198
_metatef_i `"`0'"' `"`rhscmd'"'
return local cmdline `"metatef_i `cmdline'"'
end
program define _metatef_i, rclass
version 11.0
args 0 statacmd
gettoken cmd statacmd: statacmd
if substr("`cmd'", -1, .) == "," {
local cmd = substr("`cmd'", 1, length("`cmd'") - 1)
local statacmd ,`statacmd'
}
if "`cmd'"=="stpm" | "`cmd'"=="stpm2" {
local dist 7
local glm 0
local qreg 0
local xtgee 0
local normal 0
local eqxb eq(xb)
local vminmax 1
}
else {
*xfrac_chk `cmd'
cmdchk `cmd'
if `s(bad)' {
di as error "invalid or unrecognised command, `cmd'"
exit 198
}
if "`cmd'" == "streg" local eqxb eq(_t)
/*
dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm),
5 (xtgee), 6(ereg/weibull), 7 (streg/stpm/stpm2).
*/
local dist `s(dist)'
local glm `s(isglm)'
local qreg `s(isqreg)'
local xtgee `s(isxtgee)'
local normal `s(isnorm)'
local vminmax 2
}
* `using' creates a new file of results
syntax , BY(string) FIXPowers(string) with(varname) [ ADJust(string) CENtre(string) ///
DEGree(int 0) EBayes(string) FUNction(string) GENerate(string) GENWt(string) ///
GENVariance(string) MEANzero RANdom STRata(string) STUdywise TAU(string) ]
if `dist'!=7 & "`strata'"!="" {
di as err "strata() allowed only with survival models (stcox, streg, stpm, stpm2)"
exit 198
}
if "`tau'"!="" & "`random'"=="" {
di as err "tau() valid only with random"
exit 198
}
if "`studywise'"!="" & "`fixpowers'"!="" confirm matrix `fixpowers'
if `degree'==0 & "`fixpowers'"=="" local fixpowers 1
if "`fixpowers'"!="" & `degree'!=0 {
di as err "fixpowers() invalid with degree()"
exit 198
}
if `degree'<0 {
di as err "invalid degree(), must be at least 1"
exit 198
}
if "`centre'"=="" | "`centre'"=="mean" {
local centre mean // centre on the global mean of xvar. has no influence on the fitted functions.
}
else confirm num `centre'
if "`generate'"!="" {
confirm new var `generate'
}
confirm var `by'
if "`adjust'"!="" {
unab adjust: `adjust'
}
// Parse to extract varlist from statacmd
local 0 `statacmd'
if `dist'==7 {
local vminmax 1
}
else local vminmax 2
syntax varlist(min=`vminmax' max=`vminmax' numeric) [if] [in] [aw fw pw iw] [, DEAD(str) noCONStant * ]
if !missing("`dead'") local options `options' dead(`dead')
if !missing("`constant'") local options `options' `constant'
marksample touse, novarlist
markout `touse' `varlist' `by' `with' `strata' `adjust', strok
// `with' must be binary. We recode it to 0,1 and call it `t'.
quietly {
tempvar t t_save
egen int `t' = group(`with') if `touse'
sum `t', meanonly
if r(max) > 2 {
noi di as err "`with' must have exactly 2 values"
exit 198
}
replace `t' = `t' - 1
clonevar `t_save' = `t'
// Deal with weights (includes replacing touse with 0 if zero weight)
frac_wgt `"`exp'"' `touse' `"`weight'"'
local wgt `r(wgt)'
if `dist'==7 {
local anddead "& _d==1"
local xvar `varlist'
}
else gettoken yvar xvar: varlist
tempvar x study fbag sfbag
tempname z shift scale
scalar `z'=-invnorm((100-$S_level)/200)
egen int `study'=group(`by') if `touse'
sum `study', meanonly
local nstudy=r(max)
* Extract original study codes for integer study numbers 1,...,nstudy
levelsof `by', local(levels)
local j 0
while (`"`levels'"' != "") {
local ++j
gettoken studyval`j' levels : levels
}
/*
Extract fixpowers into strings from matrix, if fixpowers used with studywise
First row of matrix is used for overall model, rest of rows are powers for each study.
*/
local fixpowmat 0
if "`studywise'"!="" & "`fixpowers'"!="" {
local fixpowmat 1
local rows=rowsof(`fixpowers')
if `rows'!=(`nstudy'+1) {
di as err "number of powers in fixpowers() not equal to number of studies plus 1"
exit 198
}
local cols=colsof(`fixpowers')
forvalues j=0/`nstudy' {
local j1=`j'+1
forvalues i=1/`cols' {
local p=`fixpowers'[`j1',`i']
if !missing(`p') local pwrs`j' `pwrs`j'' `p'
}
}
}
* Adjust xvar to its overall unweighted mean
fracgen `xvar' 0 if `touse', nogen
scalar `shift'=r(shift)
scalar `scale'=r(scale)
gen `x'=(`xvar'+`shift')/`scale' if `touse'
if "`centre'"=="mean" {
sum `xvar' if `touse'
local centre=r(mean)
}
local centx=(`centre'+`shift')/`scale'
local m `degree'
if `degree'>=1 local degree degree(`m')
else local degree
noi di as txt "[`xvar' centered at " %8.0g `centre' "]"
noi di as txt "Processing ... overall" _cont
* Pooled-data estimate of deviance
if `dist'==7 { // survival models
if ("`cmd'"=="stpm" | "`cmd'"=="stpm2") local strat stratify(`study' `strata')
else local strat strata(`study' `strata')
}
if `dist'==7 { // survival models
if ("`cmd'"=="stpm" | "`cmd'"=="stpm2") {
local strat stratify(`study' `strata')
local strat2 stratify(`strata')
}
else {
local strat strata(`study' `strata')
local strat2 strata(`strata')
}
}
if "`studywise'" == "" {
// Single set of powers to be used for all data, supplied in fixpowers().
local pwrs `fixpowers'
}
else {
// FP powers as supplied in 1st row of fixpowers() matrix
local pwrs `pwrs0'
}
local overallpowers `pwrs'
fracgen `x' `pwrs' if `touse', replace adjust(`centx') noscaling
local v `r(names)'
if `dist' == 7 {
`cmd' `yvar' `t'##c.(`v') `adjust' `wgt' if `touse', `strat' `options'
}
else {
`cmd' `yvar' `t'##c.(`v') `adjust' i.`study' `wgt' if `touse', `options'
}
local devstrat=-2*e(ll)
// Calc fixed-effects weighted mean function, fbag
noi di as txt " ... study " _cont
tempvar sumw se
local devsum 0
gen double `sumw' = 0 if `touse'
gen double `fbag' = 0 if `touse'
forvalues j=1/`nstudy' {
noi di as txt `j', _cont
tempvar f`j' v`j' w`j' // fitted values, variance, weight for study j
if "`studywise'" != "" { // use studywise fixedpowers
// Note that fracgen is calculated with centering over all studies
fracgen `x' `pwrs`j'' if `touse', replace adjust(`centx') noscaling
local v `r(names)'
}
// Refit model and make predictions for all observations
`cmd' `yvar' `t'##c.(`v') `adjust' `wgt' if `study'==`j' & `touse'==1, `options' `strat2'
local dev`j' = -2 * e(ll)
local devsum = `devsum' + `dev`j''
if "`meanzero'" != "" {
// Compute column vector of mean(s) of predictor(s) for adjustment
tempname means
local mm
tokenize `v'
while "`1'" != "" {
sum `1' if `study'==`j' & `touse'==1
local mx = r(mean)
if "`mm'" == "" local mm `mx'
else local mm `mm' \ `mx'
mac shift
}
matrix `means' = `mm'
local Mean mean(`means')
}
replace `t' = 1
xpredict2 `f`j'', with(1.`t' 1.`t'#c.(`v')) double `eqxb'
if "`meanzero'" != "" {
sum `f`j'' if `study'==`j' & `touse'==1
replace `f`j'' = `f`j'' - r(mean)
}
xpredict2 `se', with(1.`t' 1.`t'#c.(`v')) double `eqxb' stdp `Mean'
replace `t' = `t_save'
* Compute weights for each function
gen double `v`j''=`se'^2 if `touse'
gen double `w`j''=1/`v`j'' if `touse'
replace `sumw'=`sumw'+`w`j''
replace `fbag'=`fbag'+`w`j''*`f`j''
drop `se'
}
replace `fbag' = cond(`sumw'==0, 0, `fbag'/`sumw') if `touse'
* Calc SE of mean function
gen double `sfbag' = cond(`sumw'==0, 0, `sumw'^(-0.5)) if `touse'
if "`random'"!="" {
* Calculate estimate of random-effects variance, tausq
tempvar Q tausq sumw2
gen double `Q'=0 if `touse'
gen double `sumw2'=0 if `touse'
forvalues j=1/`nstudy' {
replace `Q'=`Q'+`w`j''*(`f`j''-`fbag')^2
replace `sumw2'=`sumw2'+`w`j''^2
}
gen `tausq' = cond(`sumw'==0, 0, max(0, (`Q'-(`nstudy'-1))/(`sumw'-`sumw2'/`sumw'))) if `touse'
* Calculate estimated function with random-effects weights and its SE
* This requires "wstar" random-effects weights
replace `sumw'=0 if `touse'
replace `fbag'=0 if `touse'
forvalues j=1/`nstudy' {
replace `w`j''=1/(`v`j''+`tausq')
replace `w`j'' = 0 if missing(`w`j'') & `touse'==1
replace `sumw'=`sumw'+`w`j''
replace `fbag'=`fbag'+`w`j''*`f`j''
}
replace `fbag' = cond(`sumw'==0, 0, `fbag'/`sumw') if `touse'
replace `sfbag' = cond(`sumw'==0, 0, `sumw'^(-0.5)) if `touse'
if "`ebayes'"!="" {
* Calculate empirical Bayes estimates and SE for function
* (only applicable with random-effects estimates)
forvalues j=1/`nstudy' {
cap drop `ebayes'`j'
cap drop `ebayes'se`j'
gen `ebayes'`j' = cond(`v`j''==0, 0, (`tausq'*`f`j''/`v`j''+`fbag')/(`tausq'/`v`j''+1)) if `touse'
gen `ebayes'se`j' = cond(`v`j''==0, 0, sqrt(`tausq'*`v`j''/(`tausq'+`v`j'')+(`v`j''/(`tausq'+`v`j''))^2/`sumw')) if `touse'
lab var `ebayes'`j' "empirical Bayes f(`xvar') for study `studyval`j''"
lab var `ebayes'se`j' "SE(empirical Bayes f(`xvar')) for study `studyval`j''"
}
}
}
if "`function'"!="" {
* Save individual functions and SE
forvalues j=1/`nstudy' {
cap drop `function'`j'
cap drop `function'se`j'
gen `function'`j'=`f`j''
gen `function'se`j'=sqrt(`v`j'')
lab var `function'`j' "f(`xvar') for study `studyval`j''"
lab var `function'se`j' "fixed-effects SE(f(`xvar')) for study `studyval`j''"
}
}
if "`genvariance'"!="" {
* Save within-study variance functions
forvalues j=1/`nstudy' {
cap drop `genvariance'`j'
gen `genvariance'`j'=`v`j''
lab var `genvariance'`j' "study `studyval`j''"
}
}
if "`genwt'"!="" {
* Save weight functions, standardized by sum of weights
forvalues j=1/`nstudy' {
cap drop `genwt'`j'
gen `genwt'`j' = cond(`sumw'==0, 0, `w`j''/`sumw') if `touse'
lab var `genwt'`j' "study `studyval`j''"
}
}
if "`tau'"!="" {
cap drop `tau'
rename `tausq' `tau'
lab var `tau' "Random-effects component of variance between studies"
}
drop `v'
if "`generate'"!="" {
cap drop `generate'_ll
cap drop `generate'_ul
cap drop `generate'_se
gen `generate'_se=`sfbag'
gen `generate'_ll=`fbag'-`z'*`sfbag'
gen `generate'_ul=`fbag'+`z'*`sfbag'
rename `fbag' `generate'
lab var `generate' "`frs' average function f(`xvar')"
lab var `generate'_se "standard error of f(`xvar')"
lab var `generate'_ll "lower confidence limit for f(`xvar')"
lab var `generate'_ul "upper confidence limit for f(`xvar')"
}
if "`random'"!="" {
local frs "random-effects"
}
else local frs "fixed-effects"
* Store deviance and d.f.
return scalar deviance=`devsum'
return scalar devstrat=`devstrat'
forvalues j=1/`nstudy' {
return scalar dev`j'=`dev`j''
}
}
if !`fixpowmat' {
if "`fixpowers'"!="" {
local df: word count `fixpowers'
local DF=`df'*`nstudy'
local dfhomog=`df'*(`nstudy'-1)
}
else {
if "`studywise'"=="" { // same powers for all studies and overall
local DF=`m'*`nstudy'+`m'
local dfhomog=`m'*(`nstudy'-1)
}
else {
local DF=2*`m'*`nstudy'
local dfhomog=2*`m'*(`nstudy'-1)
}
}
return scalar df=`DF'
return scalar dfhomog=`dfhomog'
return scalar aic=return(deviance)+2*return(df)
return scalar P=chi2tail(`dfhomog', return(devstrat)-return(deviance))
di as text _n(2) "Total deviance over studies: " as res _col(30) return(deviance) as text ", df = " as res `DF'
di as text "Stratified model deviance: " as res _col(30) return(devstrat) as text ", df = " as res `DF'-`dfhomog'
di as text "Heterogeneity chi-square: " as res _col(30) return(devstrat)-return(deviance) ///
as text ", df = " as res `dfhomog' as text " (P = " as res %5.3f return(P) as text ")"
}
else {
di as text _n(2) "Total deviance over studies: " as res _col(30) return(deviance) // as text ", df = " as res `DF'
di as text "Stratified model deviance: " as res _col(30) return(devstrat) // as text ", df = " as res `DF'-`dfhomog'
}
di as text _n "Powers for stratified model = " as res "`overallpowers'"
return local powers `overallpowers'
if "`studywise'"!="" {
forvalues j=1/`nstudy' {
return local powers`j' `pwrs`j''
}
}
end
* version 1.0.6 PR 15feb2013
program define xpredict2
version 8.0
syntax newvarname [if] [in], WITH(varlist fv) [ CONStant A(numlist) DOUble EQ(string) mi stdp mean(string) * ]
if "`mi'"!="" {
confirm var _mj
}
fvexpand `with'
local with `r(varlist)'
if "`mean'" != "" & "`stdp'" != "" {
/*
`mean' is a matrix containing adjustment for each variable in `with';
could be the mean of each variable; option could be called `centre'
*/
confirm matrix `mean'
marksample touse, novarlist
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
capture predict `double' `newvar' `if' `in', `stdp' `options'
local rc=_rc
if `rc' error `rc'
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'", "`mean'", "`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
* version 1.0.4 PR 30sep2005
* Based on private version 6.1.4 of frac_chk, PR 25aug2004
program define cmdchk, sclass
version 7
local cmd `1'
mac shift
local cmds `*'
sret clear
if substr("`cmd'",1,3)=="reg" {
local cmd regress
}
if "`cmds'"=="" {
tokenize clogit cnreg cox ereg fit glm logistic logit poisson probit /*
*/ qreg regress rreg weibull xtgee streg stcox stpm stpm2 /*
*/ ologit oprobit mlogit nbreg
}
else tokenize `cmds'
sret local bad 0
local done 0
while "`1'"!="" & !`done' {
if "`1'"=="`cmd'" {
local done 1
}
mac shift
}
if !`done' {
sret local bad 1
*exit
}
/*
dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm),
5 (xtgee), 6 (ereg/weibull), 7 (stcox, streg, stpm, stpmrs).
*/
if "`cmd'"=="logit" | "`cmd'"=="probit" /*
*/ |"`cmd'"=="clogit"| "`cmd'"=="logistic" /*
*/ |"`cmd'"=="mlogit"| "`cmd'"=="ologit" | "`cmd'"=="oprobit" {
sret local dist 1
}
else if "`cmd'"=="poisson" {
sret local dist 2
}
else if "`cmd'"=="cox" {
sret local dist 3
}
else if "`cmd'"=="glm" {
sret local dist 4
}
else if "`cmd'"=="xtgee" {
sret local dist 5
}
else if "`cmd'"=="cnreg" | "`cmd'"=="ereg" | "`cmd'"=="weibull" | "`cmd'"=="nbreg" {
sret local dist 6
}
else if "`cmd'"=="stcox" | "`cmd'"=="streg" | "`cmd'"=="stpm" | "`cmd'"=="stpm2" {
sret local dist 7
}
else if substr("`cmd'",1,2)=="st" {
sret local dist 7
}
else sret local dist 0
sret local isglm = (`s(dist)'==4)
sret local isqreg = ("`cmd'"=="qreg")
sret local isxtgee= (`s(dist)'==5)
sret local isnorm = ("`cmd'"=="regress"|"`cmd'"=="fit"|"`cmd'"=="rreg")
end
mata:
void zcalc(string scalar newvarname, string scalar cov_beta, string scalar xvarlist, string scalar Xbar, string scalar tousevar)
{
xbar = st_matrix(Xbar)
// Form views of data in Mata
V = st_matrix(cov_beta)
xvars = tokens(xvarlist)
st_view(X=., ., xvars, tousevar)
n = rows(X)
st_view(result=., ., st_addvar("double", newvarname), tousevar)
// Should be possibe to do this in one matrix calc!
for(i=1; i<=n; i++) {
// result[i] = x1' * V * x1 + xbar' * V * xbar - 2 * x1' * V * xbar
x1 = X[i, .]' - xbar
result[i] = sqrt(x1' * V * x1)
}
}
end