# gom.em: Discrete (Rasch) Grade of Membership Model In sirt: Supplementary Item Response Theory Models

## 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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 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{θ}_k grid (applies only for model="GOMRasch") xsi0.k Vector of ξ_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 λ 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")

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, …, 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, …, u_L \} with the restriction ∑_k g_{pk} ∑_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 λ_{ik} and membership functions g_{pk}. The membership function of person p is parameterized by a location parameter θ_p and a variability parameter ξ_p. Each class k is represented by a location parameter \tilde{θ}_k. The membership function is defined as

g_{pk} \propto \exp ≤ft[ - \frac{ (θ_p - \tilde{θ}_k)^2 }{2 ξ_p^2 } \right]

The person parameter θ_p indicates the usual 'ability', while ξ_p describes the individual tendency to change between classes 1,…,K and their corresponding locations \tilde{θ}_1, …,\tilde{θ}_K. The extremal class probabilities λ_{ik} follow the Rasch model

λ_{ik}=invlogit( \tilde{θ}_k - b_i )= \frac{ \exp( \tilde{θ}_k - b_i ) }{ 1 + \exp( \tilde{θ}_k - b_i ) }

Putting these assumptions together leads to the model equation

P(X_{pi}=1 | g_{p1}, …, g_{pK} )= P(X_{pi}=1 | θ_p, ξ_p )= ∑_k \frac{ \exp( \tilde{θ}_k - b_i ) }{ 1 + \exp(\tilde{θ}_k - b_i ) } \cdot \exp ≤ft[ - \frac{ (θ_p - \tilde{θ}_k)^2 }{2 ξ_p^2 } \right]

In the extreme case of a very small ξ_p=\varepsilon > 0 and θ_p=θ_0, the Rasch model is obtained

P(X_{pi}=1 | θ_p, ξ_p )= P(X_{pi}=1 | θ_0, \varepsilon )= \frac{ \exp( θ_0 - b_i ) }{ 1 + \exp( θ_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 λ_{ik} se.lambda Standard error for estimated response probabilities λ_{ik} mu Mean of the distribution of (θ_p, ξ_p) (only applies for model="GOMRasch") Sigma Covariance matrix of (θ_p, ξ_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.

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.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 ############################################################################# # 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)