Calculating variance of differences in Proc Glimmix on the natural scale - Proc IML the only way?

[dandar]

I am interested in being able to calculate differences in means of groups (i.e. treatment - control) on the natural scale from within Proc Glimmix. For instance the snippet below from the help document for the "lsmestimate" statement suggests only the difference on the linear predictor scale can be computed. I am in fact fitting a repeated measures model to count data so am interested in the log link, however the principle is the same.

I could work out the variance of the difference on the natural scale, and this would (I think) be a function of the covariance matrix for the fixed effects. I could I suppose therefore output this estimated covariance matrix, and fixed effects estimates, read it all into IML, and compute the variance of this difference "manually". However....if there is an easier way from within Proc Glimmix?

Any help would be greatly appreciated.

[PaigeMiller]

Does "count data" mean Poisson distribution? Or does it mean some other distribution? Poisson is built into PROC GLIMMIX and further PROC GLIMMIX allows you to define your own link or variance functions.

I wouldn't bother with LSMESTIMATE if you just want differences between means of the groups, the LSMEANS statement will do that for you, and provide variance estimates on the original scale or transformed scale.

[dandar]

Thanks for your quick response. Sorry yes my data are Poisson distributed. I know the common quantity computed are ratios of group means (so the other covariate effects cancel out under the log link), but clinical folks like to see differences. For this reason I have for a long time wanted to be able to at least have the ability to compute these differences and the variance so I can easily construct an asymptotic Normal CI on it. 

You think this can be done with a custom variance function? 

[PaigeMiller]

Poisson is built into PROC GLIMMIX, you don't need to create your own variance function.

The LSMEANS statement makes use of whatever distribution you specify and provides the proper variance estimates.

[dandar]

Ahh! So lsmeans with the "diff" option computes the difference on the natural scale with the correct variance - i.e. it calculates exp(eta1)-exp(eta2) rather than exp[eta1-eta2]? I think I have had misconceptions about lsmeans then for a long time! I always tend to use the "estimate" statement, and this calculates exp[eta1-eta2] rather like lsmestimate does. Thus I convinced myself lsmeans would do the same (or I forgot it can do better). 

[PaigeMiller]

LSMEANS makes all of this easy.

There is rarely (in my experience) a need for LSMESTIMATE unless you want to compare levels where the linear combination is something like 0.25*level 1 mean -0.6*level 2 mean + 1,4*level 3 mean

[dandar]

Thanks so much for your help.

[StatDave · Accepted solution]

No. LSMEANS/ILINK DIFF does NOT compute the difference on the mean scale. In your case, it computes exp[eta1-eta2]. To compute the difference on the mean scale, you can use the NLMeans macro. An example using count data is shown in this note. 

[dandar]

That indeed did the trick! As I stated earlier, I always thought lsmeans did not do what I needed. Thanks! I am interested in finding out if the NLmixed method is equivalent to what I was going to try in proc IML, but for now I am happy this method works.

The code was easy to implement (but time-consuming since I have 35 fixed effects). I chose to supply NLestimate macro with the fixed effects estimates and fixed effects covariance matrix estimate separately that I saved from proc glimmix (as opposed to using the "store" function). I then created a "score" dataset which seems analogous to the "OM" dataset in lsmeans which serves the purpose of supplying all the covariate values need that enter the function expression (I will have multiple rows for multiple comparisons of active doses to placebo). Thus I ended up with;

exp(b_c1*x_c1 +...b_cn_p*x_cn_p + b_trt_1*x_trt_1 + ... + b_trt_q*x_trt_q) - 

    exp(b_c1*x_c1 +...b_cn_p*x_cn_p + b_trt_1*x_trt2_1 + ... + b_trt_q*x2_trt2_q) 

= exp(b_c1*x_c1 +...b_cn_p*x_cn_p) *

    { exp(b_trt_1*x_trt_1 + ... + b_trt_q*x_trt_q )  - exp(b_trt2_1*x_trt2_1 + ... + b_trt2_q*x_trt2_q ) }

where "c" denotes effects common to both terms, and 

x_trt_1=0    x_trt_2 = 1 ...... x_trt_q = 0         (defines first treatment group of choice)

x_trt2_1=1    x_trt2_2 = 0 ...... x_trt2_q = -1         (defines second treatment group of choice)

so that the code above defines the mean for treatment group 2 - mean for treatment group q, at the values defined in the common covariates. By creating more rows in the score dataset but changing just the above treatment group dummy variables, each row generates an estimate of an active group - placebo for the same covariates.

/*proc glimixx code with manually coded estimate statements

The NLmixed macro is subsequently used to estimate differences between
the active dose means and the placebo means
*/

ods output estimates=glmm_p_est covb=glmm_p_FE_cov parameterestimates=glmm_p_FE_est ;
proc glimmix data=Model_dat_poiss order=formatted plots=all  method=laplace;
class usubjid trtp avisitn randstra;
model aval = ady trtp avisitn trtp*avisitn base base*ady randstra / dist=poisson link=log solution cl covb ;
random int exposure/ subject=usubjid type=UN  g gcorr ;
output outpred=glmm_p_out pred(blup ilink)=predp lcl(blup ilink)=lowerp ucl(blup ilink)=upperp pred(noblup ilink)=predm lcl(noblup ilink)=lowerm ucl(noblup ilink)=upperm ;
covtest glm/ cl;
lsmeans trtp / e ilink cl ;
estimate 'Severe base: PBO week 1'    intercept 1 randstra 1 0 base 21.75  ady 7 base*ady 152.25  trtp 0 0 0 0 1  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0/ilink cl  ;
estimate 'Severe base: PBO week 4'    intercept 1 randstra 1 0 base 21.75  ady 28 base*ady 609  trtp 0 0 0 0 1  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0/ilink cl ;
estimate 'Severe base: PBO week 8'    intercept 1 randstra 1 0 base 21.75  ady 56 base*ady 1218  trtp 0 0 0 0 1  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0/ilink cl ;
estimate 'Severe base: PBO week 12'    intercept 1 randstra 1 0 base 21.75  ady 84 base*ady 1827  trtp 0 0 0 0 1  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1/ilink cl ;

estimate 'Severe base: 37.5 week 1'    intercept 1 randstra 1 0 base 21.75  ady 7 base*ady 152.25  trtp 0 0 1 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 37.5 week 4'    intercept 1 randstra 1 0 base 21.75  ady 28 base*ady 609  trtp 0 0 1 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 37.5 week 8'    intercept 1 randstra 1 0 base 21.75  ady 56 base*ady 1218  trtp 0 0 1 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 37.5 week 12'    intercept 1 randstra 1 0 base 21.75  ady 84 base*ady 1827 trtp 0 0 1 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Severe base: 75 week 1'    intercept 1 randstra 1 0 base 21.75  ady 7 base*ady 152.25  trtp 0 0 0 1 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0/ ilink cl ;
estimate 'Severe base: 75 week 4'    intercept 1 randstra 1 0 base 21.75  ady 28 base*ady 609  trtp 0 0 0 1 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0/ ilink cl ;
estimate 'Severe base: 75 week 8'    intercept 1 randstra 1 0 base 21.75  ady 56 base*ady 1218  trtp 0 0 0 1 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0/ ilink cl ;
estimate 'Severe base: 75 week 12'    intercept 1 randstra 1 0 base 21.75  ady 84 base*ady 1827  trtp 0 0 0 1 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0/ ilink cl ;

estimate 'Severe base: 150 week 1'    intercept 1 randstra 1 0 base 21.75  ady 7 base*ady 152.25  trtp 1 0 0 0 0  avisitn 1 0 0 0  trtp*avisitn 1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 150 week 4'    intercept 1 randstra 1 0 base 21.75  ady 28 base*ady 609  trtp 1 0 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 150 week 8'    intercept 1 randstra 1 0 base 21.75  ady 56 base*ady 1218  trtp 1 0 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 150 week 12'    intercept 1 randstra 1 0 base 21.75  ady 84 base*ady 1827  trtp 1 0 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Severe base: 300 week 1'    intercept 1 randstra 1 0 base 21.75  ady 7 base*ady 152.25  trtp 0 1 0 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 300 week 4'    intercept 1 randstra 1 0 base 21.75  ady 28 base*ady 609  trtp 0 1 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 300 week 8'    intercept 1 randstra 1 0 base 21.75  ady 56 base*ady 1218  trtp 0 1 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Severe base: 300 week 12'    intercept 1 randstra 1 0 base 21.75  ady 84 base*ady 1827  trtp 0 1 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;


estimate 'Moderate base: PBO week 1'    intercept 1 randstra 1 0 base 15  ady 7 base*ady 105  trtp 0 0 0 0 1  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0/ilink cl ;
estimate 'Moderate base: PBO week 4'    intercept 1 randstra 1 0 base 15  ady 28 base*ady 420  trtp 0 0 0 0 1  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0/ilink cl ;
estimate 'Moderate base: PBO week 8'    intercept 1 randstra 1 0 base 15  ady 56 base*ady 840  trtp 0 0 0 0 1  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0/ilink cl ;
estimate 'Moderate base: PBO week 12'    intercept 1 randstra 1 0 base 15  ady 84 base*ady 1260  trtp 0 0 0 0 1  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1/ilink cl ;

estimate 'Moderate base: 37.5 week 1'    intercept 1 randstra 1 0 base 15  ady 7 base*ady 105  trtp 0 0 1 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 37.5 week 4'    intercept 1 randstra 1 0 base 15  ady 28 base*ady 420  trtp 0 0 1 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 37.5 week 8'    intercept 1 randstra 1 0 base 15  ady 56 base*ady 840  trtp 0 0 1 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 37.5 week 12'    intercept 1 randstra 1 0 base 15  ady 84 base*ady 1260  trtp 0 0 1 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Moderate base: 75 week 1'    intercept 1 randstra 1 0 base 15  ady 7 base*ady 105  trtp 0 0 0 1 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0/ ilink cl ;
estimate 'Moderate base: 75 week 4'    intercept 1 randstra 1 0 base 15  ady 28 base*ady 420  trtp 0 0 0 1 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0/ ilink cl ;
estimate 'Moderate base: 75 week 8'    intercept 1 randstra 1 0 base 15  ady 56 base*ady 840  trtp 0 0 0 1 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0/ ilink cl ;
estimate 'Moderate base: 75 week 12'    intercept 1 randstra 1 0 base 15  ady 84 base*ady 1260  trtp 0 0 0 1 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0/ilink cl ;

estimate 'Moderate base: 150 week 1'    intercept 1 randstra 1 0 base 15  ady 7 base*ady 105  trtp 1 0 0 0 0  avisitn 1 0 0 0  trtp*avisitn 1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 150 week 4'    intercept 1 randstra 1 0 base 15  ady 28 base*ady 420  trtp 1 0 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 150 week 8'    intercept 1 randstra 1 0 base 15  ady 56 base*ady 840  trtp 1 0 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 150 week 12'    intercept 1 randstra 1 0 base 15  ady 84 base*ady 1260  trtp 1 0 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Moderate base: 300 week 1'    intercept 1 randstra 1 0 base 15  ady 7 base*ady 105  trtp 0 1 0 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 300 week 4'    intercept 1 randstra 1 0 base 15  ady 28 base*ady 420  trtp 0 1 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 300 week 8'    intercept 1 randstra 1 0 base 15  ady 56 base*ady 840  trtp 0 1 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Moderate base: 300 week 12'    intercept 1 randstra 1 0 base 15  ady 84 base*ady 1260  trtp 0 1 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;


estimate 'Low base: PBO week 1'    intercept 1 randstra 1 0 base 6.25  ady 7 base*ady 43.75  trtp 0 0 0 0 1  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0/ilink cl ;
estimate 'Low base: PBO week 4'    intercept 1 randstra 1 0 base 6.25  ady 28 base*ady 175  trtp 0 0 0 0 1  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0/ilink cl ;
estimate 'Low base: PBO week 8'    intercept 1 randstra 1 0 base 6.25  ady 56 base*ady 350  trtp 0 0 0 0 1  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0/ilink cl ;
estimate 'Low base: PBO week 12'    intercept 1 randstra 1 0 base 6.25  ady 84 base*ady 525  trtp 0 0 0 0 1  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1/ilink cl ;

estimate 'Low base: 37.5 week 1'    intercept 1 randstra 1 0 base 6.25  ady 7 base*ady 43.76  trtp 0 0 1 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 37.5 week 4'    intercept 1 randstra 1 0 base 6.25  ady 28 base*ady 175  trtp 0 0 1 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 37.5 week 8'    intercept 1 randstra 1 0 base 6.25  ady 56 base*ady 350  trtp 0 0 1 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 37.5 week 12'    intercept 1 randstra 1 0 base 6.25  ady 84 base*ady 525  trtp 0 0 1 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Low base: 75 week 1'    intercept 1 randstra 1 0 base 6.25  ady 7 base*ady 43.76  trtp 0 0 0 1 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  1 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 75 week 4'    intercept 1 randstra 1 0 base 6.25  ady 28 base*ady 175  trtp 0 0 0 1 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 1 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 75 week 8'    intercept 1 randstra 1 0 base 6.25  ady 56 base*ady 350  trtp 0 0 0 1 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 1 0  0 0 0 0/ilink cl ;
estimate 'Low base: 75 week 12'    intercept 1 randstra 1 0 base 6.25  ady 84 base*ady 525  trtp 0 0 0 1 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 1  0 0 0 0/ilink cl ;

estimate 'Low base: 150 week 1'    intercept 1 randstra 1 0 base 6.25  ady 7 base*ady 43.76  trtp 1 0 0 0 0  avisitn 1 0 0 0  trtp*avisitn 1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 150 week 4'    intercept 1 randstra 1 0 base 6.25  ady 28 base*ady 175  trtp 1 0 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 150 week 8'    intercept 1 randstra 1 0 base 6.25  ady 56 base*ady 350  trtp 1 0 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 150 week 12'    intercept 1 randstra 1 0 base 6.25  ady 84 base*ady 525  trtp 1 0 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;

estimate 'Low base: 300 week 1'    intercept 1 randstra 1 0 base 6.25  ady 7 base*ady 43.76  trtp 0 1 0 0 0  avisitn 1 0 0 0  trtp*avisitn 0 0 0 0  1 0 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 300 week 4'    intercept 1 randstra 1 0 base 6.25  ady 28 base*ady 175  trtp 0 1 0 0 0  avisitn 0 1 0 0  trtp*avisitn 0 0 0 0  0 1 0 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 300 week 8'    intercept 1 randstra 1 0 base 6.25  ady 56 base*ady 350  trtp 0 1 0 0 0  avisitn 0 0 1 0  trtp*avisitn 0 0 0 0  0 0 1 0  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;
estimate 'Low base: 300 week 12'    intercept 1 randstra 1 0 base 6.25  ady 84 base*ady 525  trtp 0 1 0 0 0  avisitn 0 0 0 1  trtp*avisitn 0 0 0 0  0 0 0 1  0 0 0 0  0 0 0 0  0 0 0 0/ilink cl ;

run;

/* create just one row to generate a difference of treatment group dose 37.5 - placebo*/
data est_cov_dat;
length base_type $ 15. ady 4 dose 4 intercept 4 ady 4 

trt_150 4 trt_300 4 trt_375 4 trt_75 4 trt_pbo 4 avisitn_3 4 avisitn_6 4 avisitn_10 4 avisitn_14 4
trt_150_avisitn_3 4 trt_150_avisitn_6 4 trt_150_avisitn_10 4 trt_150_avisitn_14 4
trt_300_avisitn_3 4 trt_300_avisitn_6 4 trt_300_avisitn_10 4 trt_300_avisitn_14 4
trt_375_avisitn_3 4 trt_375_avisitn_6 4 trt_375_avisitn_10 4 trt_375_avisitn_14 4
trt_75_avisitn_3 4 trt_75_avisitn_6 4 trt_75_avisitn_10 4 trt_75_avisitn_14 4
trt_pbo_avisitn_3 4 trt_pbo_avisitn_6 4 trt_pbo_avisitn_10 4 trt_pbo_avisitn_14 4

trt2_150 4 trt_300 4 trt2_375 4 trt_75 4 trt2_pbo 4 
trt2_150_avisitn_3 4 trt2_150_avisitn_6 4 trt2_150_avisitn_10 4 trt2_150_avisitn_14 4
trt2_300_avisitn_3 4 trt2_300_avisitn_6 4 trt2_300_avisitn_10 4 trt2_300_avisitn_14 4
trt2_375_avisitn_3 4 trt2_375_avisitn_6 4 trt2_375_avisitn_10 4 trt2_375_avisitn_14 4
trt2_75_avisitn_3 4 trt2_75_avisitn_6 4 trt2_75_avisitn_10 4 trt2_75_avisitn_14 4
trt2_pbo_avisitn_3 4 trt2_pbo_avisitn_6 4 trt2_pbo_avisitn_10 4 trt2_pbo_avisitn_14 4

base 4  ady_base 4  

randstra_comp 4 randstra_noncomp 4   ;

base_type = "Severe";
ady = 7;
dose = 37.5;

intercept = 1;
ady = 7;

trt_150 = 0;
trt_300 = 0;
trt_375 = 1;
trt_75 = 0;
trt_pbo = 0;

trt2_150 = 0;
trt2_300 = 0;
trt2_375 = 0;
trt2_75 = 0;
trt2_pbo = -1;


avisitn_3 = 1;
avisitn_6 = 0;
avisitn_10 = 0;
avisitn_14 = 0;

trt_150_avisitn_3 = 0;
trt_150_avisitn_6 = 0;
trt_150_avisitn_10 = 0;
trt_150_avisitn_14 = 0;

trt_300_avisitn_3 = 0;
trt_300_avisitn_6 = 0;
trt_300_avisitn_10 = 0;
trt_300_avisitn_14 = 0;

trt_375_avisitn_3 = 1;
trt_375_avisitn_6 = 0;
trt_375_avisitn_10 = 0;
trt_375_avisitn_14 = 0;

trt_75_avisitn_3 = 0;
trt_75_avisitn_6 = 0;
trt_75_avisitn_10 = 0;
trt_75_avisitn_14 = 0;

trt_pbo_avisitn_3 = 0;
trt_pbo_avisitn_6 = 0;
trt_pbo_avisitn_10 = 0;
trt_pbo_avisitn_14 = 0;

trt2_150_avisitn_3 = 0;
trt2_150_avisitn_6 = 0;
trt2_150_avisitn_10 = 0;
trt2_150_avisitn_14 = 0;

trt2_300_avisitn_3 = 0;
trt2_300_avisitn_6 = 0;
trt2_300_avisitn_10 = 0;
trt2_300_avisitn_14 = 0;

trt2_375_avisitn_3 = 0;
trt2_375_avisitn_6 = 0;
trt2_375_avisitn_10 = 0;
trt2_375_avisitn_14 = 0;

trt2_75_avisitn_3 = 0;
trt2_75_avisitn_6 = 0;
trt2_75_avisitn_10 = 0;
trt2_75_avisitn_14 = 0;

trt2_pbo_avisitn_3 = -1;
trt2_pbo_avisitn_6 = 0;
trt2_pbo_avisitn_10 = 0;
trt2_pbo_avisitn_14 = 0;

base = 21.75;
ady_base = 152.25;

randstra_comp=1;
randstra_noncomp=0;


output;
run;

/*run the NLestimate macro*/

%NLEstimate(inest=glmm_p_FE_est,incovb=glmm_p_FE_cov_num, label=Rate Difference,
f=exp(b_p1*intercept + b_p2*ady + b_p8*avisitn_3 + b_p9*avisitn_6 + b_p10*avisitn_10 + b_p11*avisitn_14 + b_p32*base + b_p33*ady_base + b_p34*randstra_comp + b_p35*randstra_noncomp) *
(exp( b_p3*trt_150 + b_p4*trt_300 + b_p5*trt_375 + b_p6*trt_75 + b_p7*trt_pbo + 
      b_p12*trt_150_avisitn_3 + b_p13*trt_150_avisitn_6 + b_p14*trt_150_avisitn_10 + b_p15*trt_150_avisitn_14 +
	  b_p16*trt_300_avisitn_3 + b_p17*trt_300_avisitn_6 + b_p18*trt_300_avisitn_10 + b_p19*trt_300_avisitn_14 +
	  b_p20*trt_375_avisitn_3 + b_p21*trt_375_avisitn_6 + b_p22*trt_375_avisitn_10 + b_p23*trt_375_avisitn_14 +
	  b_p24*trt_75_avisitn_3 + b_p25*trt_75_avisitn_6 + b_p26*trt_75_avisitn_10 + b_p27*trt_75_avisitn_14 +
	  b_p28*trt_pbo_avisitn_3 + b_p29*trt_pbo_avisitn_6 + b_p30*trt_pbo_avisitn_10 + b_p31*trt_pbo_avisitn_14
	  )-
exp(b_p3*trt2_150 + b_p4*trt2_300 + b_p5*trt2_375 + b_p6*trt2_75 + b_p7*trt2_pbo +
      b_p12*trt2_150_avisitn_3 + b_p13*trt2_150_avisitn_6 + b_p14*trt2_150_avisitn_10 + b_p15*trt2_150_avisitn_14 +
	  b_p16*trt2_300_avisitn_3 + b_p17*trt2_300_avisitn_6 + b_p18*trt2_300_avisitn_10 + b_p19*trt2_300_avisitn_14 +
	  b_p20*trt2_375_avisitn_3 + b_p21*trt2_375_avisitn_6 + b_p22*trt2_375_avisitn_10 + b_p23*trt2_375_avisitn_14 +
	  b_p24*trt2_75_avisitn_3 + b_p25*trt2_75_avisitn_6 + b_p26*trt2_75_avisitn_10 + b_p27*trt2_75_avisitn_14 +
	  b_p28*trt2_pbo_avisitn_3 + b_p29*trt2_pbo_avisitn_6 + b_p30*trt2_pbo_avisitn_10 + b_p31*trt2_pbo_avisitn_14)), 
score=est_cov_dat,title=Difference means,outscore=score_out);


proc print data=score_out;
var base_type ady dose pred StdErrPred lower upper ChiSq pr;
run;

Output for dose 37.5 - placebo for severe baseline subjects at day 7;

Group estimates from glimmix - we can see the correct difference in means was  calculated