R help - Oct 2005 - finite mixture model (2-component gaussian): plotting component gaussian components?

Dear Knowledgeable R Community Members,

Please excuse my ignorance, I apologize in advance if this is an easy question,
but I am a bit stumped and could use a little guidance.

I have a finite mixture modeling problem -- for example, a 2-component gaussian
mixture -- where the components have a large overlap, and
I am trying to use the "mclust" package to solve this problem.  

I need to decompose that mixture into its 2 components which will need to be
plotted.

What I don't know how to do is:

(1) restrict the number of components to 2 in the "EMclust" function

(2) obtain and plot a component gaussian density 

Regarding (1), I think this should be the 'G' value but the
documentation is somewhat cryptic as regards the format.
Regarding (2), I think I need to use the "cdens" function, but again
the documentation is somewhat cryptic.

Here is a little test script to illustrate.  (Note: my real dataset will not
have peaks this well separated, but I needed to find a small example.)

##################
# (1) get some data
##################
data(faithful)

##################
# (2) get model
##################
library(mclust)
MyMixtureModel<-summary(EMclust(faithful$eruptions),faithful$eruptions)

##################
# (3) plot mixture model
##################
attach(MyMixtureModel)
mclust1Dplot(data=faithful$eruptions,z=z,mu=mu,sigmasq=sigmasq,pro=pro,ask=FALSE,type=c("density"))
do.call("mclust1Dplot",c(list(data=faithful$eruptions,ask=FALSE,type=c("density")),MyMixtureModel))

##################
# (4) plot components
##################
???

Any information you might be able to shed on this would be very much
appreciated.

With appreciation for your help,

Karen
---
Karen M. Green, Ph.D.
Research Investigator
Drug Design Group
Sanofi Aventis Pharmaceuticals
Tucson, AZ
USA


	[[alternative HTML version deleted]]

Dear R Community,

I was finally able to figure out how to plot the components of
the finite mixture model.

Because it was not obvious, I am including the script here --
in case anyone in the future also needs to do this.

Karen
---
Karen M. Green, Ph.D.
Research Investigator
Drug Design Group
Sanofi Aventis Pharmaceuticals
Tucson, AZ, USA

#########################################################
# name:             testMclust1Dplot
#
# date:              23 October 2005
#
# status:            validated
#
# dependencies:  package mclust
#
# purpose:          to test the mclust1Dplot function &
#                      to test the plot of components
#########################################################
# (0) preliminaries
##############################################
library(mclust)

##############################################
# (1) get some data
##############################################
data(faithful)

attach(faithful)
eruptions<-sort(eruptions)
eruptions<-data.frame(eruptions)
eruptions<-as.numeric(eruptions$eruptions)
detach(faithful)

##############################################
# (2) get model
##############################################
MyMixtureModel<-summary(EMclust(eruptions),eruptions)

##############################################
# (3) plot mixture model
##############################################

attach(MyMixtureModel)
mclust1Dplot(data=eruptions,z=z,mu=mu,sigmasq=sigmasq,pro=pro,ask=FALSE,type=c("density"))
do.call("mclust1Dplot",c(list(data=eruptions,ask=FALSE,type=c("density")),MyMixtureModel))
detach(MyMixtureModel)

##############################################
# (4) plot components
##############################################
temp <- cdens(modelName =
MyMixtureModel$modelName,data=eruptions,mu=MyMixtureModel$mu,
sigmasq=MyMixtureModel$sigmasq,decomp=MyMixtureModel$decomp)

MyMixtureComponents<-cbind(eruptions,temp)
colnames(MyMixtureComponents)<-c("eruptions","component1","component2","component3","component4")

MyMixtureComponents<-data.frame(MyMixtureComponents)
MyMixtureComponents$eruptions<-as.numeric(MyMixtureComponents$eruptions)

attach(MyMixtureComponents)
lines(x=eruptions,y=component1*MyMixtureModel$pro[1],type="l",col="red")
lines(x=eruptions,y=component2*MyMixtureModel$pro[2],type="l",col="red")
lines(x=eruptions,y=component3*MyMixtureModel$pro[3],type="l",col="red")
lines(x=eruptions,y=component4*MyMixtureModel$pro[4],type="l",col="red")

detach(MyMixtureComponents)
#########################################################
>  -----Original Message-----
> From: 	Green, Karen M. PH/US  
> Sent:	Thursday, October 20, 2005 6:42 PM
> To:	'R-help@lists.R-project.org';
'r-help-request@stat.math.ethz.ch'
> Cc:	Green, Karen M. PH/US
> Subject:	finite mixture model (2-component gaussian):  plotting component
gaussian components?
> Importance:	High
> 
> Dear Knowledgeable R Community Members,
> 
> Please excuse my ignorance, I apologize in advance if this is an easy
question, but I am a bit stumped and could use a little guidance.
> 
> I have a finite mixture modeling problem -- for example, a 2-component
gaussian mixture -- where the components have a large overlap, and
> I am trying to use the "mclust" package to solve this problem.  
> 
> I need to decompose that mixture into its 2 components which will need to
be plotted.
> 
> What I don't know how to do is:
> 
> (1) restrict the number of components to 2 in the "EMclust"
function
> 
> (2) obtain and plot a component gaussian density 
> 
> Regarding (1), I think this should be the 'G' value but the
documentation is somewhat cryptic as regards the format.
> Regarding (2), I think I need to use the "cdens" function, but
again the documentation is somewhat cryptic.
> 
> Here is a little test script to illustrate.  (Note: my real dataset will
not have peaks this well separated, but I needed to find a small example.)>
> 
> ##################
> # (1) get some data
> ##################
> data(faithful)
> 
> ##################
> # (2) get model
> ##################
> library(mclust)
> MyMixtureModel<-summary(EMclust(faithful$eruptions),faithful$eruptions)
> 
> ##################
> # (3) plot mixture model
> ##################
> attach(MyMixtureModel)
>
mclust1Dplot(data=faithful$eruptions,z=z,mu=mu,sigmasq=sigmasq,pro=pro,ask=FALSE,type=c("density"))
>
do.call("mclust1Dplot",c(list(data=faithful$eruptions,ask=FALSE,type=c("density")),MyMixtureModel))
> 
> ##################
> # (4) plot components
> ##################
> ???
> 
> Any information you might be able to shed on this would be very much
appreciated.
> 
> With appreciation for your help,
> 
> Karen
> ---
> Karen M. Green, Ph.D.
> Research Investigator
> Drug Design Group
> Sanofi Aventis Pharmaceuticals
> Tucson, AZ
> USA
> 
	[[alternative HTML version deleted]]