---
title: "Group Based Trajectory Modelling (GBTM)"
author: "Cédric NOEL - Jang SHILTZ"
email: "cedric.noel@univ-lorraine.fr"
date: "`r Sys.Date()`"
output: html_vignette
bibliography:
  - biblio.bib
vignette: >
  %\VignetteIndexEntry{Group Based Trajectory Modelling (GBTM)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

<div style="text-align: center;">
![Logo](figure/logotrajeR.png){width=20%}
</div>

# Introduction

Longitudinal data from repeated measures are often encountered in practice. For each subject, we measure a variable at several different time points. For example, in finance, we can measure the salary of employees over time [@schiltz2015]; in criminology, we can measure physical aggression over time [@10.2307/j.ctvjf9z1f]; and in medicine, we can measure electroencephalography data over time [@elmer_novel_2019]. The family of Latent Growth Models (LGM) studies inter-individual (between-subject) and intra-individual (within-subject) patterns of change over time. This pattern can be represented by time trends or latent trajectories and is influenced by one or more covariates.
The complexity of these patterns can be modeled by adding latent classes. The trajectories represent the mean development of an outcome over time for each class.

 We introduce $K$ latent classes, and the assignment of an individual to a class is based on the degree of similarity of their developmental course. These models assume no between-subject variability within a class, and the error variance is assumed to be the same in each class.
 
 Consider a time-varying variable of interest $Y_{i}$ for an individual $i$ in a population of size $N$. Let $Y_i= (y_{i_{1}},\cdots ,y_{i_{T}})$ be $T$ measures of the variable $Y$, taken at times $(t_{1},\cdots ,t_{T})$ for subject $i$ from a sample of size $n$. \\
Furthermore, we assume conditional independence for the sequential realizations of the elements $y_{i_{t}}$ over the $T$ measurement periods, i.e., a value $y_{i_{t}}$ does not depend on past values $y_{i_{t'}}$ where $t'<t$.
 
We now suppose that the population is divided into $K$ homogeneous sub-populations $C_{1},...,C_{K}$. Let $P^k(Y_i)$ be the probability of $Y_i$ given membership in group $C_k$ and
$P(Y_i)$ the unconditional probability of observing the realization $Y_i$ of $Y$. Then,
\begin{equation}
P(Y_i)=\sum_{k=1}^{K}P\left(\omega_i\in C_k\right)P^k(Y_i).
\end{equation}


The density $f$ of $Y$ is then described in the form of a mixture model [@10.2307/j.ctvjf9z1f] by:
\begin{equation}
f(y_i;\psi)=\sum_{k=1}^{K}\pi_kg_k(y_i;\theta_k).
\end{equation}
and the individuals follow a time-dependent trajectory. The different time values are stored in a vector $A_i$.

Here, the group probabilities $\pi_k>0$ denote the probability of a given subject belonging to group $k$, and they satisfy $\sum_{k=1}^{K}\pi_k=1$.
The parameters $\theta_{k}$ describe the shape of the trajectories in group $k$; hence, the model depends on the parameters $\psi=(K,\pi_1,\cdots,\pi_{K-1}, \theta_1,\cdots, \theta_K)$.

If we further suppose that the trajectories of $Y$ are influenced by a static set of risk variables $X=\left(X_{1}\cdots X_{R}\right)$, as well as by a time-dependent covariate $W$ which is independent of $X$, we can write the conditional density of $Y$ given $X$ and $W$.


\begin{equation}
f\left(y_i|x_i,w_i\right)=\sum_{k=1}^{K}\left(P\left(\omega_i\in C_k|X_i=x_i\right)\prod_{t=1}^{T}P\left(Y_{i_t}=y_{i_t}|X_i=x_i,W_i=w_i,\omega_i\in C_k \right)\right)
\end{equation}


The aim of this package is to provide tools to work with this situation with different densities for $Y$:

+ Logit
+ Censored Normal
+ Zero Inflated Poisson
+ Beta
+ Non Linear

For the first three models, we assume that the trajectories follow a polynomial shape. For the last one, we can specify a particular function for the shape of the trajectories, such as an exponential function.

This package extends the \verb|traj| software for SAS or Stata written by @jones_sas_2001.

A complete description of the model and the proof of each method used inside this package can be found in Noel's thesis, see https://orbilu.uni.lu/.

# Getting started

First, we need to load `trajeR` package.
```{r setup}
library(trajeR)
```

The main function is named `trajeR`. It fit the model and find parameters to a given degree of polynomial shape. It syntax is
```{r, eval = FALSE}
trajeR(
  Y,
  A,
  Risk = NULL,
  TCOV = NULL,
  degre,
  degre.nu = 0,
  Model,
  Method = "L",
  ssigma = FALSE,
  ymax = max(Y) + 1,
  ymin = min(Y) - 1,
  hessian = TRUE,
  itermax = 100,
  paraminit = NULL,
  ProbIRLS = TRUE,
  refgr = 1,
  fct = NULL,
  diffct = NULL,
  nbvar = NULL,
  nls.lmiter = 50
)
```
The function `trajeR` accepts R syntax using the `~` operator to define the model.

**Arguments**

+ `Y`: Matrix. A matrix containing the variables in the model.
+ `A`: Matrix. A matrix containing the time variable data.
+ `Risk`: Matrix. An optional matrix that modify the probability of belong to group. By default its value is a matrix with one column with value 1.
+ `TCOV`: Matrix. An optional matrix containing the time covariate that influence the trajectory themselves. By default its value is NULL.
+ `degre`: Vector of integer. The degree of every polynomial function.
+ `degre.nu`: Vector of integer. The degree of all Poisson part for a ZIP model.
+ `Model`: String. The model used. The value are LOGIT for a Logit Mixture model, CNORM for a Censored Normal Mixture Model or ZIP for Zero Inflated Poisson Mixture model.
+ `Method`: String. Determine the method used for find the parameters of the model. The value are L for the Maximum Likelihood Estimation, EM for Expectation Maximization method with quasi newton method inside, EMIWRLS for Expectation Maximization method with Iterative Weighted Least Square.
+ `ssigma`: Logical. By default its value is FALSE. For the CNORM model, indicate if we want the same sigma for all normal density function.
+ `ymax`: Real. For the CNORM model, indicate the maximum value of the data. It concern only the model with censored data. By default its value is the maximum value of the data plus 1.
+ `ymin`: Real. For the CNORM model, indicate the minimum value of the data. It concern only the model with censored data. By default its value is the maximum value of the data minus 1.
+ `hessian`: Logical. Indicate if we want calculate the hessian matrix. Default is FALSE. If the method use is Likelihood, the hessian is calculated by inverting the Information's Fisher Matrix. To avoid numerically singular matrix we find the pseudo inverse matrix by using the ginv function int he package MASS. If the method is EM or EMIWRLS, the hessian is calculated by using Louis method.
+ `itermax`: Integer. Indicate the maximal number of iteration for optim function or for the EM algorithm.
+ `paraminit`: Vector. The vector of initial parameters. By default trajeR calculate the initial value based of the range or the standard deviation.
+ `ProbIRLS`: Logical. Indicate the method to sue in the search of predictor's probability. If TRUE (by default) we use IRLS method and if FALSE we use optimization method.
+ `refgr`: Integer. The number of reference group. By default is 1.
+ `fct`: Function. The definition of the function f in the definition in nonlinear model.
+ `diffct`: Function. The differential of the function f in the nonlinear model.
+ `nbvar`: Integer. The number of variable in the nonlinear model.
+ `nls.lmiter`: Integer. In the case of non linear model, the maximum number of iterations allowed.

**Details**

Models for `trajeR` is, by default, a polynomial regression of the time value parameters for each groups. The number fo group is controlled by the length of the vector beta. We can specify the degree of the polynomial shape for each groups by the vector degre.

**Value**

The function `trajeR` returns an S3 object whose class depends on the specified Model: `Trajectory.LOGIT`, `Trajectory.CNORM`, `Trajectory.ZIP`, `Trajectory.NL`, or `Trajectory.BETA`.


All returned objects are lists containing the following common components:

+ `beta`: Vector of estimated trajectory parameters.
+ `delta`: Vector of parameters for time-dependent covariates (if `TCOV` is provided, otherwise NULL or NA).
+ `theta`: Vector of estimated parameters for group membership probabilities.
+ `sd`: Vector of standard errors for the estimated parameters (if `hessian = TRUE`).
+ `tab`: Data frame containing parameter estimates, standard errors, T-values, and p-values.
+ `Model`: Character string indicating the model type (`LOGIT`, `CNORM`, `ZIP`, `POIS`, `NL`, or `BETA`).
+ `groups`: Integer specifying the number of latent groups ($ng$).
+ `Names`: Character vector containing parameter names.
+ `Method`: Character string indicating the estimation method used (L, EM, or EMIRLS).
+ `Size`: Integer specifying the sample size (number of individuals $n$).
+ `Likelihood`: Real number representing the final log-likelihood value.
+ `Time`: Vector containing measurement time points.
+ `degre`: Vector of polynomial degrees for trajectory shapes.
+ `varcov`: Estimated variance-covariance matrix of the parameters.
+ `convergence`: Integer indicating optimizer convergence status ($1$ for successful convergence). See `ucmif`package.

Depending on the class, the returned object also contains model-specific components:

+ `Trajectory.CNORM`:
  - `sigma`: Vector of residual standard deviations for each group.
  - `min, max`: Real numbers indicating the lower (ymin) and upper (ymax) censoring thresholds.
+ `Trajectory.ZIP`:
  - `nu`: Vector of parameters for the zero-inflation part.
  - `degre.nu`: Vector of polynomial degrees for the zero-inflation part.
  - `period`: Integer specifying the number of measurement time points.
+ `Trajectory.BETA`:
  - `phi`: Vector of precision parameters.
  - `degre.phi`: Vector of polynomial degrees for the precision part.
  - `invH`: Inverse Hessian matrix (if computed).
+ `Trajectory.NL`:
  - `sigma`: Vector of residual standard deviations for each group.
  - `fct`: The user-defined non-linear function used for modeling trajectories.
  

# Références
<div id="refs"></div>