Find the Nearest Proper Correlation Matrix
nearcor.RdThis function “smoothes” an improper correlation matrix as it can result
from cor with use="pairwise.complete.obs" or
hetcor.
It is deprecated now, in favor of
nearPD() from package Matrix.
nearcor(R, eig.tol = 1e-6, conv.tol = 1e-07, posd.tol = 1e-8,
maxits = 100, verbose = FALSE)Arguments
- R
a square symmetric approximate correlation matrix
- eig.tol
defines relative positiveness of eigenvalues compared to largest, default=1e-6.
- conv.tol
convergence tolerance for algorithm, default=1.0e-7
- posd.tol
tolerance for enforcing positive definiteness, default=1.0e-8
- maxits
maximum number of iterations
- verbose
logical specifying if convergence monitoring should be verbose.
Details
This implements the algorithm of Higham (2002), then forces symmetry,
then forces positive definiteness using code from
posdefify. This implementation does not make
use of direct LAPACK access for tuning purposes as in the MATLAB code
of Lucas (2001). The algorithm of Knol DL and ten Berge (1989) (not
implemented here) is more general in (1) that it allows contraints to
fix some rows (and columns) of the matrix and (2) to force the
smallest eigenvalue to have a certain value.
Value
A list, with components
- cor
resulting correlation matrix
- fnorm
Froebenius norm of difference of input and output
- iterations
number of iterations used
- converged
logical
References
See those in posdefify.
See also
Examples
cat("pr is the example matrix used in Knol DL, ten Berge (1989)\n")
#> pr is the example matrix used in Knol DL, ten Berge (1989)
pr <- matrix(c(1, 0.477, 0.644, 0.478, 0.651, 0.826,
0.477, 1, 0.516, 0.233, 0.682, 0.75,
0.644, 0.516, 1, 0.599, 0.581, 0.742,
0.478, 0.233, 0.599, 1, 0.741, 0.8,
0.651, 0.682, 0.581, 0.741, 1, 0.798,
0.826, 0.75, 0.742, 0.8, 0.798, 1),
nrow = 6, ncol = 6)
ncr <- nearcor(pr)
nr <- ncr$cor
plot(pr[lower.tri(pr)],
nr[lower.tri(nr)]); abline(0,1, lty=2)
round(cbind(eigen(pr)$values, eigen(nr)$values), 8)
#> [,1] [,2]
#> [1,] 4.2132 4.20e+00
#> [2,] 0.7715 7.54e-01
#> [3,] 0.5154 5.08e-01
#> [4,] 0.3845 3.80e-01
#> [5,] 0.1780 1.61e-01
#> [6,] -0.0626 4.00e-08
cat("The following will fail:\n")
#> The following will fail:
try(factanal(cov=pr, factors=2))
#> Warning: NaNs produced
#> Error in optim(start, FAfn, FAgr, method = "L-BFGS-B", lower = lower, :
#> L-BFGS-B needs finite values of 'fn'
cat("and this should work\n")
#> and this should work
try(factanal(cov=nr, factors=2))
#>
#> Call:
#> factanal(factors = 2, covmat = nr)
#>
#> Uniquenesses:
#> [1] 0.309 0.232 0.441 0.005 0.307 0.005
#>
#> Loadings:
#> Factor1 Factor2
#> [1,] 0.735 0.388
#> [2,] 0.868 0.123
#> [3,] 0.534 0.523
#> [4,] 0.151 0.986
#> [5,] 0.508 0.660
#> [6,] 0.740 0.669
#>
#> Factor1 Factor2
#> SS loadings 2.407 2.295
#> Proportion Var 0.401 0.382
#> Cumulative Var 0.401 0.784
#>
#> The degrees of freedom for the model is 4 and the fit was 13.8129
if(require("polycor")) {
n <- 400
x <- rnorm(n)
y <- rnorm(n)
x1 <- (x + rnorm(n))/2
x2 <- (x + rnorm(n))/2
x3 <- (x + rnorm(n))/2
x4 <- (x + rnorm(n))/2
y1 <- (y + rnorm(n))/2
y2 <- (y + rnorm(n))/2
y3 <- (y + rnorm(n))/2
y4 <- (y + rnorm(n))/2
dat <- data.frame(x1, x2, x3, x4, y1, y2, y3, y4)
x1 <- ordered(as.integer(x1 > 0))
x2 <- ordered(as.integer(x2 > 0))
x3 <- ordered(as.integer(x3 > 1))
x4 <- ordered(as.integer(x4 > -1))
y1 <- ordered(as.integer(y1 > 0))
y2 <- ordered(as.integer(y2 > 0))
y3 <- ordered(as.integer(y3 > 1))
y4 <- ordered(as.integer(y4 > -1))
odat <- data.frame(x1, x2, x3, x4, y1, y2, y3, y4)
xcor <- cor(dat)
pcor <- cor(data.matrix(odat)) # cor() no longer works for factors
hcor <- hetcor(odat, ML=TRUE, std.err=FALSE)$correlations
ncor <- nearcor(hcor)$cor
try(factanal(covmat=xcor, factors=2, n.obs=n))
try(factanal(covmat=pcor, factors=2, n.obs=n))
try(factanal(covmat=hcor, factors=2, n.obs=n))
try(factanal(covmat=ncor, factors=2, n.obs=n))
}
#> Loading required package: polycor
#> Warning: there is no package called ‘polycor’