*! 3may19 minor changes in the help file *! 17mar19 change initial values for par (ml, reml) *! 16jan19 allow drmeta to return fitted cases *! N.Orsini v.1.0.0 10oct18 capture program drop drmeta program drmeta, eclass version 13 if replay() { if "`e(cmd)'"!="drmeta" { `e(cmd)' `0' error 301 } display_result `0' exit } syntax varlist (min=2) [if] [in] , [ /// Se(varname numeric) /// Data(varlist numeric min=2) /// or rr hr md smd /// Id(varname) /// Type(varname) /// Level(integer $S_level) /// eform /// vwls /// ACov(varlist max=1) MCov(string) VCov(varlist) /// Hamling /// fixed reml ml mm /// 2stage /// NORETable /// STDDEViations /// NOLRt /// ] local random "fixed" if ("`ml'"!="" | "`reml'" != "" | "`mm'" != "") local random "random" local stage "One-stage" if "`2stage'" != "" local stage "Two-stage" // check that the two-stage method is specified with multiple studies if "`id'" == "" & "`2stage'" != "" { di as err "specifies the option id()" exit 198 } // Check estimation options local mvopt "fixed" if "`ml'" != "" local mvopt "ml" if "`reml'" != "" local mvopt "reml" if "`mm'" != "" local mvopt "mm" // Check the method for estimating the covariance local covmethod "GL" if "`hamling'" != "" local covmethod "H" if "`hr'"!="" { local ptcohort = "ptcohort" } if "`or'"!="" { local casecontrol = "casecontrol" } if "`rr'"!="" { local pecohort = "pecohort" } if "`md'"!="" { local meandiff = "meandiff" } if "`smd'"!="" { local stdmeandiff = "stdmeandiff" } // check level if `level' < 10 | `level' > 99 { di in red "level() must be between 10 and 99 inclusive" exit 198 } tempname levelci global levelp = `level'/100 scalar `levelci' = `level' * 0.005 + 0.50 global levelci = `level' * 0.005 + 0.50 // to use the option if/in marksample touse, strok novarlist global touse = `touse' // check observation qui count if `touse' local nobs = r(N) if `nobs'== 0 { di in red "no observations" exit 2000 } if `nobs'== 1 { di in red "insufficient number of observations" exit 2001 } // Check options for covariances if "`data'" != "" & "`acov'" != "" { di as err "choose either data() or acov()" exit 198 } if "`data'" != "" & "`mcov'" != "" { di as err "choose either data() or mcov()" exit 198 } if "`data'" != "" & "`vcov'" != "" { di as err "choose either cov() or vcov()" exit 198 } if "`acov'" != "" & "`mcov'" != "" { di as err "choose either acov() or mcov()" exit 198 } // Parsing input variables tempvar n case control v quietly { tokenize `varlist' local depname = "`1'" mac shift local xname `*' local ncov : word count `xname' /*qui count if `se' == . if r(N) != 0 { di as err "std error cannot be missing" exit 198 } */ // qui gen double `v' = `se'^2 if `touse' tempvar order gen `order' = _n // check number of covariates and number of observations qui count if `touse' local nobstouse = r(N)-1 if `ncov' > `nobstouse' { di in r "number of covariates greater than the number of observations" exit 198 } // get the variables required to if "`data'"!="" { parse "`data'" , parse(" ") gen double `n' = `1' if `touse' gen double `case' = `2' if `touse' capture confirm numeric var `n' capture confirm numeric var `case' } local multiple = 0 if ("`id'"!="") { tempvar typestudy // id clonevar `typestudy' = `type' if `touse' capture confirm numeric var `id' tab `typestudy' if `touse' qui tab `id' if `touse' local nstudies = r(r) if `nstudies' > 1 local multiple = 1 capture confirm numeric var `typestudy' capture assert (`typestudy' == 1) | (`typestudy' == 2) | (`typestudy' == 3) | (`typestudy' == 4) | (`typestudy' == 5) if `touse' if _rc != 0 { di _n as err "variable `type' can take value 1 (odds ratio), 2 (hazard ratio), 3 (risk ratio), 4 (mean difference), or 5 (std mean difference)" } // Center the design matrix about the referent with multiple studies // The lowest std error is zero for the chosen referent foreach x of local xname { capture drop _c_`x' bysort `id' (`se'): gen _c_`x' = `x'- `x'[1] if `touse' local xnamec "`xnamec' _c_`x'" } local varlist "`depname' `xnamec'" } // sort `order' // For the 2stage approach, // check that the number of regression coefficients to be estimated is at least // as great as the number of observations (except ref) for each study. tempvar nobss tagnobss if "`2stage'" != "" { if "`id'" == "" { di as err "specify the option id()" exit 198 } bysort `id': gen `nobss' = _N-1 if `touse' gen `tagnobss' = `nobss'>=`ncov' if `touse' qui count if `tagnobss' == 0 if r(N) > 0 { di as err _n "Nr. of coefs greater than available data for one or more study" exit 198 } } // Center the design matrix with a single study if ("`id'"=="") { foreach x of local xname { capture drop _c_`x' qui su `x' if `se' == 0 & `touse' gen _c_`x' = `x'- r(mean) local xnamec "`xnamec' _c_`x'" } local varlist "`depname' `xnamec'" } } // end quietly // Multiple studies if ("`id'"!="") { local i = 1 qui levelsof `id' if `touse', local(studies) qui count if `touse' local obs = `r(N)' - `nstudies' local countstudy = 1 foreach study of local studies { if "`acov'"=="" & "`mcov'"=="" & "`vcov'"=="" { qui drmeta_cov `varlist' if `id' == `study' & `touse' , s(`se') tot(`n') c(`case') /// ts(`typestudy') `vwls' ecovmethod(`covmethod') } if "`acov'"!="" { qui drmeta_cov `varlist' if `id' == `study' & `touse' , s(`se') /// ts(`typestudy') `vwls' meancov(`acov') } // With multiple studies the user can specify the a list of covariance matrix following the order of the studies. local passm : word `countstudy++' of `mcov' if "`mcov'"!="" { qui drmeta_cov `varlist' if `id' == `study' & `touse' , s(`se') /// ts(`typestudy') `vwls' matrixcov(`passm') } if "`vcov'"!="" { qui drmeta_cov `varlist' if `id' == `study' & `touse' , s(`se') /// ts(`typestudy') `vwls' varlistcov(`vcov') } // required step for the two-stage meta-analysis tempname M_bc M_vc * get study-specific estimates of the parameters mat `M_bc' = e(BC) mat `M_vc' = e(VB) scalar sltb`study' = `M_bc'[1,1] scalar sltvar`study' = `M_vc'[1,1] mat _G`study' = e(C) scalar col`study' = colsof(_G`study') mat _VS`study' = e(VB) mat _BS`study' = e(BC) mat _ZS`study' = colsof(_BS`study') mat _X`study' = e(Z) mat _C`study' = e(C) mat _B`study' = e(B) local i = `i' + 1 } tempvar ids tempname X Y G ID sort `id', stable mkmat `xnamec' if `depname' != 0 & `touse', matrix(`X') mkmat `depname' if `depname' != 0 & `touse', matrix(`Y') mkmat `id' if `depname' != 0 & `touse', matrix(`ID') * create a block-diagonal matrix using my blockaccum mata function local i = 1 foreach study of local studies { if `i' == 1 { mat PG = _G`study' } if `i' != 1 { mata: blockaccum("PG","_G`study'") mat PG = r(CACC) } local ++i } mat `G' = PG mat drop PG // Generalized Least Squares Estimation (Pool-first approach) tempvar b tempname V BETAS b COV tau2 Q ll R2 D p_D // Pool-first approach or one-stage approach, GLS fixed-effects mata: glsest("`Y'", "`X'", "`G'") mat `G' = r(COVAR) mat `V' = r(VB) mat `b' = r(BC)' mata: id = st_matrix("`ID'") mata: _C_ = cholesky(st_matrix("`G'")) mata: ystar = luinv(_C_)*st_matrix("`Y'") mata: Xstar = luinv(_C_)*st_matrix("`X'") mata: _ydev = st_matrix("`Y'") mata: _Xdev = st_matrix("`X'") mata: _bdev = st_matrix("`b'") mata: D = (_ydev:-_Xdev*_bdev')'*invsym(st_matrix("`G'"))*(_ydev:-_Xdev*_bdev') mata: _b_ = invsym(Xstar'*Xstar)*Xstar'*ystar mata: ypred = Xstar*_b_ mata: tresidual = ystar-ypred mata: R2=1-sum((ystar-ypred):^2)/sum(ystar:^2) mata: st_numscalar("`D'", D) mata: st_numscalar("`R2'", R2) scalar `p_D' = chi2tail(rowsof(`Y')-colsof(`X'), `D') tempname ty tX tres tpred mata: st_matrix("`ty'", ystar) mata: st_matrix("`tX'", Xstar) mata: st_matrix("`tpred'", ypred) mata: st_matrix("`tres'", tresidual) local i = 1 foreach study of local studies { tempname R2_`study' mata: _bs = st_matrix("_BS`study'") mata: ypred = select(Xstar, id:==`study')*_bs mata: R2=1-sum((select(ystar,id:==`study'):-ypred):^2)/sum(select(ystar, id:==`study'):^2) mata: st_numscalar("`R2_`study''", R2) * mata: estar = invsym(st_matrix("_G`l'") )*(select(_ydev, id:==`l')-select(_Xdev, id:==`l')*_bdev') * mata: select(_ydev, id:==`l'), estar } scalar `Q' = r(Q) local Q = `Q' if ("`2stage'" == "" & "`ml'" == "" & "`reml'" == "") { mata: X = st_matrix("`X'") mata: y = st_matrix("`Y'") mata: C = st_matrix("`G'") mata: id = st_matrix("`ID'") mata: Z = X mata: nvc = cols(Z)*(cols(Z)+1)/2 mata: parf = J(1,nvc,0) quietly mata: est = glsfit_fixed(X, Z, y, C, parf, id) mata: loglik = *est[8] mata: beta = *est[6]' mata: vcov = *est[7] mata: st_matrix("r(b)", beta) mata: st_matrix("r(V)", vcov) mata: st_matrix("r(V)", vcov) mata: st_numscalar("r(loglik)", loglik) mata: st_numscalar("r(aic)", -2*loglik+2*cols(beta)) matrix `b' = r(b) matrix `V' = r(V) } if "`ml'" != "" & "`2stage'" == "" { mata: X = st_matrix("`X'") mata: y = st_matrix("`Y'") mata: C = st_matrix("`G'") mata: id = st_matrix("`ID'") mata: Z = X mata: nvc = cols(Z)*(cols(Z)+1)/2 mata: init = J(1,nvc,.5) mata: S = optimize_init() mata: optimize_init_evaluator(S, &mlprof_fn()) mata: optimize_init_evaluatortype(S, "v0") mata: optimize_init_argument(S, 1, X) mata: optimize_init_argument(S, 2, Z) mata: optimize_init_argument(S, 3, y) mata: optimize_init_argument(S, 4, C) mata: optimize_init_argument(S, 5, id) mata: optimize_init_params(S, init) mata: optimize_init_which(S,"max") quietly mata: parf = optimize(S) mata: S = . mata: est = glsfit_ml(X, Z, y, C, parf, id) mata: loglik = *est[8] mata: beta = *est[6]' mata: vcov = *est[7] mata: st_matrix("r(b)", beta) mata: st_matrix("r(V)", vcov) mata: st_numscalar("r(loglik)", loglik) mata: st_matrix("r(Psi)", par2Psi(parf, cols(Z))) mata: st_numscalar("r(aic)", -2*loglik+2*(cols(beta)+cols(parf))) mata: st_numscalar("r(nre)", cols(parf)) // get log-likelihood setting all var/cov equal to zero mata: parf = J(1, cols(parf), 0) mata: est = glsfit_ml(X, Z, y, C, parf, id) mata: st_numscalar("r(loglik0)", *est[8]) matrix `b' = r(b) matrix `V' = r(V) } if "`reml'" != "" & "`2stage'" == "" { mata: X = st_matrix("`X'") mata: y = st_matrix("`Y'") mata: C = st_matrix("`G'") mata: id = st_matrix("`ID'") mata: Z = X mata: nvc = cols(Z)*(cols(Z)+1)/2 mata: init = J(1,nvc,0.5) mata: S = optimize_init() mata: optimize_init_evaluator(S, &remlprof_fn()) mata: optimize_init_evaluatortype(S, "v0") mata: optimize_init_argument(S, 1, X) mata: optimize_init_argument(S, 2, Z) mata: optimize_init_argument(S, 3, y) mata: optimize_init_argument(S, 4, C) mata: optimize_init_argument(S, 5, id) mata: optimize_init_params(S, init) mata: optimize_init_which(S,"max") quietly mata: parf = optimize(S) mata: est = glsfit_reml(X, Z, y, C, parf, id) mata: loglik = *est[8] mata: beta = *est[6]' mata: vcov = *est[7] mata: S = . mata: st_matrix("r(b)", beta) mata: st_matrix("r(V)", vcov) mata: st_numscalar("r(loglik)", loglik) mata: st_matrix("r(Psi)", par2Psi(parf, cols(Z))) mata: st_numscalar("r(aic)", -2*loglik+2*(cols(beta)+cols(parf))) mata: st_numscalar("r(nre)", cols(parf)) // get log-likelihood setting all var/cov equal to zero mata: parf = J(1, cols(parf), 0) mata: est = glsfit_reml(X, Z, y, C, parf, id) mata: st_numscalar("r(loglik0)", *est[8]) matrix `b' = r(b) matrix `V' = r(V) } tempname loglik loglik0 Psi aic nvcp k_f k scalar `loglik' = r(loglik) scalar `loglik0' = r(loglik0) matrix `Psi' = r(Psi) scalar `aic' = r(aic) scalar `nvcp' = r(nre) scalar `k_f' = colsof(`b') scalar `k' = `k_f' + `nvcp' // two-stage approach if ("`2stage'" != "") { qui tab `id' if `touse' local nstudies = r(r) if (`nstudies'==1) { di _n as err "more than one study is needed for the two-stage approach" exit 198 } // call mvmeta command quietly capture which mvmeta if _rc != 0 qui net install st0156_2, from(http://www.stata-journal.com/software/sj15-4) if "`acov'"=="" & "`mcov'"=="" & "`vcov'"=="" { qui mvmeta_make drmeta_cov `varlist' , s(`se') tot(`n') c(`case') ts(`typestudy') `crudes' `vwls' ecovmethod(`covmethod') /// saving(_glst_mvmeta_est) replace by(`id') names(b V) } if "`acov'"!="" { qui mvmeta_make drmeta_cov `varlist' , s(`se') ts(`typestudy') `vwls' meancov(`acov') /// saving(_glst_mvmeta_est) replace by(`id') names(b V) } // With multiple studies the user can specify the a list of covariance matrix following the order of the studies. local passm : word `countstudy' of `mcov' if "`mcov'"!="" { qui mvmeta_make drmeta_cov `varlist' , s(`se') ts(`typestudy') `vwls' matrixcov(`passm') /// saving(_glst_mvmeta_est) replace by(`id') names(b V) } if "`vcov'"!="" { qui mvmeta_make drmeta_cov `varlist' , s(`se') ts(`typestudy') `vwls' varlistcov(`vcov') /// saving(_glst_mvmeta_est) replace by(`id') names(b V) } preserve qui use _glst_mvmeta_est, clear if "`mvopt'" != "fixed" capture quietly mvmeta b V , bscovariance(unstructured) `mvopt' else capture quietly mvmeta b V , `mvopt' capture estimates save _glst_mvmeta_results , replace restore tempname getb getV getSigma nvcp Q Q_df Q_p mat `getSigma' = e(Sigma) matrix `Psi' = `getSigma' scalar `nvcp' = (colsof(`getSigma')*(colsof(`getSigma')-1))/2 + colsof(`getSigma') mat `getb' = e(b) if "`mvopt'" != "fixed" local pickeq "Overall_mean:" else local pickeq "Overall:" mat `b' = `getb'[1, "`pickeq'"] mat `getV' = e(V) matrix `V' = `getV'["`pickeq'", "`pickeq'"] local df_m = e(df_m) local df = `nstudies'-1 tempname loglik loglik0 aic k_f k if "`mvopt'" != "fixed" scalar `loglik' = e(ll) else scalar `loglik' = e(ll0) scalar `loglik0' = e(ll0) scalar `k_f' = e(fixedparms) scalar `k' = `k_f' + `nvcp' scalar `aic' = -2*`loglik' + 2*`k' scalar `Q' = e(Qscalar_chi2) scalar `Q_df' = e(Qscalar_df) scalar `Q_p' = chi2tail(`Q_df', `Q') } } // end multiple studies // Single study if "`id'" == "" { if "`acov'"=="" & "`mcov'"=="" & "`vcov'"=="" { drmeta_cov `varlist' `if' `in' , s(`se') tot(`n') c(`case') /// `pecohort' `ptcohort' `casecontrol' `meandiff' `stdmeandiff' `crudes' `vwls' ecovmethod(`covmethod') tempname A B mat `A' = e(cases) mat `B' = e(noncases) } if "`acov'"!="" { drmeta_cov `varlist' `if' `in' , s(`se') /// `crudes' `vwls' meancov(`acov') } if "`mcov'"!="" { drmeta_cov `varlist' `if' `in' , s(`se') /// `crudes' `vwls' matrixcov(`mcov') } if "`vcov'"!="" { drmeta_cov `varlist' `if' `in' , s(`se') /// `crudes' `vwls' varlistcov(`vcov') } tempname b V COV R k k_f k_r R2 D p_D local obs = e(N) mat `COV' = e(C) mat `R' = e(R) mat `b' = e(BC)' mat `V' = e(VB) local loglik = r(LL) local aic = -2*r(LL)+2*colsof(`b') local nstudies = 1 scalar `k' = colsof(`b') scalar `k_f' = colsof(`b') scalar `k_r' = 0 scalar `R2' = . scalar `D' = . scalar `p_D' = . } // End Single study // Saved results mat rownames `b' = `depname' mat colnames `b' = `xname' mat rownames `V' = `xname' mat colnames `V' = `xname' if "`mvopt'" != "fixed" { mat rownames `Psi' = `xname' mat colnames `Psi' = `xname' } mat coleq `b' = "" mat roweq `V' = "" mat coleq `V' = "" tempname _bcombined matrix `_bcombined' = `b' local df = `obs' - rowsof(`V') tempvar tousecopy gen `tousecopy' = `touse' ereturn post `b' `V' , dep(`depname') obs(`obs') esample(`touse') if "`id'" == "" { ereturn matrix Sigma = `COV' ereturn matrix R = `R' if "`acov'"=="" & "`mcov'"=="" & "`vcov'"=="" { ereturn matrix FC = `A' ereturn matrix FNC = `B' } } else { ereturn matrix Sigma = `G' } ereturn scalar k = `k' ereturn scalar k_f = `k_f' if "`id'" != "" & "`mvopt'" != "fixed" ereturn scalar k_r = `nvcp' else ereturn scalar k_r = 0 ereturn scalar aic = `aic' ereturn scalar ll = `loglik' if "`id'" != "" & "`mvopt'" != "fixed" ereturn scalar ll_c = `loglik0' else ereturn scalar ll_c = . ereturn local cmd = "drmeta" ereturn local depvar = "`depname'" ereturn local se = "`se'" if "`mvopt'" == "fixed" ereturn local method "gls" else ereturn local method "`mvopt'" ereturn local cmdline "drmeta `0'" ereturn local dm "`xname'" ereturn local type `random' ereturn local proc `stage' ereturn local predict "drmeta_predict" ereturn local id "`studies'" if "`id'" != "" { ereturn scalar r2 = `R2' ereturn scalar D = `D' ereturn scalar p_D = `p_D' ereturn matrix ty = `ty' ereturn matrix tX = `tX' ereturn matrix tres = `tres' ereturn matrix tpred = `tpred' } if ("`2stage'" != "") { ereturn scalar Q = `Q' ereturn scalar Q_df = `Q_df' ereturn scalar Q_p = `Q_p' } if "`e(type)'" == "random" { foreach study of local studies { matrix S`study' = _X`study'*matrix(`Psi')*_X`study'' + _C`study' matrix blup`study' = matrix(`Psi')*_X`study''*invsym(S`study')*(_B`study'-_X`study'*`_bcombined'') matrix blup`study' = blup`study'' mat rownames blup`study' = "blup_`study'" mat colnames blup`study' = `xname' matrix xbu`study' = `_bcombined' + blup`study' mat rownames xbu`study' = "xb_blup_`study'" mat colnames xbu`study' = `xname' ereturn matrix xbu`study' = xbu`study' ereturn matrix blup`study' = blup`study' ereturn matrix Sigma`study' = _C`study' ereturn matrix X`study' = _X`study' mat _BS`study' = _BS`study'' mat rownames _BS`study' = "gls_`study'" mat colnames _BS`study' = `xname' mat rownames _VS`study' = `xname' mat colnames _VS`study' = `xname' ereturn matrix bs`study' = _BS`study' ereturn matrix vs`study' = _VS`study' ereturn scalar r2_`study' = `R2_`study'' } } if "`e(type)'" != "random" { foreach study of local studies { ereturn matrix Sigma`study' = _C`study' ereturn matrix X`study' = _X`study' ereturn matrix bs`study' = _BS`study' ereturn matrix vs`study' = _VS`study' ereturn scalar r2_`study' = `R2_`study'' } } if ("`id'"!= "") { if "`mvopt'" != "fixed" { tempname PsiC copyPsi tagPsi copyPsiC matrix `copyPsi' = `Psi' matrix `tagPsi' = J(rowsof(`Psi'), rowsof(`Psi'),0) matrix `PsiC' = corr(`Psi') matrix `copyPsiC' = `PsiC' ereturn matrix Psi = `Psi' ereturn matrix PsiC = `PsiC' } ereturn local idname "`id'" } *ereturn scalar chi2_gf = r(Q) if ("`2stage'"!= "") & (`ncov'==1) ereturn scalar chi2_gf = `Q' // Q-test heterogeneity if `df' > 0 { if "`id'" == "" eret scalar chi2_gf = r(chi2_gf2) } else { eret scalar chi2_gf = . } /* compute test of indepvars = 0 */ qui testparm `xname' eret scalar df_m = r(df) eret scalar chi2 = r(chi2) eret scalar p = r(p) eret scalar N_s = `nstudies' sort `order' display_result, level(`level') `eform' `vwls' nv(`ncov') if "`noretable'" == "" & "`mvopt'" != "fixed" { di as txt "{hline 29}{c TT}{hline 15} di as txt _col(3) "Random-effects parameters" _col(30) "{c |}" _col(34) "Estimate" di as txt "{hline 29}{c +}{hline 15} * diagonal elements local i = 1 foreach v of local xname { if "`stddeviations'" == "" di as txt "var(" abbrev("`v'", 12) "," abbrev("`v'", 12) ")" _col(30) "{c |}" as res _col(34) %9.0g `copyPsi'[`i', `i'] else di as txt "std(" abbrev("`v'", 12) "," abbrev("`v'", 12) ")" _col(30) "{c |}" as res _col(34) %9.0g sqrt(`copyPsi'[`i', `i']) local i = `i' + 1 } * out-of-diagonal elements local i = 1 foreach v of local xname { local j = 1 foreach s of local xname { if (`i' != `j') & (`tagPsi'[`j', `i']!=1) & ("`stddeviations'" == "") di as txt "cov(" abbrev("`v'", 12) "," abbrev("`s'", 12) ")" _col(30) "{c |}" as res _col(34) %9.0g `copyPsi'[`i', `j'] if (`i' != `j') & (`tagPsi'[`j', `i']!=1) & ("`stddeviations'" != "") di as txt "corr(" abbrev("`v'", 12) "," abbrev("`v'", 12) ")" _col(30) "{c |}" as res _col(34) %9.0g `copyPsiC'[`i', `j'] mat `tagPsi'[`i', `j'] = 1 local j = `j' + 1 } local i = `i' + 1 } di as txt "{hline 29}{c BT}{hline 15} } if "`nolrt'" == "" & "`mvopt'" != "fixed" { di as txt "{help j_mixedlr##|_new:LR test} vs. no random-effects model = " as res 2*abs(`e(ll)'-`e(ll_c)') _col(55) as txt "Prob >= chi2(" as res `e(k_r)' as txt ") = " _col(73) as res %6.4f chi2tail(`e(k_r)', 2*abs(`e(ll)'-`e(ll_c)')) ereturn scalar lrt_c = 2*abs(`e(ll)'-`e(ll_c)') ereturn scalar df_c = `e(k_r)' ereturn scalar p_c = chi2tail(`e(k_r)', 2*abs(`e(ll)'-`e(ll_c)')) } /* if "`undetermined'" != "" { di as txt _n "{p 0 6 2}Note: The reported degrees of freedom " di as txt "assumes the null hypothesis is not on the " di as txt "boundary of the parameter space. If this is not " di as txt "true, then the reported test is " di as txt "{help j_mixedlr##|_new:conservative}.{p_end}" } */ // drop global used macro drop loglik levelci levelp capture drop `xnamec' sort `order' if "`id'" != "" { local i = 1 qui levelsof `id' , local(studies) foreach study of local studies { capture mat drop _VS`i' capture mat drop _BS`i' capture mat drop _ZS`i' capture mat drop _X`i' capture mat drop _C`i' capture mat drop _B`i' capture mat drop _ZS`i' capture mat drop _G`i' local i = `i' + 1 } } capture mat drop R capture mat drop C capture mat drop cases capture mat drop noncases end capture program drop display_result program define display_result version 13 syntax [, Level(integer $S_level) eform type(string) proc(string) stage vwls nv(string) ] if `level' < 10 | `level' > 99 { di in red "level() must be between 10 and 99 inclusive" exit 198 } tempname pgf pm if e(df_gf) > 0 { scalar `pgf' = chiprob(e(df_gf), e(chi2_gf)) } else scalar `pgf' = . scalar `pm' = chiprob(e(df_m), e(chi2)) di _n as txt "`e(proc)' `e(type)'-effect dose-response model" _col(50) "Number of studies = " _col(72) as res %7.0g e(N_s) di as txt "Optimization = " as res "`e(method)'" _col(54) as txt "Number of obs = " _col(72) as res %7.0g e(N) di as txt " AIC = " as res %5.2f `e(aic)' _col(54) as txt "Model chi2(" as res e(df_m) as txt ") = " _col(72) as res %7.2f e(chi2) di as txt "Log likelihood = " as res e(ll) _col(56) as txt "Prob > chi2 = " _col(72) as res %7.4f `pm' if "`eform'" == "" ereturn display, level(`level') else ereturn display, level(`level') eform(exb(b)) end mata function blockaccum(string scalar a, string scalar b) { real matrix A, B, GA A = st_matrix(a) B = st_matrix(b) GA = blockdiag(A,B) st_matrix("r(CACC)", GA) } function glsest(string scalar y, string scalar x, string scalar cov) { real matrix Y, X, C, BC, VB real scalar Q, LL Y = st_matrix(y) X = st_matrix(x) C = st_matrix(cov) // Variance and betas VB = invsym(X'*invsym(C)*X) BC = invsym(X'*invsym(C)*X)*X'*invsym(C)*Y // Q goodness of fit Q = (Y-X*BC)'*invsym(C)*(Y-X*BC) C2 = diag(C) // log-likelihood LL = -.5*rows(Y)*log(2*pi())-.5*log(det(C))-.5*Q st_matrix("r(BC)", BC) st_matrix("r(VB)", VB) st_numscalar("r(Q)", Q) st_numscalar("r(LL)", LL) st_matrix("r(COVAR)", C) } real matrix par2Psi(par, q){ L = J(q,q,0) T = lowertriangle(J(q,q,.)):!=0 k = 1 for (i=1; i<=q; i++) { for (j=1; j<=q; j++) { if (T[i,j]==1) { L[i,j] = par[1,k] k = k + 1 } } } L=lowertriangle(L) return(L*L') } function glsfit_fixed(X, Z, y, C, par, id){ pointer vector info n = rows(uniqrows(id)) q = cols(Z) nall = rows(y) Psi = par2Psi(par,q) idn = uniqrows(id) for (i=1; i<=n; i++) { Zi = select(Z, id:==idn[i]) Si = select(select(C, (id:==idn[i])), (id:==idn[i])') Sigmalisti = makesymmetric(Si + Zi*Psi*Zi') Ulist = cholesky(Sigmalisti)' invUlist = luinv(Ulist) yi = select(y, id:==idn[i]) invtUylist = invUlist'*yi Xi = select(X, id:==idn[i]) invtUXlist = invUlist'*Xi tXWXlist = invtUXlist'*invtUXlist if (i==1) { invtUy = invtUylist invtUX = invtUXlist diagUlist = ln(diagonal(Ulist)) tXWXtot = tXWXlist } else { invtUy = invtUy \ invtUylist invtUX = invtUX \ invtUXlist diagUlist = diagUlist \ ln(diagonal(Ulist)) tXWXtot = tXWXtot :+ tXWXlist } } coef = qrsolve(invtUX, invtUy) Q = R = . qrd(invtUX, Q, R) R= R[|1,1\q,q|] vcov = solveupper(R, diag(J(q,q,1))) * solveupper(R, diag(J(q,q,1)))' pconst = -0.5 * nall * log(2 * pi()) pres = -0.5*(invtUy - invtUX*coef)'*(invtUy - invtUX*coef) pdet1 = -sum(diagUlist) pdet2 = -sum(log(diagonal(cholesky(tXWXtot)))) lf = (pconst + pdet1 + pres) Psip = &Psi invtUXp = &invtUX invtUyp = &invtUy diagUlistp = &diagUlist tXWXtotp = &tXWXtot coefp = &coef vcovp = &vcov lfp = &lf info = (invtUXp, invtUyp, diagUlistp, tXWXtotp, Psip, coefp, vcovp, lfp) return(info) } function glsfit_reml(X, Z, y, C, par, id){ pointer vector info n = rows(uniqrows(id)) q = cols(Z) nall = rows(y) Psi = par2Psi(par,q) idn = uniqrows(id) for (i=1; i<=n; i++) { Zi = select(Z, id:==idn[i]) Si = select(select(C, (id:==idn[i])), (id:==idn[i])') Sigmalisti = makesymmetric(Si + Zi*Psi*Zi') Ulist = cholesky(Sigmalisti)' invUlist = luinv(Ulist) yi = select(y, id:==idn[i]) invtUylist = invUlist'*yi Xi = select(X, id:==idn[i]) invtUXlist = invUlist'*Xi tXWXlist = invtUXlist'*invtUXlist if (i==1) { invtUy = invtUylist invtUX = invtUXlist diagUlist = ln(diagonal(Ulist)) tXWXtot = tXWXlist } else { invtUy = invtUy \ invtUylist invtUX = invtUX \ invtUXlist diagUlist = diagUlist \ ln(diagonal(Ulist)) tXWXtot = tXWXtot :+ tXWXlist } } coef = qrsolve(invtUX, invtUy) Q = R = . qrd(invtUX, Q, R) R= R[|1,1\q,q|] vcov = solveupper(R, diag(J(q,q,1))) * solveupper(R, diag(J(q,q,1)))' pconst = -0.5 * (nall - q) * log(2 * pi()) pres = -0.5*(invtUy - invtUX*coef)'*(invtUy - invtUX*coef) pdet1 = -colsum(diagUlist) pdet2 = -colsum(log(diagonal(cholesky(tXWXtot)))) lf = (pconst + pdet1 + pdet2 + pres) Psip = &Psi invtUXp = &invtUX invtUyp = &invtUy diagUlistp = &diagUlist tXWXtotp = &tXWXtot coefp = &coef vcovp = &vcov lfp = &lf info = (invtUXp, invtUyp, diagUlistp, tXWXtotp, Psip, coefp, vcovp, lfp) return(info) } function glsfit_ml(X, Z, y, C, par, id){ pointer vector info n = rows(uniqrows(id)) idn = uniqrows(id) q = cols(Z) nall = rows(y) Psi = par2Psi(par,q) for (i=1; i<=n; i++) { Zi = select(Z, id:==idn[i]) Si = select(select(C, (id:==idn[i])), (id:==idn[i])') Sigmalisti = makesymmetric(Si + Zi*Psi*Zi') Ulist = cholesky(Sigmalisti)' invUlist = luinv(Ulist) yi = select(y, id:==idn[i]) invtUylist = invUlist'*yi Xi = select(X, id:==idn[i]) invtUXlist = invUlist'*Xi tXWXlist = invtUXlist'*invtUXlist if (i==1) { invtUy = invtUylist invtUX = invtUXlist diagUlist = ln(diagonal(Ulist)) tXWXtot = tXWXlist } else { invtUy = invtUy \ invtUylist invtUX = invtUX \ invtUXlist diagUlist = diagUlist \ ln(diagonal(Ulist)) tXWXtot = tXWXtot :+ tXWXlist } } coef = qrsolve(invtUX, invtUy) Q = R = . qrd(invtUX, Q, R) R= R[|1,1\q,q|] vcov = solveupper(R, diag(J(q,q,1))) * solveupper(R, diag(J(q,q,1)))' pconst = -0.5 * nall * log(2 * pi()) pres = -0.5*(invtUy - invtUX*coef)'*(invtUy - invtUX*coef) pdet1 = -sum(diagUlist) pdet2 = -sum(log(diagonal(cholesky(tXWXtot)))) lf = (pconst + pdet1 + pres) Psip = &Psi invtUXp = &invtUX invtUyp = &invtUy diagUlistp = &diagUlist tXWXtotp = &tXWXtot coefp = &coef vcovp = &vcov lfp = &lf info = (invtUXp, invtUyp, diagUlistp, tXWXtotp, Psip, coefp, vcovp, lfp) return(info) } void remlprof_fn(todo, par, X, Z, y, C, id, lf, g, H){ q = cols(Z) nall = rows(y) Psi = par2Psi(par,q) get_gls = glsfit_reml(X, Z, y, C, par, id) lf = *get_gls[8] } void mlprof_fn(todo, par, X, Z, y, C, id, lf, g, H){ q = cols(Z) nall = rows(y) Psi = par2Psi(par,q) get_gls = glsfit_ml(X, Z, y, C, par, id) lf = *get_gls[8] } end