R/Yule.R

In psych: Procedures for Psychological, Psychometric, and Personality Research

"Yule" <- 
 function(x,Y=FALSE) {
  # if(!is.matrix(x)) {stop("x must be a matrix")}
   stopifnot(prod(dim(x)) == 4 || length(x) == 4)
    if (is.vector(x)) {
     x <- matrix(x, 2)}
   a  <- x[1,1]
   b <- x[1,2]
   c <- x[2,1]
   d <- x[2,2]
  if (Y) {Yule <- (sqrt(a*d) - sqrt(b*c))/(sqrt(a*d)+sqrt(b*c))} else {Yule <- (a*d- b*c)/(a*d + b*c)}
   return(Yule) #Yule Q
   }
   
 
"Yule.inv" <- 
    function(Q,m,n=NULL) {#find the cells for a particular Yule value with fixed marginals
     if (length(m) > 2) { #this old way is for one correlation at a time with all 4 marginals - 
     R1 <- m[1]
     R2 <- m[2]
     C1 <- m[3]
     C2 <- m[4] } else {   #the better way is to specify just two marginals- allows better interface with matrices
       if(is.null(n)) {n<- 1}  #do it as percentages
        
        R1 <- m[1] 
        R2 <- n -R1
        C1 <- m[2]
        C2 <- n - C1}
     
     f <- function(x) {(Q- Yule(c(x,R1-x,C1-x,R2-C1+x)))^2}
     xval <- optimize(f,c(min(R1-min(C2,R1),C1-min(R2,C1)),min(R1,C1)))
     R <- matrix(ncol=2,nrow=2)
     x <- xval$minimum
     R[1,1] <- x
     R[1,2] <- R1-x
     R[2,1] <- C1 -x
     R[2,2] <- R2-C1 + x
     return(R)
     }
     
"Yule2poly" <-
    function(Q,m,n=NULL,correct=TRUE) { #find the phi equivalent to a Yule Q with fixed marginals
    .Deprecated("Yule2tet",msg="Yule2poly has been replaced by Yule2tet, please try again") 
    t <- Yule.inv(Q,m,n=n)
    r <- tetrachoric(t,correct=correct)$rho
    return(r)
    }

"Yule2tet" <-
    function(Q,m,n=NULL,correct=TRUE) { #find the phi equivalent to a Yule Q with fixed marginals
    t <- Yule.inv(Q,m,n=n)
    r <- tetrachoric(t,correct=correct)$rho
    return(r)
    }
    
"Yule2phi" <- 
   function(Q,m,n=NULL) { #find the phi equivalent to a Yule Q with fixed marginals
    t <- Yule.inv(Q,m,n=n)
  return(phi(t,digits=15))}
  
  
  
"Yule2tetra" <-
function(Q,m,n=NULL,correct=TRUE) {
if(!is.matrix(Q) && !is.data.frame(Q)) {result <- Yule2tet(Q,c(m[1],m[2]),n=n,correct=correct) } else {
nvar <- nrow(Q)
if(nvar !=ncol(Q)) {stop('Matrix must be square')}
if (length(m) !=nvar) {stop("length of m must match the number of variables")}
result <- Q
for(i in 2:nvar) {
  for (j in 1:(i-1)) { 
  result[i,j] <- result[j,i] <- Yule2tet(Q[i,j],c(m[i],m[j]),n,correct=correct)
  }
}}
return(result) }

"YuleBonett" <- 
 function(x,c=1,bonett=FALSE,alpha=.05) {
  # if(!is.matrix(x)) {stop("x must be a matrix")}
   stopifnot(prod(dim(x)) == 4 || length(x) == 4)
    if (is.vector(x)) {
     x <- matrix(x, 2)}
     p <- x/sum(x)
    C <- c
   a  <- p[1,1]
   b <- p[1,2]
   c <- p[2,1]
   d <- p[2,2]
   Rs <- rowSums(p)
   Cs <- colSums(p)
   if(bonett) {C <- .5 - (.5 - min(Rs,Cs))^2 }
   ad <- (a*d)^C
   bc <- (b*c)^C
   Yule <- (ad-bc)/(ad+bc) 
   Ystar <- 
   #See Bonett 2007 p 434
   w <- (x[1,1] + .1)* (x[2,2] + .1)/((x[2,1] + .1)* (x[1,2] + .1)) #OR estimate
    Ystar <- (w^C -1)/(w^C+1)   #this is the small cell size adjusted value
   	vlQ <- (1/(x[1,1] + .1) + 1/ (x[2,2] + .1) + 1/(x[2,1] + .1)+ 1/ (x[1,2] + .1))
	vQ <-  (C^2/4)*(1-Ystar^2)^2 *vlQ   #equation 9  
  tanhinv <- atanh(Ystar)
  upper <- tanh(atanh(Ystar) + qnorm(1-alpha/2) *sqrt( (C^2/4) *vlQ))
   lower <- tanh(atanh(Ystar) - qnorm(1-alpha/2) *sqrt( (C^2/4) *vlQ))
   result <- list(rho=Ystar,se=sqrt(vQ), upper=upper, lower=lower)
   class(result) <- c("psych","Yule")
   return(result) #Yule Q, Y, or generalized Bonett
   }
   
   
   
#Find a correlation matrix of Yule correlations
"YuleCor" <- function(x,c=1,bonett=FALSE,alpha=.05) {
cl <- match.call()
nvar <- ncol(x)
rho <- matrix(NA,nvar,nvar)
ci <-  matrix(NA,nvar,nvar)
zp <- matrix(NA,nvar,nvar)
for(i in 1:nvar) {
  for(j in 1:i) {
    YB <- YuleBonett(table(x[,i],x[,j]),c=c,bonett=bonett,alpha=alpha)
    rho[i,j] <- rho[j,i] <- YB$rho
    ci[i,j] <- YB$lower
    ci[j,i] <- YB$upper
    zp[i,j] <- YB$rho/YB$se
    zp[j,i] <- 1-pnorm(zp[i,j])
    
    }}
colnames(rho) <- rownames(rho) <- colnames(ci) <- rownames(ci) <- colnames(x)
result <- list(rho=rho,ci=ci,zp=zp,Call=cl)
class(result) <- c("psych","yule")
return(result)
}