WLS to IRT parameters

[Mauricio Garnier-Villarreal]

Hi

I am working with ordinal indicadors with WLSMV. This works well for the estimation part. But for ordinal variables the Item Characteristic Curves, and "a" and "b" parameters are commonly used and interpretable than the WLS parametrization. 

Found the attach Mplus technical paper about the IRT parametrization from CFA. Was able to calculate the "b" parameter based on equation 22, but was unable to replicate the "a" parameter based on equations 24 and 25 for theta and delta parametrization respectively. 

For theta parametrization they refer to a "theta" parameter but I couldnt figure out which is that parameter from lavaan.

And for delta parametrization they refer to a "delta" parameter that again I could figure out from lavaan.

Here is a toy example comparing mirt (for the IRT comparison) to lavaan

Any help is appreciated

library(lavaan)

library(mirt)

data <- expand.table(LSAT7)

head(data)

(mod1 <- mirt(data, 1,type="2PL"))

coef(mod1)

coef(mod1, IRTpars=T) ### this are the target parameters

summary(mod1)

mm <- '

f1 =~ Item.1 + Item.2 + Item.3 + Item.4 + Item.5

'

fit1 <- cfa(mm, data=data, std.lv=T, 

            estimator="wlsmv",ordered =colnames(data) ,

            parameterization="theta")

summary(fit1, standardized=T)

b1 <- (-1.097 - (0.587*0))/(0.587*sqrt(1))

b1

[Hugo Harada]

Maybe you are forgetting the D=1.7 scaling? Ignore the MIRT part.  See the code below and see if it helps.

Also, check the following paper for unidimensional conversions...

@article{KamataBauer2008,

  author  = {Akihito Kamata and Daniel J. Bauer}, 

  title   = {A Note on the Relation Between Factor Analytic and Item Response Theory Models},

  journal = {Structural Equation Modeling},

  year    = 2008,

  pages   = {136–153},

  volume  = 15

}

@article{takane1987relationship,

  title={On the relationship between item response theory and factor analysis of discretized variables},

  author={Takane, Yoshio and De Leeuw, Jan},

  journal={Psychometrika},

  volume={52},

  number={3},

  pages={393--408},

  year={1987},

  publisher={Springer}

}

----------------------------------------------------------------

rm(list=ls())

library(lavaan)

library(mirt)

dat <- expand.table(LSAT6)

head(dat)

#descript(LSAT6)

##########################################################################

# Unidimensional case

##########################################################################

#estimating parameters using MIRT.

mirt.lSAT6.2PL = mirt(dat, 1, itemtype = '2PL', method = 'EM')  # 

PAR=coef(mirt.lSAT6.2PL,simplify=TRUE,IRTpars=TRUE)

PAR$items[,1]/1.7 #normal ogive scaling

PAR

#estimating parameters using lavaan.

lavaan.lsat6.2pl.model <-'

Theta =~ l1*Item_1 + l2*Item_2 + l3*Item_3 + l4*Item_4 + l5*Item_5

Item_1 | th1 *t1

Item_2 | th2 *t1

Item_3 | th3 *t1

Item_4 | th4 *t1

Item_5 | th5 *t1

'

#fitting model

lavaan.lsat6.2pl.model.fit <- cfa(lavaan.lsat6.2pl.model, data = dat , std.lv=TRUE , ordered =c("Item_1","Item_2","Item_3","Item_4","Item_5"))

#getting lambda values

lambda <- lavInspect(lavaan.lsat6.2pl.model.fit,what = 'est')$lambda

#getting tau values

tau <- lavInspect(lavaan.lsat6.2pl.model.fit,what = 'est')$tau

# Convert regression to IRT

# alpha1 := ( l1)/ sqrt (1- l1^2)

# beta1 := ( - th1 )/ sqrt (1- l1^2)

# a_IRT := alpha1*1.7

# b_IRT := -beta1/alpha1

a_sem <- (lambda/sqrt(1-lambda*lambda))

b_sem <- tau/sqrt(1-lambda*lambda)/a_sem

a_sem <- 1.7*a_sem #normal ogive scaling

#making comparison tables. a and b values are close.

a_comp <- cbind(a_sem,PAR$items[,1])

colnames(a_comp)<-c("a_sem","a_tri")

a_comp

b_comp <- cbind(b_sem,PAR$items[,2])

colnames(b_comp)<-c("b_sem","b_tri")

b_comp

[Mauricio Garnier-Villarreal]

Oh, right, the IRT model works with the logit link, while WLS works with the probit link

Thank you so much

[Balal Ezanloo]

Hi hugo

I find you post in lavaan group very useful. but i think some of your commands are Vague. for example D=1.702 is not used for b_sem parameters at all and you Unnecessarily complicatethe the bellow term. do you have any idea for it?

b_sem <- tau/sqrt(1-lambda*lambda)/a_sem

[Christie Barron]

Hi all,

In case others have a similar question, I'm sharing a function developed by Jonathan Templin that creates item response functions from lavaan's cfa() output with ordered categorical items when using DWLS. Assumes theta parameterization. 

The syntax is from the following notebook: https://jonathantemplin.com/wp-content/uploads/2017/08/EPSY906_Example6A_Ordinal_IFA_Models.nb_.html

Templin describes the function during a lecture available on his youtube channel: https://www.youtube.com/watch?v=txIEpF-rXJE 

lavaanCatItemPlot = function(lavObject, varname, sds = 3){ output = inspect(object = lavObject, what = "est") if (!varname %in% rownames(output$lambda)) stop(paste(varname, "not found in lavaan object")) if (dim(output$lambda)[2]>1) stop("plots only given for one factor models") itemloading = output$lambda[which(rownames(output$lambda) == varname),1] itemthresholds = output$tau[grep(pattern = varname, x = rownames(output$tau))] factorname = colnames(output$lambda) factormean = output$alpha[which(rownames(output$alpha) == factorname)] factorvar = output$psi[which(rownames(output$psi) == factorname)] factormin = factormean - 3*sqrt(factorvar) factormax = factormean + 3*sqrt(factorvar) factorX = seq(factormin, factormax, .01) itemloc = which(lavObject@Data@ov$name == varname) itemlevels = unlist(strsplit(x = lavObject@Data@ov$lnam[itemloc], split = "\\|")) if (length(itemthresholds)>1){ plotdata = NULL plotdata2 = NULL itemY = NULL itemY2 = NULL itemX = NULL itemText = NULL for (level in 1:length(itemthresholds)){ itemY = pnorm(q = -1*itemthresholds[level] + itemloading*factorX) itemY2 = cbind(itemY2, pnorm(q = -1*itemthresholds[level] + itemloading*factorX)) itemText = paste0("P(", varname, " > ", itemlevels[level], ")") itemText2 = paste0("P(", varname, " = ", itemlevels[level], ")") plotdata = rbind(plotdata, data.frame(factor = factorX, prob = itemY, plot = itemText)) if (level == 1){ plotdata2 = data.frame(factor = factorX, plot = itemText2, prob = matrix(1, nrow = dim(itemY2)[1], ncol=1) - itemY2[,level]) } else if (level == length(itemthresholds)){ plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = itemText2, prob = itemY2[,level-1] - itemY2[,level])) plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = paste0("P(", varname, " = ", itemlevels[level+1], ")"), prob = itemY2[,level])) } else { plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = itemText2, prob = itemY2[,level-1] - itemY2[,level])) } } names(plotdata) = c(factorname , "Probability", "Cumulative") ggplot(data = plotdata, aes_string(x = factorname, y = "Probability", colour = "Cumulative")) + geom_line(size = 2) names(plotdata2) = c(factorname, "Response", "Probability") ggplot(data = plotdata2, aes_string(x = factorname, y = "Probability", colour = "Response")) + geom_line(size = 2) } else { itemY = pnorm(q = -1*itemthresholds[1] + itemloading*factorX) itemText2 = paste0("P(", varname, " = ", itemlevels[1], ")") plotdata = data.frame(factor = factorX, prob = itemY, plot = itemText2) names(plotdata) = c(factorname , "Probability", "Response") ggplot(data = plotdata, aes_string(x = factorname, y = "Probability", colour = "Response")) + geom_line(size = 2) } } lavaanCatItemPlot(lavObject = grm2Pestimates, varname = "cia2", sds = 3)