*! version 1.1 Thursday, July 3, 2003 at 12:23 (SJ3-3: st0000) * * This program implements the algorithm for the multivariate Wald test * as described in Li, K., T. Raghunathan, and D. Rubin. 1991 program define mitestparm version 8 syntax varlist capture assert "\$mimps"~=""&"\$mi_sf"~="" if _rc { display as error "please set up your data with -{help miset}- first" exit 198 } local m = \$mimps forvalues t=1/`m' { matrix _b`t' = e(b_`t') matrix _V`t' = e(V_`t') } tokenize `varlist' * find column numbers (and hence row numbers by symmetry) for each var in `varlist' local k = 0 while "`1'"~="" { local k = `k'+1 local var`k' = colnumb(_b1,"`1'") mac shift } * build new coefficient vectors, b`t', containing just the elements corresponding to vars in `varlist' forvalues t=1/`m' { matrix b`t' = J(1,`k',0) forvalues j=1/`k' { matrix b`t'[1,`j']= _b`t'[1,`var`j''] } } * similarly build new variance-covariance matrices, V`t' forvalues t=1/`m' { matrix V`t' = J(`k',`k',0) forvalues i=1/`k' { forvalues j=1/`k' { matrix V`t'[`i',`j']= _V`t'[`var`i'',`var`j''] } } } * calculate average of coefficient vectors matrix matsum = J(1,`k',0) /*set `matsum' to 1xk zero matrix*/ forvalues t=1/`m' { matrix matsum = matsum + b`t' } matrix Qbar = 1/`m' * matsum * calculate within imputation variance, Ubar matrix matsum = J(`k',`k',0) /*set `matsum' to kxk zero matrix*/ forvalues t=1/`m' { matrix matsum = matsum + V`t' } matrix Ubar = 1/`m' * matsum * calculate between imputation variance, B matrix matsum = J(`k',`k',0) /*set `matsum' to kxk zero matrix*/ forvalues t=1/`m' { matrix matsum = matsum + (b`t'-Qbar)'*(b`t'-Qbar) } matrix B = 1/(`m'-1) * matsum * calculate total variance estimate, Ttilde matrix Ubarinv = inv(Ubar) matrix B_Ubarinv = B * Ubarinv local r = 1/(`m'-1) * trace(B_Ubarinv)/`k' matrix Ttilde = (1-`r') * Ubar * calculate test statistic, dee matrix Q_0 = J(1,`k',0) matrix Qdiff = Qbar-Q_0 matrix Ttildeinv = inv(Ttilde) matrix D = Qdiff * Ttildeinv * Qdiff'/`k' local dee = trace(D) * calculate approximation for degrees of freedom, df local a = `k'*(`m'-1) if `a'>4 { local df = 4 + (`a'-4)*(1+(1- 2/`a')/`r')^2 } else { local df = `a'*(1+1/`k')*(1+1/`r')^2/2 } * calculate p-value from F distribution local p = Ftail(`k',`df',`dee') * display results di tokenize `varlist' forvalues i=1/`k' { di as txt " ( `i')" as res " `1' = 0" mac shift } di di as txt " F(" %3.0f `k' "," %6.0f `df' ") =" as res %8.2f `dee' di as txt _col(13) "Prob > F =" as res %10.4f `p' end