Using tmvtnorm and parametric bootstrap to calculate the correlation coefficient in extreme group design.

tmvtnorm_correlation_EGD.R

## This script uses the tmvtnorm package, check it out.

if(!("tmvtnorm" %in% rownames(installed.packages()))) {
  install.packages("tmvtnorm")
}

library("tmvtnorm")

## Seed for reproducibility.
set.seed(313)

## I sample values from subject to the constrain that x > C = 1, 
## variances equal to one and a rho = 0.7, means equal to 0.

C = 1
sigma = matrix(c(1,0.7,0.7,1), nrow = 2, byrow = TRUE)
mean = c(0,0)
N = 100
xx = tmvtnorm::rtmvnorm(N, mean = mean, sigma = sigma, lower = c(C,-Inf),
                        algorithm = "rejection")

## The tmvtnorm package includes a functon for computing the MLE! 

pars = tmvtnorm::mle.tmvnorm(xx, lower = c(C,-Inf))@fullcoef

correlations = c(Naive = cor(xx)[1,2], 
                 MLE   = unname(pars[4]/(sqrt(pars[3])*sqrt(pars[5]))))
correlations

## Result:
##
## Naive       MLE 
## 0.4707635 0.8366379 

## Using parametric bootstrap. This is not super-fast, but shouldn't take more
## than a few minutes.

Nreps = 1000
sigma_boot = matrix(c(pars[3],pars[4],pars[4],pars[5]), nrow = 2, byrow = TRUE)
mean_boot = pars[1:2]

boots = replicate(Nreps, {
  xx_boots = tmvtnorm::rtmvnorm(N, mean = mean_boot, sigma = sigma_boot, 
                                lower = c(C,-Inf), algorithm = "rejection")
  
  pars_boots = tmvtnorm::mle.tmvnorm(xx_boots, lower = c(C,-Inf))@fullcoef
  c(Naive = cor(xx_boots)[1,2], 
    MLE   = pars_boots[4]/(sqrt(pars_boots[3])*sqrt(pars_boots[5])))
})

## Nice to take a look at histograms of the parametric bootstrap distribution.

par(mfrow = c(2,1))
hist(boots[1,], breaks = 40, xlab = "Correlation",
     main = "Parametric bootstrap distribution: Naive estimator")
abline(v = mean(boots[1,]), col = "grey", lty = 2)
abline(v = correlations[1], lty = 2)
hist(boots[2,], breaks = 40, xlab = "Correlation",
     main = "Parametric bootstrap distribution: MLE")
abline(v = mean(boots[2,]), col = "grey", lty = 2)
abline(v = correlations[2], lty = 2)

## And percentile confidence intervals.
CI = rbind(Naive = quantile(boots[1,], c(0.025, 0.975)),
           MLE   = quantile(boots[2,], c(0.025, 0.975)))

## Result:
##
##       2.5%     97.5%
## Naive 0.3013961 0.6239027
## MLE   0.5397033 0.9729787