NLMIXED and random effects

[Kui]

I am using proc nlmixed procedure for the zero inflated beta binomial model. It runs through without problem. Then, I tried to add random intercept and slope, because it is repeated measurement data, 22 quarters and around 4000 hospitals, in which some hospitals may only have 1 or 2 quarters data, most of hospitals have more than 20 quarters data.

In the full model, I have 20 covariates, but after error message showed optimization cannot be completed, I only keep the time variable and random effects. I tried to get the starting values by the time variable only, and tried all the convergence algorithm techniques. But, still the the same error message or Hessian eigenvalue is negative.

Now I am thinking to drop the hospitals only reported once, which may confuse the SAS to find the random effects?

Any comments or suggestions would be appriciated.

Thanks,
Kui

[Dale]

Kui,

You might want to show how you have parameterized the zero-inflated beta-binomial distribution. Have you parameterized it in terms of two shape parameters or in terms of the mean and variance? I presume that you are fitting a beta-binomial regression model, in which case I would think that you would want to parameterize the beta-binomial density function in terms of mean and variance and then compute the shape parameters as functions of the mean, variance, and number of trials on which the response (number of successes) is based.

But without seeing how you have parameterized the zero-inflated beta-binomial distribution, let's consider how many different types of overdispersion you are attempting to model here. The beta-binomial distribution represents an overdipsersed binomial distribution. Allowing for zero-inflation introduces a second type of overdispersion function. Adding random intercept and slope is yet another form of overdispersion. So, you are trying to fit a model that has three sources of overdispersion.

It comes as no surprise that you are having trouble introducing random effects on top of an zero-inflated beta-binomial. Conceivably, it is possible introduce random subject effects to a zero-inflated beta-binomial distribution. But in practice, you may find it difficult to estimate such a model.

[Kui]

Dale,

Thanks so much for your prompt reply.

Using the same data, I have run through several regressions, which included Poisson, negative binomial, zero-inflated Poisson, zero-inflated negative binomial, binomial, beta binomial, and zero-inflated beta binomial.  After that, I have the zero-inflated negative binomial with random effects, but like I said in the first thread, I cannot add the random effects with zero-inflated beta binomial and not even beta binomial.

I have two beta binomial parameterizations. (Without random, they work well).

BB1:

proc nlmixed data=dsn tech=newrap qpoints=50 maxiter=5000 maxfunc=20000;

parms b0=-2.612 b1=0.037 b2=0.669 b3=0.214 b24=-0.056 b25=0.012 b26=0.045 b27=-0.008

rho=0.035 varRin=0.5099 VarRslp=0.003 cov=-0.0164;

eta=(b0+b1*group_0+b2*group_1+b3*group_2+(b24+Rslp)*nqt+Rin+b25*groupnqt_0+b26*groupnqt_1+b27*groupnqt_2);

expeta=exp(eta);

p=expeta/(1+expeta);

N=totdenom; r=R_totnum;

A=p*(1-rho)/rho;

B=(1-p)*(1-rho)/rho;

loglike=(lgamma(n+1) - lgamma(r+1) - lgamma(n-r+1))+lgamma(A+r)+

lgamma(n+B-r)+lgamma(A+B) - lgamma(A+B+n) - lgamma(A) -lgamma(B);

model r~general(loglike);

random Rin Rslp ~ nomal([0,0], [varRin, Cov, VarRslp])  subject=hsp_ID;

run;

BB2:

proc nlmixed data=ami1_04Q1 tech=newrap maxiter=5000 maxfunc=20000;

parms b0=-2.544 b1=0.183 b2=0.66 b3=0.097 b24=-0.065 b25=0.011 b26=0.052 b27=-0.0007

A=4.27 VarRin=0.5099 VarRslp=0.002163 cov =-0.01644;

eta=(b0+b1*group_0+b2*group_1+b3*group_2+(b24+Rslp)*nqt+Rin+b25*groupnqt_0+b26*groupnqt_1+b27*groupnqt_2);

expeta=exp(eta);

p=expeta/(1+expeta);

N=totdenom; r=R_totnum;

B=A*(1-p)/p);

loglike=(lgamma(n+1)-lgamma(r+1) -lgamma(n-r+1))+lgamma(A+r)+lgamma(n+B-r)+lgamma(A+B)

- lgamma(A+B+n) - lgamma(A) -lgamma(B);

model r~general(loglike);

random Rin Rslp~ normal([0,0], [varRin, Cov, varRslp]) subject=hsp_ID;

run;

I appreciate your comment and suggestions!

Kui