R/utils-import.r

In OpenRepGrid: Tools to Analyze Repertory Grid Data

#//////////////////////////////////////////////////////////////////////////////
### Functions taken from other packages which shall not be loaded due to    ###
### too much overhead or additional dependencies. Hence they are not        ### 
### included as imports.                                                    ###
#//////////////////////////////////////////////////////////////////////////////


# function pointLabel was taken from package maptools. maptools is not imported 
# or mentioned in DESCRIPTION to reduce dependencies as maptools requires sp and 
# gpclib. Thus below is the exact maptools:::pointLabel code.
#
pointLabel <- function (x, y = NULL, labels = seq(along = x), cex = 1, 
                        method = c("SANN", "GA"), allowSmallOverlap = FALSE, 
                        trace = FALSE, doPlot = TRUE, ...) 
{
    if (!missing(y) && (is.character(y) || is.expression(y))) {
        labels <- y
        y <- NULL
    }
    labels <- as.graphicsAnnot(labels)
    boundary <- par()$usr
    xyAspect <- par()$pin[1]/par()$pin[2]
    toUnityCoords <- function(xy) {
        list(x = (xy$x - boundary[1])/(boundary[2] - boundary[1]) * 
            xyAspect, y = (xy$y - boundary[3])/(boundary[4] - 
            boundary[3])/xyAspect)
    }
    toUserCoords <- function(xy) {
        list(x = boundary[1] + xy$x/xyAspect * (boundary[2] - 
            boundary[1]), y = boundary[3] + xy$y * xyAspect * 
            (boundary[4] - boundary[3]))
    }
    z <- xy.coords(x, y, recycle = TRUE)
    z <- toUnityCoords(z)
    x <- z$x
    y <- z$y
    if (length(labels) < length(x)) 
        labels <- rep(labels, length(x))
    method <- match.arg(method)
    if (allowSmallOverlap) 
        nudgeFactor <- 0.02
    n_labels <- length(x)
    width <- (strwidth(labels, units = "figure", cex = cex) + 
        0.015) * xyAspect
    height <- (strheight(labels, units = "figure", cex = cex) + 
        0.015)/xyAspect
    gen_offset <- function(code) c(-1, -1, -1, 0, 0, 1, 1, 1)[code] * 
        (width/2) + (0+1i) * c(-1, 0, 1, -1, 1, -1, 0, 1)[code] * 
        (height/2)
    rect_intersect <- function(xy1, offset1, xy2, offset2) {
        w <- pmin(Re(xy1 + offset1/2), Re(xy2 + offset2/2)) - 
            pmax(Re(xy1 - offset1/2), Re(xy2 - offset2/2))
        h <- pmin(Im(xy1 + offset1/2), Im(xy2 + offset2/2)) - 
            pmax(Im(xy1 - offset1/2), Im(xy2 - offset2/2))
        w[w <= 0] <- 0
        h[h <= 0] <- 0
        w * h
    }
    nudge <- function(offset) {
        doesIntersect <- rect_intersect(xy[rectidx1] + offset[rectidx1], 
            rectv[rectidx1], xy[rectidx2] + offset[rectidx2], 
            rectv[rectidx2]) > 0
        pyth <- abs(xy[rectidx1] + offset[rectidx1] - xy[rectidx2] - 
            offset[rectidx2])/nudgeFactor
        eps <- 1e-10
        for (i in which(doesIntersect & pyth > eps)) {
            idx1 <- rectidx1[i]
            idx2 <- rectidx2[i]
            vect <- (xy[idx1] + offset[idx1] - xy[idx2] - offset[idx2])/pyth[idx1]
            offset[idx1] <- offset[idx1] + vect
            offset[idx2] <- offset[idx2] - vect
        }
        offset
    }
    objective <- function(gene) {
        offset <- gen_offset(gene)
        if (allowSmallOverlap) 
            offset <- nudge(offset)
        if (!is.null(rectidx1)) 
            area <- sum(rect_intersect(xy[rectidx1] + offset[rectidx1], 
                rectv[rectidx1], xy[rectidx2] + offset[rectidx2], 
                rectv[rectidx2]))
        else area <- 0
        n_outside <- sum(Re(xy + offset - rectv/2) < 0 | Re(xy + 
            offset + rectv/2) > xyAspect | Im(xy + offset - rectv/2) < 
            0 | Im(xy + offset + rectv/2) > 1/xyAspect)
        res <- 1000 * area + n_outside
        res
    }
    xy <- x + (0+1i) * y
    rectv <- width + (0+1i) * height
    rectidx1 <- rectidx2 <- array(0, (length(x)^2 - length(x))/2)
    k <- 0
    for (i in 1:length(x)) for (j in seq(len = (i - 1))) {
        k <- k + 1
        rectidx1[k] <- i
        rectidx2[k] <- j
    }
    canIntersect <- rect_intersect(xy[rectidx1], 2 * rectv[rectidx1], 
        xy[rectidx2], 2 * rectv[rectidx2]) > 0
    rectidx1 <- rectidx1[canIntersect]
    rectidx2 <- rectidx2[canIntersect]
    if (trace) 
        cat("possible intersects =", length(rectidx1), "\n")
    if (trace) 
        cat("portion covered =", sum(rect_intersect(xy, rectv, 
            xy, rectv)), "\n")
    GA <- function() {
        n_startgenes <- 1000
        n_bestgenes <- 30
        prob <- 0.2
        mutate <- function(gene) {
            offset <- gen_offset(gene)
            doesIntersect <- rect_intersect(xy[rectidx1] + offset[rectidx1], 
                rectv[rectidx1], xy[rectidx2] + offset[rectidx2], 
                rectv[rectidx2]) > 0
            for (i in which(doesIntersect)) {
                gene[rectidx1[i]] <- sample(1:8, 1)
            }
            for (i in seq(along = gene)) if (runif(1) <= prob) 
                gene[i] <- sample(1:8, 1)
            gene
        }
        crossbreed <- function(g1, g2) ifelse(sample(c(0, 1), 
            length(g1), replace = TRUE) > 0.5, g1, g2)
        genes <- matrix(sample(1:8, n_labels * n_startgenes, 
            replace = TRUE), n_startgenes, n_labels)
        for (i in 1:10) {
            scores <- array(0, NROW(genes))
            for (j in 1:NROW(genes)) scores[j] <- objective(genes[j, 
                ])
            rankings <- order(scores)
            genes <- genes[rankings, ]
            bestgenes <- genes[1:n_bestgenes, ]
            bestscore <- scores[rankings][1]
            if (bestscore == 0) {
                if (trace) 
                  cat("overlap area =", bestscore, "\n")
                break
            }
            genes <- matrix(0, n_bestgenes^2, n_labels)
            for (j in 1:n_bestgenes) for (k in 1:n_bestgenes) genes[n_bestgenes * 
                (j - 1) + k, ] <- mutate(crossbreed(bestgenes[j, 
                ], bestgenes[k, ]))
            genes <- rbind(bestgenes, genes)
            if (trace) 
                cat("overlap area =", bestscore, "\n")
        }
        nx <- Re(xy + gen_offset(bestgenes[1, ]))
        ny <- Im(xy + gen_offset(bestgenes[1, ]))
        list(x = nx, y = ny)
    }
    SANN <- function() {
        gene <- rep(8, n_labels)
        score <- objective(gene)
        bestgene <- gene
        bestscore <- score
        T <- 2.5
        for (i in 1:50) {
            k <- 1
            for (j in 1:50) {
                newgene <- gene
                newgene[sample(1:n_labels, 1)] <- sample(1:8, 
                  1)
                newscore <- objective(newgene)
                if (newscore <= score || runif(1) < exp((score - 
                  newscore)/T)) {
                  k <- k + 1
                  score <- newscore
                  gene <- newgene
                }
                if (score <= bestscore) {
                  bestscore <- score
                  bestgene <- gene
                }
                if (bestscore == 0 || k == 10) 
                  break
            }
            if (bestscore == 0) 
                break
            if (trace) 
                cat("overlap area =", bestscore, "\n")
            T <- 0.9 * T
        }
        if (trace) 
            cat("overlap area =", bestscore, "\n")
        nx <- Re(xy + gen_offset(bestgene))
        ny <- Im(xy + gen_offset(bestgene))
        list(x = nx, y = ny)
    }
    if (method == "SANN") 
        xy <- SANN()
    else xy <- GA()
    xy <- toUserCoords(xy)
    if (doPlot) 
        text(xy, labels, cex = cex, ...)
    invisible(xy)
}


# Interleave Rows of data frames or matrices
# function from package gdata written by Gregory R. Warnes #
#
interleave <- function (..., append.source = TRUE, sep = ": ", drop = FALSE) 
{
    sources <- list(...)
    sources[sapply(sources, is.null)] <- NULL
    sources <- lapply(sources, function(x) if (is.matrix(x) || 
        is.data.frame(x)) 
        x
    else t(x))
    nrows <- sapply(sources, nrow)
    mrows <- max(nrows)
    if (any(nrows != mrows & nrows != 1)) 
        stop("Arguments have differening numbers of rows.")
    sources <- lapply(sources, function(x) if (nrow(x) == 1) 
        x[rep(1, mrows), , drop = drop]
    else x)
    tmp <- do.call("rbind", sources)
    nsources <- length(sources)
    indexes <- outer((0:(nsources - 1)) * mrows, 1:mrows, "+")
    retval <- tmp[indexes, , drop = drop]
    if (append.source && !is.null(names(sources))) 
        if (!is.null(row.names(tmp))) 
            row.names(retval) <- paste(format(row.names(retval)), 
                format(names(sources)), sep = sep)
        else row.names(retval) <- rep(names(sources), mrows)
    retval
}


# function errbar form Hmisc package by Frank E Harrell Jr.
#
errbar <- function (x, y, yplus, yminus, cap = 0.015, main = NULL, sub = NULL, 
    xlab = as.character(substitute(x)), ylab = if (is.factor(x) || 
        is.character(x)) "" else as.character(substitute(y)), 
    add = FALSE, lty = 1, type = "p", ylim = NULL, lwd = 1, pch = 16, 
    Type = rep(1, length(y)), ...) 
{
    if (is.null(ylim)) 
        ylim <- range(y[Type == 1], yplus[Type == 1], yminus[Type == 
            1], na.rm = TRUE)
    if (is.factor(x) || is.character(x)) {
        x <- as.character(x)
        n <- length(x)
        t1 <- Type == 1
        t2 <- Type == 2
        n1 <- sum(t1)
        n2 <- sum(t2)
        omai <- par("mai")
        mai <- omai
        mai[2] <- max(strwidth(x, "inches")) + 0.25 * TRUE #.R.
        par(mai = mai)
        on.exit(par(mai = omai))
        plot(NA, NA, xlab = ylab, ylab = "", xlim = ylim, ylim = c(1, 
            n + 1), axes = FALSE, ...)
        axis(1)
        w <- if (any(t2)) 
            n1 + (1:n2) + 1
        else numeric(0)
        axis(2, at = c(seq.int(length.out = n1), w), labels = c(x[t1], 
            x[t2]), las = 1, adj = 1)
        points(y[t1], seq.int(length.out = n1), pch = pch, type = type, 
            ...)
        segments(yplus[t1], seq.int(length.out = n1), yminus[t1], 
            seq.int(length.out = n1), lwd = lwd, lty = lty)
        if (any(Type == 2)) {
            abline(h = n1 + 1, lty = 2, ...)
            offset <- mean(y[t1]) - mean(y[t2])
            if (min(yminus[t2]) < 0 & max(yplus[t2]) > 0) 
                lines(c(0, 0) + offset, c(n1 + 1, par("usr")[4]), 
                  lty = 2, ...)
            points(y[t2] + offset, w, pch = pch, type = type, 
                ...)
            segments(yminus[t2] + offset, w, yplus[t2] + offset, 
                w, lwd = lwd, lty = lty)
            at <- pretty(range(y[t2], yplus[t2], yminus[t2]))
            axis(side = 3, at = at + offset, labels = format(round(at, 
                6)))
        }
        return(invisible())
    }
    if (add) 
        points(x, y, pch = pch, type = type, ...)
    else plot(x, y, ylim = ylim, xlab = xlab, ylab = ylab, pch = pch, 
        type = type, ...)
    xcoord <- par()$usr[1:2]
    smidge <- cap * (xcoord[2] - xcoord[1])/2
    segments(x, yminus, x, yplus, lty = lty, lwd = lwd)
    if (par()$xlog) {
        xstart <- x * 10^(-smidge)
        xend <- x * 10^(smidge)
    }
    else {
        xstart <- x - smidge
        xend <- x + smidge
    }
    segments(xstart, yminus, xend, yminus, lwd = lwd, lty = lty)
    segments(xstart, yplus, xend, yplus, lwd = lwd, lty = lty)
    return(invisible())
}


### PMC to Fisher's Z and back from package psych written by William Revelle ###
# only needed if psych is removed from COLLATE
#
# fisherz <- function (rho) 
# {
#     0.5 * log((1 + rho)/(1 - rho))
# }
# 
# fisherz2r <- function (z) 
# {
#     (exp(2 * z) - 1)/(1 + exp(2 * z))
# }




#//////////////////////////////////////////////////////////////////////////////
## Optimal Box-Cox transformation according to
## a grid-based maximization of the correlation
## of a Normal P-P plot.
#//////////////////////////////////////////////////////////////////////////////
## Author: Ioannis Kosmidis
## Email: <I.Kosmidis at warwick dot ac dot uk>
## Latest release: 02/08/2008
## Distributed under GPL 2 or greater:
## Available at http://www.gnu.org/licenses
#//////////////////////////////////////////////////////////////////////////////
## "normal.ppplot"
## 
## Arguments:
## x: a vector of the observed values
## plot: values TRUE/FALSE depending on whether the Normal P-P
##       plot is to be shown or if only the correlation
##       coefficient of the plot should be returned
## Value:
## Either a Normal P-P plot or the correlation of the Normal
## P-P plot, depending on the values of the plot argument.
#//////////////////////////////////////////////////////////////////////////////
## "optimal.boxcox"
## 
## Arguments:
## x: a vector of the observed values
## lambda: the grid of lambda values that should be considered
##         for the grid maximization
## Value:
## A plot showing the Normal P-P plot correlations for the
## grid of values of lambda considered and the Normal P-P plot 
## for the value of lambda which resulted in higher
## correlation.
## The vector of the transformed observations is also
## returned.
#//////////////////////////////////////////////////////////////////////////////
                         
normal.ppplot <- function(x, plot=FALSE) {
  standardized.x <- (x - mean(x))/sd(x)
  obs.probs <- pnorm(sort(standardized.x))
  theor.probs <- ppoints(length(x))
  corrs.pp <- cor(obs.probs, theor.probs)
  if (plot) {
    plot(obs.probs, theor.probs,
         xlab = "Probabilities based on the standardized observed vector",
         ylab = "Theoretical probabilities")
    title(main = expression("Normal distribution P-P Plot"),
          sub = paste("Normal P-P plot correlation coefficient:",
            round(corrs.pp, 3)))
    abline(0, 1)    
  }
  else corrs.pp
}


optimal.boxcox <- function(x, lambda = seq(-2, 2, len=200), plot=FALSE) {
  ll <- length(lambda)
  correlations.pp <- numeric(ll)
  for (j in 1:ll) {
    if (lambda[j]==0) temp <- log(x)
    else temp <- (x^lambda[j]-1)/lambda[j]
    correlations.pp[j] <- normal.ppplot(temp, plot = FALSE)
  }
  m.ind <- which.max(correlations.pp)
  lambda.max <- lambda[m.ind]
  if (plot){
    par(mfrow = c(1,2))
    plot(lambda, correlations.pp, type='l',
         ylim=c(min(correlations.pp), 1.1),
         xlab = expression(lambda),
         ylab = "Normal P-P plot correlation coefficient")
    points(lambda.max, correlations.pp[m.ind], pch="+")
    text(lambda.max, correlations.pp[m.ind] + 0.05,
         bquote(lambda == .(lambda.max)))
    title(expression(paste(lambda, " versus P-P plot correlations")))
  }
  power.transformed <- (x^lambda.max - 1)/lambda.max
  normal.ppplot(power.transformed, plot = plot)
  list(x=power.transformed, lambda=lambda.max)
}