Date:         Sun, 22 Sep 2002 22:11:43 -0300
Reply-To:     Silvano Cesar da Costa <silvano@UEL.BR>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         Silvano Cesar da Costa <silvano@UEL.BR>
Subject:      Re: NLMIXED vs GENMOD
Comments: To: Dale McLerran <stringplayer_2@yahoo.com>
Content-Type: text/plain; charset="iso-8859-1"

Dear Professor Dale,

Your explications were wonderfull. Thanks a lot.

The question/answer 3 was:

> 3) If I have longitudinal data I use Genmod with Repeated option to
> fit it. How can I to do this with Nlmixed????

> A type=CS analysis can be obtained by adding a random intercept term
> to eta in the nlmixed code.

But, suppose that I want a Type=AR(1) analysis. Can I to do this with
NLMIXED????

Thanks,

Silvano.











-----Original Message-----
From: Dale McLerran
To: Silvano; SAS-L@LISTSERV.UGA.EDU
Cc: Dale McLerran
Sent: 17/09/02 20:57
Subject: Re: NLMIXED vs GENMOD

--- Silvano <silvano@uel.br> wrote:
> Dear Professor Dale and other SAS experts,
>
> I fitted a one-way design with two variables ( 3 and 4 levels) and I
> used the Genmod and Nlmixed.
> I have some questions about Genmod and Nlmixed:
>
> 1) In the Genmod I have   Log Likelihood = 3895.1283   and in the
> Nlmixed I have    -Log Likelihood = 593.048752.
> Why have I this diference???

There are some instances where Genmod does not maximize the full
likelihood, but only operates on the sufficient statistic.  Then,
Genmod produces a model likelihood which is offset from the true
likelihood by a constant amount.  Consider the binomial distribution
for a moment.  The probability model is

    p(y,mu) = c(n,y) * mu^y * (1-mu)^(n-y)

We can maximize the likelihood or we can maximize

    f(p(y,mu)) = mu^y * (1-mu)^(n-y)

and obtain the same estimate for mu.  If we use the latter function,
we are not maximizing the full likelihood function, and our model
likelihood is not the true likelihood.  Here is an example:

data test;
  n=10;
  seed=1234579;
  do i=1 to 50;
    y = ranbin(seed,n,.2);
    output;
  end;
run;

proc nlmixed data=test;
  eta = b0;
  p = exp(eta)/(1 + exp(eta));
  prob = comb(n,y) * (p**y) * (1-p)**(n-y);
  loglike = log(prob);
  model y ~ general(loglike);
run;

proc nlmixed data=test;
  eta = b0;
  p = exp(eta)/(1 + exp(eta));
  prob = (p**y) * (1-p)**(n-y);
  loglike = log(prob);
  model y ~ general(loglike);
run;

proc genmod data=test;
  model y/n = / dist=bin;
run;

Note that the log likelihood produced by Genmod is the value given
by the second invocation of Nlmixed where the full likelihood is
not operated on.  Only the sufficient statistic is employed.  But
since C(n,y) is a constant regardless of the fitted model, the
reported log likelihood is only off by a constant amount.  Therefore,
the reported log likelihood can be employed to construct likelihood
ratio tests.


>
> 2) I may to fit overdispersion in the Genmod using  SCALE = Pearson
> (p.ex.). Can I to do this in the Nlmixed??? I tried to multiplicate
> the Log Likelihood for     parms  phi=0.5,      but isn't work.

You would have to get an estimate of the overdispersion parameter
separately, and then use the overdispersion parameter in the
Nlmixed code.  Again, for the binomial model, let's generate the
overdispersion parameter employing Genmod and then show that we
can obtain the same result if we use Nlmixed.


proc genmod data=test;
  model y/n = / dist=bin pscale;
run;

proc nlmixed data=test;
  eta = b0;
  p = exp(eta)/(1+exp(eta));
  prob = comb(10,y) * (p**y) * (1-p)**(10-y);
  prob = prob**(0.8624**-2);
  loglike = log(prob);
  model y ~ general(loglike);
run;

Note that the Pearson estimate of the scale parameter was 0.8624
out of Genmod.  When we exponentiate the likelihood by 0.8624**-2,
then we obtain the same result as far as parameter estimates and
standard errors from both Genmod and Nlmixed.  We get different
likelihoods because Genmod is working only on the sufficient
statistic and Nlmixed (in the code I show above) is operating on
the full likelihood.

>
> 3) If I have longitudinal data I use Genmod with  Repeated option to
> fit it. How can I to do this with Nlmixed????

A type=CS analysis can be obtained by adding a random intercept term
to eta in the nlmixed code.  Again, let's look at a binomial response.

data test2;
  n=10;
  seed=1234579;
  beta0 = -.3;
  do cluster=1 to 50;
    gamma0 = .4*rannor(seed);
    eta = beta0 + gamma0;
    p = exp(eta)/(1+exp(eta));
    do i=1 to 5;
      y = ranbin(seed,n,p);
      output;
    end;
  end;
  keep n y cluster;
run;

proc genmod data=test2;
  class cluster;
  model y/n =;
  repeated subject=cluster / type=cs;
run;

proc nlmixed data=test2;
  eta = beta0 + gamma0;
  p = exp(eta)/(1+exp(eta));
  prob = comb(n,y) * (p**y) * (1-p)**(n-y);
  loglike = log(prob);
  model y ~ general(loglike);
  random gamma0 ~ normal([0], [vgamma0]) subject=cluster;
run;

>
> In true I have longitudinal data with random effects (latent
> variable).
>
>
> Thanks,
>
> Silvano.
>


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

__________________________________________________
Do you Yahoo!?
Yahoo! News - Today's headlines
http://news.yahoo.com