*!  Version 1.0.0         (STB-58: sg156)
program define meanscor, eclass
	version 6.0
	local varlist "required existing"     
	local if "optional"
	local in "optional"
	local options "FIRST(string) SECOND(string) ODD(integer 1)" 
	parse "`*'"

   	if "$S_E_cmd"=="meanscor" {
     		error 301
   }

if ("`first'"=="" & "`second'"==""){

	logit `varlist'
	di in red "WARNING: None of first and second stage covariates" 	
	di in red "are provided ANALYSIS USING COMPLETE CASES ONLY is produced"
	exit
	}


if "`first'"=="" {
	di in red "ERROR: Please enter the first stage covariates "
	exit
	}

if "`second'"==""{
	di in red "ERROR: Please enter the second stage covariates "
	exit
	}



******get the dependent variable **************************
      local c=1
      local depvar " "
      local temp: word `c' of `varlist'
      local depvar `temp'
***********************************************************

local n: word count `varlist'

******get the first stage sample sizes********************	
	tempvar grp cnt ntot length mincnt
	qui egen `grp'=group(`depvar' `first')
	qui egen `cnt'=count(`grp'), by(`depvar' `first')
	qui egen `mincnt'=min(`cnt')
	if (`mincnt'[1]<2) {
		di in red "ERROR:One or more groups have less than 2 subjects"
		coding `depvar' `first'
		di in red "the variance of the score can't be computed"
		di in red "algorithm terminating......"
		exit
		}
	qui egen `ntot'=count(`depvar')
	scalar `length'=`ntot'[1]
   	unabbrev `varlist'

***** get the second stage sample sizes ***************************
      tempvar grp2 cnt2
      qui egen `grp2'=group(`depvar' `first') if !missing(`second')
      qui egen `cnt2'=count(`grp2'), by(`depvar' `first'), if !missing(`second')

***** create the weight for the regression analysis
	tempvar wg pi sc
	qui gen `wg'=`cnt'/`cnt2'

****now, do weighted logistic regression**************************
	preserve
	quietly keep if !missing(`second')
	qui logit `varlist' [pweight=`wg'],score(`sc') 
	matrix b=get(_b)
      predict `pi', p


***** calculate the information matrix****************************
	tempvar info1
	qui gen `info1'=sqrt(`pi'*(1-`pi')*`wg')

	local c=2
	local model " "
	while `c'<=`n'{
		tempvar info`c' 
		local temp: word `c' of `varlist'
		qui gen `info`c''=sqrt(`pi'*(1-`pi')*`wg')*`temp'
		local model `model' `temp'
		local c=`c'+1
		}
      
	qui matrix accum Ihat=`info1'-`info`n'',noconstant 
	qui matrix Ihat=Ihat/`length'

***** get the score vectors *************************************      
	local c=2
	while `c'<=`n'{
	tempvar sc`c'
	local temp: word `c' of `varlist'
	qui gen `sc`c''=`temp'*`sc'
	local c=`c'+1
	}


***** get the variance of the score for each stratum ***********

	tempvar wgi wg2 wg3 lev levels
	qui matrix sandwich=J(`n',`n',0)
	qui egen `levels'=max(`grp2')
	scalar `lev'=`levels'[1]
	local i=1
	while `i'<=`lev'{
		qui matrix accum varsi=`sc' `sc2'-`sc`n'' if `grp2'==`i',noconstant deviations
		qui matrix varsi=varsi/(r(N)-1)
		qui gen `wgi' = (`cnt'/`cnt2')*(`cnt'-`cnt2') if `grp'==`i'
		qui egen `wg2'=max(`wgi')
		scalar `wg3'=`wg2'[1]
		drop `wgi' `wg2' 
		qui matrix varsi=varsi*`wg3'
		qui matrix varsi=syminv(Ihat)*varsi*syminv(Ihat)
		est matrix wzy`i' varsi

***** and compute the sandwhich term **************************

		qui matrix sandwich=sandwich+e(wzy`i')
		local i=`i'+1
	}

***** now get the variance of the estimates
quietly{
matrix sandwich=sandwich/`length'
matrix vce=(syminv(Ihat)+sandwich)/`length'
est matrix info Ihat
}

******* swap the coefficient vector****** 
matrix b1=b[1,`n']
matrix b2=b[1,1..`n'-1]
matrix b=(b1,b2)

matrix rownames vce = cons `model'
matrix colnames vce =cons `model'
matrix colnames b = cons `model'
est matrix covar vce
matrix vce=e(covar)
estimates post b vce

**** display*****************************
di in blue "meanscore estimates"
if `odd'==1{
	estimates display, eform("odd-ratio")}
	else{ 
	    estimates display}
	
est local cmd "meanscor"
restore,preserve
end