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,,a_{it}{n_})^t $, $W_t=(w_{i1},,w_{in_})^t $, \(\beta_k =(\beta_{k1},\cdots, \beta_{kn_\beta})\) and $k=({k1},,{kn}) $.\ 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
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Likelihood
#> Convergence : Successful (Code 4)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -1.94514 0.11179 -17.40035 0
#>
#> 2 Intercept -1.78181 1.22449 -1.45514 0.14569
#> Linear -0.426 0.64633 -0.65911 0.50986
#> Quadratic 0.1834 0.10933 1.67757 0.09349
#> Cubic -0.01211 0.00586 -2.06705 0.03878
#>
#> 3 Intercept 1.74118 0.77057 2.25959 0.02389
#> Linear 0.6527 0.78256 0.83406 0.40429
#> Quadratic -0.55852 0.25858 -2.15991 0.03083
#> Cubic 0.08247 0.03335 2.4728 0.01344
#> Quartic -0.00338 0.00147 -2.29819 0.02159
#> --------------------------------------------------------------------
#> 1 pi1 0.41814 0.03818 0 0
#> 2 pi2 0.21874 0.04173 -15.52842 0
#> 3 pi3 0.36312 0.03402 -4.14709 0
#> --------------------------------------------------------------------
#> Likelihood : -2852.746For 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
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Expectation-maximization
#> Convergence : Successful (Code 1)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -1.94766 0.09785 -19.90467 0
#>
#> 2 Intercept -1.34092 0.81565 -1.644 0.10024
#> Linear -0.61492 0.49088 -1.25268 0.21038
#> Quadratic 0.20769 0.09134 2.27389 0.02301
#> Cubic -0.01308 0.00517 -2.53057 0.01142
#>
#> 3 Intercept 1.75157 0.80598 2.17321 0.02981
#> Linear 0.7902 0.81429 0.97042 0.33189
#> Quadratic -0.61231 0.26769 -2.2874 0.02221
#> Cubic 0.08893 0.03451 2.57704 0.00999
#> Quartic -0.00362 0.00152 -2.38186 0.01726
#> --------------------------------------------------------------------
#> 1 pi1 0.41732 0.02973 14.03545 0
#> 2 pi2 0.23236 0.02801 8.29681 0
#> 3 pi3 0.35032 0.04085 8.57664 0
#> --------------------------------------------------------------------
#> Likelihood : -2852.915
solEMIRLS
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Expectation-maximization with IWRLS
#> Convergence : Successful (Code 1)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -1.94514 0.09673 -20.10856 0
#>
#> 2 Intercept -1.7818 1.01371 -1.75769 0.07886
#> Linear -0.426 0.5851 -0.72809 0.46659
#> Quadratic 0.1834 0.10483 1.74952 0.08026
#> Cubic -0.01211 0.00578 -2.09538 0.03619
#>
#> 3 Intercept 1.74118 0.76282 2.28256 0.0225
#> Linear 0.65271 0.77462 0.84262 0.39948
#> Quadratic -0.55852 0.25545 -2.18637 0.02883
#> Cubic 0.08247 0.03302 2.49736 0.01254
#> Quartic -0.00338 0.00146 -2.31562 0.02062
#> --------------------------------------------------------------------
#> 1 pi1 0.41814 0.02941 14.21838 0
#> 2 pi2 0.21874 0.02786 7.85256 0
#> 3 pi3 0.36312 0.04051 8.96429 0
#> --------------------------------------------------------------------
#> Likelihood : -2852.746We 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)
#>
#> ── Fitting the trajeR model ────────────────────────────────────────────────────
#> ℹ In progress...
#> ✔ In progress... [598ms]
#>
#> ── Result ──────────────────────────────────────────────────────────────────────
solLRisk
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Likelihood
#> Convergence : Successful (Code 4)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -1.93481 0.10745 -18.00706 0
#>
#> 2 Intercept -4.48562 4.53071 -0.99005 0.3222
#> Linear 0.88044 2.24444 0.39228 0.69487
#> Quadratic -0.0103 0.34419 -0.02993 0.97612
#> Cubic -0.0031 0.01665 -0.18607 0.8524
#>
#> 3 Intercept 1.65457 0.71553 2.31238 0.0208
#> Linear 0.49525 0.69794 0.70959 0.47799
#> Quadratic -0.492 0.23359 -2.10624 0.03523
#> Cubic 0.07471 0.03065 2.43755 0.01482
#> Quartic -0.00311 0.00137 -2.26743 0.02341
#> --------------------------------------------------------------------
#> 1 Baseline 0 NA NA NA
#>
#> 2 Intercept -0.99839 0.4274 -2.33595 0.01953
#> X1 -0.07154 0.32883 -0.21756 0.82778
#> X2 0.36804 0.34428 1.069 0.28512
#>
#> 3 Intercept 0.22836 0.2132 1.07111 0.28417
#> X1 -0.24236 0.2262 -1.07143 0.28403
#> X2 -0.38481 0.23448 -1.64113 0.10083
#> --------------------------------------------------------------------
#> Likelihood : -2849.035We 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.
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
)
#>
#> ── Fitting the trajeR model ────────────────────────────────────────────────────
#> ℹ In progress...
#> ✔ In progress... [8.6s]
#>
#> ── Result ──────────────────────────────────────────────────────────────────────
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
)
#>
#> ── Fitting the trajeR model ────────────────────────────────────────────────────
#> ℹ In progress...
#> ✔ In progress... [6.9s]
#>
#> ── Result ──────────────────────────────────────────────────────────────────────
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
)
#>
#> ── Fitting the trajeR model ────────────────────────────────────────────────────
#> ℹ In progress...
#> ✔ In progress... [1.8s]
#>
#> ── Result ──────────────────────────────────────────────────────────────────────
solLTCOV
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Likelihood
#> Convergence : Successful (Code 4)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -2.01725 0.15322 -13.16592 0
#> TCOV1 0.14188 0.17483 0.81155 0.41709
#>
#> 2 Intercept -1.84958 1.26722 -1.45956 0.14447
#> Linear -0.42661 0.66024 -0.64615 0.51822
#> Quadratic 0.18422 0.11027 1.67068 0.09485
#> Cubic -0.01217 0.00586 -2.07782 0.03778
#> TCOV1 0.11965 0.19879 0.60186 0.54729
#>
#> 3 Intercept 1.72138 0.70961 2.42582 0.01531
#> Linear 0.66141 0.73626 0.89834 0.36905
#> Quadratic -0.5608 0.24817 -2.25973 0.02388
#> Cubic 0.0827 0.03236 2.55529 0.01064
#> Quartic -0.00339 0.00144 -2.35862 0.01838
#> TCOV1 0.01868 0.11451 0.16312 0.87043
#> --------------------------------------------------------------------
#> 1 pi1 0.4187 0.04103 0 0
#> 2 pi2 0.21781 0.04439 -14.7225 0
#> 3 pi3 0.3635 0.03459 -4.0878 0
#> --------------------------------------------------------------------
#> Likelihood : -2851.937
solEMTCOV
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Expectation-maximization
#> Convergence : Successful (Code 1)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -2.01997 0.13963 -14.46618 0
#> TCOV1 0.14232 0.17354 0.82009 0.41221
#>
#> 2 Intercept -1.60997 0.90539 -1.77821 0.07543
#> Linear -0.53054 0.53194 -0.99737 0.31863
#> Quadratic 0.19795 0.09748 2.03061 0.04235
#> Cubic -0.01274 0.00546 -2.33407 0.01963
#> TCOV1 0.11604 0.1947 0.59599 0.55121
#>
#> 3 Intercept 1.72741 0.77988 2.21498 0.02681
#> Linear 0.71915 0.78989 0.91044 0.36263
#> Quadratic -0.58334 0.26019 -2.24194 0.02501
#> Cubic 0.08538 0.0336 2.54086 0.01109
#> Quartic -0.00349 0.00148 -2.34995 0.01881
#> TCOV1 0.01399 0.11608 0.12049 0.9041
#> --------------------------------------------------------------------
#> 1 pi1 0.41789 0.02964 14.09934 0
#> 2 pi2 0.22445 0.02794 8.03456 0
#> 3 pi3 0.35766 0.04073 8.78159 0
#> --------------------------------------------------------------------
#> Likelihood : -2851.971
solEMIRLSTCOV
#> Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
#> Model : Logit
#> Method : Expectation-maximization with IWRLS
#> Convergence : Successful (Code 1)
#>
#> group Parameter Estimate Std. Error T for H0: Prob>|T|
#> param.=0
#> --------------------------------------------------------------------
#> 1 Intercept -2.01725 NA NA NA
#> TCOV1 0.14188 NA NA NA
#>
#> 2 Intercept -1.84956 NA NA NA
#> Linear -0.42662 NA NA NA
#> Quadratic 0.18422 NA NA NA
#> Cubic -0.01217 NA NA NA
#> TCOV1 0.11965 NA NA NA
#>
#> 3 Intercept 1.72138 NA NA NA
#> Linear 0.66141 NA NA NA
#> Quadratic -0.5608 NA NA NA
#> Cubic 0.0827 NA NA NA
#> Quartic -0.00339 NA NA NA
#> TCOV1 0.01868 NA NA NA
#> --------------------------------------------------------------------
#> 1 pi1 0.4187 NA NA NA
#> 2 pi2 0.21781 NA NA NA
#> 3 pi3 0.3635 NA NA NA
#> --------------------------------------------------------------------
#> Likelihood : -2851.937
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.