Date:         Wed, 2 Mar 2005 16:03:15 -0800
Reply-To:     Dale McLerran <stringplayer_2@YAHOO.COM>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
Comments:     DomainKeys? See http://antispam.yahoo.com/domainkeys
From:         Dale McLerran <stringplayer_2@YAHOO.COM>
Subject:      Re: Calling R from SAS
In-Reply-To:  <0ICQ00NQAVKS2F@mailgw01.uaa.alaska.edu>
Content-Type: text/plain; charset=us-ascii

--- David Neal <afdbn@UAA.ALASKA.EDU> wrote:

> Cool.  Why don't you post your LCA code here in SAS-L, or at least
> get it into the proceedings of SUGI or a regional SUG?  SAS doesn't
> have a proc which does LCA, as far as I know, and that would be a
> nice
> little tidbit.  Alternatively, you could do like Mike Friendly and
> put
> it up on your university's website...
>
> David
> --
> David Cassell, CSC
> Cassell.David@epa.gov
> Senior computing specialist
> mathematical statistician
>
>
>
> David,
>
> Although I'm all for sharing, I'm still not happy with my LCA code.
> I
> "borrowed" large chunks of it from the R package (e1071) and it just
> doesn't
> do what I want it to do.  I may have to go out and get something like
> Latent
> Gold.  So, if I do keep going on it, you will probably see it at
> SUGI.  If I
> ever get over the fact that I'm always a bit self-conscious about my
> code
> I'll probably post some of it.
>
> David Neal
>

David,

I don't know the exact nature of your latent class problem.  I
have employed the procedure NLMIXED to estimate a type of latent
class model.  Perhaps the approach which I employ can be
modified to your particular needs.

The problem which I was addressing is one in which we have cases
(cancer patients) and controls.  Suppose that classes were not
latent and we knew the status for the K types of cases (so that
we have K+1 conditions, including controls).  Then we could employ
a generalized logits model to estimate the probability of each
case type.  For each case type, we would have

    P(Y=C{i}) = exp(eta{i}) / (1 + exp(eta1) + ... + exp(etaK))

But we do not know the classes.  We do know that the sum of the
probabilities for each of the various types of case is equal to
the total probability that an observation is a case.

    P(Y=C) = P(Y=C1) + P(Y=C2) + ... + P(Y=CK)
           = (exp(eta1) + exp(eta2) + ... + exp(etak)) /
             (1 + exp(eta1) + exp(eta2) + ... + exp(etak))

This is a latent class model that we can estimate.  The code
below lays out the approach.

  proc nlmixed;
    eta1 = b1_0 + b1_1*X1 + b1_2*X2 + ... + B1_J*XJ;
    eta2 = b2_0 + b2_1*X1 + b2_2*X2 + ... + B2_J*XJ;
    ...
    etaK = bK_0 + bK_1*X1 + bK_2*X2 + ... + BK_J*XJ;

    array eta {K};
    array exp_eta {K};
    sum_exp_eta = 0;
    do i=1 to K;
      exp_eta{i} = exp(eta{i});
      sum_exp_eta + exp_eta{i};
    end;

    pY_eq_Case = sum_exp_eta / (1 + sum_exp_eta);
    pY_ne_Case = 1 - pY_eq_Case;
    if Case then loglike = log(pY_eq_Case);
    else         loglike = log(pY_ne_Case);

    model Case ~ general(loglike);
  run;


For estimability purposes, you need to initialize the parameters
of the model.  Fixing some parameters (to zero, most often) is
another way to achieve estimability.

NLMIXED can estimate a very broad array of maximum likelihood
type problems.  Latent Class Analysis is often (usually) solved
using maximum likelihood estimation methods.  You might see if
you can write the likelihood model for your LCA employing
NLMIXED.  It might be easier to work with than IML (or not!).

HTH,

Dale


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

__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around
http://mail.yahoo.com