¿Cómo utilizar REML para estimar la correlación con datos perdidos en R?

Question

En JMP Multivariate Methods, El REML se utiliza para estimar la correlación cuando hay valores de datos que faltan (pág. 28). Sin embargo, no hay documentación que describa cómo se hace.

Estoy tratando de comparar mis resultados en R con mis resultados en JMP. Sin embargo, en R, no conozco una forma de calcular la correlación o la covarianza cuando hay NAs (sin eliminar la fila con el valor NA). En JMP parece que no eliminan la fila, por lo que JMP incluye de alguna manera los datos en la estimación de la correlación y la covarianza. ¿Hay alguna manera de replicar los resultados de JMP en R? Aquí hay un conjunto de datos de ejemplo:

Var1    Var2
0.0079  0.0046
0.0136  0.0080
0.0270  0.0108
0.0287  0.0263
0.0325  0.0400
0.0228  0.0163
0.0015  0.0014
0.1198  0.0869
0.0054  0.0046

En R:

# No NAs
df <- data.frame("Var1" = c(0.0079, 0.0136, 0.0270, 0.0287, 0.0325, 0.0228, 0.0015, 0.1198, 0.0054),
       "Var2" = c(0.0046, 0.0080, 0.0108, 0.0263, 0.0400, 0.0163, 0.0014, 0.0869, 0.0046))
df_cor <- cor(df$Var1, df$Var2) # result: 0.9685043

# With NA in first row
df_withNA <- data.frame("Var1" = c(0.0079, 0.0136, 0.0270, 0.0287, 0.0325, 0.0228, 0.0015, 0.1198, 0.0054),
             "Var2" = c(NA, 0.0080, 0.0108, 0.0263, 0.0400, 0.0163, 0.0014, 0.0869, 0.0046))
df_withNA_cor <- cor(df_withNA$Var1, df_withNA$Var2) # result: NA

# Drop row with NA
df_dropNA_cor <- cor(df_withNA$Var1, df_withNA$Var2, use = "pairwise.complete.obs")
# result: 0.9670248

Aquí está la salida en JMP donde utilizan REML para estimar la correlación cuando el primer valor de Var1 está vacío (equivalente a NA):

# JMP Correlation between Var1 and Var2 with missing value: 0.9646

¿Alguien puede explicarme cómo se calcula el valor de JMP cuando falta un valor? ¿Hay alguna forma de calcular la correlación en R cuando hay valores perdidos?

Answer 1

Primero tenemos que asegurarnos de que sabemos lo que estamos buscando. Usted menciona REML. La parte ML es "máxima verosimilitud", lo que significa que basaremos todos nuestros cálculos en un modelo de máxima verosimilitud. Es decir, hacemos una suposición que los datos provienen de un específico Distribución de probabilidad paramétrica (conocida), donde algunos de los parámetros son desconocidos.

La distribución de probabilidad en la que se basa el programa JMP es la distribución normal bivariante. Es decir, que suponga que (además de los supuestos i.i.d. habituales) que 1) la primera variable (llamémosla $X$ ) tiene una distribución normal, 2) la segunda variable ( $Y$ ) también tiene una distribución normal, y 3) que cada combinación lineal $aX+bY$ también tiene una distribución normal. Esta es una suposición bastante fuerte. (Compárese con la correlación muestral habitual, que hace no hipótesis sobre la distribución de $X$ y $Y$ excepto que sus (co)variantes existen). Obsérvese que hacemos no asumir que los medios de $X$ y $Y$ son iguales, ni que sus desviaciones estándar sean iguales.

Así que ahora tenemos un modelo, y podemos ajustarlo a nuestros datos. En R, primero necesitamos que los datos estén en formato "largo", con un fila para cada observación, para un total de 17 observaciones (una $(x,y)$ se trata como 2 observaciones):

> library(tidyr)
> df_withNA$pair = 1:nrow(df_withNA)
> d = gather(df_withNA, variable, value,
             -pair, na.rm = TRUE)
> d
   pair variable  value
1     1     Var1 0.0079
2     2     Var1 0.0136
3     3     Var1 0.0270
[…]
8     8     Var1 0.1198
9     9     Var1 0.0054
10    2     Var2 0.0080
11    3     Var2 0.0108
12    4     Var2 0.0263
[…]
16    8     Var2 0.0869
17    9     Var2 0.0046

Tenga en cuenta que hemos añadido un pair para saber qué $(x,y)$ par de cada observación originalmente vino de, y el uso na.rm = TRUE para eliminar la observación que falta.

No, ajustamos el modelo utilizando REML. Hay varias formas de hacerlo, pero la siguiente es relativamente sencilla:

library(nlme)
fit = gls(value ~ variable,
        weights = varIdent(form = ~1|variable),
        correlation = corSymm(form = ~1 | pair),
        data = d)

En el primer argumento especificamos que se permite que las medias difieran entre $X$ y $Y$ en la segunda especificamos que también se permite que las desviaciones estándar difieran, y en la tercera decimos que no hay restricciones en la correlación (excepto que debe ser una correlación válida, y que es la misma para cada par de observación) y que el pair identifica las observaciones que están correlacionadas (todos los demás pares de observaciones se suponen no correlacionados, véase la hipótesis i.i.d.).

(Si quiere suponer que las medias son idénticas, no dude en cambiar value ~ variable a value ~ 1 . Si quiere suponer que las desviaciones estándar son idénticas, no dude en eliminar el weights argumento).

Ejecutar summary(fit) da la estimación REML de la correlación (y muchos otros parámetros):

Correlation Structure: General
 Formula: ~1 | pair 
 Parameter estimate(s):
 Correlation: 
  1    
2 0.964

Es algo complicado extraer la correlación como un número para su uso posterior, pero el siguiente código funciona (técnicamente extrae la correlación estimada para el segundo par de observaciones, pero como se supone que la correlación es idéntica para cada par, esto da la respuesta correcta):

> corMatrix(fit$modelStruct$corStruct)[[2]][1,2]
[1] 0.9642673

Esto es (en una aproximación muy cercana) lo mismo que se obtiene de JMP.