topic Re: Estimating predicted probability using Proc Glimmix for set parameters in Statistical Procedures

Estimating predicted probability using Proc Glimmix for set parameters

einstein — Mon, 10 Jul 2017 14:34:45 GMT

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;

Re: Estimating predicted probability using Proc Glimmix for set parameters

StatsMan — Mon, 10 Jul 2017 18:43:39 GMT

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;

Re: Estimating predicted probability using Proc Glimmix for set parameters

einstein — Mon, 10 Jul 2017 18:54:48 GMT

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?

Re: Estimating predicted probability using Proc Glimmix for set parameters

StatsMan — Tue, 11 Jul 2017 11:36:35 GMT

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.

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;

Re: Estimating predicted probability using Proc Glimmix for set parameters

einstein — Mon, 10 Jul 2017 19:24:33 GMT

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

Re: Estimating predicted probability using Proc Glimmix for set parameters

Rick_SAS — Mon, 10 Jul 2017 22:48:20 GMT

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;

Re: Estimating predicted probability using Proc Glimmix for set parameters

einstein — Tue, 11 Jul 2017 17:13:54 GMT

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!

Re: Estimating predicted probability using Proc Glimmix for set parameters

StatsMan — Tue, 11 Jul 2017 19:12:58 GMT

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?

Re: Estimating predicted probability using Proc Glimmix for set parameters

einstein — Tue, 11 Jul 2017 20:19:39 GMT

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

Re: Estimating predicted probability using Proc Glimmix for set parameters

einstein — Fri, 11 Aug 2017 20:41:05 GMT

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;