R/stats.by.r

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

"error.crosses.by" <-
function (x,y,z,labels=NULL,main=NULL,xlim=NULL,ylim= NULL,xlab=NULL,ylab=NULL,pos=NULL,offset=1,arrow.len=.2,alpha=.05,sd=FALSE,...)  # x  and y are data frame or descriptive stats
    {if(is.null(x$mean)) {x <- describe.by(x,z,mat=TRUE)
           }
     xmin <- min(x$mean)
     xmax <- max(x$mean)
     if(sd) {max.sex <- max(x$sd,na.rm=TRUE)
                      if(is.null(xlim))  {xlim=c(xmin - max.sex,xmax + max.sex) }}  else {max.sex <- max(x$se,na.rm=TRUE)}       
     if(is.null(y$mean)) {y <- describe(y)}
     ymin <- min(y$mean)
     ymax <- max(y$mean)
     if(sd) {max.sey <- max(y$sd,na.rm=TRUE)
              if(is.null(ylim))  {ylim=c(ymin - max.sey,ymax +max.sey)}} else {   max.sey <- max(y$se,na.rm=TRUE)  } 
     
     if(is.null(xlim))  xlim=c(xmin - 2*max.sex,xmax +2*max.sex)
     if(is.null(ylim))  ylim=c(ymin - 2*max.sey,ymax +2*max.sey)
     
     if(is.null(main)) {if(!sd) { main = paste((1-alpha)*100,"% confidence limits",sep="") } else {main= paste("Means and standard deviations")} }
     if(is.null(xlab)) xlab <- "Group 1"
     if(is.null(ylab)) ylab <- "Group 2"
     plot(x$mean,y$mean,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main,...)
     
    cix <- qt(1-alpha/2,x$n)
    ciy <- qt(1-alpha/2,y$n)
     z <- dim(x)[1]
    if(sd) {x$se <- x$sd
            y$se <- y$sd
            cix <- ciy <- rep(1,z)
           }
    
     if (is.null(pos)) {locate <- rep(1,z)} else {locate <- pos}
     if (is.null(labels))  {labels <- rownames(x)}
    if (is.null(labels))  {lab <- paste("V",1:z,sep="")}  else {lab <-labels}
    
        for (i in 1:z)  
    	{xcen <- x$mean[i]
    	 ycen <- y$mean[i]
    	 xse  <- x$se[i]
    	 yse <-  y$se[i]
    	 arrows(xcen-cix[i]* xse,ycen,xcen+ cix[i]* xse,ycen,length=arrow.len, angle = 90, code=3,col = par("fg"), lty = NULL, lwd = par("lwd"), xpd = NULL)
    	 arrows(xcen,ycen-ciy[i]* yse,xcen,ycen+ ciy[i]*yse,length=arrow.len, angle = 90, code=3,col = par("fg"), lty = NULL, lwd = par("lwd"), xpd = NULL)
    	text(xcen,ycen,labels=lab[i],pos=locate[i],cex=1,offset=offset)     #puts in labels for all points
    	}	
   }
   
   
   
   "ellipse" <-    function (x,y,r1,r2,...) { 
#code adapted from John Fox
    segments=51
    angles <- (0:segments) * 2 * pi/segments
    unit.circle <- cbind(cos(angles), sin(angles))
   
    xs <- r1
    #ys <- e.size * yrange
    ellipse <- unit.circle 
    ellipse[,1] <- ellipse[,1]*r1 + x
    ellipse[,2] <- ellipse[,2]*r2+ y  #ys?
    lines(ellipse, ...)
    return(xs)
}
#    modified 18/1/15 to pass just the xvar and yvar to statsBy
#    "errorCircles" <-
# function (x,y,data,ydata=NULL,group=NULL,paired=FALSE, labels=NULL,main=NULL,xlim=NULL,ylim= NULL,xlab=NULL,ylab=NULL,add=FALSE,pos=NULL,offset=1,arrow.len=.2,alpha=.05,sd=FALSE,bars=TRUE,circles=TRUE,...) { # x  and y are data frame or descriptive stats
#      
#      xvar <- x
#      yvar <- y
#      if(!is.null(group)) {data <- statsBy(data[,c(xvar,yvar,group)],group)}
#     x <- list()
#      if(paired) {
#           	x$mean <- t(data$mean[,xvar])
#      		x$sd <- t(data$sd[,xvar])
#      		x$n <- t(data$n[,xvar]) 
#      	} else {  #the normal case
#     	 x$mean <- data$mean[,xvar]
#      	x$sd <- data$sd[,xvar]
#      	x$n <- data$n[,xvar]}  
#      
#      xmin <- min(x$mean,na.rm=TRUE)
#      xmax <- max(x$mean,na.rm=TRUE)
#      x$se <- x$sd/sqrt(x$n)
#      
#      if(sd) {max.sex <- max(x$sd,na.rm=TRUE)
#                       if(is.null(xlim))  {xlim=c(xmin - max.sex,xmax + max.sex) }}  else {max.sex <- max(x$se,na.rm=TRUE)}       
#      
#      y <- list()
#      if(!is.null(ydata)) {
#           	y$mean <- ydata$mean[,yvar]
#          	y$sd <- ydata$sd[,yvar]
#         	y$n <- ydata$n[,yvar]
#         	} else {
#        	 	y$mean <- data$mean[,yvar]
#     		 y$sd <- data$sd[,yvar]
#     	 	y$n <- data$n[,yvar]}
#      
#      ymin <- min(y$mean,na.rm=TRUE)
#      ymax <- max(y$mean,na.rm=TRUE)
#      y$se <- y$sd/sqrt(y$n)
#      if(sd) {max.sey <- max(y$sd,na.rm=TRUE)
#               if(is.null(ylim))  {ylim=c(ymin - max.sey,ymax +max.sey)}} else {   max.sey <- max(y$se,na.rm=TRUE)  } 
#      
#      if(is.null(xlim))  xlim=c(xmin - 2*max.sex,xmax +2*max.sex)
#      if(is.null(ylim))  ylim=c(ymin - 2*max.sey,ymax +2*max.sey)
#      
#      if(is.null(main)) {if(!sd) { main = paste((1-alpha)*100,"% confidence limits",sep="") } else {main= paste("Means and standard deviations")} }
#      if(paired) {if(is.null(xlab)) xlab <- "Group 1"
#             if(is.null(ylab)) ylab <- "Group 2"
#         }  else {
#         if(is.null(xlab)) xlab <- colnames(data$mean)[xvar]
#         if(is.null(ylab)) ylab <- colnames(data$mean)[yvar]
#      }
#      if(add)  {  
#       if(paired) {points(x$mean,typ="p",...) } else {points(x$mean,y$mean,typ="p",...)}
#     } else {
#      if(paired) {plot(x$mean,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main,typ="p",...) } else {plot(x$mean,y$mean,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main,typ="p",...)}
#     }
#      N <-x$n
#      Nmax <- max(N)
#     cix <- qt(1-alpha/2,x$n)
#     ciy <- qt(1-alpha/2,y$n)
#      if(paired) {z <- nrow(x$mean) } else {z <- length(x$mean)}
#     if(sd) {x$se <- x$sd
#             y$se <- y$sd
#             cix <- ciy <- rep(1,z)
#            }
#     
#      if (is.null(pos)) {locate <- rep(1,z)} else {locate <- pos}
#      if (is.null(labels))  {labels <- rownames(x$mean)}
#     if (is.null(labels))  {lab <- paste("V",1:z,sep="")}  else {lab <-labels}
#     
#         for (i in 1:z)  
#     	{ if(paired) { xcen <- x$mean[i,1]
#     	 ycen <- x$mean[i,2]
#     	 xse  <- x$se[i,1]
#     	 yse <-  x$se[i,2]
#     	  } else {
#     	xcen <- x$mean[i]
#     	 ycen <- y$mean[i]
#     	 xse  <- x$se[i]
#     	 yse <-  y$se[i]}
#     if(bars) {if(max(x$se,na.rm=TRUE) > 0) 	 arrows(xcen-cix[i]* xse,ycen,xcen+ cix[i]* xse,ycen,length=arrow.len, angle = 90, code=3,col = par("fg"), lty = NULL, lwd = par("lwd"), xpd = NULL)
#     	if(max(y$se,na.rm=TRUE) >0 )  arrows(xcen,ycen-ciy[i]* yse,xcen,ycen+ ciy[i]*yse,length=arrow.len, angle = 90, code=3,col = par("fg"), lty = NULL, lwd = par("lwd"), xpd = NULL)
#     	 }
#     	
#     	text(xcen,ycen,labels=lab[i],pos=locate[i],cex=1,offset=offset)     #puts in labels for all points
#     if(circles) { xrange <- xlim[2] - xlim[1]
#                 yrange <- ylim[2] - ylim[1]
#                 xscale <-max(x$se)
#                 yscale <-max(y$se)
#     ellipse(xcen,ycen,sqrt(xscale*x$n[i]/Nmax),sqrt( yscale*x$n[i]/Nmax))
#     	}
#     	}	
#     if(!is.null(group)) return(invisible(data))
#    }
#    
#    
   
#    "statsBy.old" <-
#    function (data,group,cors=FALSE) { #  
#   valid <- function(x) { #count the number of valid cases 
#         sum(!is.na(x))
#     }
#        gr <- which(colnames(data) == group)
#        
#        z1 <- data[,group]
#        z <- z1
#        cnames <- colnames(data)
#        for (i in 1:ncol(data)) {if(is.factor(data[,i]) || is.logical(data[,i])) {
#              data[,i] <- as.numeric(data[,i])
#             # colnames(data)[i] <- paste(cnames[i],"*",sep="")
#              }}
#        xvals <- list()
#        
#                xvals$mean <- t(matrix(unlist(by(data,z,colMeans,na.rm=TRUE)),nrow=ncol(data)))              
#                xvals$sd <-t(matrix(unlist(by(data,z,function(x) sapply(x,sd,na.rm=TRUE))),nrow=ncol(data)))
#                xvals$n <- t(matrix(unlist(by(data,z,function(x) sapply(x,valid))),nrow=ncol(data)))
#                colnames(xvals$mean) <- colnames(xvals$sd) <- colnames(xvals$n) <-  colnames(data)
#                rownames(xvals$mean) <-  rownames(xvals$sd) <- rownames(xvals$n) <- levels(z)
#                nH <- harmonic.mean(xvals$n)
#                GM <- colSums(xvals$mean*xvals$n)/colSums(xvals$n) 
#                MSb <- colSums(xvals$n*t((t(xvals$mean) - GM)^2))/(nrow(xvals$mean)-1) #weight means by n
#                MSw <- colSums(xvals$sd^2*(xvals$n-1))/(colSums(xvals$n-1)) #find the pooled sd
#                xvals$F <- MSb/MSw
#                N <- colSums(xvals$n)
#               # npr <-(N^2 - colSums(xvals$n))/(N *(nrow(xvals$n) -1))
#               # npr <- harmonic.mean(xvals$n-1)
#               npr <- (colSums(xvals$n-1)+nrow(xvals$n))/(nrow(xvals$n))
#                xvals$ICC1 <- (MSb-MSw)/(MSb + MSw*(npr-1))
#                xvals$ICC2 <- (MSb-MSw)/(MSb)
#                              if(cors) { r <- by(data,z,function(x) cor(x[-1],use="pairwise"))
#               nvars <-  ncol(r[[1]])
#               xvals$r <- r
#               lower <- lapply(r,function(x) x[lower.tri(x)])
#              xvals$within <- t(matrix(unlist(lower),nrow=nvars*(nvars-1)/2))
#              wt <- by(data,z,function(x) pairwiseCount(x[-1]))
#              lower.wt <- t(matrix(unlist(lapply(wt,function(x) x[lower.tri(x)])    )  ,nrow=nvars*(nvars-1)/2))
#              lower.wt <- t(t(lower.wt)/colSums(lower.wt,na.rm=TRUE))
#              pool  <- colSums( lower.wt * xvals$within,na.rm=TRUE)
#              pool.sd <- apply(xvals$within, 2,FUN=sd, na.rm=TRUE)
#              xvals$pooled <- matrix(NaN,nvars,nvars)
#              xvals$pooled[lower.tri(xvals$pooled)] <- pool
#              xvals$pooled[upper.tri(xvals$pooled)]  <- pool
#              diag(xvals$pooled) <- 1
#              xvals$sd.r <-  matrix(NaN,nvars,nvars)
#              xvals$sd.r[lower.tri(xvals$sd.r)] <- pool.sd
#              xvals$sd.r[upper.tri(xvals$sd.r)] <- pool.sd
#              colnames(xvals$pooled) <- rownames (xvals$pooled) <- cnames[-1]
#               }
# 
#               nvar <- ncol(data)-1
#                xvals$raw <- cor(data,use="pairwise")
#              new.data <- as.matrix( merge(xvals$mean,data,by=group,suffixes =c(".bg",""))[-1])
#              
#              diffs <- new.data[,(nvar+1):ncol(new.data)] - new.data[,1:nvar]
#              colnames(diff) <- rownames(diff) <- paste(colnames(diff),".wg",sep="")
#              xvals$rbg <- cor(new.data[,1:nvar],use="pairwise")  #the between group (means)
#              xvals$rwg <- cor(diffs,use="pairwise")  #the within group (differences)
#              colnames(xvals$rwg) <- rownames(xvals$rwg) <- paste(colnames(xvals$rwg),".wg",sep="")
#              xvals$etabg <- cor(new.data[,1:nvar],new.data[,(nvar+1):ncol(new.data)],use="pairwise") #the means with the data
#              xvals$etawg <- cor(new.data[,(nvar+1):ncol(new.data)],diffs,use="pairwise") #the deviations and the data
#             rownames(xvals$etawg)  <- paste(rownames(xvals$etawg),".wg",sep="")
#             
#     return(xvals)
#     }
# 

"cor.wt" <- function(data,vars=NULL, w=NULL,sds=NULL, cor=TRUE) {
   cl <- match.call()
   if(is.list(data) && !is.data.frame(data)) {w <- data$n   #use the output from statsBy
                   sds <- data$sd
                   x <- data$mean} else {x <- data}
    if(!is.null(vars)) {x <- x[,vars] 
                        w <- w[,vars]
                        sds <- sds[,vars] }              
   if(is.null(w)) w <- matrix(rep(rep(1/nrow(x),nrow(x)),ncol(x)),nrow=nrow(x),ncol=ncol(x))
   if(is.null(ncol(w))) {wt <- w/sum(w) } else {
   wt <- t(t(w)/colSums(w))}
   cnames <- colnames(x)
    for (i in 1:ncol(x)) {if(is.factor(x[,i]) || is.logical(x[,i])) {
             x[,i] <- as.numeric(x[,i])
             colnames(x)[i] <- paste(cnames[i],"*",sep="")
             }}
   means <- colSums(x * wt,na.rm=TRUE) 
   xc  <-  scale(x,center=means,scale=FALSE)  #these are the weighted centered data
   if(is.null(sds)) {xs <- xc /sqrt(w) } else {xs  <- xc * sds/sqrt(w)}
   xwt <- sqrt(wt) * xc
   #  added February 20, 2016 to consider missing values
   if(any(is.na(xwt))) {
   cov <- apply(xwt,2, function(x) colSums(xwt * x,na.rm=TRUE))
   
   }  else {#matrix algebra without matrices
       cov <- crossprod(xwt) }            #/(1-colSums(wt^2,na.rm=TRUE))
    if(cor) {r <- cov2cor(cov)} else {r <- cov}
    xw <- wt * xc
    result <-list(r=r,xwt = xwt,wt=wt,mean=means,xc=xc,xs=xs) 
    result$Call <- cl
    class(result) <- c("psych","cor.wt")
return(result)}