Date:         Wed, 7 May 2008 09:56:40 -0700
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: Scoring nlmixed
In-Reply-To:  <ff3a5425-fb9c-4254-a854-323f9c408b15@d1g2000hsg.googlegroups.com>
Content-Type: text/plain; charset=iso-8859-1

--- "Dr. Eli Y. Kling" <Eli.Kling@GMAIL.COM> wrote:

> On May 6, 8:48 pm, stringplaye...@YAHOO.COM (Dale McLerran) wrote:
> >
> >     proc nlmixed data=combined;
> >       eta = ...;
> >       mu = <f(eta)>;
> >       if score=1 then
> >          ll = 0;
> >       else
> >          ll = ...;
> >       model y ~ general(ll) subject=id;
> >       random ...
> >       predict mu out=predicted(where=(score=1));
> >     run;
> >
> >
>
> Hi - Sounds neet.
>
> Just to double check I understood correctly:
> My current code is:
> proc nlmixed data=Staging.SphereMartII;
>  parms b1=0.13 b2=0.08 b3=-0.24;
>  bounds b1-b2 >= 0;
>  eta = b1*ABP + b2*ABD + b3*ABN;
>  expeta = exp(eta);
>  p = expeta/(1+expeta);
>  model X ~ binomial(N,p);
>  predict eta out=eta;
>  predict X out=X;
>  predict p out=p;
> run;
>
> Following your suggestion I should change it to the folowing (using
> the indicator II)
>
> proc nlmixed data=Staging.SphereMartII;
>  parms b1=0.13 b2=0.08 b3=-0.24;
>  bounds b1-b2 >= 0;
>  eta = b1*ABP + b2*ABD + b3*ABN;
>  expeta = exp(eta);
>  p = expeta/(1+expeta);
>
>   if score=1 then
>          II = 0;
>       else
>          II =N;
>
>
>  model X ~ binomial(II,p);
>  predict eta out=eta;
>  predict X out=X;
>  predict p out=p;
> run;
>
>
> Did I understand you correctly? I was deliberating between putting a
> zero on the N parameter and 0.5 on the p baparmer of the Bin dist.
>
> Eli
>

No, that is not correct.  You need to code the log-likelihood of
the binomial response and specify to NLMIXED that you are using
your own likelihood function by using the GENERAL model specification.
For a binomial response with N observations and success probability
p, the likelihood computation is

   L(X|N,p) = (N! / ((N-X)! * X!)) * p^X * (1-p)^(N-X)

The log-likelihood value is

   ll(X|N,p) = log(N!) - log((N-X)!) - log(X!)
                       + X*log(p) + (N-X)*log(1-p)

             = lgamma(N+1) - lgamma(N-X+1) - lgamma(X+1)
                       + X*log(p) + (N-X)*log(1-p)

Thus, correct NLMIXED code would be:

  proc nlmixed data=Staging.SphereMartII;
    parms b1=0.13 b2=0.08 b3=-0.24;
    bounds b1-b2 >= 0;
    eta = b1*ABP + b2*ABD + b3*ABN;
    expeta = exp(eta);
    p = expeta/(1+expeta);

     if score=1 then
            LL = 0;
         else
            LL = lgamma(N+1) - lgamma(N-X+1) - lgamma(X+1)
                       + X*log(p) + (N-X)*log(1-p);

    model X ~ general(LL);
    predict eta out=eta;
    predict p out=p;
  run;


Note that I removed the "predict X out=X;" statement.  X is the
observed response, not a function of parameters of your model.
However, you might want to replace the statement that I removed
with a statement which predicts the expected value of X which is
N*p.  Thus, you might want to insert the statement

    predict N*p out=Xhat;

I would modify your code a bit further to eliminate the BOUNDS
statement.  Note that if the boundary condition is encountered
during estimation, proper derivatives cannot be computed.  This
can cause problems in the estimation process.  However, by modifying
your code just a little bit, you can achieve parameter estimates
which are guaranteed to be positive without the need to impose
boundary constraints.

All that you have to do in order to ensure that the parameters are
nonnegative is to write the model in terms of parameters which have
been exponentiated.  That is, let the parameters of your model be
log_b1 and log_b2 and let b1=exp(log_b1) and b2=exp(log_b2).  Thus,
you would have

  proc nlmixed data=Staging.SphereMartII;
    parms log_b1=%sysfunc(log(0.13))
          log_b2=%sysfunc(log(0.08))
          b3=-0.24;
    b1 = exp(log_b1);
    b2 = exp(log_b2);
    eta = b1*ABP + b2*ABD + b3*ABN;
    expeta = exp(eta);
    p = expeta/(1+expeta);

     if score=1 then
            LL = 0;
         else
            LL = lgamma(N+1) - lgamma(N-X+1) - lgamma(X+1)
                       + X*log(p) + (N-X)*log(1-p);

    model X ~ general(LL);
    predict eta out=eta;
    predict p out=p;
  run;


By the way, I take it that the reason for using NLMIXED here is
that you want to impose these conditions on b1 and b2.  Otherwise,
you could simply use the LOGISTIC or GENMOD procedure to fit your
model.  And, I would think that you must have fit your model
employing one of the two simpler procedures and found violation
of your expectation that the parameters b1>=0 and b2>=0.  Is that
correct?  How sure are you of these parameter constraints?  Also,
why do you not have an intercept term in your model.

You don't need to explain these things here on SAS-L (although I
would not discourage you from doing so!)  But you must be prepared
to defend your model at some point in time.  Who that defense is
made to depends on whether this model is being fit for publication
in an academic journal or whether the model is being used for some
decision purpose in a business environment or ...

Dale


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


      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile.  Try it now.  http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ