R/mixed.cor.R

#mixedCor was added April 28, 2017 to make mixed.cor easier to use
"mixed.cor" <- 
function(x=NULL,p=NULL,d=NULL,smooth=TRUE,correct=.5,global=TRUE,ncat=8,use="pairwise",method="pearson",weight=NULL)  {
cat("\nmixed.cor is deprecated, please use mixedCor.")
mixedCor(data=x,c=NULL,p=p,d=d,smooth=smooth,correct=correct,global=global,ncat=ncat,use=use,method=method,weight=weight)
}

"mixedCor" <- function(data=NULL,c=NULL,p=NULL,d=NULL,smooth=TRUE,correct=.5,global=TRUE,ncat=8,use="pairwise",method="pearson",weight=NULL) {
cl <- match.call() 
original <- colnames(data)
organize <- FALSE   #the default is to specify the continuous, the polytomous and the dichotomous
if((missing(c) | is.null(c)) && (missing(p) | is.null(p))  && (missing(d) | is.null(d))) { #figure out which kinds of variables we are using 
organize <- TRUE
nvar <- ncol(data)
data <- as.matrix(data)  #to speed things up
#first, check for data coded as factors,  these mess up.   We should quit with warning 
#also check for variables with no variance and flag them (don't remove, just quit.  Likely problem with data)

progressBar(nvar,nvar,"Preparing the data")
ans <- matrix(NA,nrow=nvar,ncol=2)
   for (i in 1:nvar) {
        if (is.numeric(data[,i])) {ans[i,2] <- 1 } else {
            if ((is.factor(data[,i])) || (is.logical(data[,i]))) {
               ans[i,2]  <- 2
            } else {
                if (is.character(data[,i])) {
               ans[i,2] <- 3
                } else {ans[i,2] <- 4}
            }
        }
       ans[i,1] <- sd(data[,i],na.rm=TRUE)
    }
    if(any(ans[,2] !=1)) {cat("\nNot all of the variables are numeric.  Please check your data. \nPotential bad items are ")
    print(which(ans[,2] !=1)) 
    stop ("\nI am stopping because of the problem with the data.")
    }
    
       
 bad <- which(ans[,1] ==0)   
 if(length(bad) > 0) {cat("\nSome of the variables have no variance.  Please remove them and try again.\nBad items are ")
      print(bad)
      stop("\nI am stopping because of the problems with the data.")
      }


tab <- apply(data,2,table)
if(is.list(tab)) {len <- lapply(tab,length)} else {len <- dim(tab)[1] }
dvars <- subset(1:nvar,len==2)   #find the dichotomous variable by number
pvars <- subset(1:nvar,((len > 2) & (len <= ncat)))  #find the polytomous variables
cvars <- subset(1:nvar,(len > ncat))  #find the continuous variables (more than ncat levels)

#This next part is not as efficient as it should be, because it is making up 3 new matrices, where we could just do it by reference -- pass the names, not the data
#if(length(dvars) > 0) {d <- as.matrix(x[,dvars],ncol=length(dvars))
#              colnames(d) <- colnames(x)[dvars]} else {d <- NULL}
if(length(pvars) > 0) {#p <- as.matrix(x[,pvars],ncol=length(pvars))
              # colnames(p) <- colnames(x)[pvars] 
              # tab <- table(p) #now check to make sure that they are all on the same scale
              # if(length(tab) > ncat) stop("I tried to figure out which were continuous and which were polytomous, but failed.  Please try again by specifying x, p, and d.")
              tab <- apply(data[,pvars,drop=FALSE],2,table)  #if all elements have the same number of categories, this is a matrix, otherwise, it is a list
             if(is.list(tab)) {
               ok <- apply(data[,pvars,drop=FALSE], 2,function (x) {if (length(table(x)) != (max(x,na.rm=TRUE) - min(x,na.rm=TRUE)+1)) {FALSE} else {TRUE}})
               if(any(!ok)) {bad <- which(!ok)
               cat("\n Some polytomous variables have fewer categories than they should.  Please check your data.  \nPotential bad items are ",colnames(p)[bad],"\n")
             
               stop("\nI am stopping because of the problem with polytomous data")
               }
               }
                         } else {p <- NULL}
# if(length(cvars) > 0) {cont <- matrix(x[,cvars],ncol=length(cvars))
#                       colnames(cont) <- colnames(x)[cvars] } else {cont <- NULL}
                       
#Rho <- mixed.cor1(cont,p, d,smooth=smooth,global=global,correct=correct,use=use,method=method,weight=weight)
progressBar(nvar,nvar,"Starting mixed.cor1")
Rho <- mixed.cor1(data,cvars,pvars, dvars,smooth=smooth,global=global,correct=correct,use=use,method=method,weight=weight)
oldorder <- c(cvars,pvars,dvars)
ord <- order(oldorder)
Rho$rho <- Rho$rho[ord,ord]
} else {# organization is specified
#if ((p+ d) == nvar)
#if(!missing(c)) x <- data[c]
#if(!missing(p)) p <- data[p]
#if(!missing(d)) d <- data[d]
progressBar(1,1,"Starting mixed.cor1")
Rho <- mixed.cor1(data,c=c,p=p,d=d,smooth=smooth,global=global,correct=correct,use=use,method=method,weight=weight)}
orig <- original %in%  colnames(Rho$rho)
orig <- original[orig]
Rho$rho <- Rho$rho[orig,orig]  #organize them in the way they came in
Rho$Call <- cl
return(Rho)
}

#modified April 29th (2014) to get around the problem of missing rows or columns in the polychoric function.
#modified 1/1/14 to add multicore capability
#deprecated as of 02/05/18
# "mixed.cor" <- 
# function(x=NULL,p=NULL,d=NULL,smooth=TRUE,correct=.5,global=TRUE,ncat=8,use="pairwise",method="pearson",weight=NULL)  {
# cl <- match.call() 
# organize <- FALSE   #the default is to specify the continuous, the polytomous and the dichotomous
# if(!is.null(x) && is.null(p)  && is.null(d)) { #figure out which kinds of variables we are using 
# organize <- TRUE
# nvar <- ncol(x)
# x <- as.matrix(x)
# tab <- apply(x,2,function(x) table(x))
# if(is.list(tab)) {len <- lapply(tab,function(x) length(x))} else {len <- dim(tab)[1] }
# dvars <- subset(1:nvar,len==2)   #find the dichotomous variables
# pvars <- subset(1:nvar,((len > 2) & (len <= ncat)))  #find the polytomous variables
# cvars <- subset(1:nvar,(len > ncat))  #find the continuous variables (more than ncat levels)
# 
# if(length(dvars) > 0) {d <- as.matrix(x[,dvars],ncol=length(dvars))
#               colnames(d) <- colnames(x)[dvars]} else {d <- NULL}
# if(length(pvars) > 0) {p <- as.matrix(x[,pvars],ncol=length(pvars))
#                colnames(p) <- colnames(x)[pvars] 
#                tab <- table(p) #now check to make sure that they are all on the same scale
#                if(length(tab) > ncat) stop("I tried to figure out which were continuous and which were polytomous, but failed.  Please try again by specifying x, p, and d.")
#                ok <- apply(p, 2,function (x) {if (length(table(x)) != (max(x,na.rm=TRUE) - min(x,na.rm=TRUE)+1)) {FALSE} else {TRUE}})
#                if(any(!ok)) {bad <- which(!ok)
#                cat("\n Some polytomous variables have fewer categories than they should.  Please check your data.  \nPotential bad items are ",colnames(p)[bad],"\n")
#              
#                stop("\nI am stopping because of the problem with polytomous data")
#                }
#                          } else {p <- NULL}
#  if(length(cvars) > 0) {cont <- matrix(x[,cvars],ncol=length(cvars))tab <- apply(x,2,table)
# if(is.list(tab)) {len <- lapply(tab,length)} else {len <- dim(tab)[1] }
# dvars <- subset(1:nvar,len==2)   #find the dichotomous variable by number
# pvars <- subset(1:nvar,((len > 2) & (len <= ncat)))  #find the polytomous variables
# cvars <- subset(1:nvar,(len > ncat))  #find the continuous variables (more than ncat levels)

#                        colnames(cont) <- colnames(x)[cvars] } else {cont <- NULL}
#                        
# Rho <- mixed.cor1(cont,p, d,smooth=smooth,global=global,correct=correct,use=use,method=method,weight=weight)
# oldorder <- c(cvars,pvars,dvars)
# ord <- order(oldorder)
# Rho$rho <- Rho$rho[ord,ord]
# } else {# organization is specified
# #if ((p+ d) == nvar) 
# Rho <- mixed.cor1(x=x,p=p,d=d,smooth=smooth,global=global,correct=correct,use=use,method=method,weight=weight)}
# 
# Rho$Call <- cl
# return(Rho)
# }

#December 22,2010
#revised July 15, 2011 to work for the various special cases
#meant to combine continuous, polytomous and dichotomous correlations
#revised October 12, 2011 to get around the sd of vectors problem
#revised Sept 10, 2013 to allow for dichotomies with different minima
"mixed.cor1" <-
function(data,c=NULL,p=NULL,d=NULL,smooth=TRUE,global=TRUE,correct=correct,use=use,method=method,weight=NULL) {
 cl <- match.call() 
 x <- c   #the continuous variables
# if(!is.null(x)) {nx <- dim(x)[2]} else {nx <- 0}
# if(!is.null(p)) {np <- dim(p)[2]} else {np <- 0}
# if(!is.null(d))  {nd <- dim(d)[2]} else {nd <- 0}
# if(is.null(nx)) nx <- 1
# if(is.null(np)) np <- 1
# if(is.null(nd)) nd <- 1
if(!is.null(x)) {nx <- length(x)} else {nx <- 0}
if(!is.null(p)) {np <- length(p)} else {np <- 0}
 if(!is.null(d))  {nd <-length(d)} else {nd <- 0}
npd  <- nx +np + nd
rpd <- NULL #in case we don't have any
#check to make sure all the data are ok for doing the appropriate analyses
#first check for polychorics
#but we did this before in mixedCor.  If the person specified  the p and d, assume they are correct
# if(np > 0) { ptable <- table(as.matrix(data[,p]))   #haven't we already done this?
# nvalues <- length(ptable)  #find the number of response alternatives 
# if(nvalues > 8) stop("You have more than 8 categories for your items, polychoric is probably not needed")}
# 
# #now test for tetrachorics
# if(nd > 0) {
# dm <- apply(data[,d,drop=FALSE],2,function(x) min(x,na.rm=TRUE))
# data[,d]  <- t(t(data[,d]) - dm)  
# #d <- d -min(d,na.rm=TRUE) #in case the numbers are not 0,1   This requires them all to be the same
# if(max(data[,d,drop=FALSE],na.rm=TRUE) > 1) {stop("Tetrachoric correlations require dichotomous data")}}

if(nx > 1) {
progressBar(nx,nx,"Finding Pearson correlations")
rx <- cor(data[,x],use=use,method=method)} else {if(nx  < 1) {rx <- NULL
   rho <- NULL} else {rx <- 1
     rho <- 1}}
     
if(np > 1) {
#cat("\n Starting to find the polychoric correlations")
progressBar(np,np,"Finding polychorics")
 rp <- polychoric(data[,p],smooth=smooth,global=global,weight=weight,correct=correct)}    else {if (np == 1) {
	rho <- 1
	names(rho) <- colnames(p)
	rp <- list(rho=rho,tau=NULL)}  else {rp <- list(rho= NULL,tau=NULL)}}
if(nd > 1) {
# cat("n Starting to find the tetrachoric correlations\n\n")
progressBar(nd,nd,"Finding tetrachorics")
  rd <- tetrachoric(data[,d],smooth=smooth,correct=correct,weight=weight)}   else {if (nd == 1) {rd <- list(rho=1,tau=NULL)}  else {rd <- list(rho=NULL,tau=NULL)}}

if(nx > 0) {if(np > 0) {rxp <- polyserial(data[,x,drop=FALSE],data[,p,drop=FALSE])   #the normal case is for all three to exist
		tmixed <- cbind(rx,t(rxp))
		lmixed <- cbind(rxp,rp$rho)
		rho <- rbind(tmixed,lmixed)} else {rho <- rx}  #we now have x and p

		if(nd > 0) { rxd <- biserial(data[,x,drop=FALSE],data[,d,drop=FALSE])
		
		    if(np > 0) {
		    rpd <- polydi(data[,p,drop=FALSE],data[,d,drop=FALSE],global=global,correct=correct)$rho   #passing the global value (July 10, 2017)
			            topright <- t(cbind(rxd,t(rpd))) 
			            } else {
			        topright <- t(rxd)}
			tmixed <- cbind(rho,topright) 
			lmixed <- cbind(t(topright),rd$rho)
			rho <- rbind(tmixed,lmixed) }

		} else {  #the case of nx =0
		   if( np > 0) { 
		      if (nd > 0 ) {
		     progressBar(nd,nd,"Starting polydi")
		      rpd <- polydi(data[,p,drop=FALSE],data[,d,drop=FALSE],global=global,correct=correct)$rho  #added global  (July 10, 2017)   But this causes problems
		     
		       tmixed <- cbind(rp$rho,rpd)
		       lmixed <- cbind(t(rpd),rd$rho)
		       rho <- rbind(tmixed,lmixed)
		       }  else {rho <- rp$rho} } else {
		    rho <- rd$rho}
		    }
colnames(rho) <- rownames(rho)
class(rho) <- c("psych","mixed")
if(!is.null(rx)) class(rx) <- c("psych","mixed")
mixed <- list(rho=rho,rx=rx,poly=rp,tetra=rd,rpd=rpd,Call=cl)
class(mixed) <- c("psych","mixed")
return(mixed)
}

Try the psych package in your browser

Any scripts or data that you put into this service are public.

psych documentation built on Sept. 26, 2023, 1:06 a.m.