/* run to create data set */ data bp_drug; input Drug Disease $ BP 3.; datalines; 1 A 119.70095982 1 A 121.36188924 1 A 119.69180851 1 A 119.60217745 1 A 120.96573695 1 A 119.18963686 1 A 120.04126616 1 A 120.6485107 1 A 121.39677755 1 A 121.2941861 2 A 134.70095982 2 A 136.36188924 2 A 134.69180851 2 A 134.60217745 2 A 135.96573695 2 A 134.18963686 2 A 135.04126616 2 A 135.6485107 2 A 136.39677755 2 A 136.2941861 3 A 139.17670991 3 A 139.7596925 3 A 139.61692792 3 A 139.47841248 3 A 140.32687732 3 A 139.19503074 3 A 139.82084261 3 A 139.33853642 3 A 139.18391356 3 A 141.16519705 4 A 151.0314732 4 A 148.45628664 4 A 150.90924814 4 A 150.36754897 4 A 150.40537109 4 A 150.27207554 4 A 148.84462813 4 A 148.98134213 4 A 148.67342047 4 A 151.22752311 1 B 159.47367321 1 B 157.8492618 1 B 159.71157707 1 B 161.08606542 1 B 157.88059214 1 B 160.70808134 1 B 158.93686029 1 B 159.59337556 1 B 161.04857471 1 B 162.4582518 2 B 149.44343695 2 B 149.51766786 2 B 151.67540302 2 B 148.43104013 2 B 150.04534684 2 B 149.53540532 2 B 149.21944044 2 B 149.29473811 2 B 149.80104295 2 B 150.2553889 3 B 140.84243239 3 B 139.80385608 3 B 140.01491297 3 B 140.84663178 3 B 139.47817427 3 B 139.34930573 3 B 140.06553297 3 B 139.90122485 3 B 140.48542888 3 B 139.48601375 4 B 129.32979505 4 B 130.1515735 4 B 130.44367458 4 B 129.75982193 4 B 129.05742185 4 B 130.10355508 4 B 130.70093917 4 B 130.88766817 4 B 131.63619676 4 B 129.8023949 1 C 126.23902904 1 C 125.75600563 1 C 125.87874742 1 C 125.04661693 1 C 122.97475788 1 C 123.59937382 1 C 125.16139122 1 C 124.73864001 1 C 123.42349457 1 C 125.38733391 2 C 125.69890209 2 C 125.56064272 2 C 123.83433943 2 C 124.6489834 2 C 123.29421225 2 C 125.24488558 2 C 123.753593 2 C 125.1478912 2 C 126.59647236 2 C 124.56548104 3 C 125.09941212 3 C 126.00783602 3 C 123.7285012 3 C 123.81351742 3 C 124.62018186 3 C 124.68909051 3 C 124.86764454 3 C 124.00137446 3 C 125.56161889 3 C 124.01149493 4 C 124.8987427 4 C 124.38445589 4 C 124.91319163 4 C 125.86522746 4 C 123.22694283 4 C 125.61937452 4 C 126.41012761 4 C 125.39699702 4 C 125.05360115 4 C 125.49297642 ; run; /* I. One-Way ANOVA */ /* explore data */ proc univariate data=bp_drug; var BP; histogram / normal; probplot; run; proc sort data=bp_drug; by drug; run; symbol color = salmon; proc boxplot data=bp_drug; plot BP*Drug / cframe = vligb cboxes = dagr cboxfill = ywh; title 'Boxplot Blood Pressure Drugs'; run; /* one-way ANOVA */ options ls=75 ps=45; /* page and line spacing options*/ proc glm data = bp_drug; class drug; model BP=drug; means drug / hovtest; output out=check r=resid p=pred; /* creating an output data set called 'check' in which 'r ' is a keyword SAS recognizes as residuals and 'p' recognized as predictors */ title 'testing for equality of mean BP by Drug with GLM'; run; quit; /* have to quit our of GLM */ /* now run gplot on the 'check' dataset created above*/ Proc gplot data=check; Plot resid*pred / haxis=axis1 vaxis=axis2 vref=0; /* can leave out if ok with defaults*/ axis1 w=2 major=(w=2) minor=none offset=(10pct); axis2 w=2 major=(w=2) minor=none; title 'plot residuals vs predictors for drugs'; run; quit; /* now run proc univariate to get histogram, normal plot, kurtosis and skewness on residuals*/ proc univariate data = check normal; var resid; histogram resid / normal; probplot;* resid / mu=est sigma=est color=blue w=1; title; run; /* compare means */ proc glm data = bp_drug; class drug; model BP=drug; lsmeans drug / pdiff=all adjust=tukey; /* pdiff=all requests all pairwaise p values */ output out=check r=resid p=pred; /* creating an output data set called 'check' in which 'r ' is a keyword SAS recognizes as residuals and 'p' recognized as predictors */ title 'testing for equality of mean BP by Drug with GLM'; run; quit; /* II. Two-Way ANOVA */ /* means for combinations of disease and drug */ proc means data=bp_drug mean var std; class disease drug; var BP; title 'Selected Descriptive Statistics for drug-disease combinations'; run; /* means plot */ proc gplot data=bp_drug; symbol c=blue w=2 interpol=std1mtj line=1;/* interpolation method gives s.e. bars */ symbol2 c=green w=2 interpol=std1mtj line=2; symbol3 c=red w=2 interpol=std1mtj line=3; plot BP*drug=disease; /* vertical by horizontal */ title 'Illustrating the Interaction Between Disease and Drug'; run; quit; /* PROC GLM for ANOVA with interaction term */ proc glm data=bp_drug; class disease drug; model BP=disease drug disease*drug;/*note interaction term*/ title 'Analyze the Effects of Drug and Disease'; title2 'Including Interaction'; run; quit; /* compare means with interaction term */ proc glm data=bp_drug; class disease drug; model BP=drug disease drug*disease; lsmeans disease*drug / adjust=tukey pdiff=all; /*note looking at combinations*/ title 'Multiple Comparisons Tests for Drug and Disease'; run; quit; /***************************************************************************************/ /* SAS SYNTAX FOR PROC LOGISTIC TEST MULTIPLICATIVE INTERACTION (ON OR SCALE) */ /* FROM DR. MELANIE WALL, COLUMBIA UNIVERSITY DEPTS. PSYCHIATRY AND BIOSTATISTICS */ /***************************************************************************************/ libname itx "X:\Dropbox\sasBoston13"; data brownharris; set itx.brownharris; run; proc contents data=brownharris; run; proc logistic data = brownharris descending; model depression (event= LAST) = stress vulnerability stress*vulnerability; oddsratio stress / at(vulnerability = 0 1); run; PROC NLMIXED DATA=brownharris; ***logistic regression model is; odds = exp(b0 +b1*stress + b2*vulnerability + b3*stress*vulnerability); pi = odds/(1+odds); MODEL depression~BINARY(pi); estimate 'p00' exp(b0)/(1+exp(b0)); estimate 'p01' exp(b0+b1)/(1+exp(b0+b1)); estimate 'p10' exp(b0+b2)/(1+exp(b0+b2)); estimate 'p11' exp(b0+b1+b2+b3)/(1+exp(b0+b1+b2+b3)); estimate'p11-p10'exp(b0+b1+b2+b3)/(1+exp(b0+b1+b2+b3))- exp(b0+b2)/(1+exp(b0+b2)); estimate 'p01-p00' exp(b0+b1)/(1+exp(b0+b1)) - exp(b0)/(1+exp(b0)); estimate 'IC= interaction contrast = p11-p10 - p01 + p00' exp(b0+b1+b2+b3)/(1+exp(b0+b1+ b2+b3)) - exp(b0+b2)/(1+exp(b0+b2)) - exp(b0+b1)/(1+exp(b0+b1)) + exp(b0)/(1+exp(b0)); run;