### Introduction to Mixture Modeling

```CROSS-SECTIONAL MIXTURE
MODELING
1
Shaunna L. Clark
Advanced Genetic Epidemiology Statistical
Workshop
October 23, 2012
OUTLINE
What is a mixture?
 Introduction to LCA (LPA)

Basic Analysis Ideas\Plan and Issues
 How to choose the number of classes

How do we implement mixtures in OpenMx?
 Factor Mixture Model
 What do classes mean for twin modeling?

2
HOMOGENEITY VS. HETEROGENEITY

Most models assume homogeneity
i.e. Individuals in a sample all follow the same model
 What have seen so far (for the most part)


But not always the case

Ex: Sex, Age, Patterns of Substance Abuse
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
3
0
Alcohol
Tobacco
Cannabis
Opiates
Heroin
WHAT IS MIXTURE MODELING
Used to model unobserved heterogeneity by
identifying different subgroups of individuals
Ex: IQ, Religiosity
4
LATENT CLASS ANALYSIS (LCA)
5
Also known as Latent Profile Analysis (LPA) if
you have continuously distributed variables
LATENT CLASS ANALYSIS
Introduced by Lazarsfeld & Henry, Goodman, Clogg, Dayton
& Mcready

Setting


Cross-sectional data
Multiple items measuring a construct


Hypothesized construct represented as latent class variable
(categorical latent variable)


12 items measuring the construct of Cannabis
Abuse/Dependence
Different categories of Cannabis Abuse\Dependence patterns
Aim
Identify items that indicate classes well
Estimate proportion of sample in each class (class
probability)
 Classify individuals into classes (posterior probabilities)


6
LATENT CLASS ANALYSIS CONT’D
1
0.9
C
0.8
0.7
0.6
Non Users
0.5
0.4
Legal Drug
Users
0.3
Multiple Illicit
Users
0.2
0.1
0
Alcohol Tobacco Cannabis Opiates
Heroin
x1
x2
x3
x4
x5
7
LATENT CLASS ANALYSIS MODEL
Dichotomous (0/1) indicators u: u1, u2, ... , ur
Categorical latent variable c: c = k ; k = 1, 2, ... , K
Marginal probability for item uj = 1,
(probability item uj =1 is the sum over all class of
the product of the probability of being in class k
and the probability of endorsing item uj given
that you are in class k)
8
JOINT PROBABILITIES


Joint probability of all u’s, assuming conditional
independence:
Probability of observing a given response pattern
is equal to the sum over all classes of the product
of being in a given class and the probability of
observing a response on item 1 given that you are
in latent class k, . . . (repeat for each item)
9
POSTERIOR PROBABILITIES
Probability of being in class k given your
response pattern
 Used to assign most likely class membership


Based on highest posterior probability
Individual
P(Class1)
P(Class2)
MLCM
A
.90
.10
1
B
.8
.2
2
10
MODEL TESTING

Log-likelihood ratio χ2 test (LLRT)


Overall test against the data with H1 being the
unrestricted multinomial
Problem: Not distributed as χ2 due to boundary
conditions


Don’t use it!!! (McLachlan & Peele, 2000)
Information Criteria
Akaike Information Criteria, AIC (Akaike,1974)
AIC = 2h-2ln(L)
 Bayesian Information Criteria, BIC (Schwartz, 1978)
BIC = -2ln(L)+h*ln(n)
 Where L = log-likelihood, h = number of parameters,
n = sample size
 Chose model with lowest value of IC

11
OTHER TESTS

Since can’t do LLRT, use test which approximate
the difference in LL values between k and k-1
class models.
Vuong-Lo-Mendell-Rubin, LMR-LRT (Lo, Mendell, &
Rubin, 2001)
 Parametric bootstrapped LRT, BLRT (McLachlan,
1987)


P-value is probability that H0 is true
H0: k-1 classes; H1: k classes
 A low p-value indicates a preference for the estimated
model (i.e. k classes)


Look for the first time the p-value is nonsignificant or greater than 0.05
12
ANALYSIS PLAN
1.
Fit model with 1-class


2.
3.
Everyone in same class
Sometimes simple is better
Fit LCA models 2-K classes
Chose best number of classes
Seems simple right???
13
NOT REALLY . . .LOTS OF KNOWN ISSUES
IN MIXTURE ANALYSIS

Global vs. Local Maximum
Log Likelihood
Log Likelihood
Local
Global
Global
Local
Parameter

Parameter
Use multiple sets of random starting values to make
sure have global solution. Make sure that best LL
value has replicated
14
DETERMINING THE NUMBER OF CLASSES:
CLASS ENUMERATION

No agreed upon way to determine the correct
number of latent classes
Statistical comparisons (i.e. ICs, LRTs)
 Interpretability and usefulness of classes

Substantive theory
 Relationship to auxiliary variables
 Predictive validity of classes
 Class size


Quality of Classifications (not my favorite)
Classification table based on posterior probabilities
 Entropy - A value close to 1 indicates good classification in
that many individuals have posterior probabilities close to 0
or 1

15
SUGGESTED STRATEGY
Nylund et al. (2007), Tofighi & Enders (2008),
among others
 Simulation studies comparing tests and
information criteria described previously
 Suggest:

Use BIC and LMR to narrow down the number of
plausible models
 Then run BLRT on those models because BLRT can
be computationally intensive

16
OPENMX:
LCA EXAMPLE SCRIPT
17
LCA_example.R
MIXTURES IN OPENMX

Specify class-specific models


Specify class probabilities


Create MxModel objects for each class
Create an MxMatrix of class
probabilities\proportions
Specify model-wide objective function
Pull everything together in a parent model with data
 Weighted sum of the class models


Estimate entire model
Note: One of potentially many ways to do this
18
CLASS SPECIFIC MODELS
nameList <- names(<dataset>)
class1 <- mxModel("Class1",
mxMatrix("Iden", name = "R", nrow = nvar, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "M", nrow = 1, ncol = nvar, free=FALSE),
mxMatrix("Full", name = "ThresholdsClass1", nrow = 1, ncol = nvar,
list("Threshold",nameList), free=TRUE),
dimnames =
mxFIMLObjective(covariance="R", means="M", dimnames=nameList,
thresholds="ThresholdsClass1",vector=TRUE))
Repeat for every class in your model
Don’t be like me, make sure to
change class numbers
19
DEFINE THE MODEL
lcamodel <- mxModel("lcamodel", class1, class2, mxData(vars,
type="raw"),
Next, specify class membership probabilities
20
CLASS MEMBERSHIP PROBABILITIES

When specifying need to remember:
1.
2.
Class probabilities must be positive
Must sum to a constant - 1
mxMatrix("Full", name = "ClassMembershipProbabilities",
nrow = nclass, ncol = 1, free=TRUE,
labels = c(paste("pclass", 1:nclass, sep=""))),
mxBounds(c(paste("pclass", 1:nclass, sep="")),0,1),
mxMatrix("Iden", nrow = 1, name = "constraintLHS"),
mxAlgebra(sum(ClassMembershipProbabilities),
name = "constraintRHS"),
21
mxConstraint(constraintLHS == constraintRHS),
MODEL-WIDE OBJECTIVE FUNCTION
Weighted sum of individual class likelihoods
 Weights are class probabilities


2LL  2*log  pk Lk


k
i1

So for two classes:
2LL  2 * log( p1L1  p2 L2 )
22
MODEL WIDE OBJECTIVE FUNCTION
CONT’D
mxAlgebra(
-2*sum(log(pclass1%x%Class1.objective
+ pclass2%x%Class2.objective)),
name="lca"),
mxAlgebraObjective("lca"))
)
Now we run the model:
model <- mxRun(lcamodel)
And we wait and wait and wait till it’s done.
23
PROFILE PLOT
One way to interpret the classes is to plot them.
 In our example we had binary items, so the
thresholds are what distinguishes between
classes



Can plot the thresholds
Or you can plot the probabilities
More intuitive
 Easier for non-statisticians to understand

24
PROFILE PLOTS IN R\OPENMX
#Pulling out thresholds
class1T <- [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */\$matrices\$Class1.ThresholdsClass1
class2T <- [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */\$matrices\$Class2.ThresholdsClass2
#Converting threshold to probabilities
class1P<-t(1/(1+exp(-class1T)))
class2P<-t(1/(1+exp(-class2T)))
25
PROFILE PLOTS CONT’D
plot(class1P, type="o", col="blue",ylim=c(0,1),axes=FALSE, ann=FALSE)
axis(1,at=1:12,lab=nameList)
axis(2,las=1,at=c(0,0.2,0.4,0.6,0.8,1))
box()
lines(class2P,type="o", pch=22, lty=2, col="red")
title(main="LCA 2 Class Profile Plot", col.main="black",font.main=4)
title(xlab="DSM Items", col.lab="black")
title(ylab="Probability", col.lab="black")
legend("bottomright",c("Class 1","Class 2"), cex=0.8,
col=c("blue","red"),pch=21:22,lty=1:2)
26
OPENMX EXERCISE
Unfortunately, it takes long time for these to run
so not feasible to do in this session
 However, I’ve run the 2-, 3-, and 4- class LCA
models for this data and (hopefully) the .Rdata
files are posted on the website


Exercise: Using the .Rdata files
1.
2.
Determine which model is better according to AIC\BIC

Want the lowest value
Make a profile plot of the best solution and interpret the
classes

What kind of substances users are there?
27
CODE TO PULL OUT LL AND COMPUTE
AIC\BIC
#Pull out LL
LL_2c <- [email protected]/* <![CDATA[ */!function(t,e,r,n,c,a,p){try{t=document.currentScript||function(){for(t=document.getElementsByTagName('script'),e=t.length;e--;)if(t[e].getAttribute('data-cfhash'))return t[e]}();if(t&&(c=t.previousSibling)){p=t.parentNode;if(a=c.getAttribute('data-cfemail')){for(e='',r='0x'+a.substr(0,2)|0,n=2;a.length-n;n+=2)e+='%'+('0'+('0x'+a.substr(n,2)^r).toString(16)).slice(-2);p.replaceChild(document.createTextNode(decodeURIComponent(e)),c)}p.removeChild(t)}}catch(u){}}()/* ]]> */\$Minus2LogLikelihood
LL_2cnsam = 1878
#parameters
npar <- (nclass-1) + (nthresh*nvar*nclass
npar
#Compute AIC & BIC
AIC_2c = 2*npar + LL_2c
AIC_2c
BIC_2c = LL_2c + (npar*log(nsam))
BIC_2c
28
TABLE OF RESULTS
# Classes
-2*LL
Npar
AIC
BIC
2
6589.96
25
6639
6778
3
6329.95
38
6405
6616
4
6308.96
51
6410
6693
29
3-CLASS PROFILE PLOT
30
FACTOR MIXTURE MODELING
31
PROBLEM WITH LCA
Once in a class, everyone “looks” the same.
 In the context of substance abuse, unlikely that
every user will have the same patterns of use

Withdrawal, tolerance, hazardous use
 There is variation within a latent class



Severity
One proposed solution is the factor mixture
model

Uses a latent class variables to classify individuals
and latent factor to model severity
32
σ2 F
FACTOR MIXTURE MODEL
F
C
λ1
x1
x2
x3
λ2
λ3
x4
λ4
λ5
x5
Classes can be indicated by item thresholds (categorical)\ item
means (continuous) or factor mean and variance
33
GENERAL FACTOR MIXTURE MODEL
yik = Λk ηik + εik ,
ηik = αk + ζik ,
where,
ζik ~ N(0, Ψk)


Similar to the FA model, except many parameters can
be class varying as indicated by the subscript k
Several variations of this model which differ in terms
of the measurement invariance

Lubke & Neale (2005), Clark et al. (2012)
34
FMM PROFILE PLOT
0.9
0.8
0.7
0.6
0.5
Non Users
0.4
Users
0.3
0.2
0.1
0
Alcohol
Tobacco
Cannabis
Opiates
Heroin
35
HOW DO WE DO THIS IN OPENMX?
You’ll have to wait till tomorrow!
 Factor Mixture Model is a generalization of the
Growth Mixture Model we’ll talk about tomorrow
afternoon.

36
MIXTURES & TWIN MODELS
37
How do we combine the ACDE model and
mixtures?
OPTION 1: ACE DIRECTLY ON THE
CLASSES
1.0 (MZ) / 0.5 (DZ)
aA
1.0
cB
cA
x1A
eA
eB
CA
CB
x2A
x3A x4A
x1B
x2B
aB
x3B x4B
38
WHAT WOULD THIS LOOK LIKE FOR THE
1.0 (MZ) / 0.5 (DZ)
FMM?
1.0
1.0 (MZ) / 0.5 (DZ)
aA
cA
eA
1.0
aA
x1A
aA
cA
cA
eA
CA
x3A x4A
cB
eB
eA
FB
FA
x2A
aB
x1B
x2B
CB
x3B x4B
39
FMM & ACE CONT’D


One of many possible ways to do FMM & ACE in the same
model
Can also have class specific ACE on the factors

Each class has own heritability
40
From Muthén et al. (2006)
ISSUE WITH OPTION 1

Model is utilizes the liability threshold model to
“covert” the latent categorical variable, C, to a
latent normal variable

This requires that classes are ordered

Ex: high, medium, low users

Don’t always have nicely ordered classes

Models are VERY time intensive

Take a vacation for a week or two
41
OPTION 2: THREE-STEP METHOD
1.
2.
3.
Estimate mixture model
Assign individuals into their most likely latent
class based on the posterior probabilities of class
membership
Use the observed, categorical variable of
assigned class membership as the phenotype in
a liability threshold model version of ACE
analysis
Note: Requires ordered classes
42
OPTION 2A

Contingency table analysis using most likely
class membership
Concordance between twins in terms of most likely
class membership
 If your classes are not ordered


Odds Ratio
 Excess
twin concordance due to stronger genetic relationship
can be represented by the OR for MZ twins compared to the OR
for DZ twins.

Place restrictions on the contingency table to test
specific hypotheses
Mendelian segregation, only shared environmental
effects
 Eaves (1993)

43
ISSUES WITH OPTION 2
Potential for biased parameter estimates and
underestimated standard errors
 Assigned membership ignores fractional class
membership suggested by posterior probabilities
 Treat the classification as not having any sampling
error
 Good option when entropy is high\ well separated
classes

Individual
P(Class1)
P(Class2)
MLCM
A
.90
.10
1
B
.8
.2
2
C
.51
.49
1
44
SELECTION OF CROSS-SECTIONAL
MIXTURE GENETIC ANALYSIS WRITINGS

Latent Class Analysis


Factor Mixture Analysis


Neale & Gillespie, 2005 (?); Clark, 2010; Clark et al.
(in preparation)
Additional References


Eaves, 1993; Muthén et al., 2006; Clark, 2010
McLachlan, Do, & Ambroise, 2004
Mixtures in Substance Abuse

Gillespie (2011, 2012)

Great cannabis examples
45
```