************************************************************; ****** CODE FOR ANALYSIS AND VISUALIZATION OF INTERACTIONS IN SAS; proc format; value prog 1 = 'jogging' 2 = 'swimming' 3 = 'reading'; value female 0 = 'male' 1 = 'female'; run; data exercise; set exercise; format prog prog.; format female female.; run; ********************SIMPLE SLOPES***********************************; *******************CONT BY CONT************************; *model; proc glm data=exercise; model loss = hours|effort / solution; run; *model; proc glm data=exercise; model loss = hours|effort / solution; store contcont; run; *example estimate to get predicted loss; proc plm restore=contcont; estimate 'pred loss, hours=2, effort=30' intercept 1 hours 2 effort 30 hours*effort 60 / e; run; *get mean and standard deviation of effort; proc means data = exercise; var effort hours; run; *mean and sd of hours 2.0024033 0.4945377; *mean and sd of effort 29.6592189 5.1427635; *get slopes; proc plm restore=contcont; estimate 'hours slope, effort=mean-sd' hours 1 hours*effort 24.52, 'hours slope, effort=mean' hours 1 hours*effort 29.66, 'hours slope, effort=mean+sd' hours 1 hours*effort 34.8 / e; run; *compare slopes; proc plm restore=contcont; estimate 'diff slopes, mean+sd - mean' hours*effort 5.14; run; *contour plot; proc plm source=contcont; effectplot contour (x=hours y=effort); run; *line plot; proc plm source=contcont; effectplot fit (x=hours) / at(effort=24.52 29.66 34.8); run; *create dataset for scoring; data scoredata; do effort = 24.52, 29.66, 34.8; do hours = 0 to 4 by 0.1; output; end; end; run; /* data scoredata; input effort hours; datalines; 24.52 0 24.52 1 24.52 2 24.52 3 24.52 4 29.66 0 29.66 1 29.66 2 29.66 3 29.66 4 34.8 0 34.8 1 34.8 2 34.8 3 34.8 4 ; run; */ *score data; proc plm source=contcont; score data=scoredata out=plotdata predicted=pred lclm=lower uclm=upper; run; *plot simple slopes; proc sgplot data=plotdata; band x=hours upper=upper lower=lower / group=effort transparency=0.5; series x=hours y=pred / group=effort; yaxis label="predicted loss"; run; *****QUADRATIC EFFECTS; proc glm data=exercise order=internal; model loss = hours|hours / solution; store quad; run; *get simple slopes; proc plm restore=quad; estimate 'hours slope, hours=mean-sd(1.5)' hours 1 hours*hours 3, 'hours slope, hours=mean(2)' hours 1 hours*hours 4, 'hours slope, hours=mean+sd(2.5)' hours 1 hours*hours 5 / e; run; *graphing what does proc plm get you; proc plm restore=quad; effectplot fit (x=hours); run; /* *get datapoints on graph; data fake; do hours = 0 to 4 by .25; output; end; run; proc glm data=exercise; model loss = hours|hours / solution; store quad; run; proc plm source=quad; score data=fake out=quadslope predicted=pred lcl=lcl ucl=ucl; run; data plotdata; set exercise quadslope; run; proc sgplot data = plotdata; band x=hours upper=ucl lower=lcl; series x=hours y = pred; scatter x=hours y=loss; run; */ *********************************CAT BY CONT*************************************; *first model; proc glm data=exercise order=internal; class prog; model loss = hours|prog / solution; store catcont; run; *simple slopes by prog; proc plm restore = catcont; estimate 'hours slope, prog=1 jogging' hours 1 hours*prog 1 0 0, 'hours slope, prog=2 swimming' hours 1 hours*prog 0 1 0, 'hours slope, prog=3 reading' hours 1 hours*prog 0 0 1 / e; run; *bad reference slope syntax; proc plm restore = catcont; estimate 'hours slope, prog=3 reading (wrong)' hours 1 / e; run; *****TEST SLOPES AGAINST EACH OTHER; proc plm restore = catcont; estimate 'diff slopes, prog=1 vs prog=2' hours*prog -1 1 0, 'diff slopes, prog=1 vs prog=3' hours*prog -1 0 1, 'diff slopes, prog=2 vs prog=3' hours*prog 0 -1 1 / e; run; *graph simple slopes; proc plm source=catcon; effectplot slicefit (x=hours sliceby=prog) / clm; run; /* ADD DATA POINTS TO CATCONT *slopes plotted on actual data; proc glm data=exercise order=internal; class prog; model loss = prog|hours / solution; output out=plotdata predicted=pred; run; quit; proc sort data = plotdata; by hours; run; proc sgplot data = plotdata; scatter x = hours y = loss / group=prog; series x = hours y = pred / group=prog; run; proc glm data=exercise order=internal; class prog; model loss = prog|hours / solution; store catcon; run; data fake; do prog = 1 to 3; do hours = 0 to 4 by .1; output; end; end; run; proc plm source=catcon; score data=fake out=slopes predicted=pred; run; data plotdata; set exercise slopes; run; proc sgplot data = plotdata; scatter x = hours y = loss / group=prog; series x = hours y = pred / group=prog; run; */ *********************************SIMPLE EFFECTS**********************************; *first get the model; proc glm data=exercise order=internal; class female prog; model loss = female|prog / solution; store catcat; run; *slice statements; proc plm restore = catcat; slice female*prog / sliceby=prog diff adj=bon plots=none nof e means; slice female*prog / sliceby=female diff adj=bon plots=none nof e means; run; *all lsmestimate simple effects; proc plm restore=catcat; lsmestimate female*prog 'male-female, prog = jogging(1)' [1, 1 1] [-1, 2 1], 'male-female, prog = swimming(2)' [1, 1 2] [-1 ,2 2], 'male-female, prog = reading(3)' [1, 1 3] [-1, 2 3], 'jogging-reading, female = male(0)' [1, 1 1] [-1, 1 3], 'jogging-reading, female = female(1)' [1, 2 1] [-1, 2 3], 'swimming-reading, female = male(0)' [1, 1 2] [-1, 1 3], 'swimming-reading, female = female(1)' [1, 2 2] [-1, 2 3], 'jogging-swimming, female = male(0)' [1, 1 1] [-1, 1 2], 'jogging-swimming, female = female(1)' [1, 2 1] [-1, 2 2] / e adj=bon; run; *example compare simple effects; proc plm restore=catcat; lsmestimate prog*female 'diff male-female, prog=1 - prog=2' [1, 1 1] [-1, 2 1] [-1, 1 2] [1, 2 2], 'diff male-female, prog=1 - prog=3' [1, 1 1] [-1, 2 1] [-1, 1 3] [1, 2 3], 'diff male-female, prog=2 - prog=3' [1, 1 2] [-1, 2 2] [-1, 1 3] [1, 2 3], 'diff jogging-reading, male - female' [1, 1 1] [-1, 1 3] [-1, 2 1] [1, 2 3], 'diff swimming-reading, male - female' [1, 1 2] [-1, 1 3] [-1, 2 2] [1, 2 3], 'diff jogging-swimming, male - female' [1, 1 1] [-1, 1 2] [-1, 2 1] [1, 2 2]/ e adj=bon; run; *graph simple effects; proc plm restore=catcat; effectplot interaction (x=female sliceby=prog) / clm connect; effectplot interaction (x=prog sliceby=female) / clm connect; run; ***********THREE WAY INTERACTION*****************; proc glm data = exercise order=internal; class female prog; model loss = female|prog|hours / solution; store catcatcon; run; *simple slopes analysis; *first calculate simple slopes; proc plm restore=catcatcon; estimate 'hours slope, male prog=jogging' hours 1 hours*female 1 0 hours*prog 1 0 0 hours*female*prog 1 0 0 0 0 0, 'hours slope, male prog=swimming' hours 1 hours*female 1 0 hours*prog 0 1 0 hours*female*prog 0 1 0 0 0 0, 'hours slope, male prog=reading' hours 1 hours*female 1 0 hours*prog 0 0 1 hours*female*prog 0 0 1 0 0 0, 'hours slope, female prog=jogging' hours 1 hours*female 0 1 hours*prog 1 0 0 hours*female*prog 0 0 0 1 0 0, 'hours slope, female prog=swimming' hours 1 hours*female 0 1 hours*prog 0 1 0 hours*female*prog 0 0 0 0 1 0, 'hours slope, female prog=reading' hours 1 hours*female 0 1 hours*prog 0 0 1 hours*female*prog 0 0 0 0 0 1 / e adj=bon; run; *within each program do slopes of genders differ - test of conditional interaction of hours*female within each prog; proc plm restore=catcatcon; estimate 'diff hours slope, male-female prog=1' hours*female 1 -1 hours*female*prog 1 0 0 -1 0 0, 'diff hours slope, male-female prog=2' hours*female 1 -1 hours*female*prog 0 1 0 0 -1 0, 'diff hours slope, male-female prog=3' hours*female 1 -1 hours*female*prog 0 0 1 0 0 -1 / e adj=bon; run; *do the difference between gender slopes DIFFER between programs?; proc plm restore=catcatcon; estimate 'diff diff hours slope, male-female prog=1-prog=2' hours*female*prog 1 -1 0 -1 1 0, 'diff diff hours slope, male-female prog=1-prog=3' hours*female*prog 1 0 -1 -1 0 1, 'diff diff hours slope, male-female prog=2-prog=3' hours*female*prog 0 1 -1 0 -1 1 / e; run; *graph slopes by gender across programs; proc plm restore=catcatcon; effectplot slicefit (x=hours sliceby=female plotby=prog) / clm; run; *within each gender, do slopes of programs differ - overall test of conditional interactions of hours*prog within each gender; proc plm restore=catcatcon; estimate 'diff hours slope, male prog=1 - prog=2' hours*prog 1 -1 0 hours*female*prog 1 -1 0 0 0 0, 'diff hours slope, male prog=1 - prog=3' hours*prog 1 0 -1 hours*female*prog 1 0 -1 0 0 0, 'diff hours slope, male prog=2 - prog=3' hours*prog 0 1 -1 hours*female*prog 0 1 -1 0 0 0 / e joint; estimate 'diff hours slope, female prog=1 - prog=2' hours*prog 1 -1 0 hours*female*prog 0 0 0 1 -1 0, 'diff hours slope, female prog=1 - prog=3' hours*prog 1 0 -1 hours*female*prog 0 0 0 1 0 -1, 'diff hours slope, female prog=2 - prog=3' hours*prog 0 1 -1 hours*female*prog 0 0 0 0 1 -1 / e joint; run; *do difference in program slopes DIFFER between genders? Notice symmetry with above; proc plm restore=catcatcon; estimate 'diff diff hours slope, prog=1-prog=2, male-female' hours*female*prog 1 -1 0 -1 1 0, 'diff diff hours slope, prog=1-prog=3, male-female' hours*female*prog 1 0 -1 -1 0 1, 'diff diff hours slope, prog=2-prog=3, male-female' hours*female*prog 0 1 -1 0 -1 1 / e joint; run; *graph slopes by program across genders; proc plm restore=catcatcon; effectplot slicefit (x=hours sliceby=prog plotby=female) / clm; run; ***************3-way SIMPLE EFFECTS************; proc univariate data=exercise; var hours; run; *what are simple effects of gender in each program at important numbers of hours; proc plm restore=catcatcon; slice female*prog / sliceby=prog diff plots=none nof e means at hours=1.51; slice female*prog / sliceby=prog diff plots=none nof e means at hours=2; slice female*prog / sliceby=prog diff plots=none nof e means at hours=2.5; run; *what are simple effects of gender in each program at important numbers of hours; proc plm restore=catcatcon; lsmestimate female*prog 'male-female, prog=jogging(1) hours=1.51' [1, 1 1] [-1, 2 1], 'male-female, prog=swimming(2) hours=1.51' [1, 1 2] [-1, 2 2], 'male-female, prog=reading(3) hours=1.51' [1, 1 3] [-1, 2 3] / e adj=bon at hours=1.51; lsmestimate female*prog 'male-female, prog=jogging(1) hours=2' [1, 1 1] [-1, 2 1], 'male-female, prog=swimming(2) hours=2' [1, 1 2] [-1, 2 2], 'male-female, prog=reading(3) hours=2' [1, 1 3] [-1, 2 3] / e adj=bon at hours=2; lsmestimate female*prog 'male-female, prog=jogging(1) hours=2.5' [1, 1 1] [-1, 2 1], 'male-female, prog=swimming(2) hours=2.5' [1, 1 2] [-1, 2 2], 'male-female, prog=reading(3) hours=2.5' [1, 1 3] [-1 ,2 3] / e adj=bon at hours=2.5; run; *is the interaction of gender and program significant at each level of hours; proc plm restore=catcatcon; lsmestimate female*prog 'diff male-female, prog=1 - prog=2, hours=1.51' [1, 1 1] [-1, 2 1] [-1, 1 2] [1, 2 2], 'diff male-female, prog=1 - prog=3, hours=1.51' [1, 1 1] [-1, 2 1] [-1, 1 3] [1, 2 3], 'diff male-female, prog=2 - prog=3, hours=1.51' [1, 1 2] [-1, 2 2] [-1, 1 3] [1, 2 3] / e at hours=1.51 joint; lsmestimate female*prog 'diff male-female, prog=1 - prog=2, hours=2' [1, 1 1] [-1, 2 1] [-1, 1 2] [1, 2 2], 'diff male-female, prog=1 - prog=3, hours=2' [1, 1 1] [-1, 2 1] [-1, 1 3] [1, 2 3], 'diff male-female, prog=2 - prog=3, hours=2' [1, 1 2] [-1, 2 2] [-1, 1 3] [1, 2 3] / e at hours=2 joint; lsmestimate female*prog 'diff male-female, prog=1 - prog=2, hours=2.5' [1, 1 1] [-1, 2 1] [-1, 1 2] [1, 2 2], 'diff male-female, prog=1 - prog=3, hours=2.5' [1, 1 1] [-1, 2 1] [-1, 1 3] [1, 2 3], 'diff male-female, prog=2 - prog=3, hours=2.5' [1, 1 2] [-1, 2 2] [-1, 1 3] [1, 2 3] / e at hours=2.5 joint; run; *plot the interaction; proc plm restore=catcatcon; effectplot interaction (x=female sliceby=prog) / at(hours = 1.51 2 2.5) clm connect; run; ***********************LOGISTIC *********************; *first model; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours / expb; store logit; run; *odds ratio statement to get simple ORs and graph them; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours / expb; oddsratio hours / at(prog=all); store logit; run; *in general, it is hard to get continuous effects estimated as differences in probabilities, so we use ORs; *we ignore the issue of linearity of continuous predictors with the logit; proc plm restore=logit; estimate 'hours OR, prog=1' hours 1 hours*prog 1 0 0, 'hours OR, prog=2' hours 1 hours*prog 0 1 0, 'hours OR, prog=3' hours 1 hours*prog 0 0 1 / e exp cl; run; *take difference of slopes and exponentiate - ratio of odds ratios; proc plm restore=logit; estimate 'ratio hours OR, prog=1/prog=2' hours*prog 1 -1 0, 'ratio hours OR, prog=1/prog=3' hours*prog 1 0 -1, 'ratio hours OR, prog=2 /prog=3' hours*prog 0 1 -1 / e exp cl; run; ****simple effects; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours / expb; oddsratio prog / at(hours = 1.51 2 2.5); store logit; run; *easier to get simple effects expressed as differences in probabilities; proc plm source = logit; lsmeans prog / at hours=1.51 ilink plots=none; lsmeans prog / at hours=2 ilink plots=none; lsmeans prog / at hours=2.5 ilink plots=none; run; *graph predicted probabilities; proc plm restore=logit; effectplot interaction (x=prog) / at(hours = 1.51 2 2.5) clm; effectplot slicefit (x=hours sliceby=prog) / clm; run; *differences in probabilities; proc plm restore=logit; lsmeans prog / at hours=2; run; *get log odds slopes; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours; estimate 'hours slope, prog=1' hours 1 hours*prog 1 0 0; estimate 'hours slope, prog=2' hours 1 hours*prog 0 1 0; estimate 'hours slope, prog=3' hours 1 hours*prog 0 0 1; run; *compare log odds slopes; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours; estimate 'hours slope, prog=1 vs prog=2' hours*prog -1 1 0; estimate 'hours slope, prog=1 vs prog=3' hours*prog -1 0 1; estimate 'hours slope, prog=2 vs prog=3' hours*prog 0 -1 1; run; *get odds ratios; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours; oddsratio hours / at(prog=all) diff=all; estimate 'hours slope, prog=1' hours 1 hours*prog 1 0 0 / exp; estimate 'hours slope, prog=2' hours 1 hours*prog 0 1 0 / exp; estimate 'hours slope, prog=3' hours 1 hours*prog 0 0 1 / exp; run; *compare odds ratios; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours; estimate 'hours slope, prog=1 vs prog=2' hours*prog -1 1 0 / exp; estimate 'hours slope, prog=1 vs prog=3' hours*prog -1 0 1 / exp; estimate 'hours slope, prog=2 vs prog=3' hours*prog 0 -1 1 / exp; run; *pred prob plot; proc logistic data = exercise descending; class prog / param=glm order=internal; model satisfied = prog|hours; store logitcatcon; run; proc plm source = logitcatcon; effectplot slicefit (x=hours sliceby=prog) / limits; run; *pred prob plot through plot option on proc logistic - VERY SLOW; proc logistic data = exercise descending plots(only)=effect(x=hours sliceby=prog clband=yes); class prog / param=glm order=internal; model satisfied = prog|hours; run; ****** END SEMINAR CODE ******************************************************************* **************************************************************************** ***************************************************************;