*! v2.6.1 Roger Harbord 4 Nov 2008
	
program define metareg, eclass byable(recall)
version 7
	if replay() {
		if "`e(cmd)'" !="metareg" {
			error 301 /* last estimates not found */
			}
		if _by() {
			error 190 /* request may not be combined with by */
			}
		Display `0'
		exit
		}

	syntax varlist(min=1 numeric) [if] [in] , [ /*
*/	  wsse(varname numeric) /*
*/   EForm                 /*
*/   NOCONStant            /*
*/   MM                    /*
*/   REML                  /*
*/   EB                    /*
*/   Knapphartung          /*
*/   z                     /* standard normal w/o K-H modification to variance
*/	  Level(passthru)       /*
*/	  PERMute(string)       /*
*/   Graph                 /* produce a bubble plot
*/   RAndomsize            /* marker size according to random-effects weights
*/   TAU2test              /* tests for tau2=0 residual het test
*/   LOg                   /* maximize_option
*/   TOLerance(real 1e-6)  /* maximize_option also used by eb method
*/   ITERate(integer 100)  /* maximize_option - default 16000 is too large!
*/   wsvar(varname numeric)/* undoc'd, for compat with v1 ( wsvar = wsse^2 )
*/   bsest(string)         /* undoc'd, for compat with v1
*/   tdist                 /* undoc'd, use t-distib w/o K-H variance mod'n
*/   ORIGinal              /* undoc'd, use version 1 (metareg_orig)
*/	  LRtau2                /* undoc'd, for compat with v.2.2 on SSC
*/   NOITer                /* passed through to metareg_orig
*/   NOTAUComp             /* no comparison of tau2 with constant-only model;
		                        for use by permute() option
*/   *                     /* other `options', parsed by -mlopts-
*/   ]

	marksample touse
	tokenize `varlist'
	local y `1'	/* study estimate */
	qui count if `touse'
	local nobs = r(N)
	if r(N) == 0 {
		error 2000  /* no observations */
		}
	if r(N) <= `: word count `*'' {
		error 2001 /* insufficient observations */
		}
	macro shift 1
	_rmcoll `*' if `touse', `noconstant'
	local xvars `r(varlist)'

	mlopts mlopts, `options' iterate(`iterate')  /* Parse  maximize_options */

	if "`wsse'"=="" & "`wsvar'"==""  {
		di as error "wsse() option required"
		exit 100
		}

	if "`wsse'"!="" & "`wsvar'"!="" {
		di as error "options wsse() and wsvar() cannot both be specified"
		exit 198
		}
	
	if "`wsse'"!="" {
		capture assert `wsse' > 0 if `touse'
		if _rc {
			di as error "zero or negative wsse() not allowed"
			exit 499
			}
		tempvar wsvar
		qui gen `wsvar' = `wsse'^2 if `touse'
		}

	if "`wsvar'"!="" { 
		capture assert `wsvar' > 0 if `touse'
		if _rc {
			di as error "zero or negative wsvar() not allowed"
			exit 499
			}
		}

	if "`noconstant'"!="" {
		/* 		local nvar : word count `varlist' */
		if `: word count `varlist'' == 1 {
			di as error "independent variables required with "  /*
*/			  "noconstant option"
			exit 102
			}
		}
	
	local wc : word count "`mm' `reml' `eb'"
	if `wc' > 1 {
		di as error "Only one of options mm, reml and eb can be specified"
		exit 198
		}
		
	local wc : word count "`mm' `reml' `eb' `bsest'"
	if `wc' > 1 {
		di as error "Option bsest() cannot be used with options  mm, reml or eb"
		exit 198
		}
		
	if "`bsest'"!="reml" & "`bsest'"!="ml" & "`bsest'"!="eb"  /*
*/	  & "`bsest'"!="mm" & "`bsest'"!="" {
		di as error "Between study variance estimation method invalid:"
		di as error  "should be either reml, ml, eb or mm"
		exit 198
		}

	if "`bsest'" == "" {
		if "`mm'" != ""  | "`permute'" != "" { /* MM is default with permute */
			local bsest "mm"
			}
		else if "`eb'" != "" {
			local bsest "eb"
			}
		else  {                               /* REML is default otherwise */
			local bsest "reml"
			}
		}

	if "`bsest'" =="ml" {
		di as text  /*
*/		  "Option bsest(ml) not supported by version 2 of metareg."
		di as text "Calling version 1.06 (sbe23):" _n(1)
		metareg_orig `0'
		exit
		}

	if "`knapphartung'" !="" & "`z'" != "" {
		di as error "Options knapphartung and z cannot both be specified"
		exit 198
		}

	if "`tdist'" != "" & ( "`z'" != "" | "`knapphartung'" !="" ) {
		di as error "Option t cannot be specified with options knapphartung or z"
		exit 198
		}

	if "`original'" =="original" {
		di as text "Calling metareg version 1.06 (sbe23):" _n(1)
		metareg_orig `0'
		exit
		}
	
	if "`log'"=="" { 
		local log "nolog"
		}

	if "`lrtau2'" != "" {
		local tau2test tau2test
		}
	
	if "`permute'" != ""  {
		if "`xvars'" == "" {
			di as error "Option permute() requires at least one covariate"
			exit 198
			}
		Permute, y(`y') xvars(`xvars') wsvar(`wsvar') permute(`permute') ///
		  nobs(`nobs') `noconstant'  ///
		  bsest(`bsest') `knapphartung' `z'  ///
		  level(`level') touse(`touse')
		exit
		}


	/*** begin main code ***/
	tempvar w hat v resid
	tempname remll_c remll sumw tau2 tau2mm b V Q df_Q tau2_0

	quietly {
		
		/** (get starting values using) method-of-moments **/
		gen `w'=1/`wsvar' if `touse'
		summ `w' if `touse', meanonly
		scalar `sumw'=r(sum)
		/* use regress not vwls so can predict hat */
		regress `y' `xvars'  [iw=`w'] if `touse', `noconstant'
		predict `hat' if `touse', hat
		/* H_ii where H = X(X'WX)^{-1}X', W=diag(w) */
		replace `hat'=`hat'*`w'^2  if `touse' /* Tr( WX(X'WX)^{-1}X'W ) */
		/* iweights mess up residual df so compute it  : */
		scalar `df_Q' = `nobs' - e(df_m) - ("`noconstant'" == "")
		scalar `Q' = e(rss)
		summ `hat' if `touse', meanonly
		scalar `tau2' = (`Q' - `df_Q') / ( `sumw' - r(sum))
		/* truncate at zero */
		scalar `tau2' = max(`tau2', 0)

		/* Fit constant-only model to compare tau2
		NB can't compare different fixed-effect models using REML log-L */
		if "`noconstant'" == "" & "`notaucomp'"=="" & "`xvars'" != "" {
			metareg `y' if `touse', wsvar(`wsvar') bsest(`bsest') notaucomp
			scalar `tau2_0' = e(tau2)
			}
		} /* end quietly */
	
	if "`bsest'" == "reml" {  /* *** REML start *** */
			
		global REML_x `xvars' /* betas not estimated by ml so pass X as global*/
		global REML_nocons `noconstant' /* have to pass as global too */

			/* fit comparison fixed-effect (tau2=0) model*/
			ml model d0 metareg_ll (`y' `wsvar' = ) /*
*/		     if `touse', maximize init(_cons=0) /*
*/		     search(off) iterate(0) nowarning nolog
			scalar `remll_c' = e(ll) /* log-L, comparison model */
		
		/* estimate tau2 using ml */
		local tau2mm = max(`tau2', 1e-4)  /* init doesn't work with scalar */
		ml model d0 metareg_ll (`y' `wsvar' = ) if `touse',/*
*/			  maximize  init(_cons=`tau2mm') search(off) nopreserve noscvars /*
*/			  `log' tolerance(`tolerance') `mlopts'

		if "`e(converged)'" != "1" {
			di as error "WARNING: REML estimation of tau2 failed to converge."
			di as error "Estimates may be wrong."
			di as error "Try adding mm option to use method-of-moments".
			}
		
		macro drop REML_x
		macro drop REML_nocons
		
		scalar `tau2' = max( _b[_cons], 0)
		if `tau2' > 0 {
			scalar `remll' = e(ll)
			}
		else {
			scalar `remll' = `remll_c'
			}
			
		}                   /*  REML end */
		
	if "`bsest'" == "eb" {  /* *** EB start *** */
		local oldtau2=-1
		local j=0

		while abs((`tau2'-`oldtau2')/(`oldtau2'+1e-8)) >= `tolerance' {
			local j = `j'+1
			tempvar wt ypred wtsq numi 
			local oldtau2=`tau2'
			qui {
				gen `wt'=(`wsvar'+`tau2')^-1 if `touse'
				regress `y' `xvars' [aw=`wt'] if `touse'
				predict `ypred' if `touse'
				gen `wtsq'=`wt'^2 if `touse'
				gen `numi'=(`wt')*(((e(N)/e(df_r))*((`y'-`ypred')^2))-`wsvar') if `touse'
				summ `numi' if `touse', meanonly
				local num=r(sum)
				summ `wt', meanonly
				local tau2=`num'/r(sum)
				}
			if "`log'" == "log" {
				di as txt  "Iteration " `j' ": tau^2 = " as res `tau2'
				}
			}
		local tau2 = max(`tau2',0)
		}                  /* EB end*/
	
	/** estimate final model by weighted LS **/
	quietly {
		gen `v' = `wsvar' + `tau2'
		regress `y' `xvars' [iw=1/`v'] if `touse', mse1 `noconstant'
		matrix `b' = e(b)
		matrix `V' = e(V)
		local df_m = e(df_m)

		/** knapp and hartung variance adjustment **/
		if "`z'" == "" & "`tdist'" == "" {
			tempname vkh
			predict `resid' if `touse', resid
			replace `resid' = `resid'*`resid' / `v' if `touse'
			summarize `resid' if `touse', meanonly
			local qkh =  r(sum) / `df_Q'
			matrix `V' = `V' * max(1, `qkh')
			}
		
		} /* end quietly */

	
	estimates post `b' `V', depname(`y') obs(`nobs') esample(`touse')
	estimates scalar tau2 = `tau2'
	estimates scalar df_m = `df_m'
	estimates scalar Q = `Q'
	estimates scalar df_Q = `df_Q'
	estimates scalar I2 = max( ( `Q' - `df_Q' ) / `Q', 0 )
	if "`z'" == ""  { /* post df_r so get t & F tests not z and chi2 */
		estimates scalar df_r = `df_Q'
		}
	if "`tdist'" =="" & "`z'" == "" {
		estimates scalar q_KH = `qkh'
		}
	if "`wsse'" == "" {
		estimates local wsvar "`wsvar'"
		}
	else {
		estimates local wsse "`wsse'"
		}

	if "`bsest'" == "mm" {
		estimates local method "Method of moments"
		}
	if "`bsest'" == "eb" {
		estimates local method "Empirical Bayes"
		}
	if "`bsest'" == "reml" {
		estimates local method "REML"
		/* return reml log-l in e(remll) NOT e(ll)
		to ensure can't do inappropriate lrtest of fixed effects */
		estimates scalar remll = `remll'
		if "`remll_c'" != "" {
			estimates scalar remll_c = `remll_c'
			if (e(remll) < e(remll_c)) | (`tau2' <= 0) {
				estimates scalar chi2_c = 0
				}
			else {
				estimates scalar chi2_c = 2*(e(remll)-e(remll_c))
				}
			}
		}

	if "`noconstant'" == "" & "`notaucomp'"=="" & "`xvars'" != "" {
		estimates scalar tau2_0 = `tau2_0'
		}

	/* overall F (or chi2) test */
	if "`xvars'" != "" {
		qui test `xvars'
		if "`r(chi2)'" != "" {
			estimates scalar chi2 = r(chi2)
			}
		else {
			estimates scalar F = r(F)
			}
		}
	
	estimates local predict "metareg_p"
	estimates local depvar "`y'"
	estimates local cmd "metareg"
	
	Display, `level' `eform' `tau2test'


	if "`graph'" != "" {
		if `: word count `xvars'' != 1 {
			di as error in smcl  ///
			  "graph option works only with a single {it:indepvar}"
			}
		else Graph `y' `xvars' `v' `w' `randomsize'
		}

end

program define Graph
	version 8.2
	args y xvars v w randomsize
	
	tempvar fit
	qui predict `fit'
	if "`randomsize'" == "" {
		local size `w'
		}
	else {
		local size 1/`v'
		}
	local yti : variable label `y'
	if "`yti'" == "" {
		local yti `y'
		}
	scatter `y' `xvars' [aw=`size'], msymbol(oh)  ///
	  || line `fit' `xvars', sort ///
	  ||, legend(off) ytitle( `yti')
	
end


program define Display
version 7
	syntax [, Level(int $S_level) EForm Tau2test q]
	tempname pval i2 
	if "`eform'" == "eform" {
		local eform eform(exp(b))
		}
	scalar `pval' = chi2tail(e(df_Q), e(Q))

	di as txt _n "Meta-regression"  /*
*/  _col(55) "Number of obs" _col(70) "= " as res %7.0f e(N)
	di as res "`e(method)'"  /*
*/	   as txt " estimate of between-study variance" /*
*/ _col(55) "tau2" _col(70) "=" as res %8.4g e(tau2)
	di as txt "% residual variation due to heterogeneity"  /*
*/	  _col(55) "I-squared_res" _col(70) "= "  as res %6.2f 100 * e(I2) "%" 


	if "`e(tau2_0)'" != "" {
		di as txt "Proportion of between-study variance explained" /*
*/    _col(55) "Adj R-squared "  _col(70) "= " _c
		di as res %6.2f 100 * ( 1 - e(tau2) / e(tau2_0)  ) "%"
		}
	
	/* Overall model fit */
	if e(df_m) > 1 {
		di as txt "Joint test for all covariates" _c
		if "`e(df_r)'" == "" {
			di as txt _col(55) "Model chi2(" as res e(df_m) as txt ")" /*
*/          _col(70) "= " as res %7.2f e(chi2)
			}
		else {
			di as txt _col(55) "Model F(" as res e(df_m) as txt "," /*
*/       as res e(df_r) as txt ")" _col(70) "= " /*
*/       as res %7.2f e(F)
			}
		}

	
	di  as res "With" _c
	if "`e(q_KH)'" == "" {
		di as res "out" _c
		}
	di as txt " Knapp-Hartung modification" _c

	if e(df_m) > 1 {
		if "`e(df_r)'" == "" {
			di as txt _col(55) "Prob > chi2" _col(70) "= " /*
*/			as res %7.4f chi2tail( e(df_m), e(chi2) )
			}
		else {
			di as txt _col(55) "Prob > F" _col(70) "= " /*
*/       as res %7.4f Ftail( e(df_m), e(df_r), e(F) )
			}
		}
	else {
		di /* end line */
		}
	

	estimates display, level(`level') `eform'

	if "`tau2test'" != "" {
		di as txt "Test for residual between-study variance (of tau2=0)"  /*
*/	    _col(55) "Q_res (" as res %1.0f e(df_Q) as txt " df)"  /*
*/	    _col(70) "= " as res %7.2f e(Q) 
		di as txt _col(55) "Prob > Q_res" _col(70) "= " as res %7.4f `pval'

		if "`e(method)'" == "REML"  {
			scalar `pval' =  chi2tail(1, e(chi2_c))*0.5
			if e(chi2_c) <=0 {
				scalar `pval'= 1
				}
			di in smcl as txt "Likelihood-ratio test of tau2=0: " /*
			*/ as txt "{help j_chibar:chibar2(01) =}" as res %6.2f e(chi2_c)  /*
			*/ as txt "  Prob > chibar2 = " as res %7.4f `pval'
			}

		}

end


program define Permute, rclass
version 7
	syntax, y(varname numeric) xvars(varlist numeric) wsvar(varname numeric) ///
	  permute(string) nobs(integer) [noconstant bsest(string) knapphartung z ///
	  level(passthru) touse(varname numeric)] 
	
	gettoken reps 0 : permute, parse(",")
	confirm integer number `reps'
	if `reps' < 1 {
		di as err "permute() must be a positive integer"
		exit 198
		}
		
	syntax [, Univariable  Detail Verbose  NOISily Joint(string) * ]

	if "`verbose'" != "" { /* verbose is undoc'd synonym for detail */
		local detail detail
		}

	if "`xvars'" == "" {
		di as error "permute() option not allowed when no covariates"
		exit 198
		}

	if "`joint'" != "" & "`univariable'" != "" {
		di as error "options univariable and joint() cannot both be specified"
		exit 198
		}
		
	if `:word count `xvars'' == 1 {
		local univariable /* ensure unset if only one xvar */
		}
	
	if "`univariable'" == "" {
		di as txt _n(1) "Monte Carlo permutation test for meta-regression"
		}
	else {
		di as txt _n(1) in smcl /*
*/			  "Monte Carlo permutation test for {it:single covariate} meta-regressions"
		}

	if "`bsest'" == "" {
		local bsest "mm" /* mm is default for permutation test */
		}
	if "`bsest'" == "mm" {
		di as res _n "Moment-based " _c 
		}
	else {
		di as res _n "REML " _c 
		}
	di as txt "estimate of between-study variance"

	if "`knapphartung'" == ""  {
		local z "z"  /* z is default for for permutation test  */
		di as res "Without " _c
		}
	else {
		di as res "With " _c
		}
	di as txt "Knapp & Hartung modification to standard errors"

		
	if "`noconstant'" != "" {
		di as error "Permutation test inappropriate with no constant"
		exit 198
		}

	local permopts `univariable' `noconstant' `bsest' `knapphartung' `z'
	
	preserve
	qui keep if `touse'

	foreach x of varlist `xvars' {
		local explist "`explist' `x'=r(z_`x')"
		local xlist "`xlist' `x'"
		}
		
	/* parse varlists in joint() option separated by \ / or | */
	local jt `joint'
	local nj 0
	while "`jt'" != "" {
		gettoken tmp jt : jt, parse("\/|")
		if strpos("\/|", "`tmp'")  {
			continue
			}
		local 0 `tmp'
		syntax varlist
		local nj = `nj' + 1
		local joint`nj' `varlist'
		di as txt "joint`nj' : " as res "`joint`nj''"
		local explistj "`explistj' joint`nj'=r(chi2`nj')"
		local jlist "`jlist' joint`nj'"
		}
		
	tempvar n
	qui gen `n' = _n if `touse'
		
	tempfile pm pmadj
	if "`noisily'" =="" {
		local quietly quietly
		/* quietly permute then replay (if detail specified)
		to suppress full -permute- header */
		}
	`quietly' permute `n' /* 
*/		  "metareg_pm `y' `wsvar' `n', xvars(`xvars') joint(`joint') `permopts'" /*
*/		  `explist' `explistj', /*
*/		reps(`reps') saving(`pm') replace `options'


	tempname R C P B
	matrix `R' = r(reps)
	matrix `C' = r(c)
	matrix `P' = `C' * inv(diag(`R'))

	local nx : word count `xvars'
	forvalues i = 1/`nx' {
		local minreps = min(`reps',`R'[1,`i'])
		}
	if `minreps' != `reps' {
		di as txt "WARNING: some permutations returned missing values" _c
		di as txt "See detail below:"
		local detail detail
		}

	/* to get adjusted p-values from -permute-,
		replace z-statistics by their largest absolute values
	for each permutation */
	quietly {
		use `pm', clear
		tempvar max
		gen `max' = 0
		foreach v of varlist `xlist' {
			replace `max' = abs(`v') if abs(`v') > `max'
			}
		foreach v of varlist `xlist' {
			replace `v' = `max'
			}
		drop `max'
		save `pmadj', replace
		permute `xlist' `jlist' using `pmadj', `level'
		} // end quietly

	matrix `B' = r(b)
	matrix `R' = `R' \ r(reps)
	matrix `C' = `C' \ r(c)
	matrix `P' = `P' \ r(c) * inv(diag(r(reps)))
	foreach M in `R' `C' `P' {
		matrix rownames `M' = unadj multadj
		if "`joint'" != "" {
			/* set matrix entries for adjusted joint tests (N/A) to missing */
			matrix `M'[2, colnumb(`M',"joint1")] = J(1,`nj',.)
			}
		}

	/* clear results from last permutation so not available to user */
	estimates clear
	restore
		
	/*** Display results of permutation test ****/

	if `nx' == 1 { /* only one xvar */
		local sep = sqrt(`P'[1,1]*(1-`P'[1,1])/`reps')
		if `sep' == 0 {
			local sep .
			}
		di as txt _n _col(9) "Number of obs =" as res %8.0f _N
		di as txt _col(9) "Permutations  =" as res %8.0f `reps'
		di in smcl as txt "{hline 13}{c TT}{hline 17}"
		di in smcl as txt  /*
*/			  %12s "`y'"  " {c |}" _s(1)  /*
*/			  %6s "p" _s(4) /*
*/			  %6s "SE(p)"
			
		di in smcl as txt "{hline 13}{c +}{hline 17}"
		di in smcl as txt /*
*/				  %12s abbrev(`"`xvars'"',10) /*
*/				  " {c |}" _s(2)  /*
*/				  as result %7.3f `P'[1,1] _s(1) /*
*/            as result %7.4f `sep'
		di in smcl as txt "{hline 13}{c BT}{hline 17}"
		}

	else {   /* more than one xvar */

		/* Table Head */
		di as txt _n "P-values unadjusted and adjusted for multiple testing"
		di as txt _n _col(14) "Number of obs =" as res %8.0f `nobs'
		di as txt _col(14) "Permutations  =" as res %8.0f `reps'
		di in smcl as txt "{hline 13}{c TT}{hline 22}"
		di in smcl as txt _s(13) "{c |}" _s(12) "P"
		di in smcl as txt  /*
*/			  %12s abbrev("`y'",10)  " {c |}"  /*
*/			  %11s "Unadjusted"  %11s "Adjusted"
		di in smcl as txt "{hline 13}{c +}{hline 22}"

		foreach x of local xvars {
			local i = `i' + 1
			di in smcl as txt /*
*/				  %12s abbrev(`"`x'"',10) /*
*/				  " {c |}" _s(4)  /*
*/				  as result %7.3f `P'[1,`i'] _s(4)  /*
*/				  %7.3f `P'[2,`i']
			}

		if "`joint'" != "" {
			di in smcl as txt "{hline 13}{c +}{hline 22}"
			forvalues j = 1/`nj' {
				di in smcl as txt /*
*/               %12s "joint`j'"  /*
*/				  " {c |}" _s(4)  /*
*/				  as result %7.3f `P'[1,`nx'+`j']
				}
			}
			
		/* Table Foot */
		di in smcl as txt "{hline 13}{c BT}{hline 22}"
		/* largest SE(p)  = sqrt( max( p(1-p)/n ) ) */
		local maxvarp 0
		forvalues i = 1/2 {
			forvalues j = 1/`=colsof(`P')'{
				local maxvarp = max( `maxvarp',  ///
				  `P'[`i',`j'] * (1 - `P'[`i',`j']) / `R'[`i',`j'] )
				}
			}
		if `maxvarp' == 0 {
			local maxvarp .
			}
		di as txt "largest Monte Carlo SE(P) ="  /*
*/			  as res %7.4f sqrt(`maxvarp')
		}

	/* Detailed results if requested, as displayed by -permute- */
	if "`detail'" != "" {
		di as txt in smcl _n(2)  "{title:Unadjusted}" 
		permute `xlist' `jlist' using `pm', `level'

		if `nx' > 1 {
			di as txt in smcl _n(2)  "{title:Adjusted for multiple testing}" 
			permute `xlist' using `pmadj', `level'
			}
		}
		
	di as txt _n "WARNING:"
	di as txt "Monte Carlo methods use random numbers, so results may differ between runs."
	di as txt "Ensure you specify enough permutations to obtain the desired precision."

	return scalar N = `nobs'
	return matrix reps = `R'
	return matrix c = `C'
	return matrix p = `P'
	return matrix b = `B'

end

	

program define metareg_orig
*! v1.06 copyright Stephen Sharp January 1998  STB-42 sbe23
version 5.0
	local varlist "req ex min(1)"
	local if "opt"
	local in "opt"
	local options "WSSe(string) WSVar(string) BSest(string) TOLeran(integer 4) Level(integer $S_level)"
	local options "`options' NOITer ORIGinal" /* original added by RMH */

	parse "`*'"
	if "`wsse'"=="" & "`wsvar'"=="" {
		di in re "Must specify a variable containing estimate of precision"
		di in re "within each trial, using either wsse() or wsvar() option"
		exit 198
		}
	if "`wsse'"~="" & "`wsvar'"~=""{
		confirm variable `wsse'
		confirm variable `wsvar'
		local i=1
		while `i'<=_N {
			if abs(`wsse[`i']'^2-`wsvar[`i']')>0.00001 {
				di in re "Within study variance should be square of within study standard error"
				exit 198
				}
			local i=`i'+1
			}
		}
	if "`wsvar'"~="" {
		confirm variable `wsvar'
		}	
	if "`wsse'"~="" {
		confirm variable `wsse'
		tempvar wsvar
		qui gen `wsvar'=`wsse'^2
		}	
	if "`bsest'"=="" {
		local bsest "reml"
		}
	if "`bsest'"~="reml" & "`bsest'"~="ml" & "`bsest'"~="eb" & "`bsest'"~="mm" & "`bsest'"~="" {
		di in re "Between study variance estimation method invalid:"
		di in re "should be either reml, ml, eb or mm"
		exit 198
		}
	if "`noiter'"~="" & "`bsest'"=="mm" {
		di in bl "Warning: mm is a non-iterative method, noiter option ignored"
		}	
	parse "`varlist'", parse(" ")

	local y "`1'"
	mac shift
	local xvars "`*'"

	tempvar touse
	preserve
	qui {
		mark `touse' `if' `in'
		markout `touse' `y' `xvars'
		keep if `touse'
		}

	qui regress `y' `xvars'
	local p=_result(3)
	local N=_result(1)

	if "`bsest'"=="mm" {
		tempvar wt ypred numi one
		qui {
			gen `wt'=`wsvar'^-1
			regress `y' `xvars' [aw=`wt']
			predict `ypred'
			gen `one'=1
			tempname X C Xt C1 XtC1 A A1 A1XtC1 C1X B
			mkmat `one' `xvars', matrix(`X')
			matrix `C'=J(`N',`N',0)
			local i=1
			while `i'<=_N {
				matrix `C'[`i',`i']=`wsvar'[`i']
				local i=`i'+1
				}
			mat `Xt'=`X''
			mat `C1'=inv(`C')
			mat `XtC1'=`Xt'*`C1'
			mat `A'=`XtC1'*`X'
			mat `A1'=inv(`A')
			mat `A1XtC1'=`A1'*`XtC1'
			mat `C1X'=`C1'*`X'
			mat `B'=`C1X'*`A1XtC1'
			local trB=trace(`B')
			summ `wt'
			local denom=_result(18)-`trB'
			gen `numi'=`wt'*((`y'-`ypred')^2)
			summ `numi'
			local num=max(_result(18)-(`N'-(`p'+1)),0)
			local newtsq=`num'/`denom'
			}	
		}

	if "`bsest'"~="mm" {

		local tsq=0.1
		local newtsq=0
		local j=1

		while abs(`tsq'-`newtsq')>=10^(-`toleran') {
			tempvar wt ypred wtsq numi 
			local tsq=`newtsq'
			if "`noiter'"=="" {
				di in gr  "Iteration " `j' ": tau^2 = " in ye `tsq'
				}
			qui {
				gen `wt'=(`wsvar'+`tsq')^-1
				regress `y' `xvars' [aw=`wt']
				predict `ypred'
				gen `wtsq'=`wt'^2
				}
			if "`bsest'"=="reml" {
				qui {
					gen `numi'=(`wtsq')*(((`N'/(`N'-(`p'+1)))*((`y'-`ypred')^2))-`wsvar')
					summ `numi'
					local num=max(_result(18),0)
					summ `wtsq'
					local denom=_result(18)
					local newtsq=`num'/`denom'
					}
				}

			if "`bsest'"=="ml" {
				qui {
					gen `numi'=(`wtsq')*(((`y'-`ypred')^2)-`wsvar')
					summ `numi'
					local num=max(_result(18),0)
					summ `wtsq'
					local denom=_result(18)
					local newtsq=`num'/`denom'
					}
				}

			if "`bsest'"=="eb" {
				qui {
					gen `numi'=(`wt')*(((`N'/(`N'-(`p'+1)))*((`y'-`ypred')^2))-`wsvar')
					summ `numi'
					local num=max(_result(18),0)
					summ `wt'
					local denom=_result(18)
					local newtsq=`num'/`denom'
					}
				}

			local j=`j'+1

			}
		}

	tempvar wt
	qui {
		gen `wt'=(`wsvar'+`newtsq')^-1
		summ `wt'
		}
	local sumwt=_result(18)

#delimit ;
	di _n
	  in gr "Meta-analysis regression"
	  _col(56) "No of studies =   " in ye  `N' _n
	  in gr _col(56) "tau^2 method      " in ye "`bsest'" _n
	  in gr _col(56) "tau^2 estimate = "  in ye %6.5g `newtsq' _n ; 
#delimit cr
	if "`bsest'"=="ml" | "`bsest'"=="reml" |  "`bsest'"=="eb" {
		di in bl "Successive values of tau^2 differ by less than 10^-"`toleran' " :convergence achieved"
		}

	qui regress `y' `xvars' [aw=`wt']

	local scpar=(`sumwt'*(_result(9)^2))/_result(1)
	local scpar1=(1/`scpar')

	matrix V=get(VCE)
	matrix b=get(_b)

	matrix v=`scpar1'*V


	mat post b v
	mat mlout, level(`level')

	global S_1 = `N'
	global S_2 = `newtsq'

	restore

end


	exit