How to use REML to estimate correlation with missing data in R?

Question

In JMP Multivariate Methods, REML is used to estimate correlation when there are missing data values (pg. 28). However, there is no documentation describing how this is done.

I'm trying to compare my results in R with my results in JMP. However, in R, I don't know of a way to calculate correlation or covariance when there are NAs (without removing the row with the NA value). In JMP it appears they don't drop the row, so JMP somehow includes the data in the estimation of correlation and covariance. Is there a way to replicate the JMP results in R? Here is an example dataset:

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

In 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

Here is the output in JMP where they use REML to estimate the correlation when the first value of Var1 is empty (equivalent to NA):

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

Can anyone explain to me how the JMP value is calculated when there is a missing value? Is there a way to calculate correlation in R when there are missing values?

Answer 1

First we need to make sure we know what we’re looking for. You mention REML. The ML part of this is ‘maximum likelihood’, which means we will base all our calculations on a complete likelihood model. That is, we make an assumption that the data come from a specific (known) parametric probability distribution, where some of the parameters are unknown.

The probability distribution JMP bases its calculations on is the bivariate normal distribution. That is, we assume (in addition to the usual i.i.d. assumptions) that 1) the first variable (let’s call it $X$) has a normal distribution, 2) the second variable ($Y$) also has a normal distribution, and 3) that every linear combination $aX+bY$ also has a normal distribution. This is a rather strong assumption. (Compare this to the usual sample correlation, which makes no assumptions on the distribution of $X$ and $Y$ except that their (co)variances exist.) Note that we do not assume that the means of $X$ and $Y$ are equal, nor that their standard deviations are equal.

So now we have a model, and can fit it to our data. In R, we first need the data to be in ‘long’ format, with one row for each observation, for a total of 17 observations (a complete $(x,y)$ pair is treated as 2 observations):

> 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

Note that we have added a pair variable to keep track of which $(x,y)$ pair each observation originally came from, and use na.rm = TRUE to remove the missing observation.

No we fit the model using REML. There are several ways of doing this, but the following is relatively simple:

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

In the first argument we specify that the means are allowed to differ between $X$ and $Y$, in the second we specify that the standard deviations too are allowed to differ, and in the third we say that there are no constrains on the correlation (except that it must be a valid correlation, and that it is the same for each obseration pair) and that the pair variable identifies observations that are correlated (all other observations pairs are assumed to be uncorrelated, cf. the i.i.d. assumption).

(If you want to assume that the means are identical, feel free to change value ~ variable to value ~ 1. If you want to assume that the standard deviations are identical, feel free to remove the weights argument.)

Running summary(fit) gives the REML estimate of the correlation (and lots of other parameters):

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

It’s somewhat complicated to extract the correlation as a number for later use, but the following code works (technically it extracts the estimated correlation for the second observation pair, but since the correlation is assumed to be identical for each pair, this gives the correct answer):

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

This is (to a very close approximation) the same as you get from JMP.