knitr::opts_chunk$set(echo = TRUE)

This is an example of multivariate Latent Process model implemented in mulLPM, estimated on the paquid dataset available in lcmm package.

We define :

Installation of the package

devtools::install_github("https://github.com/VivianePhilipps/multLPM.git")
library(lcmm)
library(multLPM)

Data creation

We change the direction of the cognitive scores so that the lower the better

paquid$invMMSE <- 30-paquid$MMSE
paquid$invIST <- 40-paquid$IST

We define age of diagnosis (the variable agedem is the estimated age of dementia)

agediag <- sapply(paquid$ID, function(i) {
  j <- which(paquid$ID==i & (paquid$age-paquid$agedem>0))
  res <- Inf
  if(length(j)) res <- min(paquid$age[j])
  return(res)})

paquid$agediag <- NA
paquid$agediag[which(paquid$dem == 0)] <- paquid$agedem[which(paquid$dem == 0)]
paquid$agediag[which(paquid$dem == 1)] <- agediag[which(paquid$dem == 1)]

Longitudinal dementia indicator (missing after diagnosis)

paquid$demlong <- NA
paquid$demlong[which((paquid$dem == 0))] <- 0
paquid$demlong[which((paquid$dem == 1) & (paquid$age < paquid$agediag))] <- 0
paquid$demlong[which((paquid$dem == 1) & (paquid$age == paquid$agediag))] <- 1

Do not use measures after diagnosis

after <- which(paquid$age>paquid$agediag)
paquid[after,c("invMMSE","invIST","CESD")] <- NA

One dimensional submodels

m01 <- lcmm::multlcmm(invMMSE+invIST~I((age-65)/10)*CEP, random=~1+I((age-65)/10), subject="ID", data=paquid, link="3-quant-splines")

m02 <- lcmm::multlcmm(CESD~I((age-65)/10)*male, random=~1+I((age-65)/10), subject="ID", data=paquid, link="3-quant-splines")

Number of subjects used in each submodel

m01$ns 
m02$ns 

We need to have the same number of subjects in JointMult function, so we restrict to subjects used for both dimensions.

subjects_m01_1 <- unique(paquid$ID[-unlist(m01$na.action[[1]])])
subjects_m01_2 <- unique(paquid$ID[-unlist(m01$na.action[[2]])])
subjects_m01 <- union(subjects_m01_1, subjects_m01_2)
subjects_m02 <- unique(paquid$ID[-unlist(m02$na.action)])

paquid498 <- paquid[which(paquid$ID %in% intersect(subjects_m01, subjects_m02)),]

Multivariate models

Without adjustment on the dementia process

load("models_vignette_multLPM.RData")
m <- JointMult(Y=list(m01,m02), D=list(longDiag(age_init, age, demlong)~1), var.time="age", RE="full", delayed=TRUE, data=paquid498, maxiter=50)
summary(m)

With adjustment on the dementia process, using estimation of model m as initial values

m_adj <- JointMult(Y=list(m01,m02), D=list(longDiag(age_init, age, demlong)~CEP+male), var.time="age", RE="full", delayed=TRUE, data=paquid498, B=c(m$bopt[1:33],0,0,m$bopt[34],0,0,m$bopt[35],0,0), maxiter=50)
summary(m_adj)

The model did not converge. We fix a transformation parameter that is closed to 0 in order to avoid numerical problems :

m_adj <- JointMult(Y=list(m01,m02), D=list(longDiag(age_init, age, demlong)~CEP+male), var.time="age", RE="full", delayed=TRUE, data=paquid498, B=m_adj$bopt, posfix=11, maxiter=30)
summary(m_adj)

Postfit computations

Comparison of subject-specific predicitons and observations

fit <- fit_ss(JointMult(Y=list(m01,m02), D=list(longDiag(age_init, age, demlong)~1), var.time="age", RE="full", delayed=TRUE, data=paquid498, B=m$bopt),
       data=paquid498,id="ID",y=c("invMMSE","invIST","CESD"),time=rep("age",3),breaks=c(seq(65,95,2),105))

par(mfrow=c(1,3))
 plot.fit_ss(fit,which="fit",type=c("p","l","l","l"),lty=c(1,1,2,2),pch=c(19,NA,NA,NA),col=1,outcome=1,xlab="age",ylab="transformed MMSE")
 plot.fit_ss(fit,which="fit",type=c("p","l","l","l"),lty=c(1,1,2,2),pch=c(19,NA,NA,NA),col=1,outcome=2,xlab="age",ylab="transformed IST")
 plot.fit_ss(fit,which="fit",type=c("p","l","l","l"),lty=c(1,1,2,2),pch=c(19,NA,NA,NA),col=1,outcome=3,xlab="age",ylab="transformed CESD")

Predicted trajectories

Predicted trajectories in each dimension for different profiles of covariate :

d00 <- data.frame(age=seq(65,100,length.out=100), CEP=0, male=0)
d10 <- data.frame(age=seq(65,100,length.out=100), CEP=1, male=0)
d01 <- data.frame(age=seq(65,100,length.out=100), CEP=0, male=1)

predL00 <- predL(m=m,newdata=d00,plot=FALSE,m1=m01,m2=m02,xlim=c(65,95),ylim=c(-2,2), confint=TRUE)
predL10 <- predL(m=m,newdata=d10,plot=FALSE,m1=m01,m2=m02,xlim=c(65,95),ylim=c(-2,2), confint=TRUE)
predL01 <- predL(m=m,newdata=d01,plot=FALSE,m1=m01,m2=m02,xlim=c(65,95),ylim=c(-2,2), confint=TRUE)

par(mfrow=c(1,2))
plot(predL00[[1]],lwd=2,xlim=c(65,100),ylim=c(-3,6),legend=NULL,main="cognition",shades=TRUE)
plot(predL10[[1]],lwd=2,xlim=c(65,100),ylim=c(-3,6),legend=NULL,add=TRUE,shades=TRUE,col=2)
legend(x="toplef",legend=c("low educated", "high educated"), lwd=2, col=1:2, box.lty=0)
plot(predL00[[2]],lwd=2,xlim=c(65,100),ylim=c(-3,6),legend=NULL,main="depression",shades=TRUE)
plot(predL01[[2]],lwd=2,xlim=c(65,100),ylim=c(-3,6),legend=NULL,add=TRUE,shades=TRUE,col=3)
legend(x="toplef",legend=c("women", "men"), lwd=2, col=c(1,3), box.lty=0)

The degradation process is computed as the linear combination of the two latent processes and plotted against time :

degr00 <- m$bopt[34] * predL00[[1]]$pred[,1] + m$bopt[35] * predL00[[2]]$pred[,1]
degr10 <- m$bopt[34] * predL10[[1]]$pred[,1] + m$bopt[35] * predL10[[2]]$pred[,1]
degr01 <- m$bopt[34] * predL01[[1]]$pred[,1] + m$bopt[35] * predL01[[2]]$pred[,1]

plot(x=seq(65,100,length.out=100), y=degr00, type="l", col=1, lwd=2, xlab="age", ylab="degradation process", ylim=c(-1,4))
lines(x=seq(65,100,length.out=100), y=degr10, type="l", col=2, lwd=2)
lines(x=seq(65,100,length.out=100), y=degr01, type="l", col=3, lwd=2)
abline(h=m$bopt[33], lty=2, col=1)

The black and green lines cross the threshold at about 93 years, suggesting that low educated subjects (women in black and men in green) will be demented at 93 years, whereas high educated women (red line) will not be diagnosed before 100 years.



VivianePhilipps/multLPM documentation built on March 13, 2021, 6:35 a.m.