gom.em: Discrete (Rasch) Grade of Membership Model
In sirt: Supplementary Item Response Theory Models
gom.em R Documentation
Discrete (Rasch) Grade of Membership Model
Description
This function estimates the grade of membership model (Erosheva, Fienberg
& Joutard, 2007; also called mixed membership model) by the EM algorithm
assuming a discrete membership score distribution. The function is restricted
to dichotomous item responses.
Usage
gom.em(dat, K=NULL, problevels=NULL, weights=NULL, model="GOM", theta0.k=seq(-5,5,len=15),
xsi0.k=exp(seq(-6, 3, len=15)), max.increment=0.3, numdiff.parm=1e-4,
maxdevchange=1e-6, globconv=1e-4, maxiter=1000, msteps=4, mstepconv=0.001,
theta_adjust=FALSE, lambda.inits=NULL, lambda.index=NULL, pi.k.inits=NULL,
newton_raphson=TRUE, optimizer="nlminb", progress=TRUE)
## S3 method for class 'gom'
summary(object, file=NULL, ...)
## S3 method for class 'gom'
anova(object,...)
## S3 method for class 'gom'
logLik(object,...)
## S3 method for class 'gom'
IRT.irfprob(object,...)
## S3 method for class 'gom'
IRT.likelihood(object,...)
## S3 method for class 'gom'
IRT.posterior(object,...)
## S3 method for class 'gom'
IRT.modelfit(object,...)
## S3 method for class 'IRT.modelfit.gom'
summary(object,...)
Arguments
dat
Data frame with dichotomous responses
K
Number of classes (only applies for model="GOM")
problevels
Vector containing probability levels for membership functions
(only applies for model="GOM"). If a specific space of probability
levels should be estimated, then a matrix can be supplied (see Example 1,
Model 2a).
weights
Optional vector of sampling weights
model
The type of grade of membership model. The default "GOM"
is the nonparametric grade of membership model. A parametric multivariate normal
representation can be requested by "GOMnormal".
The probabilities and membership
functions specifications described in Details are called via "GOMRasch".
theta0.k
Vector of \tilde{\theta}_k grid (applies only for model="GOMRasch")
xsi0.k
Vector of \xi_p grid (applies only for model="GOMRasch")
max.increment
Maximum increment
numdiff.parm
Numerical differentiation parameter
maxdevchange
Convergence criterion for change in relative deviance
globconv
Global convergence criterion for parameter change
maxiter
Maximum number of iterations
msteps
Number of iterations within a M step
mstepconv
Convergence criterion within a M step
theta_adjust
Logical indicating whether multivariate normal distribution
should be adaptively chosen during the EM algorithm.
lambda.inits
Initial values for item parameters
lambda.index
Optional integer matrix with integers indicating
equality constraints among \lambda item parameters
pi.k.inits
Initial values for distribution parameters
newton_raphson
Logical indicating whether Newton-Raphson should be
used for final iterations
optimizer
Type of optimizer. Can be "optim" or "nlminb".
progress
Display iteration progress? Default is TRUE.
object
Object of class gom
file
Optional file name for summary output
...
Further arguments to be passed
Details
The item response model of the grade of membership model
(Erosheva, Fienberg & Junker, 2002;
Erosheva, Fienberg & Joutard, 2007) with K classes
for dichotomous correct responses X_{pi}
of person p on item i is as follows (model="GOM")
P(X_{pi}=1 | g_{p1}, \ldots, g_{pK} )=\sum_k \lambda_{ik} g_{pk}
\quad, \quad \sum_{k=1}^K g_{pk}=1
\quad, \quad 0 \leq g_{pk} \leq 1
In most applications (e.g. Erosheva et al., 2007), the grade of
membership function \{g_{pk}\} is assumed to follow a Dirichlet
distribution. In our gom.em implementation
the membership function is assumed to be discretely represented
by a grid u=(u_1, \ldots, u_L) with entries between 0 and 1
(e.g. seq(0,1,length=5) with L=5).
The values g_{pk} of the membership function can then
only take values in \{ u_1, \ldots, u_L \} with the restriction
\sum_k g_{pk} \sum_l \bold{1}(g_{pk}=u_l )=1.
The grid u is specified by using the argument problevels.
The Rasch grade of membership model (model="GOMRasch") poses constraints
on probabilities \lambda_{ik} and membership functions g_{pk}.
The membership
function of person p is parameterized by a location parameter \theta_p
and a variability parameter \xi_p. Each class k is represented by
a location parameter \tilde{\theta}_k. The membership function is defined as
g_{pk} \propto
\exp \left[ - \frac{ (\theta_p - \tilde{\theta}_k)^2 }{2 \xi_p^2 } \right]
The person parameter \theta_p indicates the usual 'ability', while
\xi_p describes the individual tendency to change between classes
1,\ldots,K and their corresponding locations
\tilde{\theta}_1, \ldots,\tilde{\theta}_K.
The extremal class probabilities \lambda_{ik} follow the Rasch model
\lambda_{ik}=invlogit( \tilde{\theta}_k - b_i )=
\frac{ \exp( \tilde{\theta}_k - b_i ) }{ 1 + \exp( \tilde{\theta}_k - b_i ) }
Putting these assumptions together leads to the model equation
P(X_{pi}=1 | g_{p1}, \ldots, g_{pK} )=
P(X_{pi}=1 | \theta_p, \xi_p )=
\sum_k \frac{ \exp( \tilde{\theta}_k - b_i ) }{ 1 + \exp(\tilde{\theta}_k - b_i ) }
\cdot \exp \left[ - \frac{ (\theta_p - \tilde{\theta}_k)^2 }{2 \xi_p^2 } \right]
In the extreme case of a very small \xi_p=\varepsilon > 0 and
\theta_p=\theta_0, the Rasch model is obtained
P(X_{pi}=1 | \theta_p, \xi_p )=
P(X_{pi}=1 | \theta_0, \varepsilon )=
\frac{ \exp( \theta_0 - b_i ) }{ 1 + \exp( \theta_0 - b_i ) }
See Erosheva et al. (2002), Erosheva (2005, 2006) or Galyart (2015)
for a comparison of grade of membership models with latent trait models
and latent class models.
The grade of membership model is also published under the name
Bernoulli aspect model, see Bingham, Kaban and Fortelius (2009).
Value
A list with following entries:
deviance
Deviance
ic
Information criteria
item
Data frame with item parameters
person
Data frame with person parameters
EAP.rel
EAP reliability (only applies for model="GOMRasch")
MAP
Maximum aposteriori estimate of the membership function
EAP
EAP estimate for individual membership scores
classdesc
Descriptives for class membership
lambda
Estimated response probabilities \lambda_{ik}
se.lambda
Standard error for estimated response probabilities
\lambda_{ik}
mu
Mean of the distribution of (\theta_p, \xi_p)
(only applies for model="GOMRasch")
Sigma
Covariance matrix of (\theta_p, \xi_p)
(only applies for model="GOMRasch")
b
Estimated item difficulties (only applies for model="GOMRasch")
se.b
Standard error of estimated difficulties
(only applies for model="GOMRasch")
f.yi.qk
Individual likelihood
f.qk.yi
Individual posterior
probs
Array with response probabilities
n.ik
Expected counts
iter
Number of iterations
I
Number of items
K
Number of classes
TP
Number of discrete integration points for (g_{p1},...,g_{pK})
theta.k
Used grid of membership functions
...
Further values
References
Bingham, E., Kaban, A., & Fortelius, M. (2009).
The aspect Bernoulli model: multiple causes of presences and absences.
Pattern Analysis and Applications, 12(1), 55-78.
Erosheva, E. A. (2005). Comparing latent structures of the grade of membership,
Rasch, and latent class models. Psychometrika, 70, 619-628.
Erosheva, E. A. (2006). Latent class representation of the grade of membership
model. Seattle: University of Washington.
Erosheva, E. A., Fienberg, S. E., & Junker, B. W. (2002).
Alternative statistical models and representations for large sparse
multi-dimensional contingency tables.
Annales-Faculte Des Sciences Toulouse Mathematiques, 11,
485-505.
Erosheva, E. A., Fienberg, S. E., & Joutard, C. (2007).
Describing disability through individual-level mixture models
for multivariate binary data. Annals of Applied Statistics,
1, 502-537.
Galyardt, A. (2015).
Interpreting mixed membership models: Implications of Erosheva's representation
theorem. In E. M. Airoldi, D. Blei, E. A. Erosheva, & S. E. Fienberg (Eds.).
Handbook of Mixed Membership Models (pp. 39-65). Chapman & Hall.
See Also
For joint maximum likelihood estimation of the grade of membership model
see gom.jml.
See also the mixedMem package for estimating mixed membership
models by a variational EM algorithm.
The C code of Erosheva et al. (2007) can be downloaded from
http://projecteuclid.org/euclid.aoas/1196438029#supplemental.
Code from Manrique-Vallier can be downloaded from
http://pages.iu.edu/~dmanriqu/software.html.
See http://users.ics.aalto.fi/ella/publications/aspect_bernoulli.m
for a Matlab implementation of the algorithm in Bingham, Kaban and
Fortelius (2009).
Examples
#############################################################################
# EXAMPLE 1: PISA data mathematics
#############################################################################
data(data.pisaMath)
dat <- data.pisaMath$data
dat <- dat[, grep("M", colnames(dat)) ]
#***
# Model 1: Discrete GOM with 3 classes and 5 probability levels
problevels <- seq( 0, 1, len=5 )
mod1 <- sirt::gom.em( dat, K=3, problevels, model="GOM")
summary(mod1)
## Not run:
#-- some plots
#* multivariate scatterplot
car::scatterplotMatrix(mod1$EAP, regLine=FALSE, smooth=FALSE, pch=16, cex=.4)
#* ternary plot
vcd::ternaryplot(mod1$EAP, pch=16, col=1, cex=.3)
#***
# Model 1a: Multivariate normal distribution
problevels <- seq( 0, 1, len=5 )
mod1a <- sirt::gom.em( dat, K=3, theta0.k=seq(-15,15,len=21), model="GOMnormal" )
summary(mod1a)
#***
# Model 2: Discrete GOM with 4 classes and 5 probability levels
problevels <- seq( 0, 1, len=5 )
mod2 <- sirt::gom.em( dat, K=4, problevels, model="GOM" )
summary(mod2)
# model comparison
smod1 <- IRT.modelfit(mod1)
smod2 <- IRT.modelfit(mod2)
IRT.compareModels(smod1,smod2)
#***
# Model 2a: Estimate discrete GOM with 4 classes and restricted space of probability levels
# the 2nd, 4th and 6th class correspond to "intermediate stages"
problevels <- scan()
1 0 0 0
.5 .5 0 0
0 1 0 0
0 .5 .5 0
0 0 1 0
0 0 .5 .5
0 0 0 1
problevels <- matrix( problevels, ncol=4, byrow=TRUE)
mod2a <- sirt::gom.em( dat, K=4, problevels, model="GOM" )
# probability distribution for latent classes
cbind( mod2a$theta.k, mod2a$pi.k )
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0 0.0 0.0 0.0 0.17214630
## [2,] 0.5 0.5 0.0 0.0 0.04965676
## [3,] 0.0 1.0 0.0 0.0 0.09336660
## [4,] 0.0 0.5 0.5 0.0 0.06555719
## [5,] 0.0 0.0 1.0 0.0 0.27523678
## [6,] 0.0 0.0 0.5 0.5 0.08458620
## [7,] 0.0 0.0 0.0 1.0 0.25945016
## End(Not run)
#***
# Model 3: Rasch GOM
mod3 <- sirt::gom.em( dat, model="GOMRasch", maxiter=20 )
summary(mod3)
#***
# Model 4: 'Ordinary' Rasch model
mod4 <- sirt::rasch.mml2( dat )
summary(mod4)
## Not run:
#############################################################################
# EXAMPLE 2: Grade of membership model with 2 classes
#############################################################################
#********* DATASET 1 *************
# define an ordinary 2 latent class model
set.seed(8765)
I <- 10
prob.class1 <- stats::runif( I, 0, .35 )
prob.class2 <- stats::runif( I, .70, .95 )
probs <- cbind( prob.class1, prob.class2 )
# define classes
N <- 1000
latent.class <- c( rep( 1, 1/4*N ), rep( 2,3/4*N ) )
# simulate item responses
dat <- matrix( NA, nrow=N, ncol=I )
for (ii in 1:I){
dat[,ii] <- probs[ ii, latent.class ]
dat[,ii] <- 1 * ( stats::runif(N) < dat[,ii] )
}
colnames(dat) <- paste0( "I", 1:I)
# Model 1: estimate latent class model
mod1 <- sirt::gom.em(dat, K=2, problevels=c(0,1), model="GOM" )
summary(mod1)
# Model 2: estimate GOM
mod2 <- sirt::gom.em(dat, K=2, problevels=seq(0,1,0.5), model="GOM" )
summary(mod2)
# estimated distribution
cbind( mod2$theta.k, mod2$pi.k )
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0.243925644
## [2,] 0.5 0.5 0.006534278
## [3,] 0.0 1.0 0.749540078
#********* DATASET 2 *************
# define a 2-class model with graded membership
set.seed(8765)
I <- 10
prob.class1 <- stats::runif( I, 0, .35 )
prob.class2 <- stats::runif( I, .70, .95 )
prob.class3 <- .5*prob.class1+.5*prob.class2 # probabilities for 'fuzzy class'
probs <- cbind( prob.class1, prob.class2, prob.class3)
# define classes
N <- 1000
latent.class <- c( rep(1,round(1/3*N)),rep(2,round(1/2*N)),rep(3,round(1/6*N)))
# simulate item responses
dat <- matrix( NA, nrow=N, ncol=I )
for (ii in 1:I){
dat[,ii] <- probs[ ii, latent.class ]
dat[,ii] <- 1 * ( stats::runif(N) < dat[,ii] )
}
colnames(dat) <- paste0( "I", 1:I)
#** Model 1: estimate latent class model
mod1 <- sirt::gom.em(dat, K=2, problevels=c(0,1), model="GOM" )
summary(mod1)
#** Model 2: estimate GOM
mod2 <- sirt::gom.em(dat, K=2, problevels=seq(0,1,0.5), model="GOM" )
summary(mod2)
# inspect distribution
cbind( mod2$theta.k, mod2$pi.k )
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0.3335666
## [2,] 0.5 0.5 0.1810114
## [3,] 0.0 1.0 0.4854220
#***
# Model2m: estimate discrete GOM in mirt
# define latent classes
Theta <- scan( nlines=1)
1 0 .5 .5 0 1
Theta <- matrix( Theta, nrow=3, ncol=2,byrow=TRUE)
# define mirt model
I <- ncol(dat)
#*** create customized item response function for mirt model
name <- 'gom'
par <- c("a1"=-1, "a2"=1 )
est <- c(TRUE, TRUE)
P.gom <- function(par,Theta,ncat){
# GOM for two extremal classes
pext1 <- stats::plogis(par[1])
pext2 <- stats::plogis(par[2])
P1 <- Theta[,1]*pext1 + Theta[,2]*pext2
cbind(1-P1, P1)
}
# create item response function
icc_gom <- mirt::createItem(name, par=par, est=est, P=P.gom)
#** define prior for latent class analysis
lca_prior <- function(Theta,Etable){
# number of latent Theta classes
TP <- nrow(Theta)
# prior in initial iteration
if ( is.null(Etable) ){ prior <- rep( 1/TP, TP ) }
# process Etable (this is correct for datasets without missing data)
if ( ! is.null(Etable) ){
# sum over correct and incorrect expected responses
prior <- ( rowSums(Etable[, seq(1,2*I,2)]) + rowSums(Etable[,seq(2,2*I,2)]) )/I
}
prior <- prior / sum(prior)
return(prior)
}
#*** estimate discrete GOM in mirt package
mod2m <- mirt::mirt(dat, 1, rep( "icc_gom",I), customItems=list("icc_gom"=icc_gom),
technical=list( customTheta=Theta, customPriorFun=lca_prior) )
# correct number of estimated parameters
mod2m@nest <- as.integer(sum(mod.pars$est) + nrow(Theta)-1 )
# extract log-likelihood and compute AIC and BIC
mod2m@logLik
( AIC <- -2*mod2m@logLik+2*mod2m@nest )
( BIC <- -2*mod2m@logLik+log(mod2m@Data$N)*mod2m@nest )
# extract coefficients
( cmod2m <- sirt::mirt.wrapper.coef(mod2m) )
# compare estimated distributions
round( cbind( "sirt"=mod2$pi.k, "mirt"=mod2m@Prior[[1]] ), 5 )
## sirt mirt
## [1,] 0.33357 0.33627
## [2,] 0.18101 0.17789
## [3,] 0.48542 0.48584
# compare estimated item parameters
dfr <- data.frame( "sirt"=mod2$item[,4:5] )
dfr$mirt <- apply(cmod2m$coef[, c("a1", "a2") ], 2, stats::plogis )
round(dfr,4)
## sirt.lam.Cl1 sirt.lam.Cl2 mirt.a1 mirt.a2
## 1 0.1157 0.8935 0.1177 0.8934
## 2 0.0790 0.8360 0.0804 0.8360
## 3 0.0743 0.8165 0.0760 0.8164
## 4 0.0398 0.8093 0.0414 0.8094
## 5 0.1273 0.7244 0.1289 0.7243
## [...]
#############################################################################
# EXAMPLE 3: Lung cancer dataset; using sampling weights
#############################################################################
data(data.si08, package="sirt")
dat <- data.si08
#- Latent class model with 3 classes
problevels <- c(0,1)
mod1 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, problevels=problevels )
summary(mod1)
#- Grade of membership model with discrete distribution
problevels <- seq(0,1,length=5)
mod2 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, problevels=problevels )
summary(mod2)
#- Grade of membership model with multivariate normal distribution
mod3 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, theta0.k=10*seq(-1,1,len=11),
model="GOMnormal", optimizer="nlminb" )
summary(mod3)
## End(Not run)
sirt documentation built on May 29, 2024, 8:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| gom.em | R Documentation |
Discrete (Rasch) Grade of Membership Model
Description
This function estimates the grade of membership model (Erosheva, Fienberg & Joutard, 2007; also called mixed membership model) by the EM algorithm assuming a discrete membership score distribution. The function is restricted to dichotomous item responses.
Usage
gom.em(dat, K=NULL, problevels=NULL, weights=NULL, model="GOM", theta0.k=seq(-5,5,len=15),
xsi0.k=exp(seq(-6, 3, len=15)), max.increment=0.3, numdiff.parm=1e-4,
maxdevchange=1e-6, globconv=1e-4, maxiter=1000, msteps=4, mstepconv=0.001,
theta_adjust=FALSE, lambda.inits=NULL, lambda.index=NULL, pi.k.inits=NULL,
newton_raphson=TRUE, optimizer="nlminb", progress=TRUE)
## S3 method for class 'gom'
summary(object, file=NULL, ...)
## S3 method for class 'gom'
anova(object,...)
## S3 method for class 'gom'
logLik(object,...)
## S3 method for class 'gom'
IRT.irfprob(object,...)
## S3 method for class 'gom'
IRT.likelihood(object,...)
## S3 method for class 'gom'
IRT.posterior(object,...)
## S3 method for class 'gom'
IRT.modelfit(object,...)
## S3 method for class 'IRT.modelfit.gom'
summary(object,...)
Arguments
dat |
Data frame with dichotomous responses |
K |
Number of classes (only applies for |
problevels |
Vector containing probability levels for membership functions
(only applies for |
weights |
Optional vector of sampling weights |
model |
The type of grade of membership model. The default |
theta0.k |
Vector of |
xsi0.k |
Vector of |
max.increment |
Maximum increment |
numdiff.parm |
Numerical differentiation parameter |
maxdevchange |
Convergence criterion for change in relative deviance |
globconv |
Global convergence criterion for parameter change |
maxiter |
Maximum number of iterations |
msteps |
Number of iterations within a M step |
mstepconv |
Convergence criterion within a M step |
theta_adjust |
Logical indicating whether multivariate normal distribution should be adaptively chosen during the EM algorithm. |
lambda.inits |
Initial values for item parameters |
lambda.index |
Optional integer matrix with integers indicating
equality constraints among |
pi.k.inits |
Initial values for distribution parameters |
newton_raphson |
Logical indicating whether Newton-Raphson should be used for final iterations |
optimizer |
Type of optimizer. Can be |
progress |
Display iteration progress? Default is |
object |
Object of class |
file |
Optional file name for summary output |
... |
Further arguments to be passed |
Details
The item response model of the grade of membership model
(Erosheva, Fienberg & Junker, 2002;
Erosheva, Fienberg & Joutard, 2007) with K classes
for dichotomous correct responses X_{pi}
of person p on item i is as follows (model="GOM")
P(X_{pi}=1 | g_{p1}, \ldots, g_{pK} )=\sum_k \lambda_{ik} g_{pk}
\quad, \quad \sum_{k=1}^K g_{pk}=1
\quad, \quad 0 \leq g_{pk} \leq 1
In most applications (e.g. Erosheva et al., 2007), the grade of
membership function \{g_{pk}\} is assumed to follow a Dirichlet
distribution. In our gom.em implementation
the membership function is assumed to be discretely represented
by a grid u=(u_1, \ldots, u_L) with entries between 0 and 1
(e.g. seq(0,1,length=5) with L=5).
The values g_{pk} of the membership function can then
only take values in \{ u_1, \ldots, u_L \} with the restriction
\sum_k g_{pk} \sum_l \bold{1}(g_{pk}=u_l )=1.
The grid u is specified by using the argument problevels.
The Rasch grade of membership model (model="GOMRasch") poses constraints
on probabilities \lambda_{ik} and membership functions g_{pk}.
The membership
function of person p is parameterized by a location parameter \theta_p
and a variability parameter \xi_p. Each class k is represented by
a location parameter \tilde{\theta}_k. The membership function is defined as
g_{pk} \propto
\exp \left[ - \frac{ (\theta_p - \tilde{\theta}_k)^2 }{2 \xi_p^2 } \right]
The person parameter \theta_p indicates the usual 'ability', while
\xi_p describes the individual tendency to change between classes
1,\ldots,K and their corresponding locations
\tilde{\theta}_1, \ldots,\tilde{\theta}_K.
The extremal class probabilities \lambda_{ik} follow the Rasch model
\lambda_{ik}=invlogit( \tilde{\theta}_k - b_i )=
\frac{ \exp( \tilde{\theta}_k - b_i ) }{ 1 + \exp( \tilde{\theta}_k - b_i ) }
Putting these assumptions together leads to the model equation
P(X_{pi}=1 | g_{p1}, \ldots, g_{pK} )=
P(X_{pi}=1 | \theta_p, \xi_p )=
\sum_k \frac{ \exp( \tilde{\theta}_k - b_i ) }{ 1 + \exp(\tilde{\theta}_k - b_i ) }
\cdot \exp \left[ - \frac{ (\theta_p - \tilde{\theta}_k)^2 }{2 \xi_p^2 } \right]
In the extreme case of a very small \xi_p=\varepsilon > 0 and
\theta_p=\theta_0, the Rasch model is obtained
P(X_{pi}=1 | \theta_p, \xi_p )=
P(X_{pi}=1 | \theta_0, \varepsilon )=
\frac{ \exp( \theta_0 - b_i ) }{ 1 + \exp( \theta_0 - b_i ) }
See Erosheva et al. (2002), Erosheva (2005, 2006) or Galyart (2015) for a comparison of grade of membership models with latent trait models and latent class models.
The grade of membership model is also published under the name Bernoulli aspect model, see Bingham, Kaban and Fortelius (2009).
Value
A list with following entries:
deviance |
Deviance |
ic |
Information criteria |
item |
Data frame with item parameters |
person |
Data frame with person parameters |
EAP.rel |
EAP reliability (only applies for |
MAP |
Maximum aposteriori estimate of the membership function |
EAP |
EAP estimate for individual membership scores |
classdesc |
Descriptives for class membership |
lambda |
Estimated response probabilities |
se.lambda |
Standard error for estimated response probabilities
|
mu |
Mean of the distribution of |
Sigma |
Covariance matrix of |
b |
Estimated item difficulties (only applies for |
se.b |
Standard error of estimated difficulties
(only applies for |
f.yi.qk |
Individual likelihood |
f.qk.yi |
Individual posterior |
probs |
Array with response probabilities |
n.ik |
Expected counts |
iter |
Number of iterations |
I |
Number of items |
K |
Number of classes |
TP |
Number of discrete integration points for |
theta.k |
Used grid of membership functions |
... |
Further values |
References
Bingham, E., Kaban, A., & Fortelius, M. (2009). The aspect Bernoulli model: multiple causes of presences and absences. Pattern Analysis and Applications, 12(1), 55-78.
Erosheva, E. A. (2005). Comparing latent structures of the grade of membership, Rasch, and latent class models. Psychometrika, 70, 619-628.
Erosheva, E. A. (2006). Latent class representation of the grade of membership model. Seattle: University of Washington.
Erosheva, E. A., Fienberg, S. E., & Junker, B. W. (2002). Alternative statistical models and representations for large sparse multi-dimensional contingency tables. Annales-Faculte Des Sciences Toulouse Mathematiques, 11, 485-505.
Erosheva, E. A., Fienberg, S. E., & Joutard, C. (2007). Describing disability through individual-level mixture models for multivariate binary data. Annals of Applied Statistics, 1, 502-537.
Galyardt, A. (2015). Interpreting mixed membership models: Implications of Erosheva's representation theorem. In E. M. Airoldi, D. Blei, E. A. Erosheva, & S. E. Fienberg (Eds.). Handbook of Mixed Membership Models (pp. 39-65). Chapman & Hall.
See Also
For joint maximum likelihood estimation of the grade of membership model
see gom.jml.
See also the mixedMem package for estimating mixed membership models by a variational EM algorithm.
The C code of Erosheva et al. (2007) can be downloaded from http://projecteuclid.org/euclid.aoas/1196438029#supplemental.
Code from Manrique-Vallier can be downloaded from http://pages.iu.edu/~dmanriqu/software.html.
See http://users.ics.aalto.fi/ella/publications/aspect_bernoulli.m for a Matlab implementation of the algorithm in Bingham, Kaban and Fortelius (2009).
Examples
#############################################################################
# EXAMPLE 1: PISA data mathematics
#############################################################################
data(data.pisaMath)
dat <- data.pisaMath$data
dat <- dat[, grep("M", colnames(dat)) ]
#***
# Model 1: Discrete GOM with 3 classes and 5 probability levels
problevels <- seq( 0, 1, len=5 )
mod1 <- sirt::gom.em( dat, K=3, problevels, model="GOM")
summary(mod1)
## Not run:
#-- some plots
#* multivariate scatterplot
car::scatterplotMatrix(mod1$EAP, regLine=FALSE, smooth=FALSE, pch=16, cex=.4)
#* ternary plot
vcd::ternaryplot(mod1$EAP, pch=16, col=1, cex=.3)
#***
# Model 1a: Multivariate normal distribution
problevels <- seq( 0, 1, len=5 )
mod1a <- sirt::gom.em( dat, K=3, theta0.k=seq(-15,15,len=21), model="GOMnormal" )
summary(mod1a)
#***
# Model 2: Discrete GOM with 4 classes and 5 probability levels
problevels <- seq( 0, 1, len=5 )
mod2 <- sirt::gom.em( dat, K=4, problevels, model="GOM" )
summary(mod2)
# model comparison
smod1 <- IRT.modelfit(mod1)
smod2 <- IRT.modelfit(mod2)
IRT.compareModels(smod1,smod2)
#***
# Model 2a: Estimate discrete GOM with 4 classes and restricted space of probability levels
# the 2nd, 4th and 6th class correspond to "intermediate stages"
problevels <- scan()
1 0 0 0
.5 .5 0 0
0 1 0 0
0 .5 .5 0
0 0 1 0
0 0 .5 .5
0 0 0 1
problevels <- matrix( problevels, ncol=4, byrow=TRUE)
mod2a <- sirt::gom.em( dat, K=4, problevels, model="GOM" )
# probability distribution for latent classes
cbind( mod2a$theta.k, mod2a$pi.k )
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0 0.0 0.0 0.0 0.17214630
## [2,] 0.5 0.5 0.0 0.0 0.04965676
## [3,] 0.0 1.0 0.0 0.0 0.09336660
## [4,] 0.0 0.5 0.5 0.0 0.06555719
## [5,] 0.0 0.0 1.0 0.0 0.27523678
## [6,] 0.0 0.0 0.5 0.5 0.08458620
## [7,] 0.0 0.0 0.0 1.0 0.25945016
## End(Not run)
#***
# Model 3: Rasch GOM
mod3 <- sirt::gom.em( dat, model="GOMRasch", maxiter=20 )
summary(mod3)
#***
# Model 4: 'Ordinary' Rasch model
mod4 <- sirt::rasch.mml2( dat )
summary(mod4)
## Not run:
#############################################################################
# EXAMPLE 2: Grade of membership model with 2 classes
#############################################################################
#********* DATASET 1 *************
# define an ordinary 2 latent class model
set.seed(8765)
I <- 10
prob.class1 <- stats::runif( I, 0, .35 )
prob.class2 <- stats::runif( I, .70, .95 )
probs <- cbind( prob.class1, prob.class2 )
# define classes
N <- 1000
latent.class <- c( rep( 1, 1/4*N ), rep( 2,3/4*N ) )
# simulate item responses
dat <- matrix( NA, nrow=N, ncol=I )
for (ii in 1:I){
dat[,ii] <- probs[ ii, latent.class ]
dat[,ii] <- 1 * ( stats::runif(N) < dat[,ii] )
}
colnames(dat) <- paste0( "I", 1:I)
# Model 1: estimate latent class model
mod1 <- sirt::gom.em(dat, K=2, problevels=c(0,1), model="GOM" )
summary(mod1)
# Model 2: estimate GOM
mod2 <- sirt::gom.em(dat, K=2, problevels=seq(0,1,0.5), model="GOM" )
summary(mod2)
# estimated distribution
cbind( mod2$theta.k, mod2$pi.k )
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0.243925644
## [2,] 0.5 0.5 0.006534278
## [3,] 0.0 1.0 0.749540078
#********* DATASET 2 *************
# define a 2-class model with graded membership
set.seed(8765)
I <- 10
prob.class1 <- stats::runif( I, 0, .35 )
prob.class2 <- stats::runif( I, .70, .95 )
prob.class3 <- .5*prob.class1+.5*prob.class2 # probabilities for 'fuzzy class'
probs <- cbind( prob.class1, prob.class2, prob.class3)
# define classes
N <- 1000
latent.class <- c( rep(1,round(1/3*N)),rep(2,round(1/2*N)),rep(3,round(1/6*N)))
# simulate item responses
dat <- matrix( NA, nrow=N, ncol=I )
for (ii in 1:I){
dat[,ii] <- probs[ ii, latent.class ]
dat[,ii] <- 1 * ( stats::runif(N) < dat[,ii] )
}
colnames(dat) <- paste0( "I", 1:I)
#** Model 1: estimate latent class model
mod1 <- sirt::gom.em(dat, K=2, problevels=c(0,1), model="GOM" )
summary(mod1)
#** Model 2: estimate GOM
mod2 <- sirt::gom.em(dat, K=2, problevels=seq(0,1,0.5), model="GOM" )
summary(mod2)
# inspect distribution
cbind( mod2$theta.k, mod2$pi.k )
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0.3335666
## [2,] 0.5 0.5 0.1810114
## [3,] 0.0 1.0 0.4854220
#***
# Model2m: estimate discrete GOM in mirt
# define latent classes
Theta <- scan( nlines=1)
1 0 .5 .5 0 1
Theta <- matrix( Theta, nrow=3, ncol=2,byrow=TRUE)
# define mirt model
I <- ncol(dat)
#*** create customized item response function for mirt model
name <- 'gom'
par <- c("a1"=-1, "a2"=1 )
est <- c(TRUE, TRUE)
P.gom <- function(par,Theta,ncat){
# GOM for two extremal classes
pext1 <- stats::plogis(par[1])
pext2 <- stats::plogis(par[2])
P1 <- Theta[,1]*pext1 + Theta[,2]*pext2
cbind(1-P1, P1)
}
# create item response function
icc_gom <- mirt::createItem(name, par=par, est=est, P=P.gom)
#** define prior for latent class analysis
lca_prior <- function(Theta,Etable){
# number of latent Theta classes
TP <- nrow(Theta)
# prior in initial iteration
if ( is.null(Etable) ){ prior <- rep( 1/TP, TP ) }
# process Etable (this is correct for datasets without missing data)
if ( ! is.null(Etable) ){
# sum over correct and incorrect expected responses
prior <- ( rowSums(Etable[, seq(1,2*I,2)]) + rowSums(Etable[,seq(2,2*I,2)]) )/I
}
prior <- prior / sum(prior)
return(prior)
}
#*** estimate discrete GOM in mirt package
mod2m <- mirt::mirt(dat, 1, rep( "icc_gom",I), customItems=list("icc_gom"=icc_gom),
technical=list( customTheta=Theta, customPriorFun=lca_prior) )
# correct number of estimated parameters
mod2m@nest <- as.integer(sum(mod.pars$est) + nrow(Theta)-1 )
# extract log-likelihood and compute AIC and BIC
mod2m@logLik
( AIC <- -2*mod2m@logLik+2*mod2m@nest )
( BIC <- -2*mod2m@logLik+log(mod2m@Data$N)*mod2m@nest )
# extract coefficients
( cmod2m <- sirt::mirt.wrapper.coef(mod2m) )
# compare estimated distributions
round( cbind( "sirt"=mod2$pi.k, "mirt"=mod2m@Prior[[1]] ), 5 )
## sirt mirt
## [1,] 0.33357 0.33627
## [2,] 0.18101 0.17789
## [3,] 0.48542 0.48584
# compare estimated item parameters
dfr <- data.frame( "sirt"=mod2$item[,4:5] )
dfr$mirt <- apply(cmod2m$coef[, c("a1", "a2") ], 2, stats::plogis )
round(dfr,4)
## sirt.lam.Cl1 sirt.lam.Cl2 mirt.a1 mirt.a2
## 1 0.1157 0.8935 0.1177 0.8934
## 2 0.0790 0.8360 0.0804 0.8360
## 3 0.0743 0.8165 0.0760 0.8164
## 4 0.0398 0.8093 0.0414 0.8094
## 5 0.1273 0.7244 0.1289 0.7243
## [...]
#############################################################################
# EXAMPLE 3: Lung cancer dataset; using sampling weights
#############################################################################
data(data.si08, package="sirt")
dat <- data.si08
#- Latent class model with 3 classes
problevels <- c(0,1)
mod1 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, problevels=problevels )
summary(mod1)
#- Grade of membership model with discrete distribution
problevels <- seq(0,1,length=5)
mod2 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, problevels=problevels )
summary(mod2)
#- Grade of membership model with multivariate normal distribution
mod3 <- sirt::gom.em( dat[,1:5], weights=dat$wgt, K=3, theta0.k=10*seq(-1,1,len=11),
model="GOMnormal", optimizer="nlminb" )
summary(mod3)
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.