Code for Generalised nonlinear mixed effect model in Proc nlmixed for repeated measure count data

[sohelsayma]

Dear All,

proc nlmixed data= WORK.IMPORT;
parms beta0=1 beta1=0.0001 S2u=0 ;
eta = beta0 + beta1*year1+ u;
expeta = exp(eta);
p =1-(1+expeta)**(-interval);
mean_dead=den*p;
ll=lgamma(1/alpha + dead) - lgamma(1 + dead) - lgamma(1/alpha)
                + dead*log(alpha*mean_dead)
                - (1/alpha + dead)*log(1 + alpha*mean_dead);  * log likelihood function for negative binomial regression ;
        model dead ~ general(ll);
random u ~ normal(0,s2u)subject=plot;
predict p out=pp;
run;

For my research purpose, I need to do Generalised nonlinear mixed effect model to analyse some repeated measure count data. I have some tree mortality data from 1950 which were measure in different time interval. The same tree was measure in different census period. I have already finalized the input data structure. My data structures are Plot ID, census year, Census interval, live tree, dead tree, Soil class (Please check the attachment). I would like to model the mortality trend overtime for different soil type and aspect. The similar analysis was done in an article entitled “Tree mortality in response to climate change induced drought across Beijing, China. Climatic Change (2014) 124: 179–190”. In that article they used SAS proc nlmixed. Also another paper from Nature climate change also used SAS proc nlmixed for similar data. Nonlinear because demographic rates are compounded by the time interval, generalized because we were using count data.  I have written the following code. I am not sure whether I wrote the code correctly. This code is to answer the Question: Are mortality rates within plots changing over the observation years? Can you please check whether I wrote the code correctly? Also can you suggest me what could be the code when I will consider soil type and aspect? Any help is highly appreciated. Thanks

[lvm]

There are many parameterizations of the negative binomial, and I don't have time to work through the algebra to see if yours is correct. You should note that NLMIXED has a built-in negative binomial:

model y ~ negbin(.,.);

But you have to get the parameterization right. You should read the section in the User's Guide for NLMIXED on "Built-in Log-LIkelihood Functions" to see what I mean by the different parameterizations. I do question some of your code because the NB has two parameters, but you are only trying to estimate one of them. That is, you only have parameters for the intercept and slope, and the random plot variance. But there should be a scale parameter for the NB. This is apparently "interval" in your code, but you are using it as a constant, not a parameter to be estimated.

Your model appears to be a true GLMM. Thus, you could fit it with GLIMMIX very simply. See the text section I mentioned above for converting one parameterization to another.

proc glimmix;

class plot;

model dead = year1 / dist=negbin s link=log;

random plot;

run;

With GLIMMIX you can easily generalize things to handle more complex repeated-measures structure, and to add other variables to the model.

[sohelsayma]

Thanks for your time and suggestions. I will go through this.