Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Using tmvtnorm and parametric bootstrap to calculate the correlation coefficient in extreme group design.
## 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment