*! version 1.0.0 07oct2019 Matteo Bottai
cap program drop m_power
program m_power
	syntax , [ LEVEL(real 0.05) POWER(real 0.8) SIGMA(real 1) DIFFerence(numlist min=0 max=1) Subjects(numlist min=0 max=1) Timepoints(integer 3) D(string) R(string) ]
	if "`subjects'"!="" & "`difference'"!="" {
		di as err "The options difference snd subjects are mutually exclusive"
		exit
		}
	mata: mata clear
	if "`subjects'"!="" {
		mata: m	= `subjects'
		}
	else if "`difference'"!="" {
		mata: difference = `difference'
		}
	else {
		di as err "Either difference or subjects must be specified in the options"
		exit
		}
	if "`d'"=="" {
		local d = "I(2)"
		}
	if "`r'"=="" {
		local r = "I(n)"
		}
	mata: n = `timepoints'
	mata: sigma	= `sigma'
	mata: power	= `power'
	mata: level = `level'
	mata: D	= `d'
	mata: R	= `r'
	mata: X1 = J(n,1,1), (1::n), J(n,1,0), J(n,1,0)
	mata: X2 = J(n,1,1), (1::n), J(n,1,1), (1::n)
	mata: Z	= X1[1::n,1::2]
	mata: Vinv = invsym(Z*D*Z' :+ R*sigma^2)
	mata: se = invsym(X1'*Vinv*X1 :+ X2'*Vinv*X2)
	if "`subjects'"!="" {
		mata: difference = sqrt((invnormal(1-level/2)+invnormal(power))^2*se[4,4]/m)
		mata: st_local("difference",strofreal(difference))
		noi di _new as txt "Detectable slope difference between two groups with mixed effects model"
		di _new as txt "Number of subjects" _col(33)  "=  " as res `subjects'
		di as txt "Number of time points" _col(33)  "=  " as res `timepoints'
		di as txt "Level" _col(33)  "=  " as res `level'
		di as txt "Power" _col(33)  "=  " as res `power'
		di as txt "Residual standard deviation" _col(33)  "=  " as res `sigma'
		di as txt "Residual correlation" _col(33)  "=  " as res subinstr(subinstr("`r'","I(","Identity(",.),"(n)","(`timepoints')",.)
		di as txt "Random coefficients covariance" _col(33)  "=  " as res subinstr("`d'","I(","Identity(",.)
		di _new as txt "Estimated slope difference" _col(33) "=  " as res `difference'
		}
	else {
		mata: m = (invnormal(1-level/2)+invnormal(power))^2*se[4,4]/difference^2
		mata: st_local("subjects",strofreal(m))
		noi di _new as txt "Detectable slope difference between two groups with mixed effects model"
		di _new as txt "Slope difference" _col(33)  "=  " as res `difference'
		di as txt "Number of time points" _col(33)  "=  " as res `timepoints'
		di as txt "Level" _col(33)  "=  " as res `level'
		di as txt "Power" _col(33)  "=  " as res `power'
		di as txt "Residual standard deviation" _col(33)  "=  " as res `sigma'
		di as txt "Residual correlation" _col(33)  "=  " as res subinstr(subinstr("`r'","I(","Identity(",.),"(n)","(`timepoints')",.)
		di as txt "Random coefficients covariance" _col(33)  "=  " as res subinstr("`d'","I(","Identity(",.)
		di _new as txt "Estimated number of subjects" _col(33) "=  " as res ceil(`subjects')
		}
end