program define msopt, eclass
	version 6.0
	local varlist "required existing"     
	local if "optional"
	local in "optional"
	local options "FIRST(string) SAMPLE(string) CODING(integer 0)" 
	parse "`*'"

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


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

*******get the dependent and predictor variable************	
      local c=1
	local vars " "
	unabbrev `varlist'
	local n: word count `varlist'
	while `c'<=`n' {
		local temp: word `c' of `varlist'
		local vars `vars' `temp'
		local c=`c'+1
		}

*****determine the number of strata and second stage*******
*****sample size for each stratum**************************

      tempvar grp2 cnt2
      qui egen `grp2'=group(`depvar' `first') 
	qui egen `cnt2'=count(`grp2'), by(`depvar' `first')
	di in red "the second stage sample sizes"
	table `grp2'
	tempvar lev levels
	qui egen `levels'=max(`grp2')
	scalar `lev'=`levels'[1]

*********get the nobs and sample size for each stratum******
	tempvar ntot length sample1
	local i=1
	gen `sample1' = 0
	while `i'<=`lev' {
		qui replace `sample1'=`sample'[`i',1] in `i'
		local i=`i'+1
		}

	qui egen `ntot'=sum(`sample1')
	scalar `length'=`ntot'[1]

	tempvar cnt 
	local i=1
	gen `cnt'=1
	while `i'<=`lev'{
		qui replace `cnt'=`sample1'[`i'] if `grp2'==`i'
		local i=`i'+1
		}
	
**********checking the second stage sample sizes************
	if `coding'==1 {
		coding `depvar' `first'
		}

	tempvar grp_z grp_yz nobs
	qui egen `grp_z'=group(`first')
	qui egen `grp_yz'=group(`depvar' `first') 
	qui egen `nobs' =count(`grp_z'), by(`grp_yz')

	preserve
	if `length'<=2 {
	gen prev =`cnt'}
		else {
		gen n1=`cnt'}

	qui gen n2_pilot=`nobs'
	gen grp_z=`grp_z'
	gen grp_yz=`grp_yz'
	
	format `varlist' grp_z %5.0g

	if `length'<= 2 {
	format grp_yz prev n2_pilot %6.0g}
		else {
		format grp_yz n1 n2_pilot %6.0g}

	di in red "please check the sample sizes!"

	if `length'<= 2 {
	collapse `depvar' `first' grp_z prev n2_pilot, by (grp_yz)}
		else {
		collapse `depvar' `first' grp_z n1 n2_pilot, by (grp_yz)}

	list, nodisplay nolabel noobs
	restore, preserve
	
***********weighted logistic regression*********************
	tempvar wg pi sc
	qui gen `wg'=`cnt'/`cnt2'
	local varist `vars'
	qui logit `varist' [pweight=`wg'],score(`sc') 
	matrix b=get(_b)
      predict `pi', p

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

	local c=2
	while `c'<=`n'{
		tempvar info`c' 
		local temp: word `c' of `varist'
		qui gen `info`c''=sqrt(`pi'*(1-`pi')*`wg')*`temp'
		local c=`c'+1
		}
      
	qui matrix accum Ihat=`info1'-`info`n'',noconstant 
	qui matrix Ihat=Ihat/`length'
************************************************************
      
****get the variance of the score of each stratum***********
	local c=2
	while `c'<=`n'{
		tempvar sc`c'
		local temp: word `c' of `varist'
		qui gen `sc`c''=`temp'*`sc'
		local c=`c'+1
		}

	tempvar wgi wg2 wg3 
	qui matrix sandwich=J(`n',`n',0)
	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 `grp2'==`i'
		qui egen `wg2'=max(`wgi')
		scalar `wg3'=`wg2'[1]
		drop `wgi' `wg2' 
		qui matrix varsi=syminv(Ihat)*varsi*syminv(Ihat)
		est matrix wzy`i' varsi

		if `length'>2 {
		qui matrix varsi=e(wzy`i')*`wg3'
		qui matrix sandwich=sandwich+varsi
		}
	local i=`i'+1
	}
************************************************************

*****swap the vector of coefficient for displaying purpose*******
matrix b1=b[1,`n']
matrix b2=b[1,1..`n'-1]
matrix b=(b1,b2)
local name:colnames b

**********get the variance of meanscore estimates*********** 
*******(if first stage sample sizes was supplied!)**********

if `length'> 2 {
	quietly{
		matrix sandwich=sandwich/`length'
		matrix vce=(syminv(Ihat)+sandwich)/`length'
		matrix rownames vce =`name'
		matrix colnames vce =`name'
		est matrix covar vce
		}
	}

*****export the information matrix and number of strata ********
est scalar lev=`lev'
est matrix info Ihat 
est local cmd "msopt"
end