Nothing
##########################################################
# This function compute the asymptotic distribution-free #
# covariance matrix of correlations. #
# #
# Arguments: #
# X - matrix of predictor scores #
# y - vector of criterion scores #
# #
# Output #
# adfCovMat - Asymptotic distribution-free estimate of #
# the covariance matrix #
##########################################################
#' Asymptotic Distribution-Free Covariance Matrix of Correlations
#'
#' Function for computing an asymptotic distribution-free covariance matrix of
#' correlations.
#'
#'
#' @param X Data matrix.
#' @param y Optional vector of criterion scores.
#' @return \item{adfCorMat}{Asymptotic distribution-free estimate of the
#' covariance matrix of correlations.}
#' @author Jeff Jones and Niels Waller
#' @references Browne, M. W. (1984). Asymptotically distribution-free methods
#' for the analysis of covariance structures. \emph{British Journal of
#' Mathematical and Statistical Psychology, 37}, 62--83.
#'
#' Steiger, J. H. and Hakstian, A. R. (1982). The asymptotic distribution of
#' elements of a correlation matrix: Theory and application. \emph{British
#' Journal of Mathematical and Statistical Psychology, 35}, 208--215.
#' @keywords Statistics
#' @examples
#'
#' ## Generate non-normal data using monte1
#' set.seed(123)
#' ## we will simulate data for 1000 subjects
#' N <- 1000
#'
#' ## R = the desired population correlation matrix among predictors
#' R <- matrix(c(1, .5, .5, 1), 2, 2)
#'
#' ## Consider a regression model with coefficient of determination (Rsq):
#' Rsq <- .50
#'
#' ## and vector of standardized regression coefficients
#' Beta <- sqrt(Rsq/t(sqrt(c(.5, .5))) %*% R %*% sqrt(c(.5, .5))) * sqrt(c(.5, .5))
#'
#' ## generate non-normal data for the predictors (X)
#' ## x1 has expected skew = 1 and kurtosis = 3
#' ## x2 has expected skew = 2 and kurtosis = 5
#' X <- monte1(seed = 123, nvar = 2, nsub = N, cormat = R, skewvec = c(1, 2),
#' kurtvec = c(3, 5))$data
#'
#' ## generate criterion scores
#' y <- X %*% Beta + sqrt(1-Rsq)*rnorm(N)
#'
#' ## Create ADF Covariance Matrix of Correlations
#' adfCor(X, y)
#'
#' #> 12 13 23
#' #> 12 0.0012078454 0.0005331086 0.0004821594
#' #> 13 0.0005331086 0.0004980130 0.0002712080
#' #> 23 0.0004821594 0.0002712080 0.0005415301
#'
#' @export
adfCor <- function(X,y=NULL) {
######################## Internal Functions ########################
# row.remove: Removes rows from a symmetric transition matrix
# to create a correlation transition matrix
# (see Browne & Shapiro, (1986); Nel, 1985).
row.remove <- function(p) {
p1 <- p2 <- p
rows <- rep(1,p)
for(i in 2:p) {
rows[i] <- rows[i] + p1
p1 <- p1 + (p2-1)
p2 <- p2 - 1
}
rows
}
# Dp: Duplicator Matrix
Dp <- function(p) {
M <- matrix(nrow = p, ncol = p)
M[ lower.tri(M, diag = T) ] <- seq( p*(p + 1)/2 )
M[ upper.tri(M, diag = F) ] <- t(M)[ upper.tri(M, diag = F) ]
D <- outer(c(M), unique(c(M)),
FUN = function(x, y) as.numeric(x == y) )
D
}
###################### End Internal Functions ######################
####################################################################
Xy <- if(is.null(y)) X else cbind(X,y)
dev <- scale(Xy, scale=FALSE)
nvar <- ncol(dev)
N <- nrow(dev)
# order of the covariance matrix of covariances
ue <- nvar^2
# container for indices
s <- vector(length=ue, mode="character")
z <- 0
for(i in 1:nvar){
for(j in 1:nvar){
z<-z+1
s[z]<-paste(i,j,sep="")
}
}
# computes all possible combinations of the
# indices in s
v <- expand.grid(s, s)
# paste the index pairs togehter
V <- paste(v[,1], v[,2], sep="")
# separate the indices into their own columns
ids <- matrix(0,nrow=ue^2,4)
for(i in 1:4) ids[,i] <- as.numeric(sapply(V,substr,i,i))
covs <- matrix(0,(ue^2),1)
# compute the covariances using Steiger and Hakstian (1982) Eqn 3.4
for(i in 1:(ue^2)) {
w_ii <- cov(dev[,ids[i,1]],dev[,ids[i,1]])
w_jj <- cov(dev[,ids[i,2]],dev[,ids[i,2]])
w_kk <- cov(dev[,ids[i,3]],dev[,ids[i,3]])
w_ll <- cov(dev[,ids[i,4]],dev[,ids[i,4]])
w_ij <- cov(dev[,ids[i,1]],dev[,ids[i,2]])
w_kl <- cov(dev[,ids[i,3]],dev[,ids[i,4]])
w_ik <- cov(dev[,ids[i,1]],dev[,ids[i,3]])
w_jl <- cov(dev[,ids[i,2]],dev[,ids[i,4]])
w_ijkl <- t(dev[,ids[i,1]]*dev[,ids[i,2]])%*%(dev[,ids[i,3]]*dev[,ids[i,4]])/(N-1)
w_iikk <- t(dev[,ids[i,1]]*dev[,ids[i,1]])%*%(dev[,ids[i,3]]*dev[,ids[i,3]])/(N-1)
w_jjkk <- t(dev[,ids[i,2]]*dev[,ids[i,2]])%*%(dev[,ids[i,3]]*dev[,ids[i,3]])/(N-1)
w_iill <- t(dev[,ids[i,1]]*dev[,ids[i,1]])%*%(dev[,ids[i,4]]*dev[,ids[i,4]])/(N-1)
w_jjll <- t(dev[,ids[i,2]]*dev[,ids[i,2]])%*%(dev[,ids[i,4]]*dev[,ids[i,4]])/(N-1)
w_iikl <- t(dev[,ids[i,1]]*dev[,ids[i,1]])%*%(dev[,ids[i,3]]*dev[,ids[i,4]])/(N-1)
w_jjkl <- t(dev[,ids[i,2]]*dev[,ids[i,2]])%*%(dev[,ids[i,3]]*dev[,ids[i,4]])/(N-1)
w_ijkk <- t(dev[,ids[i,1]]*dev[,ids[i,2]])%*%(dev[,ids[i,3]]*dev[,ids[i,3]])/(N-1)
w_ijll <- t(dev[,ids[i,1]]*dev[,ids[i,2]])%*%(dev[,ids[i,4]]*dev[,ids[i,4]])/(N-1)
r_ij <- w_ij/sqrt(w_ii*w_jj)
r_kl <- w_kl/sqrt(w_kk*w_ll)
r_ik <- w_ik/sqrt(w_ii*w_kk)
r_jl <- w_jl/sqrt(w_jj*w_ll)
r_ijkl <- w_ijkl/sqrt(w_ii*w_jj*w_kk*w_ll)
r_iikk <- w_iikk/sqrt(w_ii*w_ii*w_kk*w_kk)
r_jjkk <- w_jjkk/sqrt(w_jj*w_jj*w_kk*w_kk)
r_iill <- w_iill/sqrt(w_ii*w_ii*w_ll*w_ll)
r_jjll <- w_jjll/sqrt(w_jj*w_jj*w_ll*w_ll)
r_iikl <- w_iikl/sqrt(w_ii*w_ii*w_kk*w_ll)
r_jjkl <- w_jjkl/sqrt(w_jj*w_jj*w_kk*w_ll)
r_ijkk <- w_ijkk/sqrt(w_ii*w_jj*w_kk*w_kk)
r_ijll <- w_ijll/sqrt(w_ii*w_jj*w_ll*w_ll)
covs[i] <- (r_ijkl + .25*r_ij*r_kl*(r_iikk + r_jjkk + r_iill + r_jjll) -
.5*r_ij*(r_iikl + r_jjkl) - .5*r_kl*(r_ijkk + r_ijll))
}
adfFull <- matrix(covs,ue,ue,byrow=T)
# Create symmetric transition matrix (see Nel, 1985)
Kp <- solve(t(Dp(nvar)) %*% Dp(nvar)) %*% t(Dp(nvar))
# Create correlation transition matrix (see Browne & Shapiro, 1986).
Kpc <- as.matrix(Kp[-row.remove(nvar),] )
if(ncol(Kpc) == 1) Kpc <- t(Kpc)
adfCovMat <- Kpc %*% adfFull %*% t(Kpc)/N
## add row and column labels
rc<-expand.grid(1:nvar,1:nvar)
rc<-rc[!(rc[,1]==rc[,2]),]
rc<-unique(apply(t(apply(rc,1,sort)),1,paste, collapse=""))
dimnames(adfCovMat)<-list(rc,rc)
adfCovMat
}
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.