knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
In this section we consider a latent variable $y^_{it}$ such that \begin{equation} y^{it} = f(a{it}; \beta_k, \delta_k)+\epsilon_{it} =\beta_k A_{it}+\delta_k W_t+\epsilon_{it} \end{equation} where $\epsilon_{it}\sim \mathcal{N}\left(0;\ \sigma_k\right)$, $A_{it}=(1,a_{it},a_{it}^2,\cdots,a_{it}^{n_\beta-1})^t $, $W_t=(w_{i1},\cdots,w_{in_\delta})^t $, $\beta_k =(\beta_{k1},\cdots, \beta_{kn_\beta})$ and $\delta_k=(\delta_{k1},\cdots,\delta_{kn_\delta}) $.\ It is then classically assumed that the binary variable $y_{it}=1$ if $y^{it}>0$ and $y{it}=0$ if $y^_{it}\leq 0$.\
If $\epsilon_{it}$ is assumed to follow a logistic distribution, we obtain the logit model. Let $\rho_{ikt}=P(Y_{it}=1|W_i=w_i,C_i=k)$ be the probability of $y_{it}=1$ given membership in group $k$. \begin{equation} \rho_{ikt}=\dfrac{e^{\beta_k A_{it}+\delta_k W_{it}}}{1+e^{\beta_k A_{it}+\delta_k W_{it}}} \end{equation}
We use artificial data contained in the library. These are 500 longitudinal data points of 0 and 1 values generated by the parameters below:
The data contains the time-dependent variable $Y_i$ in data[,2:11], the time variable in data[,12:21], a time-dependent covariate in data[,24:33], and a time-invariant covariate that influences the group membership probability in data[,48:49].
data[,2:11] is a matrix with 0 and 1 values.data[,12:21] is a matrix with time 1 to 10.data[,24:33] is a time-dependent covariate that influences the shape of the trajectory. It is a matrix with 0 and 1 values, representing the presence or absence of a characteristic for the individual.data[,48:49] is a matrix with 0 and 1 values, representing a time-invariant covariate.library(trajeR) data(data_LOGIT) matplot( t(data_LOGIT[, 12:21]),t(data_LOGIT[, 2:11]), pch = 1, type='b', col="black", lty=1, xlab = "Times", ylab = "Values", main = "Plot of the individual's trajectories")
We use the each method to fit the model. For all method, we specify the number of group of our model, ng=3, the degree of the polynomial shape of the trajectories. Here we choice a line parallel to abscissa axis, a cubic polynomial and a quadric polynomial. So degree is vector $(0,3,4)$.
We specify hessian=TRUE to ask the calculus of the hessian matrix.
Itermax is set to 300 to assure a good approximation of the parameters.
For the Likelihood method we call trajeR with option Method ="L".
solL <- trajeR( Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "L", hessian = TRUE ) solL
For the EM method, we use the same parameters except Method = EM or EMIRLS. We specifiy itermax = 300 to ensure convergence.
#EM solEM <- trajeR(Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "EM", hessian = TRUE, itermax = 300) #EMIRLS solEMIRLS <- trajeR(Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], ng = 3, degre=c(0,3,4), Model = "LOGIT", Method = "EMIRLS", hessian = TRUE, itermax = 300) solEM solEMIRLS
We can add risk covariate that influence the belonging probability. By default, the effect of risk variable are compared to the first group that is reference group but we can change this setting by the option refgr.
solLRisk <- trajeR(Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], Risk = data_LOGIT[,22:23], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "L", hessian = TRUE, itermax = 300) solLRisk
We plot only the trajectory on a graph. We have just to use the trajeR plot function. It use the class of the object to draw the polynomial shape. The values in abscissa are the first row of the time variable.
plotrajeR(solLRisk)
We can add longitudinal data to this graph.
If we want them on the plot we have to specify Y and A in the function plot().
For more visibility, we have enlarge the 0 and 1 value on the graph and plot the data with a little and random shift, control by the index dec in the function plot.
By default colors are gray scale, but we can specify colors we want.
# colour's defintion trans <- "70" col1 <- "#034569" col1.1 <- paste0("#64AAD0", trans) col2 <- "#750062" col2.1 <- paste0("#D962C7", trans) col3 <- "#A68900" col3.1 <- paste0("#FFE773", trans) cols1 <- c(col1.1, col2.1, col3.1) cols2 <- c(col1, col2, col3) vcol <- c(cols1, cols2) plotrajeR(solLRisk, Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], dec = 5, col = vcol)
A time-dependent covariate can directly influence the shape of the trajectories. In the data, we have a time-dependent covariate composed of 0 and 1 values. We can account for its effects by using the TCOV option in the trajeR command. Any of the three methods above can be used.
solLTCOV <- trajeR( Y = data_LOGIT[,2:11], A =data_LOGIT[,12:21], TCOV = data_LOGIT[,24:33], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "L", hessian = TRUE, itermax = 300 ) solEMTCOV <- trajeR( Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], TCOV = data_LOGIT[,24:33], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "EM", hessian = TRUE, itermax = 300 ) solEMIRLSTCOV <- trajeR( Y = data_LOGIT[,2:11], A = data_LOGIT[,12:21], TCOV = data_LOGIT[,24:33], ng = 3, degre = c(0,3,4), Model = "LOGIT", Method = "EMIRLS", hessian = FALSE, itermax = 300 ) solLTCOV solEMTCOV solEMIRLSTCOV plotrajeR(solLTCOV, col = vcol)
If we want show the impact of a particular value of the time covariate in the trajectory, we can add this to the plot by plotcov option.
The fill line is the trajectory with the time covariate matrix and the dashed one show the impact on this trajectory of a particular value.
plotrajeR(solLTCOV, col = vcol, plotcov = c(0,0,0,0,0,1,1,1,1,1))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.