--- Tim Wade <wadetj@GMAIL.COM> wrote: > Hello: I have the following Proc NLMIXED code which works perfectly: > > proc nlmixed data=body1; > parms a1=-2.80 beta1=0.2539 logs2u=-0.1; > eta=a1+beta1*m_l10count_dy+u; > expeta=exp(eta); > p=expeta/(1+expeta); > s2u = exp(2*logs2u); > model hcgi~binary(p); > random u~normal(0, s2u) subject=beachnum; > run; > > I would like to try the same model with an identity link. I tried to > use the following, with initial parameter estimates from Proc Genmod > (dist=family, link=identity) > > > proc nlmixed data=body1; > parms a1=0.0244 beta1=0.0290 beta2=0.1202 beta3=0.0607 s2u=-0.1; > eta=a1+beta1*m_l10count_dy+beta2*gicontact_any+beta3*venfest+u; > *expeta=exp(eta); > p=eta; > s2u = exp(2*logs2u); > model hcgi~binary(p); > random u~normal(0, s2u) subject=beachnum; > run; > > but it won't converge. Am I on the right track here? Any suggestions? > Thanks much. Tim >

Tim,

Yes, you are on the right track (although I can't say that I really like using the identity link with a binary distribution). Below, I show code for a fixed effect model where there is a binary response and we assume an identity link. Parameter estimates are first obtained from PROC GENMOD and then from PROC NLMIXED. Point estimates and standard errors are the same for both procedures.

/* generate some data where the response is binary and the */ /* probability of the response depends on a single var X */ data test; seed=1234579; do i=1 to 100; x = (ranuni(seed)-0.5); /* X centered around zero */ eta = x; p = exp(eta) / (1 + exp(eta)); y = ranbin(seed,1,p); output; end; keep x y; run;

proc genmod data=test descending; model y = x / dist=bin link=identity; run;

proc nlmixed data=test; parms b0=0.5 b1=0; eta = b0 + b1*x; p = eta; model y ~ binary(p); run;

When I first attempted to fit the NLMIXED version, I did not initialize parameters with a PARMS statement. Note that the effect of excluding the PARMS statement is that both b0 and b1 were initialized to be 1. That means that the probability of the response could be as high as 1+max(x)=1.5. For all X>0, the probability has initial value >1. Employing the default parameters, the model did not iterate. NLMIXED returned the error message in the log that "No valid parameter points were found."

In general, the binary response model with identity link appears to fail if the probability estimates are frequently <=0 or >=1. If the population probability of the response is very low or very high (I might guess very low from your intercept term), then the added random intercept may push many probability estimates beyond the [0,1] limits.

I would also note that there is an inconsistency in your code such that you express the VARIABLE s2u=exp(logs2u) in a programming statement, but you initialize the PARAMETER s2u to -0.1 on the PARMS statement. How NLMIXED deals with this inconsistency, I do not know. However, I would modify the PARMS statement to give logs2u=-10 or some other very small value. That will initialize the model such that the random effects will, initially anyway, be quite small. That is, given a large negative value for logs2u, you should have initial model which is essentially the fixed effect binary response with identity link. You can allow the random effects to increase from there.

HTH,

Dale

--------------------------------------- Dale McLerran Fred Hutchinson Cancer Research Center mailto: dmclerra@NO_SPAMfhcrc.org Ph: (206) 667-2926 Fax: (206) 667-5977 ---------------------------------------

__________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com