Introducción

En el artículo (http://www.biomedcentral.com/1471-2288/13/95) se simularon tiempos Weibull correlacionados a partir de variables normales mediante el siguiente algoritmo:

rho <- 0.3     # Correlacion inicial
k   <- 1.5
c   <- exp(3)  # parametro de la weibull para el grupo control
c1  <- exp(2)  # parametro de la weibull para el grupo tratamiento

id <- 2000  # Numero sujetos
n  <- 3     # Numero maximo de eventos

w <- (pi*log(2)^2)/(pi*log(2)^2 +(1-log(2))^2)
rho_0 <- (-w+sqrt(w^2+2*rho*(1-w)))/(1-w)

y el valor de correlación calculado para \(\rho_0\) es,

# Correlacion final
print(rho_0)
## [1] 0.316

Una vez calculada la correlación \(\rho_0\) generamos las variables aleatorias normales:

# Generamos n variables aleatorias normales con media 0 y 
# varianza 1-(rho_0)^2
# -------------------------------------------------------
set.seed(123456)
Z <- cbind(rnorm(id,0,1), matrix(rnorm(id*(n-1),0, 
           sqrt(1-rho_0*rho_0)), nrow = id))
colnames(Z) <- c("Z1", "Z2", "Z3")
# Generamos normales bivariantes
# ------------------------------
for (j in 1:(n-1)) {
  Z[,j+1] = Z[,j+1] + rho_0 * Z[,j]
}

y obtenemos la siguiente matriz de correlación,

Z1 Z2 Z3
Z1 1.000 0.354 0.140
Z2 0.354 1.000 0.319
Z3 0.140 0.319 1.000

que podemos observar en la siguiente figura ,

luego transformamos estas Normales en variables Uniformes,

#uniformes
U <- pnorm(Z)

y finalmente estas Uniformes pasan a convertirse en variables Weibull,

trat <- c(rep(0,id/2),rep(1,id/2)) 
# weibull con c si control o c1 si tratamiento
W <- ((-log(U)/c)^(1/k))*(1-trat)+((-log(U)/c1)^(1/k))*trat

las que se describen en el siguiente listado,

##        W1              W2              W3       
##  Min.   :0.000   Min.   :0.000   Min.   :0.001  
##  1st Qu.:0.075   1st Qu.:0.075   1st Qu.:0.081  
##  Median :0.140   Median :0.145   Median :0.148  
##  Mean   :0.179   Mean   :0.179   Mean   :0.179  
##  3rd Qu.:0.249   3rd Qu.:0.244   3rd Qu.:0.237  
##  Max.   :0.883   Max.   :0.923   Max.   :1.058

junto a la matriz de correlación,

Cópulas

Otra manera de crear estas variables Weibull correlacionadas es a través del uso de cópulas. Una cópula se define como una distribución multivariante cuyas marginales son todas uniformes (0,1). Sea U un vector aleatorio \(p\)-dimensional definido en un cubo unitario, se define la cópula C por la función de distribución, \[C(u_1, \ldots, u_p)=Pr(U_1 \leq u_1, \ldots, U_p \leq u_p)\] Dado que cualquier variable aleatoria continua puede transformarse en Uniforme (0,1) mediante su función de distribución, las cópulas permiten construir estructuras de dependencia multivariantes con distribuciones marginales fijas. Sea F una función de distribución de dimensión p cuyas marginales son \(F_1, \ldots, F_p\). Sklar(1965) demostró que existe una cópula C de dimensión \(p\) tal que para todo x en el dominio de F, \[F(x_1, \ldots, x_p)=C[F_1(x_1), \ldots, F_p(x_p)]\] Con la librería copula podemos ilustrar la manera de generar variables weibull correlacionadas a partir de variables normales.

En este ejemplo usaremos los siguientes valores de correlación: \(\rho_{12}=0.3\), \(\rho_{13}=0.1\) y \(\rho_{23}=0.3\) para las normales \(Z_1\), \(Z_2\) y \(Z_3\) y como marginales usaremos distribuciones Weibull cada una con parámetros de forma igual a \(1.5\) y de escala \(.25\):

library(copula)
# Creacion normal multivariante con parametros mu y
# correlacion=(0.3,0.1,0.3)
myCop <- normalCopula(param = c(0.3,0.1,0.3), dim = 3, 
                      dispstr = 'un')
myMvd <- mvdc(copula = myCop, 
              margins = c('weibull', 'weibull', 'weibull'),
              paramMargins = list(list(shape = 1.5, scale = .25),
                                  list(shape = 1.5, scale = .25), 
                                  list(shape = 1.5, scale = .25)))
Weib <- rMvdc(2000, myMvd)

en el siguiente listado observamos algunas medidas descriptivas para cada una de las Weibull anteriores,

##        W1              W2              W3       
##  Min.   :0.002   Min.   :0.002   Min.   :0.000  
##  1st Qu.:0.113   1st Qu.:0.115   1st Qu.:0.112  
##  Median :0.201   Median :0.203   Median :0.191  
##  Mean   :0.227   Mean   :0.234   Mean   :0.224  
##  3rd Qu.:0.316   3rd Qu.:0.324   3rd Qu.:0.305  
##  Max.   :0.976   Max.   :0.987   Max.   :0.901

y la matriz de correlaciones,

Cópulas elípticas

Dentro de las clases de cópulas existentes más comunes están las cópulas elípticas y las arquimedianas. Entre las cópulas elípticas mas conocidas tenemos la Gaussiana y la basada en la distribución \(t\). Una cópula elíptica es la cópula correspondiente a una distribución elíptica según el teorema de Sklar. Sea \(F\) la función de distribución acumulada multivariante de una distribución elíptica. Sea \(F_i\) la función de distribución acumulada de la marginal \(i\) y \(F_i^{-1}\) su función inversa, \(i=1,\ldots,p\). La cópula elíptica para \(F\) es,

\[C(u_1, \ldots, u_p) = F[F_1^{-1}(u_1), \ldots, F_p^{-1}(u_p)]\]

Las cópulas elípticas que trae implementada la librería copula son la Gaussiana a través de la función normalCopula basada en la distribución normal multivariante y la \(t\) a través de la función tCopula basada en la distribución t multivariante. Ambas cópulas poseen una matriz de dispersión que es herencia de las distribuciones elípticas. Como las cópulas son invariantes a las transformaciones monótonas de las marginales la matriz de dispersión o matriz de correlación define la estructura de dependencia. Esta librería permite las siguientes estructuras de independencia: autoregresiva de orden 1 (ar1), intercambiable (ex), Toeplitz (toep) y la no estructurada (uns). Cada una de estas estructuras de dependencia tienen la siguiente forma matricial cuando la dimensión es \(p=3\),

\[ %\begin{displaymath} \begin{pmatrix} 1 & \rho_1 & \rho_1^2 \\ \rho_1 & 1 & \rho_1 \\ \rho_1^2 & \rho_1 & 1 \end{pmatrix} , \: \begin{pmatrix} 1 & \rho_1 & \rho_1 \\ \rho_1 & 1 & \rho_1 \\ \rho_1 & \rho_1 & 1 \end{pmatrix} , \: \begin{pmatrix} 1 & \rho_1 & \rho_2 \\ \rho_1 & 1 & \rho_1 \\ \rho_2 & \rho_1 & 1 \end{pmatrix} \: y \: \begin{pmatrix} 1 & \rho_1 & \rho_2 \\ \rho_1 & 1 & \rho_3 \\ \rho_2 & \rho_3 & 1 \end{pmatrix} %\end{displaymath} \]

Referencias

Diaz, W. (2013). Contribuciones a la dependencia y dimensionalidad en cópulas, Universitat de Barcelona. Departament d’Estadística. (http://www.tdx.cat/handle/10803/104204)

Sklar AW (1959). Fonctions de répartition à n dimension et leurs marges. Publications de l’Institut de Statistique de l’Université de Paris, 8, 229–231.

Yan, J. (2007). Enjoy the Joy of Copulas: With a Package copula, Journal of Statistical Software, 21(4),1-21. (http://www.jstatsoft.org/v21/i04/paper)