Date: Mon, 23 Feb 2004 15:25:10 -0800
Reply-To: Dale McLerran <stringplayer_2@YAHOO.COM>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: Dale McLerran <stringplayer_2@YAHOO.COM>
Subject: Re: ZI_NOR NLMIXED_Problem with 2 Codes
In-Reply-To: <ED92B7C78985D611BB3B00B0D0FC32BA19B50C@SH0CX691>
Content-Type: text/plain; charset=us-ascii
Hamani,
First, I would note that the PDF function does not specify the
variance of the response, but the standard deviation of the
response. In order to have equivalent models, you would want
to write
prob = pdf('NORMAl', &depvar, mu, sqrt(s2e))
That said, the only effect that should have is that your estimate
of the variance would differ between the two models. What happens
in your second invocation if you change
prob = pdf('NORMAl', &depvar, mu, s2e)
loglike = log(prob);
to
resid = y - mu;
loglike = -0.5*(log(2*constant('PI')) + log(s2e) + (resid**2)/s2e);
Perhaps there is some execution error which is arising in the
use of the PDF function which is avoided if you write the density
directly. Let me also point out that you might gain some insight
into a more robust specification of your likelihood by looking at
the log file for the following code:
%macro normal;
%let data=projline2;
%let depvar=habitud;
proc nlmixed data=&data ;
parms b0=0.6 b1=1 s2e=14 s2u1=37 ;
mu =b0+ b1*time ;
mu =mu+u1 ;
model &depvar ~ normal(mu,s2e);
random u1 ~ normal(0,s2u1 ) subject=id ;
run;
%mend;
options mprint;
%normal
When you examine the log file after executing the macro, you
will observe that the procedure NLMIXED includes some exception
handling for the log likelihood. You might want to code your
own likelihood function specification with some exception
handling.
Dale
--- "Elmaache, Hamani" <Hamani.Elmaache@CCRA-ADRC.GC.CA> wrote:
> Hi there.
>
> I have hard time to understand Code 1 works very well, but NOT Code
> 2.
> For Code 2, I get the following error message:
>
> ERROR: Quadrature accuracy of 0.000100 could not be achieved with 31
> points.
> The achieved
> accuracy was 1.000000.
>
> In fact, the CODE 2 is first step to set Zero-infated Normal
>
> Thanks for any help.
>
>
> /********** Code 1 ********/
> %let data=projline2;
> %let depvar=habitud;
> proc nlmixed data=&data ;
> parms b0=0.6 b1=1 s2e=14 s2u1=37 ;
> mu =b0+ b1*time ;
> mu =mu+u1 ;
> model &depvar ~ normal(mu,s2e);
> random u1 ~ normal(0,s2u1 ) subject=id ;
> run;
>
> /********** Code 2 ********/
>
> proc nlmixed data=&data ;
> parms b0=0.6 b1=1 s2e=14 s2u1=37 ;
> mu =b0+ b1*time ;
> mu =mu+u1 ;
> prob = pdf('NORMAL',&depvar,mu,s2e);
> loglike = log(prob);
> model &depvar ~ general(loglike);
> random u1 ~ normal(0,s2u1 ) subject=id ;
> run;
=====
---------------------------------------
Dale McLerran
Fred Hutchinson Cancer Research Center
mailto: dmclerra@fhcrc.org
Ph: (206) 667-2926
Fax: (206) 667-5977
---------------------------------------
__________________________________
Do you Yahoo!?
Yahoo! Mail SpamGuard - Read only the mail you want.
http://antispam.yahoo.com/tools