/* Program to accompany "Simultaneous confidence intervals for multinomial proportions" by Rick Wicklin. Published on The DO Loop blog http://blogs.sas.com/content/iml/2017/02/15/confidence-intervals-multinomial-proportions.html */ /*******************************************************/ /* The following functions will construct simultaneous */ /* confidence intervals for multinomial proportions */ /* based on several methods of construction. */ /* Values for the variable 'Method' are: */ /* 1 = Quesenberry and Hurst (1964) */ /* 2 = Goodman (1965) */ /* 3 = naive binomial */ /* 4 = Fitzpatrick and Scott */ /* (alpha <= 0.15 only) */ /* 5 = Q & H type and estimated variance */ /* 6 = Goodman type and estimated variance */ /* */ /* Original code: May and Johnson (1997). "A SAS Macro */ /* for Constructing Simultaneous Confidence Intervals*/ /* for Multinomial Proportions," Computer Methods and*/ /* Programs in Biomedicine, p. 153-162. */ /* */ /* Rewritten by Rick Wicklin, Feb 2016 */ /* */ /*******************************************************/ /* Define library of SAS/IML function */ proc iml; /* helper function that computes interval based on ChiSq value */ start MultCIFromChiSquare(Count, chiSq); n = Count[+]; p = Count / n; a = chiSq + n; b = chiSq + 2*n*p; c = n*p##2; root = sqrt(b##2 - 4*a*c); Lower = (b-root) / (2*a); Upper = (b+root) / (2*a); return (p || Lower || Upper); finish; /* helper function that computes interval based on ChiSq value and estimated variance factor */ start MultCIFromChiSquareVar(Count, chiSq); n = Count[+]; p = Count / n; a = sqrt( chiSq * p#(1-p) / n ); Lower = p - a; Upper = p + a; return (p || Lower || Upper); finish; /* Quesenberry and Hurst (1964) */ start MultCI_QH(Count, alpha); k = nrow(Count); chiSq = quantile("ChiSquare", 1-alpha, k-1); return( MultCIFromChiSquare(Count, chiSq) ); finish; /* Quesenberry and Hurst with estimated variance factor */ start MultCI_QHVar(Count, alpha); k = nrow(Count); chiSq = quantile("ChiSquare", 1-alpha, k-1); return( MultCIFromChiSquareVar(Count, chiSq) ); finish; /* Goodman (1965) */ start MultCI_Goodman(Count, alpha); k = nrow(Count); chiSq = quantile("ChiSquare", 1-alpha/k, 1); /*Bonferroni adjustment */ return( MultCIFromChiSquare(Count, chiSq) ); finish; /* Goodman with estimated variance factor */ start MultCI_GoodmanVar(Count, alpha); k = nrow(Count); chiSq = quantile("ChiSquare", 1-alpha/k, 1); /*Bonferroni adjustment */ return( MultCIFromChiSquareVar(Count, chiSq) ); finish; /* naive binomial method */ start MultCI_Binomial(Count, alpha); chiSq = quantile("chisquare", 1-alpha, 1); n = Count[+]; p = Count / n; a = sqrt(chiSq / (4*n)); Lower = p - a; Upper = p + a; return (p || Lower || Upper); finish; /* Fitzpatrick and Scott (1987) */ start MultCI_FS(Count, alpha); if alpha <= 0.016 then c2 = quantile("chisquare", 1-alpha/2, 1); else if alpha > 0.016 && alpha <= 0.15 then c2 = (8/9)*quantile("chisquare", 1-alpha/3, 1); else do; print 'alpha must be less than 0.15 to use the Fitzpatrick and Scott method'; return( j(nrow(Count), 3, .) ); end; n = Count[+]; p = Count / n; a = sqrt(c2 / (4*n)); Lower = p - a; Upper = p + a; return (p || Lower || Upper); finish; /*******************************************************/ /* Driver function for May and Johnson (1997) functions*/ /* Print the CI results for various tests. */ /* Values for the variable 'type' are: */ /* 1 = Quesenberry and Hurst (1964) */ /* 2 = Goodman (1965) [default] */ /* 3 = naive binomial */ /* 4 = Fitzpatrick and Scott */ /* (alpha<=0.15 only) */ /* 5 = Q & H type using observed */ /* proportions as variance */ /* 6 = Goodman type using observed */ /* proportions as variance */ /*******************************************************/ start MultCI(Count, alpha=0.05, Method=2); if Method=1 then return MultCI_QH(Count, alpha); else if Method=2 then return MultCI_Goodman(Count, alpha); else if Method=3 then return MultCI_Binomial(Count, alpha); else if Method=4 then return MultCI_FS(Count, alpha); else if Method=5 then return MultCI_QHVar(Count, alpha); else if Method=6 then return MultCI_GoodmanVar(Count, alpha); finish; start MultCIPrint(Category, Count, alpha=0.05, Method=2); title = { "Method of Quesenberry and Hurst (1964)" "Method of Goodman (1965)" "Naive binomial method" "Method of Fitzpatrick and Scott (1987)" "Method of Quesenberry and Hurst (1964), using Var = p(1-p)" "Method of Goodman (1965), using Var = p(1-p) " }; print (title[Method]); if type(Category)='C' then Cat = Category; else Cat = putn(Category, "Best6."); alphaText = putn(alpha, "Best6."); labl = "Simultaneous Confidence Intervals"; CI = MultCI(Count, alpha, Method); result = Count || CI; onemalpha = strip(putn(1-alpha, "Percent7.")); varNames = {"Count" "Proportion"} || (onemalpha + " Lower") || (onemalpha + " Upper"); print result[Label=labl rowname=Cat colname=varNames format=Best6.]; finish; store module=_all_; quit; /* Example calling syntax: data Psych; input Category \$21. Count; datalines; Neurotic 91 Depressed 49 Schizophrenic 37 Personality disorder 43 ; proc iml; load module=(MultCI MultCIPrint); use Psych; read all var {"Category" "Count"}; close; alpha = 0.05; call MultCIPrint(Category, Count, alpha, 2); * Goodman = 2; call MultCIPrint(Category, Count, alpha, 4); * Fitzpatrick and Scott; CI = MultCI(Count, alpha, 2); */