*! version 1.0.0 07jan2014

vers 12.1
m :

// mata clear

void mPwmc(real matrix stats, 
		string rowvector prc, 
		real scalar cilev,
		real colvector lvls,
		| string scalar over)
{	
	real colvector n, mu, sd, df
	real scalar k, kstar, alpha, col, c
	real rowvector diff, Var, nuhat, pvals
	real matrix A, CI
	
	string colvector olvls
	string matrix vsnams, prcnams, cinams
	string scalar onam

	n = stats[., 1]
	mu = stats[., 2]
	sd = stats[., 3]	
	
	k = rows(mu)
	kstar = k * (k - 1)/2
	df = (n :- 1)
	alpha = cilev/100
	
	diff = J(1, kstar, .)
	nuhat = J(1, kstar, .)
	Var = J(1, kstar, .)
	A = J(kstar, cols(prc), .)
	CI = J((cols(prc)*kstar), 2, .)
	pvals = J(kstar, cols(prc), .)
	
	vsnams = J(kstar, 2, "")
	prcnams = J(cols(prc), 1, ""), prc'
	olvls = strofreal(lvls)
	if (over != "") onam = "." + over
	
	// Dunnett's C
	if (anyof(prc, "c")) {
		sr = invtukeyprob(k, df, alpha) :* (sd :^ 2 :/ n)
	}
	
	// calculate diff, Var, nuhat and A
	c = 0
	for (j = 1; j <= kstar; ++j) {
		for (i = (j + 1); i <= k; ++i) {
			++c
			
			diff[1, c] = mu[i, 1] - mu[j, 1]
			
			Var[1, c] = (sd[i, 1]^2/n[i, 1]) + (sd[j, 1]^2/n[j, 1])
			
			nuhat[1, c] = Var[1, c]^2 / ///
			((sd[i, 1]^4/(n[i, 1]^2 * df[i, 1])) ///
			+ (sd[j, 1]^4/(n[j, 1]^2 * df[j, 1])))
			
			vsnams[c, 2] = olvls[i, 1] + "vs" + olvls[j, 1] + onam
			
			// Dunnett's C
			if (anyof(prc, "c")) {
				col = select(J(1, 1, (1..cols(prc))), (prc :== "c"))
				A[c, col] = ((sr[i, 1] + sr[j, 1])/Var[1, c]) / sqrt(2)
			}
		}
	}	
	
	// Games and Howell
	if (anyof(prc, "gh")) {
		col = select(J(1, 1, (1..cols(prc))), (prc :== "gh"))
		A[., col] = (invtukeyprob(k, nuhat, alpha) :/ sqrt(2))'
		pvals[., col] = ///
		(1 :- tukeyprob(k, nuhat, abs(diff :/ sqrt(Var)) * sqrt(2)))'
	}
	
	// Tamhane's T2
	if (anyof(prc, "t2")) {
		col = select(J(1, 1, (1..cols(prc))), (prc :== "t2"))
		A[., col] = (invttail(nuhat, (1 - alpha^(1/kstar)) / 2))'
		pvals[., col] = ///
		(1 :- (1 :- 2*ttail(nuhat, abs(diff :/ sqrt(Var)))) :^ kstar)'
	}

	// calculate CI
	cinams = vec(J(kstar, 1, prc)), J(cols(prc), 1, vsnams[., 2])
	CI = vec(J(1, cols(prc), diff') :- A :* sqrt(Var')), ///
	vec(J(1, cols(prc), diff') :+ A :* sqrt(Var'))

	// return in r()
	st_rclear()
	
	st_numscalar("r(k)", k)
	st_numscalar("r(ks)", kstar)
	st_numscalar("r(level)", cilev)
	
	st_global("r(procedure)", st_local("procedure"))
	st_global("r(over)", over)
	st_global("r(depvar)", st_local("varlist"))
	st_global("r(cmd)", "pwmc")
	
	st_matrix("r(A)", A)
	st_matrixcolstripe("r(A)", prcnams)
	st_matrixrowstripe("r(A)", vsnams)
	
	st_matrix("r(p_adj)", pvals)
	st_matrixcolstripe("r(p_adj)", prcnams)
	st_matrixrowstripe("r(p_adj)", vsnams)
	
	st_matrix("r(t)", diff :/ sqrt(Var))
	st_matrixcolstripe("r(t)", vsnams)
	
	st_matrix("r(ci)", CI)
	st_matrixcolstripe("r(ci)", ((""\ ""), ("ll"\ "ul")))
	st_matrixrowstripe("r(ci)", cinams)
	
	st_matrix("r(nuhat)", nuhat)
	st_matrixcolstripe("r(nuhat)", vsnams)
	
	st_matrix("r(Var)", Var)
	st_matrixcolstripe("r(Var)", vsnams)
	
	st_matrix("r(diff)", diff)
	st_matrixcolstripe("r(diff)", vsnams)
	
	st_matrix("r(levels_over)", lvls, "hidden")
	
	// return old results
	for (i = 1; i <= cols(A); ++i) {
		onam = "r(A_" + prc[1, i] + ")"
		st_matrix(onam, A[., i]', "hidden")
		st_matrixcolstripe(onam, vsnams)
	}
	
	for (i = 1; i <= rows(CI); i = i + kstar) {
		onam = "r(ll_" + cinams[i, 1] + ")"		
		st_matrix(onam, CI[(i::i + kstar - 1), 1]', "hidden")
		st_matrixcolstripe(onam, vsnams)
		
		onam = "r(ul_" + cinams[i, 1] + ")"		
		st_matrix(onam, CI[(i::i + kstar - 1), 2]', "hidden")
		st_matrixcolstripe(onam, vsnams)
	}
}

// mata mosave mPwmc() ,dir(PERSONAL) replace

end