How do I calculate LS Means in Proc NLMixed?

[elsie5355]

Hello - 

I am estimating a zero truncated Poisson model in PROC NLMIXED I would like to get the LS Means of the interaction of group*time, but am unsure how to do this or if it is possible. If the model was not zero truncated, I could use PROC GLIMMIX. For example:

proc glimmix data=mydata;

class id time group var1;

model y=time group var1 group*time group*var1 group*time*var1/link=log dist=poisson;

random _residual_/subject=id;

lsmeans group*time/ilink pdiff;

lsmeans group*time*var1/ilink pdiff;

run;

Is it possible to get similar output in NLMIXED using the Estimate statements? I've also seen the NLMeans macro, but I had a hard time deciphering how to output the data needed for the macro. 

proc nlmixed data=mydata;
log_lambda = intercept +b_var1*var1 +  b_time1*time1+ b_time2*time2+b_group*group+ b_grouptime1*grouptime1+ b_grouptime2*grouptime2+ b_grouptime1var1*grouptime1var1+ b_grouptime2var1*grouptime2var1+u;
lambda = exp(log_lambda);
ll = y*log_lambda - lambda - log(1-exp(-lambda)) - lgamma(y+1);
model y ~ general(ll);
random u ~ normal(0,s2u) subject=id;
run;

Thanks.

[StatDave]

Use the E option in your LSMEANS statements in GLIMMIX to display the linear combination of model parameters that it estimates. Then use the ESTIMATE statement in NLMIXED to estimate the same linear combination in that procedure. But make sure that you parameterize the categorical variables the same way as in GLIMMIX. 

[elsie5355]

Do I use the coefficients in the Row1, Row2, etc. from the table? The estimates are much lower than I was expecting. 

For example, would the mean for group =1 be

Estimate 'Group=1' exp(intercept*1+b_Var1*.5+b_group*1+b_Time1*.5+b_Time2*.5+b_Time1Group*.5+b_Time2Group*.5+b_Var1Group*.5)

Coefficients for Group Least Squares Means
Effect	Time1	Time2	Group	Time1Group	Time2Group	Var1	Var1Group	Row1	Row2
Intercept								1	1
Var1						0		0.5	0.5
Var1						1		0.5	0.5
Group			0					1
Group			1						1
Time1	0							0.5	0.5
Time1	1							0.5	0.5
Time2		0						0.5	0.5
Time2		1						0.5	0.5
Time1Group				0				0.5	0.5
Time1Group				1				0.5	0.5
Time2Group					0			0.5	0.5
Time2Group					1			0.5	0.5
Var1Group							0	0.5	0.5
Var1Group							1	0.5	0.5

[StatDave]

If you just add the E option in the LSMEANS statements in the GLIMMIX code you showed, there would not be TIME1 and TIME2. There should just be a (presumably) two-level TIME variable. Also, you seem to have changed the interactions into binary variables. The only change to the GLIMMIX code you sent should be to add the E options, not any change to the MODEL. However, if your VAR1, GROUP, and TIME variables are all 0,1-coded binary predictors in the CLASS and MODEL statements, then all of the non-zero parameters from GLIMMIX will be associated with the 0 levels of those variables. It will be easiest if you add the REF=FIRST option following a slash in the CLASS statement. This will cause the nonzero parameters to be associated with the 1 levels. Then, your approach to writing the ESTIMATE statement from the E option tables seems correct. Note that for the GROUP=0 LS-mean you will simply leave b_group out, or equivalently, multiply it by zero.