#
# Analysis of Pearson crab ratio data
# Mixture distribution, 2 normal components
# MC Neale 8 March 2013
# Load libraries & data
require(OpenMx)
load("pearson.Rda")
summary(pearson)
head(pearson)

# Knock out the one with infinite value
pearson <- pearson[-29,]

# Set number of variables & classes
nvar <- 1
nClass <-2 
selVars<-c("ratio")

g1Model <-  mxModel("group1",
                mxMatrix("Symm", nvar, nvar, values=.1, name="expCov", lbound=1.E-4, free=T),
                mxMatrix("Full", 1, nvar, free=T, name="expMean", values=.4),
                mxData(pearson, type="raw"), 
                mxFIMLObjective("expCov", "expMean", dimnames="ratio", vector=T)
                )
g2Model <-  mxModel(g1Model, name="group2", 
                mxMatrix("Full", 1, nvar, values=.6, free=T, name="expMean")
                )
# can repeat above “duplication” step for more classes, in a loop if needed, and stick them in a list[]

mixtureModel <- mxModel("mixture", g1Model, g2Model,
                    mxMatrix(type="Full", nrow=nClass, ncol=1, values=runif(nClass), 
						free=c(rep(T,nClass-1),F), lbound=1.E-1, name="praw"),
					mxAlgebra(praw %x% (1/sum(praw)), name="p"),
                    mxMatrix(type="Full", nrow=dim(pearson)[1], ncol=1, 
						values=as.matrix(pearson$freq), name="freq"),
                    mxAlgebra( -2*sum (freq * (log(cbind(group1.objective, group2.objective) %*% p ))), 
						name="min2LL"), 
                    mxAlgebraObjective("min2LL")
                )
# Run twice to check solution - could randomize starting values also
summary(mixtureModelFit <- mxRun(mixtureModel, unsafe=T))
summary(mixtureModelFit <- mxRun(mixtureModelFit))



# Copy Model
nonMixtureModel <- mixtureModel

# Fix proportion parameter to 1 and fix mean & variance parameters of group 2
nonMixtureModel$praw@free[]<-F
nonMixtureModel$praw@values[]<-c(1,0)

nonMixtureModel$group2.expCov@free[]<-F
nonMixtureModel$group2.expMean@free[]<-F

summary(nonMixtureModelFit <- mxRun(nonMixtureModel))

mxCompare(mixtureModelFit,nonMixtureModelFit)