Proc nlmixed: about bernoulli distribution

[Jie111]

proc nlmixed data=lab5;
parms mu 0.9 sigma 1 gamma0 0.5 gamma1 -0.5 ;  
bounds 0 < sigma; 
sqrt2pi = sqrt(2*constant('pi'));
eta = exp((gamma0 + gamma1*y)) / (1 + exp((gamma0 + gamma1 * y)));

f1 = -log(sigma) - log(sqrt2pi*y) - (log(y)-mu)**2  / (2*sigma**2);
f2 = (eta**(d=1) )* ((1 - eta)**(d=2)); /*d is binary variable, 1 or 2*/

ll= log (f1 * f2);
model y~general(ll);
run;

i want to calculate the parameters for joint distribution.

f(y,d) = f(y)f(d|y)

f(y) follows a normal distribution, f(d|y) follow Bernoulli distribution. Then I wrote the syntax above.

however, there was no results and NOTE: Execution error for observation 1.

is there anyone know the problem and how to deal with it?

[SteveDenham]

I don't see anything that jumps out at me about the syntax, so it might be that your data has something unusual in it.  Since the error shows at the first observation, could you post that first record so we could see if something goes to zero, such that the log does not exist.

SteveDenham

[Jie111]

Thanks for the reply. I attached the dataset as follows:

[SteveDenham]

You have some observations with d=0.  I think those would result in f2=1, rather than either eta or 1-eta.  I don't know if this is what causes the error message, but it seems to me that it results in an incorrect value for f(d|y), which is a problem in itself.  What happens if you exclude the observations where d=0?  That may not be interesting, but it would let you know if that is a problem.

SteveDenham

[Jie111]

I did the analyses with deleting those with d=0. however, the Note Execution error for observation 1 is still there.

data nn;
set lab5;
if d=0 then delete;
run;

proc nlmixed data=nn;
parms mu 0.9 sigma 1 gamma0 0.5 gamma1 -0.5 ;  
bounds 0 < sigma; 
sqrt2pi = sqrt(2*constant('pi'));
eta = exp((gamma0 + gamma1*y)) / (1 + exp((gamma0 + gamma1 * y)));

f1 = -log(sigma) - log(sqrt2pi*y) - (log(y)-mu)**2  / (2*sigma**2);
f2 = (eta**(d=1) )* ((1 - eta)**(d=2)); /*d is binary variable, 1 or 2*/

ll= log (f1 * f2);
model y~general(ll);
run;

The new dataset:

[SteveDenham]

This is going to sound amazingly obtuse, but just for sanity's sake, could you define d1 as d1 = d-1 and run this

proc nlmixed data=nn;
parms mu 0.9 sigma 1 gamma0 0.5 gamma1 -0.5 ;  
bounds 0 < sigma;
d1=d - 1: 
sqrt2pi = sqrt(2*constant('pi'));
eta = exp((gamma0 + gamma1*y)) / (1 + exp((gamma0 + gamma1 * y)));

f1 = -log(sigma) - log(sqrt2pi*y) - (log(y)-mu)**2  / (2*sigma**2);
f2 = (eta**d1 )* ((1 - eta)**d1); /*d is binary variable, 1 or 2*/

ll= log (f1 * f2);
model y~general(ll);
run;

Plugging in your starting values, I don't get anything odd to cause this to fail, so maybe the indicator variable is the source.

This is mystifying me.

SteveDenham

[SteveDenham]

Pay no attention to the previous post. Instead I tried to see what the assignment (d=1) or (d=2) was doing.  I got this after simplifying terms.

proc nlmixed data=nn;
parms mu 0.9 sigma 1 gamma0 0.5 gamma1 -0.5 ;  
bounds 0 < sigma;
d1 = d -1
sqrt2pi = sqrt(2*constant('pi'));
eta = exp((gamma0 + gamma1*y)) / (1 + exp((gamma0 + gamma1 * y)));

f1 = -log(sigma) - log(sqrt2pi*y) - (log(y)-mu)**2  / (2*sigma**2);
if d1=0 then
f2 = eta: 
else if d1=1 then
f2= ( 1 - eta);


ll= log (f1 * f2);
model y~general(ll);
run

Since eta is a logit, I'm not real sure that the product is the log likelihood of the conditional probability.  I am sure I am missing something here.

SteveDenham

[Jie111]

Really thanks for the reply. The syntax for f1 was wrong. I revised it and now it workes well.

[StatDave]

You might want to take a look at the joint modeling capabilities available in PROC GLIMMIX to see if this does what you want. See the example titled "Joint Modeling of Binary and Count Data" in the GLIMMIX documentation which illustrates this capability for data with a binary response and a count response.

[Jie111]

Thanks for the reply. I would try it.

[Rick_SAS · Accepted solution]

Check your loglikelihood. You are using 

ll= log (f1 * f2);

but f1 is negative.

[Jie111]

Yeah, you are right. I just noticed that the syntax for the log-normal distribution is wrong.

Thanks all of you a lot.

attached with the right syntax.

data nn;
  set lab5;
  if d ne 0 then output;
run;

proc nlmixed data=lab5;
parms mu 0.9 sigma 1 gamma0 0.5 gamma1 -0.5 ;  
bounds 0 < sigma; 
sqrt2pi = sqrt(2*constant('pi'));

eta = exp((gamma0 + gamma1*y)) / (1 + exp((gamma0 + gamma1 * y)));

f1 = exp (- (log (y)-mu)**2 / (2*sigma**2) )/ sigma / sqrt2pi / y;
f2 = (eta**(d=1) )* ((1 - eta)**(d=2));

ll= log (f1 * f2);
model y~general(ll);
where d ne 0;
run;