*! version 4.0.1 PR 09Nov95. program define aov1cont /* contrast in oneway anova from means, SDs and Ns. */ version 4.0 local varlist "req ex min(3) max(3)" local if "opt" local in "opt" local options "Coeffs(string) TRend LEvel(integer $S_level)" parse "`*'" if "`coeffs'"=="" & "`trend'"=="" { di in red "no contrast coefficients given" exit 198 } if "`coeffs'"!="" { if "`trend'"!="" { di in red "cannot have both trend and coeffs" exit 198 } confirm var `coeffs' } parse "`varlist'", parse(" ") local n `1' local mean `2' local sd `3' tempvar touse tmp1 tmp2 a quietly { mark `touse' `if' `in' markout `touse' `n' `mean' `sd' count if `touse' local ngrp = _result(1) if `ngrp'<2 { di in red "at least two groups are required" exit 2001 } if "`trend'"!="" { gen `a' = sum(1) if `touse' } else gen `a' = `coeffs' if `touse' count if `touse' & `a'!=. local nc = _result(1) if `nc'<2 { di in red /* */ "at least two numeric contrast coefficients are required" exit 2001 } /* Standardize contrast coeffs */ sum `a' replace `a' = `a'-_result(3) /* Calc pooled SD (i.e. residual SD from oneway anova) */ gen `tmp1' = `sd'^2 if `touse' gen `tmp2' = `n'-1 if `touse' sum `tmp1' [weight=`tmp2'] local rsd = sqrt(_result(3)) sum `tmp2' local rdf = `ngrp'*_result(3) if `rdf'==0 | `rdf'==. { di in red "zero residual degrees of freedom" exit 2001 } local obs = `rdf'+`ngrp' /* Weighted mean, SSQ between groups, ANOVA F test */ sum `mean' [aw=`n'] if `touse' local xbar = _result(3) replace `tmp1' = `n'*(`mean'-`xbar')^2 if `touse' sum `tmp1' local ssgrp = _result(1)*_result(3) local vratio = (`ssgrp'/(`ngrp'-1))/`rsd'^2 local Pgrp = fprob(`ngrp'-1, `rdf', `vratio') /* Contrast and its SE (by literal method) */ replace `tmp1' = `mean'*`a' if `touse' sum `tmp1' local L = _result(1)*_result(3) replace `tmp1' = `a'^2/`n' if `touse' sum `tmp1' local seL = `rsd'*sqrt(_result(1)*_result(3)) local t = `L'/`seL' } *!! Following code is a horrible kludge to deal with -list-'s inability *!! to cope with user-supplied var labels. Don't need __a to be permanent. cap drop __a rename `a' __a list `n' `mean' `sd' `coeffs' __a if `touse' drop __a local invt = invt(`rdf', `level'/100) local lower = `L' - `invt'*`seL' local upper = `L' + `invt'*`seL' local P = fprob(1, `rdf', `t'^2) /* Output results */ _jprcout `mean' group `obs' `L' `seL' `lower' `upper' `t' `P' /* */ `vratio' `Pgrp' `rdf' `rsd' `level' global S_1 `obs' global S_2 `ngrp' global S_3 `L' global S_4 `seL' global S_5 `lower' global S_6 `upper' global S_7 `t' global S_8 `P' end