topic Re: Proc nlmixed: about bernoulli distribution in Statistical Procedures

Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 13:02:13 GMT

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?

Re: Proc nlmixed: about bernoulli distribution

SteveDenham — Fri, 29 May 2020 13:58:38 GMT

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

Re: Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 14:03:25 GMT

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

Re: Proc nlmixed: about bernoulli distribution

SteveDenham — Fri, 29 May 2020 14:18:31 GMT

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

Re: Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 14:22:23 GMT

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:

Re: Proc nlmixed: about bernoulli distribution

StatDave — Fri, 29 May 2020 15:25:25 GMT

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.

Re: Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 15:41:42 GMT

Thanks for the reply. I would try it.

Re: Proc nlmixed: about bernoulli distribution

SteveDenham — Fri, 29 May 2020 17:26:59 GMT

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

Re: Proc nlmixed: about bernoulli distribution

SteveDenham — Fri, 29 May 2020 17:56:48 GMT

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

Re: Proc nlmixed: about bernoulli distribution

Rick_SAS — Fri, 29 May 2020 18:03:12 GMT

Check your loglikelihood. You are using

ll= log (f1 * f2);

but f1 is negative.

Re: Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 18:16:30 GMT

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;

Re: Proc nlmixed: about bernoulli distribution

Jie111 — Fri, 29 May 2020 18:17:30 GMT

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