Nothing
#' Generate Random NPD R matrices from a user-supplied population R
#'
#' Generate a list of Random NPD (pseudo) R matrices with a user-defined fixed
#' minimum eigenvalue from a user-supplied population R using the method of
#' alternating projections.
#'
#'
#' @param Rpop input (PD or PSD) p x p Population correlation matrix.
#' @param Lp desired minimum eigenvalue in the NPD matrices.
#' @param NNegEigs number of eigenvalues < 0 in Rnpd.
#' @param NSmoothPosEigs number of eigenvalues > 0 to smooth: the smallest
#' NSmoothPosEigs > 0 will be smoothed toward 0.
#' @param NSubjects sample size (required when NSamples > 0) parameter used to
#' generate sample correlation matrices. Default = NULL.
#' @param NSamples generate NSamples sample R matrices. If NSamples = 0 the
#' program will attempt to find Rnpd such that ||Rpop - Rnpd||_2 is minimized.
#' @param MaxIts maximum number of projection iterations.
#' @param PRINT (logical) If TRUE the program will print the iteration history
#' for Lp. Default = NULL.
#' @param Seed Optional seed for random number generation.
#' @return \item{Rpop}{population (PD) correlation matrix.} \item{R}{sample
#' correlation matrix.} \item{Rnpd}{NPD improper (pseudo) correlation matrix.}
#' \item{Lp}{desired value of minimum eigenvalue.} \item{minEig}{observed value
#' of minimum eigenvalue of Rnpd.} \item{convergence}{0 = converged; 1 = not
#' converged in MaxIts iterations of the alternating projections algorithm.}
#' \item{feasible}{logical) TRUE if max(abs(r_ij)) <= 1. If FALSE then one or
#' more values in Rnpd > 1 in absolute value.} \item{Seed}{saved seed for
#' random number generator.} \item{prbs1}{vector probabilities used to generate
#' eigenvalues < 0.} \item{prbs2}{vector of probabilities used to smooth the
#' smallest NSmoothPosEigs towards zero.}
#' @author Niels G. Waller
#' @keywords fungible
#' @export
#' @import MASS
#' @examples
#'
#' library(MASS)
#'
#' Nvar = 20
#' Nfac = 4
#' NSubj = 2000
#' Seed = 123
#'
#' set.seed(Seed)
#'
#' ## Generate a vector of classical item difficulties
#' p <- runif(Nvar)
#'
#' cat("\nClassical Item Difficulties:\n")
#'
#' print(rbind(1:Nvar,round(p,2)) )
#'
#' summary(p)
#'
#'
#' ## Convert item difficulties to quantiles
#' b <- qnorm(p)
#'
#'
#' ## fnc to compute root mean squared standard deviation
#' RMSD <- function(A, B){
#' sqrt(mean( ( A[lower.tri(A, diag = FALSE)] - B[lower.tri(B, diag = FALSE)] )^2))
#' }
#'
#'
#' ## Generate vector of eigenvalues with clear factor structure
#' L <- eigGen(nDimensions = Nvar,
#' nMajorFactors = Nfac,
#' PrcntMajor = .60,
#' threshold = .50)
#'
#'
#' ## Generate a population R matrix with the eigenvalues in L
#' Rpop <- rGivens(eigs = L)$R
#'
#' ## Generate continuous data that will reproduce Rpop (exactly)
#' X <- mvrnorm(n = NSubj, mu = rep(0, Nvar),
#' Sigma = Rpop, empirical = TRUE)
#'
#' while( any(colSums(X) == 0) ){
#' warning("One or more variables have zero variance. Generating a new data set.")
#' X <- mvrnorm(n = NSubj, mu = rep(0, Nvar),
#' Sigma = Rpop, empirical = TRUE)
#' }
#'
#' ## Cut X at thresholds given in b to produce binary data U
#' U <- matrix(0, nrow(X), ncol(X))
#' for(j in 1:Nvar){
#' U[X[,j] <= b[j],j] <- 1
#' }
#'
#' ## Compute tetrachoric correlations
#' Rtet <- tetcor(U, Smooth = FALSE, PRINT = TRUE)$r
#' # Calculate eigenvalues of tetrachoric R matrix
#' Ltet <- eigen(Rtet)$values
#'
#' if(Ltet[Nvar] >= 0) stop("Rtet is P(S)D")
#'
#' ## Simulate NPD R matrix with minimum eigenvalue equal to
#' # min(Ltet)
#' out <- RnpdMAP(Rpop,
#' Lp = Ltet[Nvar],
#' NNegEigs = Nvar/5,
#' NSmoothPosEigs = Nvar/5,
#' NSubjects = 150,
#' NSamples = 1,
#' MaxIts = 15000,
#' PRINT = FALSE,
#' Seed = Seed)
#'
#' ## RLp is a NPD pseudo R matrix with min eigenvalue = min(Ltet)
#' RLp <- out[[1]]$Rnpd
#'
#' ## Calculate eigenvalues of simulated NPD R matrix (Rnpd)
#' Lnpd <- eigen(RLp, only.values = TRUE)$values
#'
#' ## Scree plots for observed and simulated NPD R matrices.
#' ytop <- max(c(L,Lnpd,Ltet))
#' pointSize = .8
#' plot(1:Nvar, L, typ = "b", col = "darkgrey", lwd=3,
#' lty=1,
#' main =
#' "Eigenvalues of Rpop, Tet R, and Sim Tet R:
#' \nSimulated vs Observed npd Tetrachoric R Matrices",
#' ylim = c(-1, ytop),
#' xlab = "Dimensions",
#' ylab = "Eigenvalues",
#' cex = pointSize,cex.main = 1.2)
#' points(1:Nvar, Lnpd, typ="b",
#' col = "red", lwd = 3, lty=2, cex=pointSize)
#' points(1:Nvar, Ltet, typ="b",
#' col = "darkgreen", lwd = 3, lty = 3, cex= pointSize)
#'
#' legend("topright",
#' legend = c("eigs Rpop", "eigs Sim Rnpd", "eigs Emp Rnpd"),
#' col = c("darkgrey", "red","darkgreen"),
#' lty = c(1,2,3),
#' lwd = c(4,4,4), cex = 1.5)
#'
#' abline(h = 0, col = "grey", lty = 2, lwd = 4)
#'
#' cat("\nRMSD(Rpop, Rtet) = ", round(rmsd(Rpop, Rtet), 3))
#' cat("\nRMSD(Rpop, RLp) = ", round(rmsd(Rpop, RLp), 3))
#'
RnpdMAP <- function (Rpop,
Lp = NULL,
NNegEigs = 1,
NSmoothPosEigs = 4,
NSubjects = NULL,
NSamples = 0,
MaxIts = 15000,
PRINT=FALSE,
Seed = NULL)
{
#############################################################################
## Sample NPD R matrices from a user-supplied
## population R using the alternating projections algorithm
##
## Niels Waller
## Dec 12 2017 3:10pm
##
## Arguments:
## Rpop Input (PD or PSD) p x p Population correlation matrix
##
## Lp Desired minimum eigenvalue in the NPD matrices
##
## NNegEigs Number of eigenvalues < 0 in Rnpd
##
##
## NSmoothPosEigs Number of eigenvalues > 0 to smooth: the smallest
## NSmoothPosEigs > 0 be smoothed toward 0.
##
## NSubjects Sample size (required when NSamples > 0) parameter used to
## generate sample correlation matrices. Default = NULL.
##
## NSamples Generate NSamples sample R matrices. If NSamples = 0
## the program will attempt to find Rnpd such that
## ||Rpop - Rnpd||_2 is minimized.
##
## MaxIts maximum number of projection interations.
##
## PRINT (Logical) If TURE the program will print the iteration history for Lp.
## Default = NULL.
## Seed Optional seed for random number generation.
##
## Values:
## Rpop Population correlation matrix.
##
## R Sample correlation matrix.
##
## Rnpd NPD improper (pseudo) correlation matrix.
##
## Lp desired value of minimum eigenvalue.
##
## minEig observed value of minimum eigenvalue of Rnpd.
##
## convergence 0 = converged
## 1 = not converged in MaxIts.
##
## feasible (Logical) TRUE if max(abs(r_ij)) <=1
##
## Seed saved seed for random number generator.
##
## prbs1
##
## prbs2
##########################################################################
if(NSmoothPosEigs == 0) NSmoothPosEigs <- NULL
prbs2 <- NULL
if(is.null(Seed)){
Seed<- as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31)
}
set.seed(Seed)
## Generate monotonically decreasing negative eigenvalues with
## min eig = Lp
prbs1Fnc <- function(){
prbs <- 1
if(NNegEigs > 1){
prbs <- 1- c(log(sort(runif(NNegEigs-1,1,100), decreasing = TRUE))/log(100),0)
} #end if else
prbs
} ## END prbsFnc
if( NSamples > 0 ){
## If NSamples > 0 then Generate sample correlation matrices from Wishart distribution
WishartArray <- rWishart(n = NSamples, df = NSubjects-1, Sigma = Rpop)
WishartList <-lapply( seq(dim(WishartArray)[3]), function(x) WishartArray[ , , x] )
## scale to corr matrices with cov2cor
corList<-lapply(WishartList,cov2cor)
}
else if (NSamples == 0)
corList <- list(Rpop)
## Lpop = eigenvalues of population PD R
Lpop <- eigen(Rpop, only.values = TRUE)$values
## Project matrix onto Symmetric A with min eig Lp
## and (NNegEigs) smoothed eigs < 0
proj.S.Lp <- function(A, Lp, Nvar, prbs1, prbs2){
VLV <- eigen(A)
V <- VLV$vectors
L <- VLV$values
L[L < Lp] <- Lp
if(NNegEigs == 1){
L[Nvar] <- Lp
} else{
# Smooth the smaller eigevalues < 0
L[ (Nvar - (NNegEigs - 1)):Nvar] <- prbs1 * Lp
# Smooth the smaller eigevalues > 0
if(!is.null(NSmoothPosEigs)){
L[(Nvar - (NNegEigs + NSmoothPosEigs )):(Nvar - (NNegEigs + 1 ))] <-
prbs2 * Lpop[(Nvar - (NNegEigs + NSmoothPosEigs )):(Nvar - (NNegEigs + 1 ))]
}
} # End else
# reconstruct symmetric Cov matrix
A <- V %*% diag(L) %*% t(V)
(A + t(A))/2
} #End proj.S.Lp
## Project onto symmetric matrix set with unit diagonals in bounding cube
proj.U <- function(A){
diag(A) <- 1
# find row col indices for |rij| > 1
# and bring elements into cube
rc.ind <- which(abs(A) > 1, arr.ind = TRUE)
A[rc.ind] <- sign(A[rc.ind])
A
}
## Find NPD pseudo R matrix
RAPA <- function(R, Lp, MaxIts = 1000){
Nvar <- ncol(R)
EigTest <- FALSE
tries <- 0
converged <- TRUE
EigHx <- rep(0, MaxIts)
Rk <- R
# prbs1 and prbs2 are use to smooth the rate of
# descrease in the smallest eigenvalues.
# prbs1 smooths NNegEigs eigs < 0
prbs1 <- prbs1Fnc()
# prbs2 smooths the smallest NSmoothPosEigs eigs > 0
if(!is.null(NSmoothPosEigs)){
prbs2 <- log(sort(runif(NSmoothPosEigs,1,100), decreasing = TRUE))/log(100)
}
while (!EigTest)
{
tries <- tries + 1
if (tries > MaxIts) {
converged = FALSE
break
}
## Project onto S with min Eig Lp
Sk <- proj.S.Lp(Rk, Lp, Nvar, prbs1, prbs2)
## project onto space of Symmetric matrices, Uk, with unit diagonal (in bounding hypercube)
Uk <- proj.U(Sk)
# Rk+1 = Uk
Rk <- Uk
minEig <- eigen(Rk, symmetric = TRUE, only.values = TRUE)$val[Nvar]
EigHx[tries] <- minEig
if(PRINT){
cat("\nAt iter ", tries, " min eig = ", eigen(Rk)$val[Nvar])
}
## test convergence
if(1e+10*(minEig - Lp)^2 < 1e-12) #
EigTest <- TRUE
} ## END while loop
EigHx <- EigHx[1:tries]
colnames(Rk) <- rownames(Rk) <- colnames(R)
feasible <- TRUE
## if any |rxy| > 1 then reject solution
if(max(abs(Rk)) > 1) feasible <- FALSE
list(Rpop = Rpop,
R=R,
Rnpd = Rk,
Lp = Lp,
ObsLp = minEig,
EigHx = EigHx,
converged = converged,
feasible = feasible,
Seed = Seed,
prbs1 = prbs1,
prbs2 = prbs2)
} #End internal function RAPA
## Call Rnpd and return output as a list
lapply(corList, RAPA, Lp, MaxIts)
} ## END FUNCTION
#out <- RnpdMAP(Rpop=Rpop, Lp = -.5, NSubjects = 200, NSamples = 10, MaxIts = 5000, PRINT=FALSE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.