*! fwer1stage v0.02 cap program drop fwer1stage cap mata: mata drop mvnpb() program def fwer1stage, rclass version 10 /* Calculate FWER for multi-arm 1-stage trials algebraically v0.02 - can specify 1 alpha for each arm or a common alpha */ syntax, arms(int) alpha(string) aratio(real) local K = `arms'-1 // # E arms local A = `aratio' local nopts: word count `alpha' if `nopts'!=1 & `nopts'!=`K' { di as err "1 or `K' options should be specified for alpha" exit 198 } forvalues k = 1/`K' { if `nopts'==1 local z`k' = invnormal(1-`alpha') else { local a`k': word `k' of `alpha' local z`k' = invnormal(1-`a`k'') } } // Correlation matrix local r = `A'/(`A'+1) tempname R matrix def `R' = I(`K') forvalues k = 1/`K' { forvalues l = 1/`K' { if `k'!=`l' matrix def `R'[`k',`l'] = `r' else matrix def `R'[`k',`l'] = 1 } } // Vector of z values local z1ma `z1' forvalues k = 2/`K' { local z1ma `z1ma', `z`k'' } tempname Z matrix `Z' = (`z1ma') local rep = 5000 mata: mvnpb("`Z'", "`R'", `rep') local fwer = 1-r(p) return scalar fwer = `fwer' di as text "FWER = " as res %5.4f `fwer' end mata: mata clear void mvnpb(string scalar xx, string scalar vv, real scalar reps) { /* Assumes Hammersley sequences are to be generated, without antithetics. */ real vector opt x = st_matrix(xx) // row or col vector of args at which probability is required V = st_matrix(vv) // variance-covariance matrix opt = (2, reps, 1, 0) // 2 for Hammersley p = ghk( x, V, opt, rank=.) st_numscalar("r(p)", p) } end