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

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?

## 1 Answer

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. • Thank you for the clear answer! This is what I was looking for. I did see that gls() used REML, but most examples for gls() were mixed effect models (which I'm not familiar with), so didn't attempt to use it. It looks like you may have "created" a mixed effect model to use gls()? – ken Sep 24 '15 at 17:13 • @ken Both yes and no. It’s a model fitted using generalised least squares (which becomes a likelihood-based model if you put in an extra assumption that things are normally distributed), and is not per se a mixed-effects model (there are no random effects). But it’s equivalent to a random effects model (one with a random effect for each observation pair, and with heteroscedastic error terms). So you can fit a mixed-effects model (using nls()) and get the same correlation estimate, but it’s a bit more complicated to do so. – Karl Ove Hufthammer Sep 24 '15 at 17:22 • This seems interesting, but I don’t get what gls does, and how the pair #1 can be informative. Can you write a model? and the restricted likelihood which is maximized? – Elvis Sep 24 '15 at 20:29 • @Elvis Regarding how the first sample can be informative, this is rather fascinating. One would expect that having just one observation in a pair couldn’t possible give any information on the correlation. But it does! Remember, we have a full likelihood (Gaussian) model. Simply put, if this one value (of$X$, say) is far away from the rest of$X$s, this implies that the variance of$X$is large(r), which again implies that the correlation is large(r) (which, BTW, is also the reason correlation should never be used in a regression setting, i.e. where one marginal is fixed by design.) – Karl Ove Hufthammer Sep 24 '15 at 20:55 • @Elvis I’ll try to explain the model in more detail / other words. You have a bunch of$X$s and$Y$s. The$(X_i, Y_i)$pairs are independent and have a bivariate normal distribution (ℱ). (This also imples that$X_i$and$Y_j$for$i\neq j$are independent).$X_1\$ has a distribution equal to the first marginal distribution of ℱ. You can write out the likelihood for this, and maximise it to get ML estimates for the five parameters (means, standard deviations and correlation). – Karl Ove Hufthammer Sep 24 '15 at 21:08