Nothing
Normal and Censored Normal model
In trajeR: Group Based Modeling Trajectory
Censored Normal model
We suppose that the variable $Y_{it}$ is a censored variable, i.e., its values are bounded by two numbers, $y_{min}$ and $y_{max}$.
We consider a variable $Y^_{it}$ that is normally distributed 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\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}) $
Furthermore, we can link $y^{it}$ to the observed and censored variable $y{it}$ as follows:
\begin{align}
& y_{it}= y_{min} \text{ if } y^{it}< y{min}\
& y_{it}= y^{it} \text{ if } y{min}\leq y^{it}\leq y{max}\
&y_{it}= y_{max} \text{ if } y^*{it}> y{max}\
\end{align}
By setting $\mu_{ikt}=\beta_k A_{it}+\delta_k W_t$, we can write:
\begin{align}
P(Y_{it}=y_{it}|W_i=w_i,C_i=k)=\left\lbrace
\begin{array}{l}
\Phi\left(\dfrac{y_{min}-\mu_{ikt}}{\sigma_k}\right) \text{ if } y^{it}< y{min}\
\dfrac{1}{\sigma_k}\phi\left(\dfrac{y_{it}-\mu_{ikt}}{\sigma_k}\right)\text{ if } y_{min}\leq y^{it}\leq y{max}\
1-\Phi\left(\dfrac{y_{max}-\mu_{ikt}}{\sigma_k}\right) \text{ if } y^*{it}> y{max}
\end{array}\right.
\end{align}
Parameters
Example
We use the artificial data_CNORM dataset contained in the library. It contains 500 longitudinal observations generated with the parameters below:
- 3 groups;
- probabilities are $\pi_1=0.32$, $\pi_2=0.54$, and $\pi_3=0.14$;
- period is 10;
- we use 3 polynomial shapes:
- degree 4 and $\beta_1=(2.797,8.809,-3.201,0.463,-0.021)$;
- degree 0 and $\beta_2=7$;
- degree 3 and $\beta_{3}=(19.545,-0.297,-0.407,0.026)$.
- $\sigma$ is the same for the 3 groups, $\sigma= 4$.
The data_CNORM dataset contains the time-dependent variable $Y_i$ in data_CNORM[,2:11], the time variable in data_CNORM[,12:21], a time-dependent covariate in data_CNORM[,22:41], and a covariate that influences the group membership probability in data_CNORM[,42:43].
data_CNORM[,2:11] is a matrix of real values.data_CNORM[,12:21] is a matrix containing time points from 1 to 10.data_CNORM[,22:41] is a time-dependent covariate that influences the shape of the trajectories. It is a matrix with 0 and 1 values, representing the presence or absence of a characteristic for the individual.data_CNORM[,42:43] is a matrix of real values.
library(trajeR)
data("data_CNORM")
matplot(
t(data_CNORM[, 12:21]),
t(data_CNORM[, 2:11]),
pch = 1,
type = 'b',
col = "black",
lty = 1,
xlab = "Times",
ylab = "Values",
main = "Plot of the individual's trajectories"
)
We use each method to fit the model. For all methods, we specify the number of group of our model, ng=3, the degree of the polynomial shape of the trajectory. Here we choice a line parallel to abscissa axis, a cubic polynomial and a quadric polynomial. So degree is a vector $(0,3,4)$.
We specify hessian=TRUE to ask the calculus of the hessian matrix.
Likelihood method
For the Likelihood method we call trajeR with option Method ="L". We specify if we want the same sigma in each group with the parameters ssigma. If it is TRUE then the algorithm search the same sigma in all group. To ensure convergence, given that the time matrix contains large values such as 10, we set stepmax to 1e-2, consistent with the ucminf package.
# Likelihood different sigma
solL <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
hessian = TRUE,
ssigma = FALSE,
control = list(stepmax = 1e-2)
)
solL
If we want force the use unique sigma for each groups we write ssigma = TRUE.
# Likelihood same sigma
solLs <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
control = list(stepmax = 1e-2)
)
solLs
EM method
For use EM method we write the same syntax but change the variable Method to EM.
# EM
solEM <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = FALSE,
hessian = TRUE
)
solEMs <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE
)
solEM
solEMs
t1 <- solL$tab[, 1]
t1[14:16] <- exp(t1[14:16]) / sum(exp(t1[14:16]))
t2 <- solLs$tab[, 1]
t2[14:16] <- exp(t2[14:16]) / sum(exp(t2[14:16]))
tab <- round(cbind(t1, solEM$tab[, 1], t2, solEMs$tab[, 1]), 5)
colnames(tab) <- c("SolL", "SolEM", "SolLs", "SolEMs")
library(kableExtra)
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Pi", 14, 16)
Plot
We plot only the trajectories on a graph. We have just to use the native 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(solL)
We can add longitudinal data_CNORM to this graph. If we want this on the plot we have to specify Y| and A in the function plotrajeR().
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(solEM, Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], col = vcol)
Adding risk covariate
We can add covariate in the calculus of the parameters with or without option same sigma.
If we add risk variable we use option Risk =.
solLRisk <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
Risk = data_CNORM[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = FALSE,
hessian = TRUE,
control = list(stepmax = 1e-2)
)
solLRisks <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
Risk = data_CNORM[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
control = list(stepmax = 1e-2)
)
solEMRisk <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
Risk = data_CNORM[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = FALSE,
hessian = TRUE
)
solEMRisks <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
Risk = data_CNORM[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE
)
tab <- cbind(
solLRisk$tab[, 1],
solEMRisk$tab[, 1],
solLRisks$tab[, 1],
solEMRisks$tab[, 1]
)
for (i in 1:4) {
tab[14:22, i] <- tab[14:22, i] - tab[14:16, i]
}
tab <- tab[-c(14:16), ]
tab <- round(tab, 5)
colnames(tab) <- c("SolLRisk", "SolEMRisk", "SolLRisks", "SolEMRisks")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Theta - First group 0", 14, 19)
Adding time covariate
If we add time dependent covariate we use option TCOV.
solLTCOV2 <- trajeR(
Y = data_CNORM[, 2:11],
A = data_CNORM[, 12:21],
TCOV = data_CNORM[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = FALSE,
hessian = TRUE,
control = list(stepmax = 1e-2)
)
solLTCOV2
solLTCOV2s <- trajeR(
data_CNORM[, 2:11],
data_CNORM[, 12:21],
TCOV = data_CNORM[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
control = list(stepmax = 1e-2)
)
solEMTCOV2 <- trajeR(
data_CNORM[, 2:11],
data_CNORM[, 12:21],
TCOV = data_CNORM[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = FALSE,
hessian = TRUE
)
solEMTCOV2s <- trajeR(
data_CNORM[, 2:11],
data_CNORM[, 12:21],
TCOV = data_CNORM[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE
)
Here we can use all method and with or without same sigma.
We have
t1 <- solLTCOV2$tab[, 1]
t1[20:22] <- exp(t1[20:22]) / sum(exp(t1[20:22]))
t2 <- solLTCOV2s$tab[, 1]
t2[20:22] <- exp(t2[20:22]) / sum(exp(t2[20:22]))
tab <- round(cbind(t1, solEMTCOV2$tab[, 1], t2, solEMTCOV2s$tab[, 1]), 5)
colnames(tab) <- c("SolLTCOV2", "SolEMCOV2", "SolLTCOV2s", "SolEMCOV2s")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Delta 1", 14, 15) %>%
pack_rows("Delta 2", 16, 17) %>%
pack_rows("Delta 3", 18, 19) %>%
pack_rows("Pi", 20, 22)
Censored Normal distribution
To illustrate this party we use artificial data contained in this package CNORM_data_Censored.
We have modified the previous data in such way that they become censored data. All values upper to 23 become 23 and those smaller to 2 become 2.
So we obtained data that follow a censored normal distribution. If we want take into consideration this fact,
library(trajeR)
data("data_CNORM_Censored")
matplot(
t(data_CNORM_Censored[, 12:21]),
t(data_CNORM_Censored[, 2:11]),
pch = 1,
type = 'b',
col = "black",
lty = 1,
xlab = "Times",
ylab = "Values",
main = "Plot of the individual's trajectories"
)
We use each method to fit the model. For all methods, we specify the number of group of our model, ng=3, the degree of the polynomial shape of the trajectory. Here we choice a line parallel to abscissa axis, a cubic polynomial and a quadric polynomial. So degree is a vector $(0,3,4)$.
We specify hessian=TRUE to ask the calculus of the hessian matrix.
For the Likelihood method we call trajeR with option Method ="L". We specify if we want the same sigma in each group with the parameters ssigma. If it is TRUE then the algorithm search the same sigma in all group.
# Likelihood different sigma
solLC <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
hessian = TRUE,
ssigma = FALSE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solLC
# 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(
solLC,
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
col = vcol
)
If we search the solution without considering the fact that the data are not censored, we find parameters that are wrong. This amounts to considering the maximum and minimum are really maximum and minimum while the real values are omitted. To illustrate this fact we compute the parameters without censored.
solLnC <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
hessian = TRUE,
ssigma = FALSE,
control = list(stepmax = 1e-2)
)
tab <- round(cbind(solLC$tab[, 1], solLnC$tab[, 1]), 5)
tab[14:16, 1] <- exp(tab[14:16, 1]) / sum(exp(tab[14:16, 1]))
tab[14:16, 2] <- exp(tab[14:16, 2]) / sum(exp(tab[14:16, 2]))
colnames(tab) <- c("Censored", "Not Censored")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Pi", 14, 16)
par(mfrow = c(1, 2))
plotrajeR(
solLC,
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
col = vcol,
main = "Censored"
)
plotrajeR(
solLnC,
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
col = vcol,
main = "Not Censored"
)
For all the method we have
solLCs <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solEMC <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = FALSE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
solEMCs <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
tab <- cbind(solLC$tab[, 1], solEMC$tab[, 1], solLCs$tab[, 1], solEMCs$tab[, 1])
tab[14:16, 1] <- exp(tab[14:16, 1]) / sum(exp(tab[14:16, 1]))
tab[14:16, 3] <- exp(tab[14:16, 3]) / sum(exp(tab[14:16, 3]))
tab <- round(tab, 5)
colnames(tab) <- c("SolLC", "SolEMC", "SolLCs", "SolEMCs")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Pi", 14, 16)
We can add time covariate to the censored data.
solLCRisk <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
Risk = data_CNORM_Censored[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solLCRisk
solLCRisks <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
Risk = data_CNORM_Censored[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solEMCRisk <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
Risk = data_CNORM_Censored[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
solEMCRisks <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
Risk = data_CNORM_Censored[, 42:43],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
library(kableExtra)
tab <- cbind(
solLCRisk$tab[, 1],
solEMCRisk$tab[, 1],
solLCRisks$tab[, 1],
solEMCRisks$tab[, 1]
)
for (i in 1:4) {
tab[14:22, i] <- tab[14:22, i] - tab[14:16, i]
}
tab <- tab[-c(14:16), ]
tab <- round(tab, 5)
colnames(tab) <- c("SolLRisk", "SolEMRisk", "SolLRisks", "SolEMRisks")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Theta - First group 0", 14, 19)
If we add time dependent covariate we use option TCOV.
solLCTCOV2 <- trajeR(
Y = data_CNORM_Censored[, 2:11],
A = data_CNORM_Censored[, 12:21],
TCOV = data_CNORM_Censored[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = FALSE,
hessian = TRUE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solLTCOV2
Here we can use all method and with or without same sigma.
solLCTCOV2s <- trajeR(
data_CNORM_Censored[, 2:11],
data_CNORM_Censored[, 12:21],
TCOV = data_CNORM_Censored[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "L",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
control = list(stepmax = 1e-2)
)
solEMCTCOV2 <- trajeR(
data_CNORM_Censored[, 2:11],
data_CNORM_Censored[, 12:21],
TCOV = data_CNORM_Censored[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = FALSE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
solEMCTCOV2s <- trajeR(
data_CNORM_Censored[, 2:11],
data_CNORM_Censored[, 12:21],
TCOV = data_CNORM_Censored[, 22:41],
ng = 3,
degre = c(0, 3, 4),
Model = "CNORM",
Method = "EM",
ssigma = TRUE,
hessian = TRUE,
ymin = 2,
ymax = 23,
)
We have
t1 <- solLCTCOV2$tab[, 1]
t1[20:22] <- exp(t1[20:22]) / sum(exp(t1[20:22]))
t2 <- solLCTCOV2s$tab[, 1]
t2[20:22] <- exp(t2[20:22]) / sum(exp(t2[20:22]))
tab <- round(cbind(t1, solEMCTCOV2$tab[, 1], t2, solEMCTCOV2s$tab[, 1]), 5)
colnames(tab) <- c("SolLTCOV2", "SolEMCOV2", "SolLTCOV2s", "SolEMCOV2s")
kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>%
kable_styling() %>%
add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>%
pack_rows("Beta 1", 1, 1) %>%
pack_rows("Beta 2", 2, 5) %>%
pack_rows("Beta 3", 6, 10) %>%
pack_rows("Sigma", 11, 13) %>%
pack_rows("Delta 1", 14, 15) %>%
pack_rows("Delta 2", 16, 17) %>%
pack_rows("Delta 3", 18, 19) %>%
pack_rows("Pi", 20, 22)
Try the trajeR package in your browser
Any scripts or data that you put into this service are public.
trajeR documentation built on Aug. 4, 2026, 1:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Censored Normal model
We suppose that the variable $Y_{it}$ is a censored variable, i.e., its values are bounded by two numbers, $y_{min}$ and $y_{max}$. We consider a variable $Y^_{it}$ that is normally distributed 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\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}) $ Furthermore, we can link $y^{it}$ to the observed and censored variable $y{it}$ as follows: \begin{align} & y_{it}= y_{min} \text{ if } y^{it}< y{min}\ & y_{it}= y^{it} \text{ if } y{min}\leq y^{it}\leq y{max}\ &y_{it}= y_{max} \text{ if } y^*{it}> y{max}\ \end{align}
By setting $\mu_{ikt}=\beta_k A_{it}+\delta_k W_t$, we can write: \begin{align} P(Y_{it}=y_{it}|W_i=w_i,C_i=k)=\left\lbrace \begin{array}{l} \Phi\left(\dfrac{y_{min}-\mu_{ikt}}{\sigma_k}\right) \text{ if } y^{it}< y{min}\ \dfrac{1}{\sigma_k}\phi\left(\dfrac{y_{it}-\mu_{ikt}}{\sigma_k}\right)\text{ if } y_{min}\leq y^{it}\leq y{max}\ 1-\Phi\left(\dfrac{y_{max}-\mu_{ikt}}{\sigma_k}\right) \text{ if } y^*{it}> y{max} \end{array}\right. \end{align}
Parameters
Example
We use the artificial data_CNORM dataset contained in the library. It contains 500 longitudinal observations generated with the parameters below:
- 3 groups;
- probabilities are $\pi_1=0.32$, $\pi_2=0.54$, and $\pi_3=0.14$;
- period is 10;
- we use 3 polynomial shapes:
- degree 4 and $\beta_1=(2.797,8.809,-3.201,0.463,-0.021)$;
- degree 0 and $\beta_2=7$;
- degree 3 and $\beta_{3}=(19.545,-0.297,-0.407,0.026)$.
- $\sigma$ is the same for the 3 groups, $\sigma= 4$.
The data_CNORM dataset contains the time-dependent variable $Y_i$ in data_CNORM[,2:11], the time variable in data_CNORM[,12:21], a time-dependent covariate in data_CNORM[,22:41], and a covariate that influences the group membership probability in data_CNORM[,42:43].
data_CNORM[,2:11]is a matrix of real values.data_CNORM[,12:21]is a matrix containing time points from 1 to 10.data_CNORM[,22:41]is a time-dependent covariate that influences the shape of the trajectories. It is a matrix with 0 and 1 values, representing the presence or absence of a characteristic for the individual.data_CNORM[,42:43]is a matrix of real values.
library(trajeR) data("data_CNORM") matplot( t(data_CNORM[, 12:21]), t(data_CNORM[, 2:11]), pch = 1, type = 'b', col = "black", lty = 1, xlab = "Times", ylab = "Values", main = "Plot of the individual's trajectories" )
We use each method to fit the model. For all methods, we specify the number of group of our model, ng=3, the degree of the polynomial shape of the trajectory. Here we choice a line parallel to abscissa axis, a cubic polynomial and a quadric polynomial. So degree is a vector $(0,3,4)$.
We specify hessian=TRUE to ask the calculus of the hessian matrix.
Likelihood method
For the Likelihood method we call trajeR with option Method ="L". We specify if we want the same sigma in each group with the parameters ssigma. If it is TRUE then the algorithm search the same sigma in all group. To ensure convergence, given that the time matrix contains large values such as 10, we set stepmax to 1e-2, consistent with the ucminf package.
# Likelihood different sigma solL <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", hessian = TRUE, ssigma = FALSE, control = list(stepmax = 1e-2) ) solL
If we want force the use unique sigma for each groups we write ssigma = TRUE.
# Likelihood same sigma solLs <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, control = list(stepmax = 1e-2) ) solLs
EM method
For use EM method we write the same syntax but change the variable Method to EM.
# EM solEM <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = FALSE, hessian = TRUE ) solEMs <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE ) solEM solEMs
t1 <- solL$tab[, 1] t1[14:16] <- exp(t1[14:16]) / sum(exp(t1[14:16])) t2 <- solLs$tab[, 1] t2[14:16] <- exp(t2[14:16]) / sum(exp(t2[14:16])) tab <- round(cbind(t1, solEM$tab[, 1], t2, solEMs$tab[, 1]), 5) colnames(tab) <- c("SolL", "SolEM", "SolLs", "SolEMs") library(kableExtra) kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Pi", 14, 16)
Plot
We plot only the trajectories on a graph. We have just to use the native 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(solL)
We can add longitudinal data_CNORM to this graph. If we want this on the plot we have to specify Y| and A in the function plotrajeR().
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(solEM, Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], col = vcol)
Adding risk covariate
We can add covariate in the calculus of the parameters with or without option same sigma.
If we add risk variable we use option Risk =.
solLRisk <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], Risk = data_CNORM[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = FALSE, hessian = TRUE, control = list(stepmax = 1e-2) ) solLRisks <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], Risk = data_CNORM[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, control = list(stepmax = 1e-2) ) solEMRisk <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], Risk = data_CNORM[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = FALSE, hessian = TRUE ) solEMRisks <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], Risk = data_CNORM[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE )
tab <- cbind( solLRisk$tab[, 1], solEMRisk$tab[, 1], solLRisks$tab[, 1], solEMRisks$tab[, 1] ) for (i in 1:4) { tab[14:22, i] <- tab[14:22, i] - tab[14:16, i] } tab <- tab[-c(14:16), ] tab <- round(tab, 5) colnames(tab) <- c("SolLRisk", "SolEMRisk", "SolLRisks", "SolEMRisks") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Theta - First group 0", 14, 19)
Adding time covariate
If we add time dependent covariate we use option TCOV.
solLTCOV2 <- trajeR( Y = data_CNORM[, 2:11], A = data_CNORM[, 12:21], TCOV = data_CNORM[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = FALSE, hessian = TRUE, control = list(stepmax = 1e-2) ) solLTCOV2
solLTCOV2s <- trajeR( data_CNORM[, 2:11], data_CNORM[, 12:21], TCOV = data_CNORM[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, control = list(stepmax = 1e-2) ) solEMTCOV2 <- trajeR( data_CNORM[, 2:11], data_CNORM[, 12:21], TCOV = data_CNORM[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = FALSE, hessian = TRUE ) solEMTCOV2s <- trajeR( data_CNORM[, 2:11], data_CNORM[, 12:21], TCOV = data_CNORM[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE )
Here we can use all method and with or without same sigma.
We have
t1 <- solLTCOV2$tab[, 1] t1[20:22] <- exp(t1[20:22]) / sum(exp(t1[20:22])) t2 <- solLTCOV2s$tab[, 1] t2[20:22] <- exp(t2[20:22]) / sum(exp(t2[20:22])) tab <- round(cbind(t1, solEMTCOV2$tab[, 1], t2, solEMTCOV2s$tab[, 1]), 5) colnames(tab) <- c("SolLTCOV2", "SolEMCOV2", "SolLTCOV2s", "SolEMCOV2s") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Delta 1", 14, 15) %>% pack_rows("Delta 2", 16, 17) %>% pack_rows("Delta 3", 18, 19) %>% pack_rows("Pi", 20, 22)
Censored Normal distribution
To illustrate this party we use artificial data contained in this package CNORM_data_Censored.
We have modified the previous data in such way that they become censored data. All values upper to 23 become 23 and those smaller to 2 become 2. So we obtained data that follow a censored normal distribution. If we want take into consideration this fact,
library(trajeR) data("data_CNORM_Censored") matplot( t(data_CNORM_Censored[, 12:21]), t(data_CNORM_Censored[, 2:11]), pch = 1, type = 'b', col = "black", lty = 1, xlab = "Times", ylab = "Values", main = "Plot of the individual's trajectories" )
We use each method to fit the model. For all methods, we specify the number of group of our model, ng=3, the degree of the polynomial shape of the trajectory. Here we choice a line parallel to abscissa axis, a cubic polynomial and a quadric polynomial. So degree is a vector $(0,3,4)$.
We specify hessian=TRUE to ask the calculus of the hessian matrix.
For the Likelihood method we call trajeR with option Method ="L". We specify if we want the same sigma in each group with the parameters ssigma. If it is TRUE then the algorithm search the same sigma in all group.
# Likelihood different sigma solLC <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", hessian = TRUE, ssigma = FALSE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solLC # 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( solLC, Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], col = vcol )
If we search the solution without considering the fact that the data are not censored, we find parameters that are wrong. This amounts to considering the maximum and minimum are really maximum and minimum while the real values are omitted. To illustrate this fact we compute the parameters without censored.
solLnC <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", hessian = TRUE, ssigma = FALSE, control = list(stepmax = 1e-2) )
tab <- round(cbind(solLC$tab[, 1], solLnC$tab[, 1]), 5) tab[14:16, 1] <- exp(tab[14:16, 1]) / sum(exp(tab[14:16, 1])) tab[14:16, 2] <- exp(tab[14:16, 2]) / sum(exp(tab[14:16, 2])) colnames(tab) <- c("Censored", "Not Censored") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Pi", 14, 16) par(mfrow = c(1, 2)) plotrajeR( solLC, Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], col = vcol, main = "Censored" ) plotrajeR( solLnC, Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], col = vcol, main = "Not Censored" )
For all the method we have
solLCs <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solEMC <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = FALSE, hessian = TRUE, ymin = 2, ymax = 23, ) solEMCs <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, )
tab <- cbind(solLC$tab[, 1], solEMC$tab[, 1], solLCs$tab[, 1], solEMCs$tab[, 1]) tab[14:16, 1] <- exp(tab[14:16, 1]) / sum(exp(tab[14:16, 1])) tab[14:16, 3] <- exp(tab[14:16, 3]) / sum(exp(tab[14:16, 3])) tab <- round(tab, 5) colnames(tab) <- c("SolLC", "SolEMC", "SolLCs", "SolEMCs") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Pi", 14, 16)
We can add time covariate to the censored data.
solLCRisk <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], Risk = data_CNORM_Censored[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solLCRisk
solLCRisks <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], Risk = data_CNORM_Censored[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solEMCRisk <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], Risk = data_CNORM_Censored[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, ) solEMCRisks <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], Risk = data_CNORM_Censored[, 42:43], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, )
library(kableExtra) tab <- cbind( solLCRisk$tab[, 1], solEMCRisk$tab[, 1], solLCRisks$tab[, 1], solEMCRisks$tab[, 1] ) for (i in 1:4) { tab[14:22, i] <- tab[14:22, i] - tab[14:16, i] } tab <- tab[-c(14:16), ] tab <- round(tab, 5) colnames(tab) <- c("SolLRisk", "SolEMRisk", "SolLRisks", "SolEMRisks") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Theta - First group 0", 14, 19)
If we add time dependent covariate we use option TCOV.
solLCTCOV2 <- trajeR( Y = data_CNORM_Censored[, 2:11], A = data_CNORM_Censored[, 12:21], TCOV = data_CNORM_Censored[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = FALSE, hessian = TRUE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solLTCOV2
Here we can use all method and with or without same sigma.
solLCTCOV2s <- trajeR( data_CNORM_Censored[, 2:11], data_CNORM_Censored[, 12:21], TCOV = data_CNORM_Censored[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "L", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, control = list(stepmax = 1e-2) ) solEMCTCOV2 <- trajeR( data_CNORM_Censored[, 2:11], data_CNORM_Censored[, 12:21], TCOV = data_CNORM_Censored[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = FALSE, hessian = TRUE, ymin = 2, ymax = 23, ) solEMCTCOV2s <- trajeR( data_CNORM_Censored[, 2:11], data_CNORM_Censored[, 12:21], TCOV = data_CNORM_Censored[, 22:41], ng = 3, degre = c(0, 3, 4), Model = "CNORM", Method = "EM", ssigma = TRUE, hessian = TRUE, ymin = 2, ymax = 23, )
We have
t1 <- solLCTCOV2$tab[, 1] t1[20:22] <- exp(t1[20:22]) / sum(exp(t1[20:22])) t2 <- solLCTCOV2s$tab[, 1] t2[20:22] <- exp(t2[20:22]) / sum(exp(t2[20:22])) tab <- round(cbind(t1, solEMCTCOV2$tab[, 1], t2, solEMCTCOV2s$tab[, 1]), 5) colnames(tab) <- c("SolLTCOV2", "SolEMCOV2", "SolLTCOV2s", "SolEMCOV2s") kable(tab, "html", booktabs = T, escape = FALSE, align = "r") %>% kable_styling() %>% add_header_above(c("Different sigma" = 2, "Same sigma" = 2)) %>% pack_rows("Beta 1", 1, 1) %>% pack_rows("Beta 2", 2, 5) %>% pack_rows("Beta 3", 6, 10) %>% pack_rows("Sigma", 11, 13) %>% pack_rows("Delta 1", 14, 15) %>% pack_rows("Delta 2", 16, 17) %>% pack_rows("Delta 3", 18, 19) %>% pack_rows("Pi", 20, 22)
Try the trajeR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.