EVERYTHING INSIDE /* AND */ IS TAKEN AS A COMMENT in SAS programs These are shown in green here /*--------------------------------------------------- ex5wt.sas G Raab /Helen Storkey Nov/Dec 04 Some SAS graphs are plotted by the code here but are not shown in the output file First set up libname to tell sas where data can be found Use yout own directory name --------------------------------------------------------*/ libname exemp5 'C:\Documents and Settings\gillian raab\My Documents\aprojects\peas\ex5datafiles\data'; proc means data=exemp5.ex5wt min max mean sum; * get population and respondent totals as a check; var npop ayra_num ; run; back to top
/*------------------------------- non-response modelling-----------------------
commented out lines show earlier stages of model building
-------------------------------------------------------------------------------*/
proc logistic data=exemp5.ex5wt;
class agegrp gender;
*model ayra_num/npop = agegrp gender agegrp*gender sacc sinc semp shlth sed shse;
***************** only sacc sinc shse drop others one by one;
*model ayra_num/npop = agegrp gender agegrp*gender sacc sinc semp shlth shse;
*model ayra_num/npop = agegrp gender agegrp*gender sacc sinc shlth shse;
*model ayra_num/npop = agegrp gender agegrp*gender sacc sinc shse;
**************** now try interactions of sinc (most important predictor) with age and gender;
*model ayra_num/npop = agegrp gender agegrp*gender*sinc sacc sinc shse;
*model ayra_num/npop = agegrp gender agegrp*gender agegrp*sinc gender*sinc sacc sinc shse;
*model ayra_num/npop = agegrp gender agegrp*gender agegrp*sinc sacc sinc shse;
****************** housing can now be dropped to give final model;
model ayra_num/npop = agegrp gender agegrp*gender agegrp*sinc sacc sinc ;
run;
/*--------------------------------------------------------------
now check linearity of important term sinc by setting it as a grouped variable
and also putting in a quadratic term
------------------------------------------------------------------*/
data new; set exemp5.ex5wt;
sincl=sinc;
sinc2=sinc**2;
run;
proc univariate data=aapop;* get its percentiles in data zone distribution;
output out=temp pctlpre=P_ pctlpts=10 to 100 by 10;
var sinc ; run;
* this divides sinc into 5 groups with same no of data zones in each;
proc format; value sinc low-<5.8=1 5.8-<12.4=2 12.4-<18.3=3 18.3-<26.1=4 26.1-high=5;
proc logistic data=new;
class agegrp gender sinc;
format sinc sinc.;
model ayra_num/npop = agegrp gender agegrp*gender agegrp*sincl sacc sinc sincl;
* grouped model did not improve on linear one;
run;
proc logistic data=new;
class agegrp gender ;
format sinc sinc.;
model ayra_num/npop = agegrp gender agegrp*gender agegrp*sincl sacc sinc sinc2;
* neither did quadratic term;
run;
/*----------------------------------------------------------------------
now make data set to add in to examine predictions at all possible combinations
of agegroup sex and a point in each of the 5 groups for sinc
sacc is held at its mean value because it is of minor importance
------------------------------------------------------------------------*/
data predict (drop=i j);* set predictions at p10..p30...tp p90 (mid points of 5 groups);
do agegrp=1 to 13;
do i= 1 to 2;
if i=1 then gender='M';
if i=2 then gender='F';
do j=1 to 5;
if j=1 then sinc=32.3;
if j=2 then sinc=21.8;
if j=3 then sinc=15.7;
if j=4 then sinc=9;
if j=5 then sinc=3.3;
sacc=0;
shse=15;
output;
end;
end;
end;
run;
data new;
set exemp5.ex5wt predict;
run;
proc logistic data=new; * rerun to get predictions from final model above;
class agegrp gender ;
*rerunning model on data set with prediction cases;
model ayra_num/npop = agegrp gender agegrp*gender agegrp*sinc sacc sinc ;
output out=pred p=phat;
run;
data pred;* just keep the ones for prediction with missing outcome;
set pred;
if npop=.;
keep agegrp gender sinc phat ;
format gender 1.0 agegrp 2.0;
run;
* now make plots of predicted values;
proc sort data=pred; by gender;run;
proc gplot data=pred;
goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate
hpos=0 vpos=0 htext= ftext= ctext= target= gaccess= gsfmode= ;
goptions device=WIN ctext=blue graphrc interpol=join
colors=(black red green blue cyan gold brown purple olive orange grey violet);
axis1 order=(0 to 0.035 by 0.001) color=blue width=2.0 ;
axis2 color=blue order=(0 to 13 by 1) width=2.0 ;* excludes outside this range;
symbol1 value=dot colour=;
plot phat*agegrp=sinc/haxis=axis2 vaxis=axis1;
by gender;
run;
run;
back to top
/*--------------- now try quadratic function of age instead of groups---------------*/
proc logistic data=new;
class gender ;
format sinc sinc.;
model ayra_num/npop = agegrp gender agegrp*gender agegrp*sinc
agegrp*agegrp agegrp*agegrp*gender agegrp*agegrp*sinc sacc sinc ;
output out=pred p=phat;
run;
data pred2;
set pred;
if npop=.; * just the ones for prediction;
keep agegrp gender sinc phat ;
format gender 1.0 agegrp 2.0;
run;
*----------------- and plot these too;
proc sort data=pred2; by gender; run;
proc gplot data=pred2;
goptions reset=(axis, legend, pattern, symbol, title, footnote) norotate
hpos=0 vpos=0 htext= ftext= ctext= target= gaccess= gsfmode= ;
goptions device=WIN ctext=blue graphrc interpol=join
colors=(black red green blue cyan gold brown purple olive orange grey violet);
axis1 order=(0 to 0.035 by 0.001) color=blue width=2.0 ;
axis2 color=blue order=(0 to 13 by 1) width=2.0 ;* excludes outside this range;
symbol1 value=dot colour=;
plot phat*agegrp=sinc/haxis=axis2 vaxis=axis1;
by gender;
run;
run;
back to top
*----------------------- these look OK will go with this;
/*---------------------- now calculate weights only need those with any in sample
ready to merge with data file-----------------------------*/
data predx;
set pred;
if ayra_num>0 and npop^=.; * just keep those in sample;
gweight=1/phat; *grossing weight;
avweight=gweight/67.76*4165/4147; * rescale to mean 1;
keep dz ayra_num npop dz agegrp gender avweight gweight;
run;
proc means data=predx mean max min sum;* get mean weight for sample;
weight ayra_num ;
var avweight gweight ;
run;
/*-------------------------------------------------
next step would be to merge weights back into
the data file
------------------------------------------------*/