Introduction to Mixture Modeling
Download
Report
Transcript 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
i1
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 <- model@output$matrices$Class1.ThresholdsClass1
class2T <- model@output$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 <- model@output$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