clear
capture program drop paired
program define paired, rclass
        version 14.0
        drop _all
        syntax [, nobs(integer 25) rho(real 0.70) mu1(real 50) mu2(real 70) sigma1(real 20) sigma2(real 30) beta(real 5) ]
		set obs `nobs'
		gen n1 = int(`nobs'*0.50)
		gen z=0
		replace z=1 if _n > n1
		gen z1=invnormal(uniform())
		gen z2=invnormal(uniform())
		gen z3=invnormal(uniform())
*
		gen time1=sqrt(`rho')*z1+sqrt(1-`rho')*z2
		replace time1 = time1*`sigma1' + `mu1'
		gen time2=sqrt(`rho')*z1+sqrt(1-`rho')*z3
		replace time2 = time2*`sigma2' + `mu2'
		replace time2 = time2 + `beta'*z
*
		reg time2 z
		scalar ttest=_b[z]
		scalar zttest = _b[z]/_se[z]
		scalar pttest = 2*(1-normal(abs(zttest)))
*
*
		gen diff=time2-time1
		ttest diff, by(z)
		scalar paired = r(mu_2) - r(mu_1)
		scalar ppaired= r(p)

*
		gen id=_n
		drop z1 z2 z2
		reshape long time , i(id)
		xtset id
		gen jj=_j-1
	    xtmixed time i.z#i.jj || id:, 
		lincom (1.z#1.jj - 1.z#0.jj) - (0.z#1.jj - 0.z#0.jj)
		scalar mxt = r(estimate)
		scalar zmxt = r(estimate)/r(se)
		scalar pmxt = 2*(1-normal(abs(zmxt)))
        return scalar ttest = ttest
        return scalar pttest = pttest
        return scalar paired=paired
        return scalar ppaired=ppaired
        return scalar mxt = mxt
        return scalar pmxt = pmxt
end
/* 
     The block below allows you to specify the parameters you will use in the
     simulation. Simply change the values to simulate different settings

	nsim					the number of simulations to be performed
	nobs 					the number of people in the trial (1/2 placebo, 1/2 minocycline)
	rho						the correlation between global deficit scores at baseline and wk 24
	mu1						the mean global deficit score at baseline 
	mu2						the mean global deficit score at week 24 (placebo group)
	sigma1					the SD of global deficit scores at baseline
	sigma2					the SD of global deficit scores at week 24
	beta: 					the effect of minocycline at week 24 (mean difference in GDS betw' minocycline and placebo)

													*/
													
local nsim=2000
local nobs=100
local rho=0.7
local mu1=1.0
local mu2=0.90
local sigma1=0.35
local sigma2=0.35
local beta=-0.15
set more off
set seed 966
capture log using simresults.log, replace
simulate ttest=r(ttest) pttest=r(pttest) mxt=r(mxt) pmxt=r(pmxt) paired=r(paired) ppaired=r(ppaired), reps(`nsim'):paired, nobs(`nobs') rho(`rho') mu1(`mu1') mu2(`mu2') sigma1(`sigma1') sigma2(`sigma2') beta(`beta')
foreach x in ttest mxt paired {
gen power_`x'=0
replace power_`x'=1 if p`x' < 0.05
}
summ
log close