This vignette illustrates how to estimate a longitudinal IRT following the methodology described in Proust-Lima et al. (2021 - https://arxiv.org/abs/2109.13064). The model combines a graded response measurement model to link a series of ordinal items measured repeatedly over time with their underlying latent process. Simultaneously, a linear mixed model describes the trajectory of the underlying latent process over time.

Continuous-time longitudinal IRT model estimation

We consider here a longitudinal IRT model with natural cubic splines on time on the waiting list to account for a possible nonlinear trajectory over time, and we adjust the trajectory for the group. We consider 2 internal knots placed at 7 and 15, and shift the right boundary at 60 due to the long tail of the distribution. In the main analysis, we estimate a model with no differential item functioning (DIF) on the group and no response shift (RS) on time.

Estimation of the model with no DIF/RS

To reduce the computation time, we first estimate a model with only 1 random-effect, a random intercept. And then, we use the estimates as initial values for the estimation of the model that also includes random-effects on the time functions. Estimation involves a numerical integration approximated by quasi Monte-Carlo (QMC) with 1000 points. This induces very intensive and long computations but was shown to give accurate results in simulations.

modIRT_i <- multlcmm(hads_2 + hads_4 +hads_6 + hads_8 +hads_10+hads_12 + hads_14 ~ ns(time,knots=c(7,15),Boundary.knots = c(0,60))*grp,random=~1,data=simdataHADS,subject="ID",link="thresholds",methInteg="QMC",nMC=1000)
# use the estimates as initial values - the vector c(0,1,0,0,1,0,0,0,1) initializes the cholesky matrix of the random-effects
Binit <- c(modIRT_i$best[1:7],c(0,1,0,0,1,0,0,0,1),modIRT_i$best[8:length(modIRT_i$best)])

modIRT <- multlcmm(hads_2 + hads_4 +hads_6 + hads_8 +hads_10+hads_12 + hads_14 ~ ns(time,knots=c(7,15),Boundary.knots = c(0,60))*grp,random=~1+ns(time,knots=c(7,15),Boundary.knots = c(0,60)),data=simdataHADS,subject="ID",link="thresholds",methInteg="QMC",nMC=1000, B=Binit)

The summary of the estimation:

summary(modIRT)
#> General latent class mixed model 
#>      fitted by maximum likelihood method 
#>  
#> multlcmm(fixed = hads_2 + hads_4 + hads_6 + hads_8 + hads_10 + 
#>     hads_12 + hads_14 ~ ns(time, knots = c(7, 15), Boundary.knots = c(0, 
#>     60)) * grp, random = ~1 + ns(time, knots = c(7, 15), Boundary.knots = c(0, 
#>     60)), subject = "ID", link = "thresholds", data = simdataHADS, 
#>     methInteg = "QMC", nMC = 1000)
#>  
#> Statistical Model: 
#>      Dataset: simdataHADS 
#>      Number of subjects: 561 
#>      Number of observations: 7980 
#>      Number of latent classes: 1 
#>      Number of parameters: 44  
#>      Link functions: Thresholds for hads_2 
#>                      Thresholds for hads_4 
#>                      Thresholds for hads_6 
#>                      Thresholds for hads_8 
#>                      Thresholds for hads_10 
#>                      Thresholds for hads_12 
#>                      Thresholds for hads_14 
#>  
#> Iteration process: 
#>      Convergence criteria satisfied 
#>      Number of iterations:  30 
#>      Convergence criteria: parameters= 2e-08 
#>                          : likelihood= 3.2e-06 
#>                          : second derivatives= 1.8e-06 
#>  
#> Goodness-of-fit statistics: 
#>      maximum log-likelihood: -7218.96  
#>      AIC: 14525.91  
#>      BIC: 14716.42  
#>  
#> Maximum Likelihood Estimates: 
#>  
#> Fixed effects in the longitudinal model:
#> 
#>                                                                coef      Se
#> intercept (not estimated)                                   0.00000        
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1     -0.13037 0.14971
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2      0.20745 0.35059
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3      0.55413 0.20665
#> grp                                                        -0.60195 0.19524
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1:grp  0.23105 0.20314
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2:grp  0.46499 0.44974
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3:grp -0.06219 0.32068
#>                                                              Wald p-value
#> intercept (not estimated)                                                
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1     -0.871 0.38387
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2      0.592 0.55404
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3      2.681 0.00733
#> grp                                                        -3.083 0.00205
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1:grp  1.137 0.25537
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2:grp  1.034 0.30118
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3:grp -0.194 0.84622
#> 
#> 
#> Variance-covariance matrix of the random-effects:
#> (the variance of the first random effect is not estimated)
#>                                                        intercept
#> intercept                                                1.00000
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1  -0.34385
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2  -1.00748
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3  -0.09015
#>                                                        ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1
#> intercept                                                                                                    
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1                                                0.65616
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2                                                0.79507
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3                                                0.40139
#>                                                        ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2
#> intercept                                                                                                    
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1                                                       
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2                                                2.13117
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3                                                0.30796
#>                                                        ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3
#> intercept                                                                                                    
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))1                                                       
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))2                                                       
#> ns(time, knots = c(7, 15), Boundary.knots = c(0, 60))3                                                 0.2747
#> 
#>                            hads_2   hads_4   hads_6   hads_8  hads_10  hads_12
#> Residual standard error:  0.78685  0.63468  1.03439  1.01073  1.05546  0.67910
#>                           hads_14
#> Residual standard error:  1.63853
#> 
#> Parameters of the link functions:
#> 
#>                     coef      Se   Wald p-value
#> hads_2-Thresh1  -0.54296 0.16481 -3.294 0.00099
#> hads_2-Thresh2   1.12581 0.05349 21.046 0.00000
#> hads_2-Thresh3   0.88949 0.05032 17.678 0.00000
#> hads_4-Thresh1  -0.23896 0.15916 -1.501 0.13325
#> hads_4-Thresh2   1.01287 0.04844 20.909 0.00000
#> hads_4-Thresh3   1.04562 0.05856 17.856 0.00000
#> hads_6-Thresh1  -0.38855 0.16408 -2.368 0.01788
#> hads_6-Thresh2   1.35285 0.06836 19.790 0.00000
#> hads_6-Thresh3   1.32439 0.09443 14.025 0.00000
#> hads_8-Thresh1  -1.39105 0.20071 -6.931 0.00000
#> hads_8-Thresh2   1.34308 0.06446 20.836 0.00000
#> hads_8-Thresh3   1.10197 0.05825 18.918 0.00000
#> hads_10-Thresh1 -0.02481 0.16367 -0.152 0.87949
#> hads_10-Thresh2  1.00550 0.05464 18.403 0.00000
#> hads_10-Thresh3  1.10645 0.06561 16.863 0.00000
#> hads_12-Thresh1 -0.20523 0.15916 -1.289 0.19724
#> hads_12-Thresh2  0.99398 0.04808 20.675 0.00000
#> hads_12-Thresh3  1.02543 0.05776 17.753 0.00000
#> hads_14-Thresh1  0.91742 0.20988  4.371 0.00001
#> hads_14-Thresh2  1.49255 0.09869 15.124 0.00000
#> hads_14-Thresh3  0.87180 0.10143  8.595 0.00000

Predicted underlying depressive symptomatolgy trajectory

plot of predicted mean trajectories

The predicted trajectory of the underlying process can be obtained with predictL function and the associated plot function.

datnew <- data.frame(time = seq(0,75,by=1))
datnew$grp <- 0
pIRT0 <- predictL(modIRT,datnew,var.time="time",confint = T)
datnew$grp <- 1
pIRT1 <- predictL(modIRT,datnew,var.time="time",confint = T)
plot(pIRT0,col=1,lwd=2,ylim=c(-1.5,1.5),legend=NULL,main="",ylab="latent depressive symptomatology",xlab="months since entry on the waiting list",type="l",bty="l",shades=T)
plot(pIRT1,add=T,col=4,lwd=2,shades=T)
legend(x="topleft",legend=c("dialysed","preemptive"),lty=c(1,1),col=c(1,4),lwd=2,bty="n")

Posteriori distribution

To better appreciate the range of the underlying depressive symptomatology, the empirical posterior distribution can be computed

beta <- modIRT$best
t <- 0:72
Z <- cbind(rep(1,length(t)),ns(t,knots=c(7,15),Boundary.knots = c(0,60)))
chol <- matrix(0,ncol=4,nrow=4)
chol[upper.tri(chol, diag = T)] <- c(1,beta[7:15])
library(mvtnorm)
Lambda0 <- rmvnorm(10000,mean = Z%*%c(0,beta[1:3]),Z%*%t(chol)%*%chol%*%t(Z))
Lambda1 <- rmvnorm(10000,mean = Z%*%beta[4:7],Z%*%t(chol)%*%chol%*%t(Z))
Lambda <- data.frame(Lambda = as.vector(rbind(Lambda0,Lambda1)))
phist <- ggplot(Lambda,aes(x=Lambda))+ geom_density(color="grey", fill="grey") + theme_bw() +
  xlab("underlying depressive symptomatology") +xlim(-7,7)
phist
Warning: Removed 14464 rows containing non-finite values (stat_density).

The 95% range of the underlying distribution is approximately:

quantile(Lambda$Lambda,p=c(0.025,0.975))
#>      2.5%     97.5% 
#> -4.589935  4.663775

Location and discrimination of the items (with SE by Delta-Method)

The location and discrimination parameters are transformations of the estimated parameters. They are retrieved with the following code:

isolate the corresponding parameters within the vector of estimated parameters (and the corresponding variance covariance matrix)

## Parameters
nlevel <- 4
nitems <- 7
levels <- rep(nlevel,nitems)
npm <- length(modIRT$best)
seuils <- modIRT$best[(npm-(nlevel-1)*(nitems)+1):(npm)]
err <- abs(modIRT$best[(npm-(nlevel-1)*(nitems)-(nitems-1)):(npm-(nlevel-1)*(nitems))])
seuils
#>     Thresh1     Thresh2     Thresh3     Thresh1     Thresh2     Thresh3 
#> -0.54295998  1.12581248  0.88948647 -0.23896189  1.01286765  1.04561833 
#>     Thresh1     Thresh2     Thresh3     Thresh1     Thresh2     Thresh3 
#> -0.38854575  1.35284533  1.32439060 -1.39104501  1.34308298  1.10196668 
#>     Thresh1     Thresh2     Thresh3     Thresh1     Thresh2     Thresh3 
#> -0.02481498  1.00550201  1.10644595 -0.20523446  0.99397714  1.02543122 
#>     Thresh1     Thresh2     Thresh3 
#>  0.91741673  1.49254985  0.87180230
err
#> std.err 1 std.err 2 std.err 3 std.err 4 std.err 5 std.err 6 std.err 7 
#> 0.7868467 0.6346774 1.0343869 1.0107349 1.0554643 0.6790975 1.6385266
# Variance
Vseuils <- VarCov(modIRT)[(npm-(nlevel-1)*(nitems)+1):(npm),(npm-(nlevel-1)*(nitems)+1):(npm)]
Verr <- VarCov(modIRT)[(npm-(nlevel-1)*(nitems)-(nitems-1)):(npm-(nlevel-1)*(nitems)),(npm-(nlevel-1)*(nitems)-(nitems-1)):(npm-(nlevel-1)*(nitems))]

retransform these parameters to obtain the locations and the discrimination of each item with corresponding standard deviation (computed via Delta-method)

# generic function
location <- function(min,max,param,Vparam){
  loc <- param[1]
  se <- sqrt(Vparam[1,1])
  param2 <- rep(0,length(param))
  param2[1] <- 1
  if ((max-min)>1) {
    for (l in 1:(max-min-1)) {
      param2[l+1] <- 2*param[l+1]
      loc[l+1] <- loc[l] + param[1+l]^2
      se[l+1] <- sqrt(t(param2) %*%Vparam %*%param2)
    }
  }
  return(c(loc,se))
}
# application
ItemLoc <- sapply(1:nitems,function(k){location(min=0,max=nlevel-1,param=seuils[((nlevel-1)*(k-1)+1):((nlevel-1)*k)],Vparam=Vseuils[((nlevel-1)*(k-1)+1):((nlevel-1)*k),((nlevel-1)*(k-1)+1):((nlevel-1)*k)])})
colnames(ItemLoc) <- paste("Item",(1:nitems)*2)
ItemLoc <- ItemLoc[c(1,4,2,5,3,6),]
rownames(ItemLoc) <- rep(c("Threshold","SE"),nlevel-1)
discrimination <- 1/abs(err)
sediscr <- diag(err^(-2))%*%Verr%*%diag(err^(-2))

The 3 thresholds and discrimination estimates of each item are:

t(rbind(ItemLoc,discrimination,Se=sqrt(diag(sediscr))))
#>           Threshold        SE Threshold        SE Threshold        SE
#> Item 2  -0.54295998 0.1648116 0.7244938 0.1761869  1.515680 0.2184180
#> Item 4  -0.23896189 0.1591583 0.7869390 0.1767165  1.880257 0.2422033
#> Item 6  -0.38854575 0.1640761 1.4416447 0.2220543  3.195655 0.3954593
#> Item 8  -1.39104501 0.2007092 0.4128269 0.1685202  1.627157 0.2307004
#> Item 10 -0.02481498 0.1636731 0.9862193 0.1968426  2.210442 0.2870308
#> Item 12 -0.20523446 0.1591645 0.7827561 0.1782252  1.834265 0.2411234
#> Item 14  0.91741673 0.2098824 3.1451218 0.4376183  3.905161 0.5382206
#>         discrimination         Se
#> Item 2       1.2708956 0.12409358
#> Item 4       1.5756036 0.15681998
#> Item 6       0.9667562 0.10083591
#> Item 8       0.9893791 0.09850672
#> Item 10      0.9474504 0.10031786
#> Item 12      1.4725427 0.14449255
#> Item 14      0.6103044 0.07789280

item category probability curve

The curve of each item category probability according to the underlying level of depressive symptomatology can be obtain usinf the ItemInfo function.

## computations
info_modIRT <- ItemInfo(modIRT, lprocess=seq(-6,6,0.1))

meaning <- c("Enjoy","Laugh","Cheerful" ,"Slow" ,"Appearance" ,"Looking Forward" ,"Leisure")
items <- paste("hads", seq(2,14,2), sep="_")

## automatic graph
par(mfrow=c(2,4), mar=c(3,2,2,1), mgp=c(2,0.5,0))
for(k in 1:7)
{     
 plot(info_modIRT, which="LevelProb", outcome=items[k],
      main=paste("Item",2*k,"-",meaning[k]), lwd=2, legend=NULL)
}
plot(0,axes=FALSE, xlab="", ylab="", col="white")
legend("center", legend=paste("Level",0:3), col=1:4, lwd=2, box.lty=0)

## graph with ggplot
p <- NULL
for (k in 1:7){
ICC  <- info_modIRT$LevelInfo[which(info_modIRT$LevelInfo[,1]==items[k]),]
p[[k]] <- (ggplot(ICC)
      + geom_line(aes(x = Lprocess, y = Prob, group = Level,color=Level), show.legend = F,alpha = 1,size=1.2)
      # + stat_smooth(aes(x = time, y = hads_scorea), method = "loess", size = 0.75)
      + theme_bw()
      + labs(title=paste("Item",2*k,"-",meaning[k]))
      + xlab("construct")
      + ylim(0,1)
      + ylab("Probability of item level"))
}
p[[8]] <- (ggplot(ICC)
      + geom_line(aes(x = Lprocess, y = Prob, group = Level,color=Level),alpha = 1,size=1.2)
      + theme_bw()
)
legend <- get_legend(p[[8]])
grid.arrange(p[[1]],p[[2]],p[[3]],p[[4]],p[[5]],p[[6]],p[[7]],as_ggplot(legend),ncol=4)

Item characteristic curves

The following code computes the expectation of each item according to the underlying level of depressive symptomatology. This is achieved with predictYcond function with two plot possibilities: direct plot function or ggplot

expe <- predictYcond(modIRT,lprocess = seq(-6,6,by=0.1))
# via the internal plot function 
plot(expe, xlab="underlying depressive symptomatology", main="Item Expectation Curves",
     legend=paste("Item",(1:nitems)*2), lwd=2)

# via ggplot
j <- table(expe$pred$Yname)[1]
expe$pred$item <- as.factor(c(rep(2,j),rep(4,j),rep(6,j),rep(8,j),rep(10,j),rep(12,j),rep(14,j)))
p <- (ggplot(expe$pred)
      + geom_line(aes(x = Lprocess, y = Ypred, group=item,color=item),alpha = 1,size=1.2)
      + theme_bw()
      + xlab("underlying depressive symptomatology")
      + ylim(0,3)
      + ylab("Item Expectation"))
p

Item Information Function

The level of information brought by each item category (information share) and brought in total by each item is also computed by the ItemInfo function. The curves can be again plotted directly with options which=“LevelInfo” and which=“ItemInfo” respectively.

by Category

par(mfrow=c(2,4), mar=c(3,2,2,1), mgp=c(2,0.5,0))
for(k in 1:7)
{     
 plot(info_modIRT, which="LevelInfo", outcome=items[k],
 main=paste("Item",2*k,"-",meaning[k]), lwd=2, legend=NULL, ylim=c(0,1.3))
}
plot(0,axes=FALSE, xlab="", ylab="", col="white")
legend("center", legend=paste("Level",0:3), col=1:4, lwd=2, box.lty=0)

by Item

plot(info_modIRT, which="ItemInfo", lwd=2, legend.loc="topleft")

Predicted item trajectory according to time

Item predicted trajectories according to a specific profile of covariates can be computed using predictY function:

head(datnew)
datnew$grp <- 0
ns0 <- predictY(modIRT,var.time = "time",newdata=datnew,methInteg = 1,nsim=2000,draws=T)
datnew$grp <- 1
ns1 <- predictY(modIRT,var.time = "time",newdata=datnew,methInteg = 1,nsim=2000,draws=T)

par(mfrow=c(2,4), mar=c(3,2,2,1), mgp=c(2,0.5,0))
for(k in 1:7){
plot(ns0,outcome = k,shades = T,ylim=c(0,3),bty="l",legend=NULL,main=paste("Item",2*k,"-",meaning[k]),ylab="Item level",xlab="months on the waiting list")
plot(ns1,outcome=k,shades=T,add=T,col=2)
}

Continuous-time longitudinal IRT model of PREDIALA study

Importation

Dataset

Description of the sample

Timescale of interest: time on the waiting list