Programming the statistical procedures from SAS

Estimating predicted probability using Proc Glimmix for set parameters

Accepted Solution Solved
Reply
Contributor
Posts: 70
Accepted Solution

Estimating predicted probability using Proc Glimmix for set parameters

Hello, I'm new at this, so not sure if I'm explaining this correctly, but I'm trying to estimate predicted probably of a doctor having enough interaction with a patient based on predictors having a set value.  So probabilty Doctor 1, Doctor 2, Doctor 3, etc  has enough interaction with the patient if the patient is age 18-34, female, with a college degree or higher.  I tried using this code below using "estimate", but it doesn't seem to be working.  Anyone know what I might be doing wrong?

 

proc glimmix data=tables2 method=laplace noclprint;
class attending;
model doc_interaction (descending) = age18_34 age35_54 gender_m edu_hs_m edu_some_m edu_col_M
/CL DIST=Binary link=logit solution;
estimate "Intercept Age 18-34, female, college grad or higher" 
int 1 age18_34 1 gender_m 1 edu_col_m 1 / ILINK; 
random intercept /  subject=attending S CL Type=VC ;
COVTEST / wald;
run;

 


Accepted Solutions
Solution
‎07-10-2017 03:24 PM
SAS Employee
Posts: 11

Re: Estimating predicted probability using Proc Glimmix for set parameters

[ Edited ]

You will need one ESTIMATE statement for each attending you have in your data.  For the second level of attending, the ESTIMATE statement might look something like

 

estimate '2' int 1 x1 1 x2 1 | int 1 / subject 0 1 ilink e; 

 

You include the random effects you need in your ESTIMATE statement after the | operator (the vertical bar operator).  Since we also have a SUBJECT= effect in the RANDOM statement, you specify the subject of interest for the RANDOM statement with the SUBJECT option.  Treat this effect like you would any other CLASS effect in your model.

 

If you have a lot of levels to the attending effect, this could get rather tedious.  You will need one ESTIMATE statement for each of those levels.  However, SAS Macro code can make writing all of these ESTIMATE statements a snap.  Try something like

 

%do i=1 %to 10;
estimate "&i" int 1 x1 1 x2 1| int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;

 

That macro loop will generate the estimates you need for the first ten attending levels.  Just change the 10 on that %do line to generate more ESTIMATE statements.  You will need to put your PROC GLIMMIX in a macro to make this work.

 

I have ammended the code below to include this suggestion and to include Rick's suggestion just below my answer.

 

 

 

data test;
call streaminit(6354153); /* initialize random number stream */
do a=1 to 10; /* 10 levels to the random effect */
ra=rand('normal',.7); /* generate random intercept adjustment for each level */
do i=1 to 20;
x1=rand('uniform'); /* generate random covariates */
x2=rand('uniform');
x3=rand('normal')*.5;
x4=rand('normal')*.5;
trt=ceil(rand('uniform')*3); /* generate trt and group levels */
grp=ceil(rand('uniform')*5);

logit=-2 + 2*x1 + x2 + x3 - x4 + (trt-2) + (grp-3) + ra; /* calculate the logit */
p=exp(-logit)/(1+exp(-logit)); /* and then the predicted probability */
if rand('uniform')>p then y=1; else y=0; /* assign the response */
output;
end;
end;
run;

proc glimmix data=test method=laplace;
class a trt grp;
model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
random intercept / subject=a;
estimate '1' int 1 x1 1 x2 1 | int 1 / subject 1 ilink e;
estimate '2' int 1 x1 1 x2 1 | int 1 / subject 0 1 ilink e;
run;

 

%macro doests;

%do i=1 %to 10;
estimate "&i" int 1 x1 1 x2 1| int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;

%mend;


proc glimmix data=test method=laplace;
class a trt grp;
model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
random intercept / subject=a;
%doests;
run;

View solution in original post


All Replies
SAS Employee
Posts: 11

Re: Estimating predicted probability using Proc Glimmix for set parameters

Are you getting an error or warning message from the GLIMMIX procedure?  The simulation below does something similar to your model, with the addition of some fixed effects.  The MEAN column in the ESTIMATE statement output gives you the predicted probability for the variables you have specified (and averaged over any additional fixed effects in the model).

 

data test;
call streaminit(6354153); /* initialize random number stream */
do a=1 to 10; /* 10 levels to the random effect */
ra=rand('normal',.7); /* generate random intercept adjustment for each level */
do i=1 to 20;
x1=rand('uniform'); /* generate random covariates */
x2=rand('uniform');
x3=rand('normal')*.5;
x4=rand('normal')*.5;
trt=ceil(rand('uniform')*3); /* generate trt and group levels */
grp=ceil(rand('uniform')*5);

logit=-2 + 2*x1 + x2 + x3 - x4 + (trt-2) + (grp-3) + ra; /* calculate the logit */
p=exp(-logit)/(1+exp(-logit)); /* and then the predicted probability */
if rand('uniform')>p then y=1; else y=0; /* assign the response */
output;
end;
end;
run;

proc glimmix data=test method=laplace;
class a trt grp;
model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
random intercept / subject=a;
estimate '1' int 1 x1 1 x2 1 / ilink e;
run;

Contributor
Posts: 70

Re: Estimating predicted probability using Proc Glimmix for set parameters

Nope, I don't get an error.  It seems to run fine, but I'm not getting the output I expect.  

 

Is there a way to get at an estimate for each physician ("attending" in my dataset) as opposed to just one estimate in the ESTIMATE statement output, so that I can get at a predicted probability for each physician for the variables I have specified?  

Solution
‎07-10-2017 03:24 PM
SAS Employee
Posts: 11

Re: Estimating predicted probability using Proc Glimmix for set parameters

[ Edited ]

You will need one ESTIMATE statement for each attending you have in your data.  For the second level of attending, the ESTIMATE statement might look something like

 

estimate '2' int 1 x1 1 x2 1 | int 1 / subject 0 1 ilink e; 

 

You include the random effects you need in your ESTIMATE statement after the | operator (the vertical bar operator).  Since we also have a SUBJECT= effect in the RANDOM statement, you specify the subject of interest for the RANDOM statement with the SUBJECT option.  Treat this effect like you would any other CLASS effect in your model.

 

If you have a lot of levels to the attending effect, this could get rather tedious.  You will need one ESTIMATE statement for each of those levels.  However, SAS Macro code can make writing all of these ESTIMATE statements a snap.  Try something like

 

%do i=1 %to 10;
estimate "&i" int 1 x1 1 x2 1| int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;

 

That macro loop will generate the estimates you need for the first ten attending levels.  Just change the 10 on that %do line to generate more ESTIMATE statements.  You will need to put your PROC GLIMMIX in a macro to make this work.

 

I have ammended the code below to include this suggestion and to include Rick's suggestion just below my answer.

 

 

 

data test;
call streaminit(6354153); /* initialize random number stream */
do a=1 to 10; /* 10 levels to the random effect */
ra=rand('normal',.7); /* generate random intercept adjustment for each level */
do i=1 to 20;
x1=rand('uniform'); /* generate random covariates */
x2=rand('uniform');
x3=rand('normal')*.5;
x4=rand('normal')*.5;
trt=ceil(rand('uniform')*3); /* generate trt and group levels */
grp=ceil(rand('uniform')*5);

logit=-2 + 2*x1 + x2 + x3 - x4 + (trt-2) + (grp-3) + ra; /* calculate the logit */
p=exp(-logit)/(1+exp(-logit)); /* and then the predicted probability */
if rand('uniform')>p then y=1; else y=0; /* assign the response */
output;
end;
end;
run;

proc glimmix data=test method=laplace;
class a trt grp;
model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
random intercept / subject=a;
estimate '1' int 1 x1 1 x2 1 | int 1 / subject 1 ilink e;
estimate '2' int 1 x1 1 x2 1 | int 1 / subject 0 1 ilink e;
run;

 

%macro doests;

%do i=1 %to 10;
estimate "&i" int 1 x1 1 x2 1| int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;

%mend;


proc glimmix data=test method=laplace;
class a trt grp;
model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
random intercept / subject=a;
%doests;
run;

Contributor
Posts: 70

Re: Estimating predicted probability using Proc Glimmix for set parameters

Oh my goodness - it worked!  Thank you so much!

SAS Super FREQ
Posts: 3,753

Re: Estimating predicted probability using Proc Glimmix for set parameters

Great answer by StatsMan. I'll merely add that you don't need the whole GLIMMIX call in a macro, only the ESTIMATE statements:

 

%macro doEst;
%do i=1 %to 10;
 estimate "&i" int 1 x1 1 x2 1| int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;
%mend;

proc glimmix data=test method=laplace;
 class a trt grp;
 model y=trt grp x1 x2 x3 x4 / s dist=bin link=logit;
 random intercept / subject=a;
%doEst;
run;
Contributor
Posts: 70

Re: Estimating predicted probability using Proc Glimmix for set parameters

Just tried this and it worked great!  Thank you!

 

But I also just noticed that the labels of the providers don't actually match the attending codes I have.  Is there a way to fix this?  So my Labels are Attending 1, Attending 2, etc. but the attending codes I have starts with 2 and skips around.  Its starting to get confusing.

 

Also, is there a way to see the confidence intervals without me having to manually calculate them in the estimates section?  Is there an option for that?

 

Thanks! 

SAS Employee
Posts: 11

Re: Estimating predicted probability using Proc Glimmix for set parameters

The CL option on the ESTIMATE statement does not provide a confidence interval on the inverse link scale.  You would have to manually calculate that result, and those confidence bands tend to be pretty wide.

 

Are you talking about the labels that are written to the ESTIMATE statement results themselves?  Using the macro code, yes, that just gives the numbers 1,2,3, etc. to the estimates in order.

 

The ordering of the subjects is either given in the output of the CLASS levels table (if the SUBJECT= effect is a CLASS variable) or is the order the subjects appear in your input data set if the SUBJECT= effect does not appear on the CLASS statement.  You can create a SAS data set that has the ordered subject levels in either scenario and create macro variables from this set of subject values.  The macro variables can then be used in the label for the ESTIMATE statements.  I do not have code handy for this task, however.  Perhaps someone else has that available?

Contributor
Posts: 70

Re: Estimating predicted probability using Proc Glimmix for set parameters

Thanks for the reply. I guess I'll just have to live with manually
calculating the CIs and the labels that are written to the estimate
statement. Thanks again for all your help!


##- Please type your reply above this line. Simple formatting, no
attachments. -##
Contributor
Posts: 70

Re: Estimating predicted probability using Proc Glimmix for set parameters

I just tried to get at the predicted means using the output statement, but it doesn't appear to take into account the "estimate" statements I added in.  Is there a way to do this?  Or do I need to somehow create an output statement for each attending using a where statement to get at the predicted means for each attending?

 

%macro doEst;
%do i=1 %to 34;
estimate "Attending &i" int 1 age18_34 1 gender_m 1 | int 1 / subject %do k=1 %to %eval(&i-1); 0 %end; 1 ilink e;
%end;
%mend;

proc glimmix data=tables method=quad(qpoints=31) noclprint;
class attending;
model contact (descending) = age18_34 age35_54 gender_m 

/CL DIST=Binary link=logit solution;

random intercept /  subject=attending S CL Type=VC ;

COVTEST / wald;

output out=glm_out pred(ilink) lcl(ilink) ucl(ilink) stderr(ilink);

%doest;
run;
☑ This topic is solved.

Need further help from the community? Please ask a new question.

Discussion stats
  • 9 replies
  • 341 views
  • 7 likes
  • 3 in conversation