Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

- Home
- /
- Analytics
- /
- Stat Procs
- /
- Estimating predicted probability using Proc Glimm...

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

07-10-2017 10:34 AM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to einstein

07-10-2017 03:13 PM - edited 07-11-2017 07:36 AM

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;

All Replies

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to einstein

07-10-2017 02:43 PM

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;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to StatsMan

07-10-2017 02:54 PM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to einstein

07-10-2017 03:13 PM - edited 07-11-2017 07:36 AM

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;

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;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to StatsMan

07-10-2017 03:24 PM

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

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to StatsMan

07-10-2017 06:48 PM

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;
```

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

07-11-2017 01:13 PM

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!

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to einstein

07-11-2017 03:12 PM

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?

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to StatsMan

07-11-2017 04:19 PM

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. -##

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. -##

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

08-11-2017 04:41 PM

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;