[S] Summary: max likelihood question

[VALLIAR1@westat.com (VALLIAR1)]

I received the following replies to an earlier query on MLEs.  Thanks to Albyn
Jones, Brad Biggerstaff, and Andrzej Galecki.

The question:
> 
> This is a question on how to compute MLEs in a special type of multivariate
> normal distribution.  Suppose the sample is clustered.  Units in different
> clusters are independent.  Within a cluster, units have a common variance,
> sigma.sq, and different units have a common correlation, rho.  Plus the
> within-cluster structure is the same for every cluster (same sigma.sq , same
> rho).
> 
> Is there an easy way to get MLEs of sigma.sq and rho in SPLUS, assuming
> multivariate normality?
> 


________________________________________________________________________
Albyn Jones:
Sounds like a standard mixture distribution.

The EM algorithm works pretty well for this problem.  Here is a simple
univariate version that someone posted a while back for two
component mixtures:

"EM" <- 
function(y, mu0, mu1, var0, var1, prob, eps = 1/100000)
{
        new.params <- c(mu0, mu1, var0, var1, prob)
        err <- 1
        while(err > eps) {
                # E step
                bayes <- (prob * dnorm(y, mu1, sqrt(var1)))/((prob * dnorm(y,
                        mu1, sqrt(var1))) + ((1 - prob) * dnorm(y, mu0, sqrt(
                        var0))))  
                # M step
                mu1 <- sum(bayes * y)/sum(bayes)
                mu0 <- sum((1 - bayes) * y)/sum(1 - bayes)
                var1 <- sum(bayes * (y - mu1)^2)/sum(bayes)
                var0 <- sum((1 - bayes) * (y - mu0)^2)/sum(1 - bayes)
                prob <- mean(bayes)
                old.params <- new.params
                new.params <- c(mu0, mu1, var0, var1, prob)
                err <- max(abs((old.params - new.params)/new.params))
         }
        list(mu=c(mu0,mu1),v=c(var0,var1),p=prob)
}

I have several other mixture functions, but none that do exactly what you
want.  The multivariate version is pretty straight forward I think:
compute weighted means and covariance matrices, given the "bayes" factor,
which you could think of as estimating the posterior probability for each
observation of being in the clusters, given the current parameter estimates.
I have been meaning to write such a function, so if you can't figure it
out from this, let me know.  See also the book by Titterington, Smith, 
and Makov on Finite Mixture distributions.


__________________________________________________________________________
Brad Biggerstaff:
Sounds like lme (linear mixed effects) would do you just fine.  You can
specify estimation to be maximum likelihood, where the default is restricted
maximum likelihood, I think.

I imagine varcomp might help you, too, but that's if you're interested only
in the variance components.  For fixed effects, too, use lme.



__________________________________________________________________________
Andrzej Galecki:
Your description implies compound symmetry model for covariance
matrix. To obtain sigma and rho estimates I think you can use lme
function contributed by Jose Pinheiro and Doug Bates.  



Richard Valliant            email: valliar1@westat.com
Westat              Phone: 301-610-5118
1650 Research Blvd, RE442       FAX: 301-294-2034
Rockville MD 20850
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu.  To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message:  unsubscribe s-news