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

# 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.

```{r}
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. 


```{r, message = FALSE}
# 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`.


```{r, message = FALSE}
# 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`.


```{r, message = FALSE}
# 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
```


```{r, echo = FALSE}
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.

```{r}
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.


```{r}
# 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 =`.


```{r, message = FALSE}
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
)
```


```{r, echo = FALSE}
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`.

```{r, message = FALSE}
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
```

```{r, message = FALSE}
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

```{r, echo = FALSE}
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,

```{r}
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.


```{r, message = FALSE}
# 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.

```{r, echo = FALSE, message = FALSE}
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)
)
```


```{r, echo = FALSE}
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

```{r, message = FALSE}
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,
)
```

```{r, echo = FALSE}
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.

```{r, message = FALSE}
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
```


```{r, message = FALSE}
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,
)
```


```{r, echo = FALSE}
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`.


```{r, message = FALSE}
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.


```{r, message = FALSE}
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


```{r, echo = FALSE}
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)
```