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,,a_{it}{n_})^t $, $W_t=(w_{i1},,w_{in_})^t $, \(\beta_k =(\beta_{k1},\cdots, \beta_{kn_\beta})\) and $k=({k1},,{kn}) $ 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}\]
Example
We use the artificial data_CNORM dataset contained in
the library. It contains 500 longitudinal observations generated with
the parameters below:
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.
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 4)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.0494 0.08383 84.08865 0
##
## 2 Intercept 19.30454 0.60769 31.76723 0
## Linear -0.09315 0.45407 -0.20514 0.83747
## Quadratic -0.45614 0.0939 -4.85767 0
## Cubic 0.02919 0.00565 5.16293 0
##
## 3 Intercept 1.6695 1.36105 1.22663 0.22002
## Linear 10.11827 1.53861 6.57626 0
## Quadratic -3.70726 0.53348 -6.94915 0
## Cubic 0.53764 0.07137 7.53355 0
## Quartic -0.02459 0.00323 -7.6119 0
## --------------------------------------------------------------------
## 1 sigma1 3.95795 0.05967 66.33372 0
## 2 sigma2 4.11085 0.06874 59.79864 0
## 3 sigma3 4.00173 0.0927 43.16709 0
## --------------------------------------------------------------------
## 1 pi1 0.45891 0.02797 0 0
## 2 pi2 0.34901 0.02185 -12.52982 0
## 3 pi3 0.19208 0.01746 -49.88924 0
## --------------------------------------------------------------------
## Likelihood : -14564.35
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 2)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.05316 0.08769 80.42978 0
##
## 2 Intercept 19.31467 0.60044 32.16748 0
## Linear -0.08935 0.44818 -0.19936 0.84199
## Quadratic -0.45753 0.09238 -4.9527 0
## Cubic 0.02927 0.00554 5.28434 0
##
## 3 Intercept 1.6695 1.35745 1.22988 0.2188
## Linear 10.11834 1.5372 6.58231 0
## Quadratic -3.70726 0.53298 -6.95578 0
## Cubic 0.53764 0.07115 7.55611 0
## Quartic -0.02459 0.00321 -7.65495 0
## --------------------------------------------------------------------
## 1 sigma1 4.02028 0.0401 100.26225 0
## 2 sigma2 4.02028 0.0401 100.26225 0
## 3 sigma3 4.02028 0.0401 100.26225 0
## --------------------------------------------------------------------
## 1 pi1 0.45978 0.02765 0 0
## 2 pi2 0.34814 0.0214 -12.9987 0
## 3 pi3 0.19207 0.01752 -49.83043 0
## --------------------------------------------------------------------
## Likelihood : -14565.75
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Expectation-maximization
## Convergence : Successful (Code 1)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.0494 0.08356 84.36189 0
##
## 2 Intercept 19.30454 0.60682 31.81271 0
## Linear -0.09315 0.45447 -0.20496 0.83761
## Quadratic -0.45614 0.09383 -4.8614 0
## Cubic 0.02919 0.00563 5.18733 0
##
## 3 Intercept 1.6695 1.34917 1.23743 0.21599
## Linear 10.11827 1.52068 6.65379 0
## Quadratic -3.70726 0.52667 -7.03905 0
## Cubic 0.53764 0.07038 7.6394 0
## Quartic -0.02459 0.00318 -7.72826 0
## --------------------------------------------------------------------
## 1 sigma1 3.95795 0.05919 66.86856 0
## 2 sigma2 4.11085 0.07048 58.32633 0
## 3 sigma3 4.00173 0.09153 43.72174 0
## --------------------------------------------------------------------
## 1 pi1 0.45891 0.02247 20.42129 0
## 2 pi2 0.34901 0.02151 16.22397 0
## 3 pi3 0.19208 0.03111 6.17446 0
## --------------------------------------------------------------------
## Likelihood : -14564.35
## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Expectation-maximization
## Convergence : Successful (Code 1)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.05316 0.08479 83.18379 0
##
## 2 Intercept 19.31467 0.5943 32.50011 0
## Linear -0.08935 0.44501 -0.20078 0.84088
## Quadratic -0.45753 0.09187 -4.98033 0
## Cubic 0.02927 0.00551 5.31318 0
##
## 3 Intercept 1.6695 1.35545 1.2317 0.21812
## Linear 10.11834 1.52775 6.62302 0
## Quadratic -3.70726 0.52912 -7.00645 0
## Cubic 0.53764 0.07071 7.60399 0
## Quartic -0.02459 0.0032 -7.69242 0
## --------------------------------------------------------------------
## 1 sigma1 4.02028 0.06143 65.44769 0
## 2 sigma2 4.02028 0.06687 60.12312 0
## 3 sigma3 4.02028 0.09258 43.42646 0
## --------------------------------------------------------------------
## 1 pi1 0.45978 0.02247 20.45779 0
## 2 pi2 0.34814 0.0215 16.19353 0
## 3 pi3 0.19207 0.0311 6.17567 0
## --------------------------------------------------------------------
## Likelihood : -14565.75
| SolL | SolEM | SolLs | SolEMs |
|---|---|---|---|
| Beta 1 | |||
| 7.04940 | 7.04940 | 7.05316 | 7.05316 |
| Beta 2 | |||
| 19.30454 | 19.30454 | 19.31467 | 19.31467 |
| -0.09315 | -0.09315 | -0.08935 | -0.08935 |
| -0.45614 | -0.45614 | -0.45753 | -0.45753 |
| 0.02919 | 0.02919 | 0.02927 | 0.02927 |
| Beta 3 | |||
| 1.66950 | 1.66950 | 1.66950 | 1.66950 |
| 10.11827 | 10.11827 | 10.11834 | 10.11834 |
| -3.70726 | -3.70726 | -3.70726 | -3.70726 |
| 0.53764 | 0.53764 | 0.53764 | 0.53764 |
| -0.02459 | -0.02459 | -0.02459 | -0.02459 |
| Sigma | |||
| 3.95795 | 3.95795 | 4.02028 | 4.02028 |
| 4.11085 | 4.11085 | 4.02028 | 4.02028 |
| 4.00173 | 4.00173 | 4.02028 | 4.02028 |
| Pi | |||
| 0.45891 | 0.45891 | 0.45978 | 0.45978 |
| 0.34901 | 0.34901 | 0.34814 | 0.34814 |
| 0.19208 | 0.19208 | 0.19207 | 0.19207 |
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.
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)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
)| SolLRisk | SolEMRisk | SolLRisks | SolEMRisks |
|---|---|---|---|
| Beta 1 | |||
| 7.04923 | 7.04923 | 7.05306 | 7.05306 |
| Beta 2 | |||
| 19.30497 | 19.30497 | 19.31517 | 19.31517 |
| -0.09451 | -0.09451 | -0.09075 | -0.09075 |
| -0.45584 | -0.45584 | -0.45723 | -0.45723 |
| 0.02917 | 0.02917 | 0.02926 | 0.02926 |
| Beta 3 | |||
| 1.66955 | 1.66955 | 1.66955 | 1.66955 |
| 10.11828 | 10.11828 | 10.11834 | 10.11834 |
| -3.70723 | -3.70722 | -3.70723 | -3.70723 |
| 0.53763 | 0.53763 | 0.53763 | 0.53763 |
| -0.02459 | -0.02459 | -0.02459 | -0.02459 |
| Sigma | |||
| 3.95767 | 3.95767 | 4.02025 | 4.02025 |
| 4.11119 | 4.11119 | 4.02025 | 4.02025 |
| 4.00162 | 4.00162 | 4.02025 | 4.02025 |
| Theta - First group 0 | |||
| -0.65570 | -0.65570 | -0.65904 | -0.65904 |
| 0.19567 | 0.19567 | 0.19531 | 0.19531 |
| 0.57547 | 0.57547 | 0.57361 | 0.57361 |
| -1.18319 | -1.18319 | -1.18428 | -1.18428 |
| 0.10679 | 0.10679 | 0.10630 | 0.10630 |
| 0.52683 | 0.52683 | 0.52553 | 0.52553 |
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 4)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 6.75767 0.20105 33.61208 0
## TCOV1 -0.00356 0.17728 -0.02009 0.98398
## TCOV2 0.58131 0.30223 1.92339 0.05449
##
## 2 Intercept 19.45638 0.66801 29.12601 0
## Linear -0.09667 0.47145 -0.20504 0.83755
## Quadratic -0.4549 0.09717 -4.68171 0
## Cubic 0.02909 0.00581 5.00416 0
## TCOV1 -0.17069 0.20426 -0.83565 0.40339
## TCOV2 -0.13292 0.33975 -0.39122 0.69565
##
## 3 Intercept 1.81446 1.41199 1.28504 0.19884
## Linear 10.11476 1.56079 6.48055 0
## Quadratic -3.70808 0.53846 -6.88646 0
## Cubic 0.53803 0.07175 7.49863 0
## Quartic -0.02462 0.00324 -7.60361 0
## TCOV1 0.0805 0.2734 0.29444 0.76844
## TCOV2 -0.35937 0.46975 -0.76502 0.4443
## --------------------------------------------------------------------
## 1 sigma1 3.95414 0.06166 64.12813 0
## 2 sigma2 4.11016 0.07184 57.21556 0
## 3 sigma3 4.00019 0.09175 43.60022 0
## --------------------------------------------------------------------
## 1 pi1 0.45879 0.02339 0 0
## 2 pi2 0.34913 0.02184 -12.50831 0
## 3 pi3 0.19208 0.01716 -50.73231 0
## --------------------------------------------------------------------
## Likelihood : -14561.48
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
| SolLTCOV2 | SolEMCOV2 | SolLTCOV2s | SolEMCOV2s |
|---|---|---|---|
| Beta 1 | |||
| 6.75767 | 6.75763 | 6.76379 | 6.76375 |
| Beta 2 | |||
| 19.45638 | 19.45636 | 19.46854 | 19.46852 |
| -0.09667 | -0.09667 | -0.09328 | -0.09328 |
| -0.45490 | -0.45490 | -0.45624 | -0.45624 |
| 0.02909 | 0.02909 | 0.02917 | 0.02917 |
| Beta 3 | |||
| 1.81446 | 1.81447 | 1.81445 | 1.81446 |
| 10.11476 | 10.11476 | 10.11485 | 10.11485 |
| -3.70808 | -3.70808 | -3.70809 | -3.70809 |
| 0.53803 | 0.53803 | 0.53803 | 0.53803 |
| -0.02462 | -0.02462 | -0.02462 | -0.02462 |
| Sigma | |||
| 3.95414 | 3.95414 | 4.01801 | 4.01801 |
| 4.11016 | 4.11016 | 4.01801 | 4.01801 |
| 4.00019 | 4.00019 | 4.01801 | 4.01801 |
| Delta 1 | |||
| -0.00356 | -0.00354 | -0.00763 | -0.00762 |
| 0.58131 | 0.58136 | 0.58090 | 0.58095 |
| Delta 2 | |||
| -0.17069 | -0.17068 | -0.16725 | -0.16724 |
| -0.13292 | -0.13289 | -0.13780 | -0.13777 |
| Delta 3 | |||
| 0.08050 | 0.08049 | 0.08045 | 0.08044 |
| -0.35937 | -0.35939 | -0.35932 | -0.35934 |
| Pi | |||
| 0.45879 | 0.45879 | 0.45970 | 0.45970 |
| 0.34913 | 0.34913 | 0.34823 | 0.34823 |
| 0.19208 | 0.19208 | 0.19207 | 0.19207 |
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 4)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.06352 0.08341 84.68581 0
##
## 2 Intercept 19.23116 0.6122 31.41334 0
## Linear -0.08916 0.45751 -0.19488 0.84549
## Quadratic -0.44929 0.09514 -4.72232 0
## Cubic 0.02841 0.00576 4.93387 0
##
## 3 Intercept 1.47075 1.35333 1.08677 0.27719
## Linear 10.30213 1.52528 6.75425 0
## Quadratic -3.7672 0.53406 -7.05396 0
## Cubic 0.54581 0.07223 7.55644 0
## Quartic -0.02499 0.0033 -7.56716 0
## --------------------------------------------------------------------
## 1 sigma1 3.92236 0.06322 62.04722 0
## 2 sigma2 4.10491 0.08035 51.0897 0
## 3 sigma3 4.00573 0.10231 39.15176 0
## --------------------------------------------------------------------
## 1 pi1 0.45829 0.02799 0 0
## 2 pi2 0.34961 0.02137 -12.66693 0
## 3 pi3 0.1921 0.01808 -48.07821 0
## --------------------------------------------------------------------
## Likelihood : -13317.83
# 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.
| Censored | Not Censored |
|---|---|
| Beta 1 | |
| 7.0635200 | 7.2426600 |
| Beta 2 | |
| 19.2311600 | 18.8462400 |
| -0.0891600 | 0.0646300 |
| -0.4492900 | -0.4836800 |
| 0.0284100 | 0.0324600 |
| Beta 3 | |
| 1.4707500 | 2.3257000 |
| 10.3021300 | 9.3122900 |
| -3.7672000 | -3.4175300 |
| 0.5458100 | 0.5012500 |
| -0.0249900 | -0.0233200 |
| Sigma | |
| 3.9223600 | 3.5984100 |
| 4.1049100 | 3.5483900 |
| 4.0057300 | 3.6011100 |
| Pi | |
| 0.4582880 | 0.4588387 |
| 0.3496098 | 0.3489012 |
| 0.1921022 | 0.1922600 |
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,
)| SolLC | SolEMC | SolLCs | SolEMCs |
|---|---|---|---|
| Beta 1 | |||
| 7.06352 | 7.06352 | 7.06261 | 7.06261 |
| Beta 2 | |||
| 19.23116 | 19.23116 | 19.23109 | 19.23109 |
| -0.08916 | -0.08916 | -0.08164 | -0.08164 |
| -0.44929 | -0.44929 | -0.45135 | -0.45135 |
| 0.02841 | 0.02841 | 0.02859 | 0.02859 |
| Beta 3 | |||
| 1.47075 | 1.47075 | 1.47257 | 1.47257 |
| 10.30213 | 10.30213 | 10.30028 | 10.30028 |
| -3.76720 | -3.76720 | -3.76648 | -3.76648 |
| 0.54581 | 0.54581 | 0.54571 | 0.54571 |
| -0.02499 | -0.02499 | -0.02498 | -0.02498 |
| Sigma | |||
| 3.92236 | 3.92236 | 3.99933 | 3.99933 |
| 4.10491 | 4.10491 | 3.99933 | 3.99933 |
| 4.00573 | 4.00573 | 3.99933 | 3.99933 |
| Pi | |||
| 0.45829 | 0.45829 | 0.45937 | 0.45937 |
| 0.34961 | 0.34961 | 0.34855 | 0.34855 |
| 0.19210 | 0.19210 | 0.19208 | 0.19208 |
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 2)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 7.06254 0.0882 80.07425 0
##
## 2 Intercept 19.23158 0.61684 31.17764 0
## Linear -0.08304 0.46168 -0.17987 0.85726
## Quadratic -0.45105 0.09566 -4.71536 0
## Cubic 0.02857 0.00577 4.95372 0
##
## 3 Intercept 1.47263 1.31226 1.12221 0.26183
## Linear 10.30026 1.4842 6.93994 0
## Quadratic -3.76644 0.521 -7.2293 0
## Cubic 0.5457 0.07058 7.73141 0
## Quartic -0.02498 0.00323 -7.73433 0
## --------------------------------------------------------------------
## 1 sigma1 3.99935 0.04462 89.63184 0
## 2 sigma2 3.99935 0.04462 89.63184 0
## 3 sigma3 3.99935 0.04462 89.63184 0
## --------------------------------------------------------------------
## 1 Baseline 0 NA NA NA
##
## 2 Intercept -0.65849 0.27284 -2.41342 0.01584
## V42 0.19917 0.36437 0.54661 0.58467
## V43 0.57275 0.35488 1.6139 0.10661
##
## 3 Intercept -1.18416 0.32452 -3.64898 0.00027
## V42 0.10819 0.43516 0.24863 0.80366
## V43 0.52529 0.4286 1.2256 0.22041
## --------------------------------------------------------------------
## Likelihood : -13317.65
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,
)| SolLRisk | SolEMRisk | SolLRisks | SolEMRisks |
|---|---|---|---|
| Beta 1 | |||
| 7.06254 | 7.06254 | 7.06254 | 7.06254 |
| Beta 2 | |||
| 19.23158 | 19.23158 | 19.23158 | 19.23158 |
| -0.08304 | -0.08304 | -0.08304 | -0.08304 |
| -0.45105 | -0.45105 | -0.45105 | -0.45105 |
| 0.02857 | 0.02857 | 0.02857 | 0.02857 |
| Beta 3 | |||
| 1.47263 | 1.47263 | 1.47263 | 1.47263 |
| 10.30026 | 10.30026 | 10.30026 | 10.30026 |
| -3.76644 | -3.76644 | -3.76644 | -3.76644 |
| 0.54570 | 0.54570 | 0.54570 | 0.54570 |
| -0.02498 | -0.02498 | -0.02498 | -0.02498 |
| Sigma | |||
| 3.99935 | 3.99935 | 3.99935 | 3.99935 |
| 3.99935 | 3.99935 | 3.99935 | 3.99935 |
| 3.99935 | 3.99935 | 3.99935 | 3.99935 |
| Theta - First group 0 | |||
| -0.65849 | -0.65849 | -0.65849 | -0.65849 |
| 0.19917 | 0.19917 | 0.19917 | 0.19917 |
| 0.57275 | 0.57275 | 0.57275 | 0.57275 |
| -1.18416 | -1.18416 | -1.18416 | -1.18416 |
| 0.10819 | 0.10819 | 0.10819 | 0.10819 |
| 0.52529 | 0.52529 | 0.52529 | 0.52529 |
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## Call TrajeR with 3 groups and a 0,3,4 degrees of polynomial shape of trajectory.
## Model : Censored Normal
## Method : Likelihood
## Convergence : Successful (Code 4)
##
## group Parameter Estimate Std. Error T for H0: Prob>|T|
## param.=0
## --------------------------------------------------------------------
## 1 Intercept 6.75767 0.20105 33.61208 0
## TCOV1 -0.00356 0.17728 -0.02009 0.98398
## TCOV2 0.58131 0.30223 1.92339 0.05449
##
## 2 Intercept 19.45638 0.66801 29.12601 0
## Linear -0.09667 0.47145 -0.20504 0.83755
## Quadratic -0.4549 0.09717 -4.68171 0
## Cubic 0.02909 0.00581 5.00416 0
## TCOV1 -0.17069 0.20426 -0.83565 0.40339
## TCOV2 -0.13292 0.33975 -0.39122 0.69565
##
## 3 Intercept 1.81446 1.41199 1.28504 0.19884
## Linear 10.11476 1.56079 6.48055 0
## Quadratic -3.70808 0.53846 -6.88646 0
## Cubic 0.53803 0.07175 7.49863 0
## Quartic -0.02462 0.00324 -7.60361 0
## TCOV1 0.0805 0.2734 0.29444 0.76844
## TCOV2 -0.35937 0.46975 -0.76502 0.4443
## --------------------------------------------------------------------
## 1 sigma1 3.95414 0.06166 64.12813 0
## 2 sigma2 4.11016 0.07184 57.21556 0
## 3 sigma3 4.00019 0.09175 43.60022 0
## --------------------------------------------------------------------
## 1 pi1 0.45879 0.02339 0 0
## 2 pi2 0.34913 0.02184 -12.50831 0
## 3 pi3 0.19208 0.01716 -50.73231 0
## --------------------------------------------------------------------
## Likelihood : -14561.48
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
| SolLTCOV2 | SolEMCOV2 | SolLTCOV2s | SolEMCOV2s |
|---|---|---|---|
| Beta 1 | |||
| 6.77065 | 6.77061 | 6.77210 | 6.77206 |
| Beta 2 | |||
| 19.48164 | 19.48160 | 19.48197 | 19.48193 |
| -0.10155 | -0.10155 | -0.09435 | -0.09435 |
| -0.44601 | -0.44601 | -0.44805 | -0.44805 |
| 0.02818 | 0.02818 | 0.02836 | 0.02836 |
| Beta 3 | |||
| 1.55721 | 1.55721 | 1.55922 | 1.55923 |
| 10.28483 | 10.28483 | 10.28281 | 10.28281 |
| -3.76279 | -3.76279 | -3.76200 | -3.76200 |
| 0.54545 | 0.54545 | 0.54534 | 0.54534 |
| -0.02498 | -0.02498 | -0.02498 | -0.02498 |
| Sigma | |||
| 3.91882 | 3.91882 | 3.99687 | 3.99687 |
| 4.10362 | 4.10362 | 3.99687 | 3.99687 |
| 4.00396 | 4.00396 | 3.99687 | 3.99687 |
| Delta 1 | |||
| 0.01436 | 0.01437 | 0.00826 | 0.00828 |
| 0.56591 | 0.56596 | 0.56734 | 0.56739 |
| Delta 2 | |||
| -0.20902 | -0.20900 | -0.20407 | -0.20406 |
| -0.26964 | -0.26960 | -0.27281 | -0.27277 |
| Delta 3 | |||
| 0.13478 | 0.13478 | 0.13465 | 0.13465 |
| -0.27838 | -0.27839 | -0.27831 | -0.27832 |
| Pi | |||
| 0.45819 | 0.45819 | 0.45931 | 0.45931 |
| 0.34971 | 0.34971 | 0.34861 | 0.34861 |
| 0.19210 | 0.19210 | 0.19208 | 0.19208 |