mreggersimexonce.ado

*! 0.1.0 23jul2017 Tom Palmer
program mreggersimexonce, rclass
syntax varlist(min=2 max=2) [aweight] [if] [in] , gxse(varname numeric) ///
	[ivw fe re ///
	reslope recons HETerogi noRESCale PENWeighted Level(cilevel) ///
	tdist lambda(numlist) seed(string) simreps(int 100) *]

tokenize `"`varlist'"'
/*
2 variables:
1: gd beta
2: gp beta
aw: gd SE
*/
tempvar mreggery mreggerx
qui gen double `mreggery' = `1'*sign(`2') `if' `in'
qui gen double `mreggerx' = abs(`2') `if' `in'        

** values of lambda to loop through
if "`lambda'" == "" {
	local lambda ".5 1 1.5 2"
}

** number of genotypes
qui count `if' `in'
local n = r(N)

** id variable
tempvar id
qui gen `id' = _n
    
tempvar xunknown
qui putmata `mreggerx' `gxse' `id' `if' `in', replace
tempvar res
tempname slope cons
local lambdalength = wordcount("`lambda'")
local nres = `simreps'*`lambdalength'
mata `res' = J(`nres', 3, .)

** lambda = 0: implies original data
tempname slopelambda0 conslambda0
qui mregger `mreggery' `mreggerx' [aw `exp'] `if' `in', ///
	`ivw' `fe' `re' ///
	`reslope' `recons' `heterogi' `rescale' ///
	`penweighted' level(`level') ///
	`tdist' `options'
scalar `slopelambda0' = _b[slope]
scalar `conslambda0' = _b[_cons]
local j = 0
foreach lam of numlist `lambda' {
	local j = `j' + 1
	forvalues i=1/`simreps' {
		** the new variable with measurement error
		mata `xunknown' = `mreggerx' + sqrt(`lam'):*`gxse':*rnormal(`n', 1, 0, 1)
		qui getmata `xunknown', double replace id(`id')
		qui mregger `mreggery' `xunknown' [aw `exp'] `if' `in', ///
			`ivw' `fe' `re' ///
			`reslope' `recons' `heterogi' `rescale' ///
			`penweighted' level(`level') ///
			`tdist' `options'
		scalar `slope' = _b[slope]
		scalar `cons' = _b[_cons]
		local resrow = `i' + (`j' - 1)*`simreps'
		mata `res'[`resrow',] = (`lam', st_numscalar("`slope'"), st_numscalar("`cons'"))
	}
}

preserve
drop _all
qui set obs `nres'
getmata (lambda slope cons)=`res', double
collapse (mean) slope cons, by(lambda)
local colln = _N
qui set obs `=_N+1'
gsort -lambda
qui replace lambda = 0 in `=_N'
qui replace slope = `slopelambda0' in `=_N'
qui replace cons = `conslambda0' in `=_N'
qui gen double lambda2 = lambda*lambda
qui putmata simexests=(lambda slope cons), replace

mata simexests2 = J(6,3,.)
mata simexests2[1..5,] = simexests
qui regress slope lambda lambda2
local ext -1
local pred = _b[lambda]*`ext' + _b[lambda2]*`ext'*`ext' + _b[_cons]
mata simexests2[6,1..2] = (-1, `pred')
ret sca slope = `pred'

qui regress cons lambda lambda2
local pred = _b[lambda]*`ext' + _b[lambda2]*`ext'*`ext' + _b[_cons]
mata simexests2[6,3] = (`pred')
ret sca cons = `pred'

mata st_matrix("simexests2", simexests2)
local names "lambda slope cons"
mat colnames simexests2 = `names'
return matrix simexests = simexests2

qui mata mata drop simexests simexests2 `gxse' `mreggerx' `id' `res' `xunknown'
end