/*--------------------------------------------------------------------*/ /* */ /* This macro is based on the article "The Analysis of Two-Factor */ /* Interactions in Fixed Effects Linear Models," by R. J. Boik, */ /* Journal of Educational Statistics, 1993, 18, 1--40. */ /* */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /* This iml program computes the F ratio for the maximal product */ /* interaction contrast as well as F ratios for partial interactions */ /* and product interaction contrasts that correspond to coefficient */ /* vectors input by the user. */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /*-----------------Required Arguments on First Use--------------------*/ /*--------------------------------------------------------------------*/ /* */ /* response = var, where var is a variable or a list of variables. */ /* */ /* data = dataset, where dataset is the name of the sas data set that */ /* contains the data. */ /* */ /* class = var, where var is a variable or a list of variables */ /* */ /* model_matrix = var, where var is a variable or a list */ /* of variables. If model_matrix is omitted, then */ /* model_matrix = ; is assumed. */ /* */ /* FactorA = classification variable that codes for Factor A */ /* */ /* FactorB = classification variable that codes for Factor B */ /* */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /*-----------Required Arguments on Subsequent Use---------------------*/ /*--------------------------------------------------------------------*/ /* */ /* None */ /* */ /*--------------------------------------------------------------------*/ /*--------------------------------------------------------------------*/ /*-------------------Optional Arguments-------------------------------*/ /*--------------------------------------------------------------------*/ /* */ /* smr = value, where value is either yes or no. If smr = yes, then */ /* maximal F for product contrats is computed. Default = yes. */ /* */ /* partialA = value, where value is an a x 1 vector of coefficients */ /* for Factor A. The partial interaction F test is computed */ /* and the b x 1 vector of simple effects of A at each level of */ /* B is computed. */ /* */ /* labelA = value, where value is a label (title) to be printed for */ /* the partial interaction test that corresponds to partialA. */ /* */ /* partialB = value, where value is a b x 1 vector of coefficients */ /* for Factor B. The partial interaction F test is computed */ /* and the a x 1 vector of simple effects of B at each level of */ /* A is computed. */ /* */ /* labelB = value, where value is a label (title) to be printed for */ /* the partial interaction test that corresponds to partialB. */ /* */ /* productA = value, where value is an a x 1 vector of coefficients */ /* for Factor A. The product interaction corresponding to */ /* productA and productB is computed. */ /* */ /* productB = value, where value is a b x 1 vector of coefficients */ /* for Factor B. The product interaction corresponding to */ /* productA and productB is computed. */ /* */ /* labelAB = value, where value is a label for the product contrast. */ /* */ /* glmoptions = value, where value = noprint. Use this argument if */ /* you do not want printed output from proc glm as called in the */ /* macro. The default is print. */ /* */ /* rerun = value, where value = yes or no. If value = yes, then */ /* lsmeans will be recomputed even though a file containing */ /* lsmeans already exists. Default is no. */ /* Caution: if the model has changed, then rerun=yes is */ /* necessary in order to obtain correct analyses. */ /* */ /* nameA = label for Factor A, default is "Factor A" */ /* */ /* nameB = label for Factor B, default is "Factor B" */ /* */ /*--------------------------------------------------------------------*/ %macro twoway_interaction(response = , model_matrix = , class = , factorA = , factorB = , data = ,smr = yes, partialA = 1, labelA = , partialB = 1, labelB = , productA = 1, productB = 1, labelAB = , glmoptions = ,rerun = no, nameA = Factor A, nameB = Factor B); %if %sysfunc(exist(means_smr)) = 0 | &rerun = yes %then %do; proc glm data = &data &glmoptions outstat=temp_stat_out; class &factorB &factorA &class; model &response = &model_matrix; lsmeans &factorA*&factorB/ cov out=means_smr; run; %end; proc IML; Start als; wp=inv((wq @ I(p))`*Phi*(wq @ I(p)))*(wq @ I(p))`*psi; wq=inv((I(q) @ wp)`*Phi*(I(q) @ wp))*(I(q) @ wp)`*psi; epsi=psi`*(wq @ wp)-R; R = R + epsi; finish; start iterate; do while(epsi >=.00001); run als; end; finish; use temp_stat_out; reset noname; read all var _num_ into Z; dfe=Z[1,1]; SSE=Z[1,2]; MSE=SSE/dfe; use means_smr; reset noname; read all var _num_ into X; a=ncol(design(X[,2])); b=ncol(design(X[,1])); p=a-1; q=b-1; Sigma = X[,6:5+a*b]; mu=X[,3]; Ha=I(a)-J(a,a,1/a); Ka=Ha[,1:p]; Hb=I(b)-J(b,b,1/b); Kb=Hb[,1:q]; *; * Caution: transpose[shape(vector,b,a)] = dvec(vector,a,b); * Caution: shape(Matrix, a*b,1) = vec[transpose(Matrix)]; *; M=shape(mu,b,a)`; nlength=%length(&partialA); vecCa = J(nlength,1,0); %let j1=1; %do %while (%qscan(&partialA,%eval(&j1),' ') ^= ); jj=%eval(&j1); vecCa[jj]=%qscan(&partialA,%eval(&j1),' '); %let j1 = %eval(&j1+1); %end; nA = %eval(&j1)-1; if nA <=1 then do; Ca={0}; nA=0; end; else Ca=shape(vecCa[1:nA],nA/a,a)`; nlength=%length(&partialB); vecCb = J(nlength,1,0); %let j1=1; %do %while (%qscan(&partialB,%eval(&j1),' ') ^= ); jj=%eval(&j1); vecCb[jj]=%qscan(&partialB,%eval(&j1),' '); %let j1 = %eval(&j1+1); %end; nB = %eval(&j1)-1; if nB <=1 then do; Cb={0}; nB=0; end; else Cb=shape(vecCb[1:nB],nB/b,b)`; nlength=%length(&productA); vecproductA = J(nlength,1,0); %let j1=1; %do %while (%qscan(&productA,%eval(&j1),' ') ^= ); jj=%eval(&j1); vecproductA[jj]=%qscan(&productA,%eval(&j1),' '); %let j1 = %eval(&j1+1); %end; nABA = %eval(&j1)-1; if nABA <=1 then do; productA = {0}; nABA = 0; end; else productA=shape(vecproductA[1:nABA],nABA/a,a)`; nlength=%length(&productB); vecproductB = J(nlength,1,0); %let j1=1; %do %while (%qscan(&productB,%eval(&j1),' ') ^= ); jj=%eval(&j1); vecproductB[jj]=%qscan(&productB,%eval(&j1),' '); %let j1 = %eval(&j1+1); %end; nABB = %eval(&j1)-1; if nABB <=1 then do; productB={0}; nABB = 0; end; else productB=shape(vecproductB[1:nABB],nABB/b,b)`; %if &smr = yes %then %do; Phi = (Kb @ Ka)`*Sigma*(Kb @ Ka); Psi = Ka`*M*Kb; call svd(U,D,V,Psi); wq=V[,1]; psi=shape(Psi`,p*q,1); epsi = 1; R = 0; run iterate; Va=(J(b,1,1) @ Ka)`*Sigma*(J(b,1,1) @ Ka); cca=Ka*inv(Va)*(J(b,1,1) @ Ka)`*mu; SS=mu`*(J(b,1,1) @ Ka)*inv(Va)*(J(b,1,1) @ Ka)`*mu; print "Maximal F Ratio for &nameA Contrast",, 'SS for Main Effect Contrast =' (SS*MSE),, 'F = SS/MSE =' SS,; print "Contrast Coeff of Maximal &nameA Contrast" (cca/sqrt(cca`*cca)); print '---------------------------------------------------'; Vb=(Kb @ J(a,1,1))`*Sigma*(Kb @ J(a,1,1)); ccb=Kb*inv(Vb)*(Kb @ J(a,1,1))`*mu; SS=mu`*(Kb @ J(a,1,1))*inv(Vb)*(Kb @ J(a,1,1))`*mu; print "Maximal F Ratio for &nameB Contrast",, 'SS for Main Effect Contrast =' (SS*MSE),, 'F = SS/MSE =' SS,; print "Contrast Coeff of Maximal &nameB Contrast" (ccb/sqrt(ccb`*ccb)); print '---------------------------------------------------'; print "Maximal F Ratio for Product Contrast",, 'SS for Product Interaction Contrast =' (R*MSE),, 'F = SS/MSE =' R,, 'p, q, dfe:' p q dfe,; print "&nameA Coeff of Maximal Product Contrast:" (Ka*wp/sqrt(wp`*Ka`*Ka*wp)),; print "&nameB Coeff of Maximal Product Contrast:" (Kb*wq/sqrt(wq`*Kb`*Kb*wq)); print '---------------------------------------------------'; %end; if nA > 0 then do; Phi=(Kb @ Ca)`*Sigma*(Kb @ Ca); psi=(Kb @ Ca)`*mu; T=psi`*inv(Phi)*psi; F=T/q; tmpa=Ca`; oneb=J(b,1,1)/b; MpCa= M`*Ca; main_effect=oneb`*MpCa; se_main_effect=sqrt((oneb @ Ca)`*Sigma*(oneb @ Ca)); ta=(main_effect/se_main_effect); ssa=MSE*ta##2; print "&labelA",, "c_A' =" tmpa,, 'Estimate of Main Effect Contrast =' main_effect, 'SS for Main Effect Contrast =' ssa, 'SE of Main Effect Contrast Estimator =' se_main_effect, 't = Estimate/SE for Main Effect Contrast =' ta,, 'SS for Partial Interaction = ' (T*MSE), 'Q = SS/MSE for Partial Interaction = ' T, 'F = MS/MSE for Partial Interaction = ' F, 'q, dfe:' q dfe; one=char(1,3); namerow= concat("&nameB"," level",one); seMpCa=(I(b) @ Ca)`*Sigma*(I(b) @ Ca); outM=J(b,2,0); do i=1 to b; b_level=char(i,3); if i > 1 then namerow=namerow//concat("&nameB"," level",b_level); outM[i,1]=MpCa[i]; outM[i,2]=sqrt(seMpCa[i,i]); end; print "Simple Effects of &nameA at Levels of &nameB"; mattrib outM rowname = namerow colname = {'Estimate','St. Error'} label = ' '; print outM; print '---------------------------------------------------'; end; if nB > 0 then do; Phi=(Cb @ Ka)`*Sigma*(Cb @ Ka); psi=(Cb @ Ka)`*mu; T=psi`*inv(Phi)*psi; F=T/p; tmpb=Cb`; onea=J(a,1,1)/a; MCb= M*Cb; main_effect=onea`*MCb; se_main_effect=sqrt((Cb @ onea)`*Sigma*(Cb @ onea)); tb=(main_effect/se_main_effect); ssb=MSE*tb##2; print "&labelB",, "c_B' =" tmpb,, 'Estimate of Main Effect Contrast =' main_effect, 'SS for Main Effect Contrast =' ssb, 'SE of Main Effect Contrast Estimator =' se_main_effect, 't = Estimate/SE for Main Effect Contrast =' tb,, 'SS for Partial Interaction = ' (T*MSE), 'Q = SS/MSE for Partial Interaction = ' T, 'F = MS/MSE for Partial Interaction = ' F, 'p, dfe:' p dfe; one=char(1,3); namerow= concat("&nameA"," level",one); seMCb=(Cb @ I(a))`*Sigma*(Cb @ I(a)); outM=J(a,2,0); do i=1 to a; a_level=char(i,3); if i > 1 then namerow=namerow//concat("&nameA"," level",a_level); outM[i,1]=MCb[i]; outM[i,2]=sqrt(seMCb[i,i]); end; print "Simple Effects of &nameB at Levels of &nameA"; mattrib outM rowname = namerow colname = {'Estimate','St. Error'} label = ' '; print outM; print '---------------------------------------------------'; end; if nABA*nABB > 0 then do; Phi=(productB @ productA)`*Sigma*(productB @ productA); psi=(productB @ productA)`*mu; T=psi`*inv(Phi)*psi; se=sqrt(Phi); tmpa=productA`; tmpb=productB`; print "&labelAB",, 'c_A` =' tmpa,, 'c_B` =' tmpb,, 'Estimate of Product Interaction Contrast = ' psi,, 'Standard Error = ' se,, 'SS for Product Contrast =' (T*MSE),, 'F = SS/MSE for Product Interaction = ' T,; print '---------------------------------------------------'; end; quit; %mend twoway_interaction;