*! 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