Gems.R
From Organic Design wiki
"BiplotC" <-
structure(function(xx, comps = c(1, 2), cex = c(.8, .4), pchs = 0:6,
colour.c = c("black"), colour.v = c("blue"), scale = 1,
clust = 1, no.clust = 3, lab.cex = 1.2, ax.cex = 0.9,
xylab = "PC", expand = 1, xlim = NULL, ylim = NULL,
main = NULL, xlab = NULL, ylab = NULL)
{
## Purpose: Adds more options to BiplotB
## ----------------------------------------------------------------------
## Modified from: BiplotB
## ----------------------------------------------------------------------
## Arguments:
## colour.c: colour vector for cases
## colour.v: colour vector for variables
## clust: are we clustering cases (1) or variables (2)?
## If no clustering is wanted, use BiplotB2
## no.clust number of clusters
## xylab: Prefix to x- and y-labels
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 16 Aug 2006, 09:49
if(class(xx) == "princomp"){
scores <- xx$scores[, paste("Comp", comps, sep = ".")]
loadings <- xx$loadings[, paste("Comp", comps, sep = ".")]
}
if(class(xx) == "prcomp"){
scores <- xx$x[, paste("PC", comps, sep = "")]
loadings <- xx$rotation[, paste("PC", comps, sep = "")]
}
importance <- round(100*xx$sdev/max(cumsum(xx$sdev)), 1) # Variation attributed
import.lab <- importance[comps]
x <- scores
y <- loadings
n <- nrow(scores)
p <- nrow(loadings)
## Copy in lam stuff from the default (omit n.obs stuff)
lam <- xx$sdev[comps]
lam <- lam * sqrt(n)
if (scale < 0 || scale > 1)
warning("'scale' is outside [0, 1]")
if (scale != 0)
lam <- lam^scale
else lam <- 1
x <- t(t(scores)/lam)
y <- t(t(loadings)*lam)
unsigned.range <- function(x) c(-abs(min(x)), abs(max(x)))
rangx1 <- unsigned.range(x[, 1])
rangx2 <- unsigned.range(x[, 2])
rangy1 <- unsigned.range(y[, 1])
rangy2 <- unsigned.range(y[, 2])
### fix up xlim and ylim
if (missing(xlim) && missing(ylim))
xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2)
else if (missing(xlim))
xlim <- rangx1
else if (missing(ylim))
ylim <- rangx2
### Ratio of variable sclae to cases scale
ratio <- max(rangy1/rangx1, rangy2/rangx2)/expand
blank.plot(range(x), range(x), xlim = xlim, ylim = ylim)
## work out scaling factor for the loadings
xax.tik <- pretty(range(xlim))
xax.tik.at <- xax.tik
yax.tik <- pretty(range(ylim))
yax.tik.at <- yax.tik
box()
abline(v = 0, col = "grey")
abline(h = 0, col = "grey")
### Decide on what clustering is to happen
if(clust == 1 | clust == "case"){ # We're clustering cases
if(length(colour.c) < no.clust)
colour.pal.c <- brewer.pal(8, "Dark2")[1:no.clust]
else colour.pal.c <- colour.c
case.cluster <- pam(x, no.clust)
colour.c <- colour.pal.c[case.cluster$clustering]
}
if(clust == 2 | clust == "variable"){ # We're clustering variables
if(length(colour.v) < no.clust)
colour.pal.v <- brewer.pal(8, "Dark2")[1:no.clust]
else colour.pal.v <- colour.v
variable.cluster <- pam(y, no.clust)
colour.v <- colour.pal.v[variable.cluster$clustering]
}
items <- rownames(scores) # can stay non-generic for now
## Get same letters into groups
## regions <- gsub("[0-9]", "", items)
text(x[,1], x[,2], items, cex = cex[1], col = colour.c)
axis(2, cex.axis = ax.cex, mgp = c(1, .8, 0), at = yax.tik.at,
labels = yax.tik)
axis(1, cex.axis = ax.cex, mgp = c(1, .8, 0), at = xax.tik.at,
labels = xax.tik)
if(is.null(xlab))
xlab <- ppaste(xylab, comps[1], " (",import.lab[1], "%)")
if(is.null(ylab))
ylab <- ppaste(xylab, comps[2], " (",import.lab[2], "%)")
mtext(xlab, side = 1, line = 2.2, cex = lab.cex)
mtext(ylab, side = 2, line = 2.5, cex = lab.cex)
text(y[, 1]/ratio, y[, 2]/ratio,
rownames(loadings), cex = cex[2], col = colour.v)
}
, comment = "18/08/2006")
"JIM" <-
structure(function()
{
## Purpose: Loads Jim Lindsay's libraries for gnlm
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 4 Feb 2002, 07:57
library(rmutil, lib.loc = "/home/hrapgc/Rstuff/library/")
library(gnlm, lib.loc = "/home/hrapgc/Rstuff/library/")
}
, comment = "01/07/2002")
"MAP" <-
structure(function()
{
### Installs maps library
library(maps, lib.loc = "/home/pat/Rstuff/library/")
}
, comment = "06/06/2001")
"MASS" <-
structure(function()
{
# Installs MASS library
library(MASS, lib.loc = "/home/hrapgc/Rstuff/library")
}
, comment = "03/10/2001")
"PLS" <-
structure(function()
{
### Installs Partial Least Squares library
library(pls, lib.loc = "/home/hrapgc/patrick/library")
}
, comment = "09/05/2001")
"Plot.rpart" <-
structure(function(x, uniform=FALSE, branch=1, compress=FALSE,
nspace, margin=0, minbranch=.3, col = "turquoise",...)
{
## Purpose: Allows different colours for lines in plot.rpart
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 26 Jan 2005, 09:54
if(!inherits(x, "rpart"))
stop("Not an rpart object")
if (!is.null(x$frame$splits)) x <- rpconvert(x) #help for old objects
if (compress & missing(nspace)) nspace <- branch
if (!compress) nspace <- -1 #means no compression
dev <- dev.cur()
if (dev == 1) dev <- 2
assign(paste(".rpart.parms", dev, sep = "."),
list(uniform=uniform, branch=branch, nspace=nspace,
minbranch=minbranch), envir=.GlobalEnv)
#define the plot region
temp <- rpart:::rpartco(x)
xx <- temp$x
yy <- temp$y
temp1 <- range(xx) + diff(range(xx))*c(-margin, margin)
temp2 <- range(yy) + diff(range(yy))*c(-margin, margin)
plot(temp1, temp2, type='n', axes=FALSE, xlab='', ylab='', ...)
# Draw a series of horseshoes or V's, left son, up, down to right son
# NA's in the vector cause lines() to "lift the pen"
node <- as.numeric(row.names(x$frame))
temp <- rpart:::rpart.branch(xx, yy, node, branch)
if (branch>0) text(xx[1], yy[1], '|')
lines(c(temp$x), c(temp$y), col = col)
invisible(list(x=xx, y=yy))
}
, comment = "26/01/2005")
"SURV" <-
structure(function() library(survival5, lib.loc = "/home/hrapgc/patrick/library")
, comment = "27/04/2001")
"Text.rpart" <-
structure(function(x, splits = TRUE, label = "yval", FUN = text, all=FALSE,
pretty = NULL, digits = getOption("digits") - 3,
use.n=FALSE, fancy=FALSE, fwidth=.8, fheight =.8,
leaf.col = "darkorange",...)
{
## Purpose: Modifies text.rpart to position text better
## ----------------------------------------------------------------------
## Modified from: text.rpart in the basic distribution
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 21 Jan 2005, 16:39
if(!inherits(x, "rpart")) stop("Not legitimate rpart")
if(!is.null(x$frame$splits)) x <- rpconvert(x)#Backwards compatability
frame <- x$frame
col <- names(frame)
ylevels <- attr(x,'ylevels')
if(!is.null(ylevels <- attr(x, "ylevels")))
col <- c(col, ylevels)
if(is.na(match(label, col)))
stop("Label must be a column label of the frame component of the tree"
)
cxy <- par("cxy") #character width and height
if(!is.null(srt <- list(...)$srt) && srt == 90)
cxy <- rev(cxy)
xy <- rpart:::rpartco(x)
node <- as.numeric(row.names(x$frame))
is.left <- (node%%2 ==0) #left hand sons
node.left <- node[is.left]
parent <- match(node.left/2, node)
##Put left splits at the parent node
if(splits) {
left.child <- match(2 * node, node)
right.child <- match(node * 2 + 1, node)
rows <- labels(x, pretty = pretty)
if(fancy) {
## put split labels on branches instead of nodes
xytmp <- rpart.branch(x=xy$x,y=xy$y,node=node)
leftptx <- (xytmp$x[2,]+xytmp$x[1,])/2
leftpty <- (xytmp$y[2,]+xytmp$y[1,])/2
rightptx <- (xytmp$x[3,]+xytmp$x[4,])/2
rightpty <- (xytmp$y[3,]+xytmp$y[4,])/2
FUN(leftptx,leftpty+.52*cxy[2],
rows[left.child[!is.na(left.child)]],...)
FUN(rightptx,rightpty-.52*cxy[2],
rows[right.child[!is.na(right.child)]],...)
}
else FUN(xy$x, xy$y + 0.5 * cxy[2], rows[left.child], ...)
}
leaves <- if(all) rep(TRUE, nrow(frame)) else frame$var == "<leaf>"
if (is.null(frame$yval2))
stat <- x$functions$text(yval=frame$yval[leaves],
dev=frame$dev[leaves],
wt=frame$wt[leaves],
ylevel=ylevels, digits=digits,
n=frame$n[leaves], use.n=use.n)
else
stat <- x$functions$text(yval=frame$yval2[leaves,],
dev=frame$dev[leaves],
wt=frame$wt[leaves],
ylevel=ylevels, digits=digits,
n=frame$n[leaves], use.n=use.n)
oval <- function(middlex,middley,a,b) {
theta <- seq(0,2*pi,pi/30)
newx <- middlex + a*cos(theta)
newy <- middley + b*sin(theta)
polygon(newx,newy,border=TRUE,col=0)
## polygon(newx,newy,border=T)
}
rectangle <- function(middlex, middley,a,b) {
newx <- middlex + c(a,a,-a,-a)
newy <- middley + c(b,-b,-b,b)
polygon(newx,newy,border=TRUE,col=0)
## polygon(newx,newy,border=T)
}
if(fancy) {
## find maximum length of stat
maxlen <- max(string.bounding.box(stat)$columns) + 1
maxht <- max(string.bounding.box(stat)$rows) +1
if(fwidth<1) a.length <- fwidth*cxy[1]*maxlen
else a.length <- fwidth*cxy[1]
if(fheight<1) b.length <- fheight*cxy[2]*maxht
else b.length <- fheight*cxy[2]
### create ovals and rectangles here
## sqrt(2) creates the smallest oval that fits around the
## best fitting rectangle
for(i in parent) oval(xy$x[i],xy$y[i],
a=sqrt(2)*a.length/2, b=sqrt(2)*b.length/2)
child <- match(node[frame$var=="<leaf>"],node)
for(i in child) rectangle(xy$x[i],xy$y[i],
a=a.length/2,b=b.length/2)
}
##if FUN=text then adj=1 puts the split label to the left of the
## split rather than centered
##Allow labels at all or just leaf nodes
## stick values on nodes
if(fancy) FUN(xy$x[leaves], xy$y[leaves] + .5 * cxy[2], stat, ...)
else FUN(xy$x[leaves], xy$y[leaves] -0.9 * cxy[2], stat, adj=.5,
col = leaf.col, ...)
invisible()
}
, comment = "09/02/2005")
"Updown" <-
structure(function(cv)
{
## Bill Venables's function (adjusted) to make uppercase first
## letters, and lower the rest.
cv <- casefold(cv)
first <- substring(cv, 1, 1)
FIRST <- casefold(first, upper = T)
paste(FIRST, substring(cv, 2), sep = "")
}
, comment = "05/04/2002")
"add.text" <-
structure(function(xp = 0.5, yp = 0.5, tex = NULL, ygap = 0.05, cex = 0.8, adj = 0,
adj.y = 0, font = 1, col = "black", srt = 0, ...)
{
### Cut down from bar.legs: used only for text; can also be rotated
### labels is character vector of legend text
### the length of tex sets the length of other vectors
###
###
### Establish positions in page co-ordinates [0,1] then use regular text function
### From Splus a/c began 11/07/1997
### Added vertical adj bit 17/6/04
uu <- par()$usr
xd <- diff(uu[1:2])
yd <- diff(uu[3:4])
yp <- yp - (seq(along = tex) - 1) * ygap
if(length(xp) == 1)
xp <- rep(xp, length(tex))
if(length(cex) == 1)
cex <- rep(cex, length(xp))
if(length(adj) == 1)
adj <- rep(adj, length(xp))
if(length(font) == 1)
font <- rep(font, length(xp))
if(length(col) == 1)
col <- rep(col, length(xp))
xh <- uu[1] + xd * xp
yh <- uu[3] + yd * yp
for(i in 1:length(tex)) {
text(xh[i], yh[i], tex[i], adj = c(adj[i], adj.y), cex = cex[i],
font = font[i], srt = srt, col = col[i], ...)
}
}
, comment = "25/10/2006")
"agg" <-
structure(function(data, by, FUNS)
{
#
# DATE WRITTEN: 14 May 1998
# AUTHOR: John R. Wallace (jrw@fish.washington.edu)
#
out <- aggregate.data.frame(data, by, eval(parse(text = FUNS[1])))
for(i in 2:length(FUNS)) {
tmp <- aggregate.data.frame(data, by, eval(parse(text = FUNS[i])))
tmp <- tmp[, ncol(tmp)]
out <- cbind(out, tmp)
}
dimnames(out)[[2]][ - (1:length(by))] <- FUNS
out
}
, comment = "15/09/1998")
"allfit" <-
structure(function(choice = 1, mn = NULL, xfun = function(x)
x, xfun.inv = NULL, datafun = glean.ramp, link = "cloglog", interval = T, pc =
c(50, 95, 99), pc.monotone = c(loess = 50, loess = 99), plot.monotone =
NULL, cutoff.mono = 0, span.mono = 1, pc.spline = NULL, plot.spline = F,
p.min = NULL, printdat = F, save.resid = F, cm.strategy = "adjust.later",
ext = "", cm.code = NULL, print.plot = F, full.robust = F, ...)
{
## Purpose:
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: John Maindonald, Creation date: 7 Sep 1994 (about)
#Changed May 1999 to make pc = 95 part of the default
xtra <- list(...)
new.page <- T
if(!is.null(names(xtra))) {
h <- xtra$h
ri <- xtra$ri
cj <- xtra$cj
if(any(is.na(match(names(xtra), c("h", "ri", "cj")))))
cat("\n*** There are one or more illegal parameters to allfit()\n",
"from among the names:", names(xtra), "***\n")
}
else {
ri <- 3
cj <- 4
h <- NULL
}
if(is.character(choice))
cat(choice)
if((!is.null(plot.monotone)) | plot.spline) {
oldpar <- par(mfrow = c(ri, cj), oma = c(0, 1, 5, 0))
on.exit(par(oldpar))
}
cat(" ", date(), fill = T)
cat("** Link =", link, "** x-transformation =", as.character(substitute(
xfun)), "**", fill = T)
if(!is.null(h))
u <- datafun(choice = choice, h = h)
else u <- datafun(choice = choice)
if(is.null(cm.code))
cm.code <- u$cm.code
if(is.null(cm.code))
cm.code <- 0
cm.allcodes <- u$cm.allcodes
temp <- u$temp
stage <- u$stage
if(all(xfun(exp(1:5)) == (1:5)))
takelog <- T
else takelog <- F
idset <- u$id
uset <- u$trial.id.order
if(is.null(uset))
uset <- unique(idset)
if(is.null(mn))
mn <- 1:length(uset)
if(max(mn) > length(uset)) {
maxj <- length(uset)
cat("You have specified a final dataset no.", max(mn), "\nwhereas", maxj,
"datasets only are available\n")
keep <- mn <= maxj
if(sum(keep) > 0)
mn <- mn[keep]
else stop(paste("Dataset(s)", mn, "not found."))
}
alltimes <- u$times
alltot <- u$total
alldead <- u$dead
cm <- u$cm
cmtot <- u$cmtot
cutx <- u$cutx
offset <- u$offset
xaxtitle <- u$xaxtitle
treat <- u$treat
maint <- u$maint
add.main <- u$add.main
legend <- u$legend
leg.brief <- u$leg.brief
xxlabels <- u$xlabels
here <- !is.na(alltimes) & !is.na(alldead) & !is.na(alltot)
print(maint)
if(!is.null(pc)) {
ltmat <- matrix(, nrow = length(mn), ncol = length(pc))
semat <- matrix(, nrow = length(mn), ncol = length(pc))
dimnames(ltmat) <- list(NULL, paste(pc))
}
else {
ltmat <- NULL
semat <- NULL
}
if(!is.null(pc.spline)) {
ltmat.spline <- matrix(, nrow = length(mn), ncol = length(pc.spline))
dimnames(ltmat.spline) <- list(NULL, paste(pc.spline))
}
else ltmat.spline <- NULL
if(!is.null(pc.monotone)) {
ltmat.monotone <- matrix(, nrow = length(mn), ncol = length(pc.monotone)
)
dimnames(ltmat.monotone) <- list(NULL, paste(pc.monotone))
}
else ltmat.monotone <- NULL
dimnames(ltmat) <- list(NULL, paste(pc))
dimnames(semat) <- list(NULL, paste(pc))
b0 <- array(, length(mn))
b1 <- array(, length(mn))
cm.est <- array(, length(mn))
dev <- array(, length(mn))
df <- array(, length(mn))
dev.robust <- array(, length(mn))
df.robust <- array(, length(mn))
times.list <- vector("list", length(mn))
total.list <- vector("list", length(mn))
resid.list <- vector("list", length(mn))
resid.dev.list <- vector("list", length(mn))
leg.allfit <- array(, length(mn))
i <- 0
for(j in mn) {
i <- i + 1
idj <- uset[j]
here.j <- idset == idj & here
leg <- legend[j]
leg.allfit[i] <- leg
dead <- alldead[here.j]
tot <- alltot[here.j]
mins <- alltimes[here.j]
if(!is.null(cutx))
cutoff <- cutx[j]
else cutoff <- NULL
if(!is.null(cm))
cmj <- cm[j]
else cmj <- NULL
uu <- fitconf(link = link, mins, dead, tot, cutoff, offset, j, leg = leg,
interval = interval, pc = pc, pc.spline = pc.spline, pc.monotone =
pc.monotone, plot.spline = plot.spline, plot.monotone = plot.monotone,
cutoff.mono = cutoff.mono, span.mono = span.mono, xfun = xfun, cm =
cmj, numcm = cmtot[j], p.min = p.min, printdat = F, save.resid =
save.resid, cm.strategy = cm.strategy, cm.code = cm.code, cm.allcodes
= cm.allcodes, xfun.inv = xfun.inv, full.robust = full.robust)
if((!is.null(plot.monotone)) && prod(par()$mfg[1:2]) == 1) {
new.page <- T
mtext(paste(maint, add.main), 3, line = 2.75, outer = T, cex = 1)
mtext("Curves & LT's are from smooths following monotone regression",
3, line = -0.75, outer = T, cex = 0.80000000000000004)
}
if((!is.null(plot.monotone)) && prod(par()$mfg[1:2]) == prod(par()$mfg[3:
4]) && print.plot && new.page) {
dev.print(device = pscript, horiz = T)
new.page <- F
}
if(printdat) {
print(legend[j])
print(cbind(mins, dead, tot))
}
ltmat[i, ] <- uu$lt
if(!is.null(ltmat.spline))
ltmat.spline[i, ] <- uu$lt.spline
if(!is.null(ltmat.monotone))
ltmat.monotone[i, ] <- uu$lt.monotone
semat[i, ] <- uu$se
b0[i] <- uu$b0
b1[i] <- uu$b1
dev[i] <- uu$dev
df[i] <- uu$df
cm.est[i] <- uu$cm
dev.robust[i] <- uu$dev.robust
df.robust[i] <- uu$df.robust
resid.list[[i]] <- uu$resid
resid.dev.list[[i]] <- uu$resid.dev
times.list[[i]] <- uu$times
total.list[[i]] <- uu$total
}
if((!is.null(plot.monotone)) && print.plot && new.page)
dev.print(device = pscript, horiz = T) #
invisible(list(lt50 = NULL, lt99 = NULL, se50 = NULL, se99 = NULL,
intercept = b0, slope = b1, lt = ltmat, se = semat, lt.spline =
ltmat.spline, lt.monotone = ltmat.monotone, span.mono = span.mono, cutx
= cutx, dev = dev, df = df, cm.strategy = cm.strategy, cm = cm.est,
cm.code = cm.code, dev.robust = dev.robust, df.robust = df.robust,
takelog = takelog, total = total.list, resid = resid.list, resid.dev =
resid.dev.list, times = times.list, datafun = as.character(substitute(
datafun)), choice = choice, temp = temp, stage = stage, add.main =
add.main, legend = leg.allfit, leg.brief = leg.brief))
}
, comment = "08/09/2004")
"ang" <-
structure(function(u)
{
(180 * asin(sqrt(u)))/atan(1)/4
}
, comment = "26/09/2001")
"ang.bt" <-
structure(function(y)
{
100 * (sin(pi/180 * y))^2
}
, comment = "16/02/2005")
"append.cols" <-
structure(function(data = natural.df, repeat.cols = 1:4, collapse.cols = 5:6, resp.id =
"Resp.id", out.col = "Response")
{
# By default repeats the first four columns twice and
# makes columns 5 and 6 into one
data <- as.data.frame(data)
if(is.numeric(repeat.cols))
repeat.cols <- names(data)[repeat.cols]
if(is.numeric(collapse.cols))
collapse.cols <- names(data)[collapse.cols]
repeat.df <- as.data.frame(data[, repeat.cols])
collapse.df <- data[, collapse.cols]
ind <- 1:nrow(repeat.df)
repeat.times <- ncol(collapse.df)
repeated.df <- as.data.frame(repeat.df[rep(ind, repeat.times), ])
names(repeated.df) <- repeat.cols # if only one column
repeated.df[, resp.id] <- rep(collapse.cols, rep(max(ind), repeat.times))
right <- unlist(collapse.df)
repeated.df[[out.col]] <- right
unfactor(repeated.df)
}
, comment = "30/04/2001")
"append.rows" <-
structure(function(x, id.cols = 1, spread.col = 3, spread.id = 2,
new.cols = NULL)
{
### Purpose: Unstacks dataframe -- converse of append.cols
### ----------------------------------------------------------------------
### Modified from: append.cols (a bit)
### ----------------------------------------------------------------------
### Arguments: new.cols: names of new columns
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 20 Feb 2003, 09:28
x <- as.data.frame(x)
### more useful in character form
if(is.numeric(spread.col))
spread.col <- names(x)[spread.col]
if(is.numeric(spread.id))
spread.id <- names(x)[spread.id] # column identifier to get new names
if(is.numeric(id.cols))
id.cols <- names(x)[id.cols] # rows that identify general groups
if(is.null(new.cols))
new.cols <- unique(x[, spread.id])
nrow.out <- nrow(x)/length(new.cols)
### set up dataframe and then add columns to it from rows of original
x <- df.sort(x, spread.id) # make sure it's ordered
out.df <- as.data.frame(x[seq(nrow.out), id.cols])
row.names(out.df) <- seq(nrow(out.df))
names(out.df) <- id.cols
spread.mat <- matrix(x[, spread.col], ncol = length(new.cols))
dimnames(spread.mat) <- list(NULL,new.cols)
cbind(out.df, spread.mat)
}
, comment = "20/02/2003")
"average.down" <-
structure(function(x)
{
# uses fill.down() to average the missing values with ones present
y <- fill.down(x)
z <- rev(fill.down(rev(x)))
xx <- (y + z)/2
}
, comment = "10/08/1999")
"axis.fudge" <-
structure(function(ax.cex = .8, las = 2, y.tix, mgp = c(2, .7, 0), fudge = 0)
{
## Purpose: Fudges y-axis labels (if the default is too low)
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments: y.ax.fudge is the amount to bump up the position of
## ytick labels
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 17 Jun 2004, 13:44
axis(2, cex.axis = ax.cex, mgp = mgp, at = y.tix, labels = FALSE)
axis(2, cex.axis = ax.cex, las = 2, mgp = mgp,
at = y.tix + fudge, labels = paste(y.tix), tick = FALSE)
}
, comment = "18/06/2004")
"bar.legs" <-
structure(function(xp = 0.5, yp = 0.5, ebv = NULL, labs = NULL, colour = F, pchs = 0:10,
ltys = 1:9, line.break = 0,
cols = rep(1, 20), wid = 0.007, eq.tex.gap = T,
ygap = ifelse(eq.tex.gap, 0.05, 0.02), xgap = 0.01, leg.cex = 0.6,
point.cex = 0.8, adj = 0, hershey = F, font = 1, line.leng = ifelse(is.null(pchs),
0.03, 0.05), mixed = F, bar.lwd = par()$lwd, line.lwd = par()$lwd,
pt.lwd = bar.lwd, pt.lty = 1)
{
### For placing legends with error bars and adjusting the spaces between lines
###
### xp and yp are x & y positions to begin
### ebv is vector of errorbars widths: or a list of two vectors (max and min)
### labels is character vector of legend text
### wid is "width" of errorbars: adjusts eps in errorbars function called from this one
### ygap space between elements of legend
### pchs is vector of plotting characters
### mixed is T when mixed.text is to be used for mixing characters in labs
###
### font is becomes vfont if herschey fonts are to be used, in which case, it is a
### vector of two elements: otherwise a numeric single value.
uu <- par()$usr
xd <- diff(uu[1:2])
yd <- diff(uu[3:4])
xh <- uu[1] + xd * xp
yh <- uu[3] + yd * yp
if(is.vector(labs))
labs <- as.list(labs)
hershey <- is.character(font)
if(is.null(ebv)) xh <- xh - xgap * xd
if(length(point.cex) < length(pchs))
point.cex <- rep(point.cex, length(pchs))
if(length(line.lwd) < length(ltys))
line.lwd <- rep(line.lwd, length.out = length(ltys))
if(length(cols) < length(ltys))
cols <- rep(cols, length.out = length(ltys))
### puts lines() back to xp
if(!is.list(ebv)) {
### If ebv is a single vector, only one lot of error bars:
for(i in 1:length(labs)) {
### If space between errorbars is to be made equal, start top of first error bar
### at beginning of legend. Otherwise, the legend text is at that position.
### If errorbar spaces are to be equal, half adjustment is done before drawing one,
### otherwise, all adjustment happens after drawing each one.
if(!eq.tex.gap) yh <- yh - ebv[i] ###
if(!is.null(ebv))
errorbars(xh, yh, ebv[i], eps = wid * xd, lwd = bar.lwd)
if(!is.null(pchs))
points(xh + (line.leng/2 + xgap) * xd, yh, pch = pchs[i], cex =
point.cex[i], col = cols[i], lwd = pt.lwd, lty = pt.lty)
if(!is.null(ltys))
if(line.break > 0) {
{
lines(xh + (c(0, line.leng/2 - line.break) + xgap) * xd,
rep(yh, 2), lty = ltys[i], col = cols[i], lwd = line.lwd[i])
lines(xh + (c(line.leng/2 + line.break, line.leng) + xgap) *
xd, rep(yh, 2), lty = ltys[i], col = cols[i], lwd =
line.lwd[i])
}
}
else lines(xh + (c(0, line.leng) + xgap) * xd, rep(yh, 2), lty =
ltys[i], col = cols[i], lwd = line.lwd[i])
if(!is.null(labs)) {
if(is.expression(labs[[i]]) & hershey) stop(
"Hershey fonts cannot be used with expressions\n")
### Use herschey fonts, font becomes vfont: adj is a bit different too
if(hershey)
text(xh + (line.leng + xgap * 2) * xd, yh, labels = labs[[i]],
adj = c(adj, 0.5), cex = leg.cex, vfont = font)
else text(xh + (line.leng + xgap * 2) * xd, yh, labels = labs[[i]],
adj = c(adj, 0.5), cex = leg.cex, font = font)
}
### Move down the appropriate amount for the next one:
yh <- yh - ifelse(eq.tex.gap, 0, ebv[i]) - ygap * yd
}
}
else {
### Two error bars, maximums and minimums are drawn
### Move starting point to allow for extra errorbar
xh <- xh + (xgap + wid) * xd
ebv1 <- ebv$max
ebv2 <- ebv$min
for(i in 1:length(labs)) {
if(!eq.tex.gap)
yh <- yh - ebv1[i]
errorbars(xh - (xgap + wid) * xd, yh, ebv1[i], eps = wid * xd, lwd =
bar.lwd)
errorbars(xh, yh, ebv2[i], eps = wid * xd, lwd = bar.lwd)
if(!is.null(pchs))
points(xh + (line.leng/2 + xgap) * xd, yh, pch = pchs[i], cex =
point.cex[i], col = cols[i], lwd = pt.lwd)
if(line.break > 0) {
{
lines(xh + (c(0, line.leng/2 - line.break) + xgap) * xd,
rep(yh, 2), lty = ltys[i], col = cols[i], lwd = line.lwd[i])
lines(xh + (c(line.leng/2 + line.break, line.leng) + xgap) * xd,
rep(yh, 2), lty = ltys[i], col = cols[i], lwd = line.lwd[i])
}
}
else lines(xh + (c(0, line.leng) + xgap) * xd, rep(yh, 2),
lty = ltys[i], col = cols[i], lwd = line.lwd[i])
if(!is.null(labels)) {
if(mixed)
mixed.text(xh + (line.leng + xgap * 2) * xd, yh, labs[[i]], adj
= adj, cex = leg.cex, vfont = vfont)
else text(xh + (line.leng + xgap * 2) * xd, yh, labs[[i]], adj =
adj, cex = leg.cex, vfont = vfont)
}
yh <- yh - ebv1[i] - ygap * yd
}
}
}
, comment = "13/08/2003")
"blank.plot" <-
structure(function(x = 0, y = 0, ...)
{
# Does a blank plot.
#
# Useful if you want to change the axes on the current plot (in which case specify
# par(new=F))
# Sometimes you might want to prevent plotting in one part of a multiplot page.
plot(x, y, xlab = "", ylab = "", xaxt = "n", yaxt = "n", pch = " ", bty =
"n", ...)
}
, comment = "06/05/1997")
"boxplot.fn" <-
structure(function(x = delbbd.datf, floop = NULL, rloop = NULL, trans = function(z)
z, mod1 = 1, mod2 = 0, notch = F, subset = NULL, paste.cols = NULL)
{
# 1a) Adding Paste Columns
if(!is.null(paste.cols) & !is.logical(paste.cols)) {
indc <- apply(x[paste.cols], 1, paste, collapse = ":")
#ind <- match(indc, unique(indc))
x <- cbind(indc, x)
dimnames(x)[[2]][1] <- paste("treat", paste(paste.cols, collapse = ""),
sep = "")
if(is.null(floop)) {
floop <- 1
}
else {
floop <- c(1, floop + 1)
}
}
# 1b) For loop for Responses
for(j in rloop) {
response.names <- dimnames(x)[[2]][j]
response <- x[, j]
if(!is.null(subset)) {
response[response > subset] <- NA
}
yval <- trans(response/mod1 + mod2) # 2) For loop for factors
for(i in floop) {
factors <- factor(x[, i])
factor.names <- dimnames(x)[[2]][i] # 3) Boxplot function
plot.factor(factors, yval, xlab = factor.names, ylab = response.names,
notch = notch, boxmeans = T, style.bxp = "old")
if(i == floop[1] & j == rloop[1])
mtext(paste(as.character(substitute(trans)), " ", as.character(
substitute(x))), 3, cex = 1, line = 1.3999999999999999)
dev.ask()
}
}
}
, comment = "06/12/1996")
"catalogue" <-
structure(function(x = NULL, outfile = "add.dates.sc", new = TRUE, remove = FALSE)
{
### by default, adds today's date as a comment on any object that has none.
### new: delete any previous src file?
### src.now: delete any previous src file?
### remove: delete this src file?
###
del.comm <- paste("rm", outfile) # deleting previous src file
src.comm <- ppaste("source('", outfile, "')") # to use src file
if(new)
try(system(del.comm))
if(is.null(x))
x <- ls(pos = 1)
new.obj <- 0
today <- system("date +%d/%m/%Y", TRUE)
for(i in x) {
source.line <- ppaste("comment(", i, ") <- '", today, "'")
if(is.null(comment(get(i)))) {
write(source.line, file = outfile, append = TRUE)
new.obj <- new.obj + 1
}
}
if(remove)
system(del.comm)
if(new.obj > 0)
cat(paste("Now you can add the date comments to", new.obj, ifelse(
new.obj > 1, "objects", "object"), "using:\n", src.comm, "\n"))
else cat(paste("No objects need updating:\n"))
}
, comment = "02/04/2001")
"catalogue.gems" <-
structure(function(x = NULL, outfile = "gems.date.sc", new = TRUE, remove = FALSE)
{
###
### modifies catalogue() to incorporate dates in Gems functions from Splus
###
### x is character string of object names of interest. All if not specified
###
month.name <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug",
"Sep", "Oct", "Nov", "Dec")
month.no <- 1:12
names(month.no) <- month.name
attach(gems.df)
on.exit(detach("gems.df"))
del.comm <- paste("rm", outfile) # deleting previous src file
src.comm <- ppaste("source('", outfile, "')") # to use src file
if(new)
try(system(del.comm))
if(is.null(x))
x <- ls(pos = 1)
today <- system("date +%d/%m/%Y", TRUE)
new.obj <- 0
for(i in x) {
yr.i <- Year[Object == i]
dy.i <- Dy[Object == i]
mo.i <- month.no[Mo[Object == i]]
if(length(yr.i) > 0) {
dy.j <- substring(ppaste("0", dy.i), nchar(dy.i))
mo.j <- substring(ppaste("0", mo.i), nchar(mo.i))
today.i <- paste(dy.j, mo.j, yr.i, sep = "/")
source.line <- ppaste("comment(", i, ") <- '", ifelse(today.i == "",
today, today.i), "'")
if(is.null(comment(get(i)))) {
write(source.line, file = outfile, append = TRUE)
new.obj <- new.obj + 1
}
}
}
if(new.obj > 0)
cat(paste("Now you can add the date comments to", new.obj, ifelse(
new.obj > 1, "object", "objects"), "using:\n", src.comm, "\n"))
else cat(paste("No objects need updating:\n"))
}
, comment = "04/04/2001")
"catalogue.save.gems" <-
structure(function(x = NULL, outfile = "gems.date.save.sc", new = TRUE)
{
###
### modifies catalogue() to incorporate dates in Gems functions from Splus
###
### x is character string of object names of interest. All if not specified
###
attach(gem.df)
on.exit(detach("gem.df"))
del.comm <- paste("rm", outfile) # deleting previous src file
src.comm <- ppaste("source('", outfile, "')") # to use src file
if(new)
try(system(del.comm))
if(is.null(x))
x <- ls(pos = 1)
today <- system("date +%d/%m/%Y", TRUE)
new.obj <- 0
for(i in x) {
date.i <- Date[Object == i]
source.line <- ppaste("comment(", i, ") <- '", date.i, "'")
if(is.null(comment(get(i)))) {
write(source.line, file = outfile, append = TRUE)
new.obj <- new.obj + 1
}
}
if(new.obj > 0)
cat(paste("Now you can add the date comments to", new.obj, ifelse(
new.obj > 1, "object", "objects"), "using:\n", src.comm, "\n"))
else cat(paste("No objects need updating:\n"))
}
, comment = "31/05/2001")
"check.fact" <-
structure(function(df)
{
#Checks if columns of dataframes are factors:
sapply(df, is.factor)
}
, comment = "30/04/2001")
"check.na" <-
structure(function(df)
{
# Returns the number of NAs in any column of a dataframe
sapply(df, function(x)
sum(is.na(x)))
}
, comment = "21/09/1997")
"circ.boxplot" <-
structure(function(x, rad, fat = .1, col = "blue", rad.lim = 1,
pch = 16, cex = .8, seg = 600, test = FALSE)
{
## Purpose: Curved boxplots
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments: x- vector of values (in degrees),
## rad- how far from centre to draw
## fat- how thick the boxes are to be
## rad.lim- max radial length on plot
## test(mark individual values to test the box position)
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 29 Mar 2006, 13:20
deg <- boxplot(x, plot = FALSE)
box.stats <- deg$stats[,1]
names(box.stats) <- c("beginWhisker", "beginBox", "Median", "endBox", "endWhisker")
rad <- rad/rad.lim
fat <- fat/rad.lim # get to unit scale if not already
### draw inner and outer sides of boxes and the whiskers
circle(,, rad + fat/2, box.stats["beginBox"], box.stats["endBox"], col = col, seg = seg)
circle(,, rad - fat/2, box.stats["beginBox"], box.stats["endBox"], col = col, seg = seg)
circle(,, rad , box.stats["beginWhisker"], box.stats["beginBox"], col = col, seg = seg)
circle(,, rad , box.stats["endBox"], box.stats["endWhisker"], col = col, seg = seg)
### co-ordinates for box edges
x.box.i <- (rad - fat/2)* cos(box.stats * pi/180)[c("beginBox", "endBox")]
y.box.i <- (rad - fat/2) * sin(box.stats * pi/180)[c("beginBox", "endBox")]
x.box.o <- (rad + fat/2)* cos(box.stats * pi/180)[c("beginBox", "endBox")]
y.box.o <- (rad + fat/2) * sin(box.stats * pi/180)[c("beginBox", "endBox")]
## draw two radial lines for ends of boxes
lines(c(x.box.i["beginBox"], x.box.o["beginBox"]),
c(y.box.i["beginBox"], y.box.o["beginBox"]), col = col)
lines(c(x.box.i["endBox"], x.box.o["endBox"]),
c(y.box.i["endBox"], y.box.o["endBox"]), col = col)
### co-ordinates for whisker edges
x.whisker.i <- (rad - fat/4)*
cos(box.stats * pi/180)[c("beginWhisker", "endWhisker")]
y.whisker.i <- (rad - fat/4) *
sin(box.stats * pi/180)[c("beginWhisker", "endWhisker")]
x.whisker.o <- (rad + fat/4)*
cos(box.stats * pi/180)[c("beginWhisker", "endWhisker")]
y.whisker.o <- (rad + fat/4) *
sin(box.stats * pi/180)[c("beginWhisker", "endWhisker")]
## draw two radial lines for ends of whiskers
lines(c(x.whisker.i["beginWhisker"], x.whisker.o["beginWhisker"]),
c(y.whisker.i["beginWhisker"], y.whisker.o["beginWhisker"]), col = col)
lines(c(x.whisker.i["endWhisker"], x.whisker.o["endWhisker"]),
c(y.whisker.i["endWhisker"], y.whisker.o["endWhisker"]), col = col)
### co-ordinates for median edges
x.median.i <- (rad - fat/2) * cos(box.stats[c("Median")] * pi/180)
y.median.i <- (rad - fat/2) * sin(box.stats[c("Median")] * pi/180)
x.median.o <- (rad + fat/2) * cos(box.stats[c("Median")] * pi/180)
y.median.o <- (rad + fat/2) * sin(box.stats[c("Median")] * pi/180)
## draw radial line for median
lines(c(x.median.i["Median"], x.median.o["Median"]),
c(y.median.i["Median"], y.median.o["Median"]), lwd = 3, col = col)
### Mark in outliers
x.out <- rad * cos(deg$out * pi/180)
y.out <- rad * sin(deg$out * pi/180)
points(x.out, y.out, pch = pch, cex = cex, col = col)
if(test){ # if we're testing to see if the box is in the right quadrant
points(rad * cos(x * pi/180), rad * sin(x * pi/180),
pch = 1, col = col, cex = cex )
}
}
, comment = "02/11/2006")
"circle" <-
structure(function(cen.x = 0, cen.y = 0, rad = 1, from = 0, to = 360, seg = 400, ...)
{
## Purpose: draws a circle (or part thereof) at cen.x, cen.y, with radius rad
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 3 Feb 2006, 16:08
bits <- seq(from, to, length = 400)
x.coords <- cen.x + rad * cos (bits * pi/180)
y.coords <- cen.y + rad * sin (bits * pi/180)
lines(x.coords, y.coords, ...)
}
, comment = "29/03/2006")
"clg" <-
structure(function()
{
# Clears graphics window
plot(0, 0, new = T, xaxt = "n", pch = " ", yaxt = "n", bty = "n", xlab = "",
ylab = "")
}
, comment = "09/05/1997")
"clipline" <-
structure(function(xlim, ylim, a, b, funlog = function(x)
x, unlog = function(x)
x, lty = 1, lwd = 1, loud = T, ...)
{
x0 <- xlim[1]
y0 <- a + b * funlog(x0)
x1 <- xlim[2]
y1 <- a + b * funlog(x1)
if(y0 < ylim[1]) {
y0 <- ylim[1]
x0 <- unlog((ylim[1] - a)/b)
}
else if(y0 > ylim[2]) {
y0 <- ylim[2]
x0 <- unlog((ylim[2] - a)/b)
}
if(y1 < ylim[1]) {
y1 <- ylim[1]
x1 <- unlog((ylim[1] - a)/b)
}
else if(y1 > ylim[2]) {
y1 <- ylim[2]
x1 <- unlog((ylim[2] - a)/b)
}
if(loud)
cat(" Join (", round(x0, 3), round(y0, 3), ") to (", round(x1, 3),
round(y1, 3), ")", fill = T)
lines(c(x0, x1), c(y0, y1), lty = lty, lwd = lwd, ...)
}
, comment = "13/09/2002")
"cloglog" <-
structure(function(pr)
{
log( - log(1 - pr))
}
, comment = "16/04/1996")
"cloglog.bt" <-
structure(function(x)
{
# To calculate the proportion that gives the known complememtary log log value
1 - exp( - exp(x))
}
, comment = "24/11/1997")
"cols.for.genstat" <-
structure(function(x, num = F)
{
# Returns a string of the column names to use in genstat
# need only to change the opening and closing " and to '
y <- names(x)
paste(y, collapse = ", ")
}
, comment = "12/04/1999")
"comb" <-
structure(function(N, x)
{
gamma(N + 1)/gamma(x + 1)/gamma(N - x + 1)
}
, comment = "23/07/1997")
"complot" <-
structure(function(data, xtitle = "Time (Minutes)", main = "", yfun = "cloglog", logx = F,
plimits = c(0.04, 0.9995), xtix = NULL, ytit = NULL, labp
= NULL, lty = 1, colours = F, dot.marks = T, xlab.line = 2, x.mgp = 0.2,
y.mgp = ifelse(post, 0.5, 0.6), ylab.line = 1.6, pchs = c(0, 1, 2, 5, 6:16),
id = NULL, ab.lty = 2, lab.cex = 1, ax.cex = 0.8, main.cex = 1.1, point.cex
= 0.9, mark.control = F, ...)
{
# Modified (very) version of simplot
#
# data is a dataframe of the sets to be plotted on one plot.
# with names of columns: Time, Dead, Total, Id
#
# xfun is the function to be applied to the x axis data before plotting it.
# It ia a list which will be extended to the same length as the number of sets
# of points to be plotted (unless already given as that length).
#
# y.mgp is the second value in an mgp parameter
# x.mgp is the amount by which the corresponding x axis parameter is decreased
#
#
# dot.marks: Mark 100% and 0% points with dots?
post <- names(dev.cur()) == "postscript"
par(...)
if(logx)
f <- log
else f <- function(x)
x
if(logx)
stop(
"\ncomplot function needs fixing to do log transformation on x-axis")
u.id <- unique(data$Id)
g <- link.function(yfun)
# Identify Control and Treatment data:
if(logx)
data <- data[data$Time > 0, ]
data$Mort <- data$Dead/data$Total
if(is.null(xtix))
x.range <- range(data$Time)
else x.range <- range(xtix)
p.range <- range(data$Mort) #
# Calculate adjustments for the space on the y-axis:
ylow <- g(plimits[1])
yhigh <- g(plimits[2])
eps <- (yhigh - ylow) * 0.01
ylow <- ylow - eps
yhigh <- yhigh + eps
if(is.null(labp))
labp <- c(1, 5, 10, 25, 50, 75, 95, 99)
labp <- labp[labp/100 >= plimits[1] & labp/100 <= plimits[2]] #
# Set up blank plot: (uses blank.plot() in /home/hrapgc/Gems/.Data/ )
blank.plot(f(x.range), c(ylow, yhigh), ...)
usr <- par()$usr
epsx <- 0.0125 * diff(usr[1:2])
epsy <- 0.0125 * diff(usr[3:4])
midhigh <- 0.5 * (yhigh + usr[4])
midlow <- 0.5 * (ylow + usr[3])
hundred.line <- F
zero.line <- F
colour <- is.character(colours)
for(i in seq(u.id)) {
df.i <- data[data$Id == u.id[i], ]
attach(df.i) #
# Plot "normal" points:
time <- Time[Mort > plimits[1] & Mort < plimits[2]]
mort <- Mort[Mort > plimits[1] & Mort < plimits[2]]
points(f(time), g(mort), pch = pchs[i], col = ifelse(colour, colours[i],
1), cex = point.cex) #
# Plot "high" points (i.e. ones near 100%)
high.mort <- Mort[Mort > plimits[2]]
high.time <- Time[Mort > plimits[2]]
points(f(high.time), rep(midhigh, length(high.time)), pch = pchs[i], cex
= point.cex, col = ifelse(colour, colours[i], 1))
if(any(Mort == 1)) {
hundred.line <- T
full.mort <- Mort[Mort == 1]
full.time <- Time[Mort == 1]
if(dot.marks)
points(f(full.time), rep(midhigh, length(full.time)), pch = ".",
cex = point.cex, col = ifelse(colour, colours[i], 1))
}
# Plot "low" points (i.e. ones near 0%)
low.mort <- Mort[Mort < plimits[1]]
low.time <- Time[Mort < plimits[1]]
points(f(low.time), rep(midlow, length(low.time)), pch = pchs[i], cex =
point.cex, col = ifelse(colour, colours[i], 1))
if(any(Mort == 0)) {
zero.line <- T
zero.mort <- Mort[Mort == 0]
zero.time <- Time[Mort == 0]
if(logx) {
zero.time <- zero.time[zero.time > 0]
zero.mort <- zero.mort[zero.time > 0]
}
if(dot.marks)
points(f(zero.time), rep(midlow, length(zero.time)), pch = ".",
cex = point.cex)
}
cont.mort <- Mort[Time == 0]
if(mark.control) {
# Mark control mortalities:
if(any(cont.mort < plimits[1])) text(0, midlow, "C")
else text(0, g(cont.mort), "C", cex = point.cex,
col = ifelse(colour, colours[i], 1))
}
detach("df.i")
}
# Do tricky things to y-axis:
mgp <- c(3, y.mgp, 0)
axis(2, at = g(labp/100), labels = paste(labp), adj = 1, mgp = mgp, cex.axis =
ax.cex, las = 1)
axis(2, at = c(ylow + 2 * epsy, yhigh - 2 * epsy), labels = F, tck = 0, las = 1)
for(k in 1:2) {
y.end <- c(ylow, yhigh)[k]
axis(2, at = c(y.end, usr[2 + k]), labels = F, tck = 0, las = 1)
}
if(zero.line)
abline(h = ylow, lty = ab.lty)
if(hundred.line)
abline(h = yhigh, lty = ab.lty)
axis(2, at = midlow, tck = 0, labels = "0", adj = 1, mgp = mgp, cex.axis =
ax.cex, las = 1)
lines(usr[1] + c( - epsx, epsx), c(ylow - 1.5 * epsy, ylow + 1.5 * epsy),
xpd = T)
lines(usr[1] + c( - epsx, epsx), c(ylow + 0.5 * epsy, ylow + 3.5 * epsy),
lty = lty, xpd = T)
axis(2, at = midhigh, tck = 0, labels = "100", adj = 1, mgp = mgp, cex.axis =
ax.cex, las = 1)
lines(usr[1] + c( - epsx, epsx), c(yhigh - 1.5 * epsy, yhigh + 1.5 * epsy),
xpd = T)
lines(usr[1] + c( - epsx, epsx), c(yhigh - 3.5 * epsy, yhigh - 0.5 * epsy),
lty = lty, xpd = T) #
# Make x-axis and labels
if(is.null(xtix))
xlab.pos <- pretty(range(data$Time))
else xlab.pos <- xtix
if(logx)
xlab.pos <- xlab.pos[xlab.pos > 0]
axis(1, at = f(xlab.pos), labels = paste(xlab.pos), cex.axis = ax.cex, mgp = mgp -
c(0, x.mgp, 0))
cat(xlab.pos, " \n") # xlab.pos specifies numeric labels
# Add break in x-axis if necessary: fill in left space
if(logx) {
x.cut1 <- diff(usr[1:2]) * 0.03 + usr[1]
x.cut2 <- diff(usr[1:2]) * 0.5 + usr[1]
axis(1, at = c(usr[1], x.cut1), labels = F, tck = 0)
axis(1, at = c(x.cut2, usr[2]), labels = F, tck = 0)
lines(x.cut1 + c( - epsx, epsx), usr[3] + 2 * c( - epsy, epsy), xpd = T)
lines(x.cut2 + c( - epsx, epsx), usr[3] + 2 * c( - epsy, epsy), xpd = T)
}
else lines(usr[1:2], usr[c(3, 3)])
mtext(side = 1, line = xlab.line, text = xtitle, cex = lab.cex)
mtext(side = 2, line = ylab.line, text = ifelse(is.null(ytit), paste(
"% Mortality,", yfun, "scale"), ytit), cex = lab.cex)
mtext(side = 3, text = main, cex = main.cex)
box(bty = "7")
}
, comment = "13/09/2002")
"count.nas" <-
structure(function(df)
{
# counts nas in columns of a dataframe have nas.
#
sapply(df, function(x)
sum(is.na(x)))
}
, comment = "01/12/1997")
"data.spot" <-
structure(function(ab, use, id = NULL, separate = T)
{
### To show table of times, dead and total to correspond to selected
### parts of ab type list of the form ab$use
### If id is not specified, all values will be shown
### separate will group each section and names individually, othewise
### only one table will be made
###
### WARNING: Will crash if the legends are not unique
names.ab <- names(ab)
if(length(use) > 1)
stop("use can be of length 1 only")
if(is.character(use))
ab.sub <- (1:length(names.ab))[!is.na(match(names.ab, use))]
else ab.sub <- use
aab <- ab[[ab.sub]]
glean <- aab$datafun
choice <- aab$choice
glean.fn <- get(glean)
glean.list <- glean.fn(choice)
df <- data.frame(id = glean.list$id, time = glean.list$times, dead =
glean.list$dead, total = glean.list$total)
df$live <- df$total - df$dead
df <- df[, c("id", "time", "live", "dead", "total")]
df$mort <- round(df$dead/df$total * 100)
if(is.null(id))
out.df <- df
else out.df <- df[!is.na(match(df[["id"]], id)), ]
print(glean.list$maint)
out.list <- list()
browser()
if(separate) {
if(is.null(id))
id <- unique(df$id)
for(i in id)
out.list[[glean.list$legend[i]]] <- df[df$id == i, ]
out.list
}
else {
if(!is.null(id))
print(glean.list$legend[id])
out.df
}
}
, comment = "18/04/1999")
"ddif" <-
structure(function(d1, d2)
{
### Purpose: Finds the number of days between two different dates
### ----------------------------------------------------------------------
### Modified from:
### ----------------------------------------------------------------------
### Arguments: d1 and d2 must be POSIXct objects with tz = "GMT"
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 20 Feb 2004, 09:09
d1s <- as.numeric(d1) # number of seconds since origin
d2s <- as.numeric(d2) #
(d1s - d2s)/24/3600
}
, comment = "20/02/2004")
"df.mat" <-
structure(function(df)
{
# Coerces a dataframe to a matrix, retaining levels instead of codes
# for any factors: see unfactor()
mat <- as.matrix(unfactor(df))
mat
}
, comment = "18/10/1996")
"df.sort" <-
structure(function(dfm, sort.by = names(dfm)[1], backwards = F)
{
### To sort the rows of a dataframe into an order set by vector "sort.by"
### If sort.by is numeric it is the column nos; if character, the column names.
if(is.character(sort.by)) {
sort.by.n <- match(sort.by, names(dfm))
if(any(is.na(sort.by.n))) {
stop(paste("\n", sort.by[is.na(sort.by.n)],
"is not a recognised name"))
}
else sort.by <- sort.by.n
}
for(i in sort.by) {
dfm.i <- dfm[order(dfm[, i]), ]
if(backwards)
dfm.i <- dfm[(rev(order(dfm[, i]))), ]
dfm <- dfm.i
}
dfm
}
, comment = "13/11/2002")
"df.summary" <-
structure(function(data, y.cols = NULL, split.cols = NULL, rep.col = "Rep", decimals = 3,
stats = "mean", help = F)
{
### cut down and modified from summarize()
### Keep Stephen's beginning:
### Requires text.bp() and df.to.file()
### [anything else will be in /home/hrapgc/Gems/.Data]
summary.names <- c("N", "Missing", "Mean", "Median", "Trimmed Mean",
"Std. Dev.", "SEM", "Min", "Max", "1st Quartile", "3rd Quartile", "Sum",
"Boxplots")
summary.types <- c("num", "missing", "mean", "median", "trm", "stdev",
"sem", "min", "max", "q1", "q3", "sum", "boxplots")
help.df <- (data.frame("Select one of these" = summary.types, "to get this"
= summary.names))
if(help) {
cat("\n In the argument \"stats\"\n")
print(help.df)
return()
}
### A few error messages:
if(length(stats) > 1)
stop(
"\n\tOnly one statistic can be calculated at a time:\n\tIf you require multiple stats with only one y.cols, \n\tyou can use summarize() with df.out set to TRUE"
)
### Make special functions to use in apply calls
q1 <- function(x)
quantile(x, na.rm = T)[2]
q3 <- function(x)
quantile(x, na.rm = T)[4]
no.missing <- function(x)
length(x[is.na(x)])
round.mean <- function(x, y = decimals)
round(mean(x, na.rm = T), y)
round.sem <- function(x, y = decimals)
round(sem(x, na.rm = T), y)
round.sum <- function(x, y = decimals)
round(sum(x, na.rm = T), y)
add.boxplot <- function(x, y.range = range.j, width = bp.width)
text.bp(x, y.range, width) #
### Functions are now set up:
###
### Make named vector of functions to be used:
use.functions <- c("length", "no.missing", "round.mean", "median", "trm",
"stdev", "round.sem", "min", "max", "q1", "q3", "round.sum",
"add.boxplot")
names(use.functions) <- summary.types
names(summary.names) <- summary.types #
### Use dataframes instead of matrices:
if(is.matrix(data))
data <- unfactor(as.data.frame(data))
else data <- unfactor(data) # avoid pesky factors
data <- df.sort(data, rev(split.cols)) #
### Make character vectors of y.cols and split.cols if necessary:
if(is.numeric(y.cols))
y.cols <- names(data)[y.cols]
if(is.numeric(split.cols)) split.cols <- names(data)[split.cols] #
### Create index and adjust for single split columns if necessary:
indx.df <- data.frame(data[, split.cols])
names(indx.df) <- split.cols # Necessary for single split.cols
attach(indx.df)
indx <- NULL
for(i in split.cols)
indx <- paste(indx, get(i), sep = ":")
detach("indx.df")
out.df <- data.frame(indx.df[match(unique(indx), indx), split.cols])
names(out.df) <- split.cols # Necessary for single split.cols
dimnames(out.df) <- list(paste(1:nrow(out.df)), names(out.df))
for(j in y.cols)
out.df[, j] <- s.tapply(data[, j], indx, get(use.functions[stats]))
out.df
}
, comment = "30/04/2001")
"df.to.file" <-
structure(function(df = unfactor(junk.df), file = "", pad = "mid", append = F,
row.nos = F, col.names = T)
{
### To print a data frame to a text file and having names line up with values
### underneath, and no need to use tab spacing
###
### Set up some padding functions: Need columns with same number of characters
### all the way down.
pad.char.cols <- function(y)
{
# For padding character columns
biggest <- max(nchar(y))
col.width <- nchar(y)
padding <- paste(rep(" ", biggest - min(col.width)), collapse = "")
out <- substring(paste(y, padding, sep = ""), 1, biggest)
out
}
pad.left <- function(z, pads)
{
padding <- paste(rep(" ", pads), collapse = "")
paste(padding, z, sep = "")
}
pad.mid <- function(z, pads)
{
# Centres text in available space
padding.right <- paste(rep(" ", pads %/% 2), collapse = "")
padding.left <- paste(rep(" ", pads - pads %/% 2), collapse = "")
paste(padding.left, z, padding.right, sep = "")
}
pad.right <- function(z, pads)
{
padding <- paste(rep(" ", pads), collapse = "")
paste(z, padding, sep = "")
}
### Character columns need different treatment from numeric columns
#browser()
char.cols <- sapply(df, is.character)
if(any(char.cols))
df[char.cols] <- sapply(df[char.cols], pad.char.cols)
if(any(!char.cols))
df[!char.cols] <- sapply(df[!char.cols], format) #
### Sometimes the names of columns are wider than the columns they name,
### sometimes vice versa.
names.width <- nchar(names(df))
row.width <- sapply(df, function(x)
max(nchar(x))) #
#browser()
### (the width of the characters in the columns (and rows) as distinct
### from their names
name.pads <- row.width - names.width
row.pads <- name.pads * -1
name.pads[name.pads < 0] <- 0
row.pads[row.pads < 0] <- 0
pad.names <- name.pads > 0
pad.rows <- row.pads > 0 #
### Pad out the column names if necessary:
pad.funct <- get(paste("pad.", pad, sep = ""))
if(any(pad.names)) {
stretch.names <- names(df)[pad.names]
for(i in stretch.names) {
names(df)[names(df) == i] <- pad.funct(i, name.pads[i])
}
}
### likewise for the rows and columns
if(any(pad.rows)) {
stretch.rows <- names(df)[pad.rows]
for(j in stretch.rows)
df[, j] <- pad.funct(df[, j], row.pads[j])
}
### the write() function doesn't use the names of the matrix, so put them
### into it.
df2 <- df
if(row.nos) {
df2 <- cbind(seq(nrow(df)), df)
no.wid <- nchar(nrow(df))
no.col.name <- paste(rep(" ", no.wid), collapse = "")
names(df2) <- c(no.col.name, names(df))
}
if(col.names)
df2 <- rbind(names(df2), as.matrix(df2))
else df2 <- as.matrix(df2)
write(t(df2), file, ncol(df2), append = append)
}
, comment = "24/07/2001")
"disp" <-
structure(function(x, r = 4)
{
### Purpose: Estimates dispersion of a glm
### ----------------------------------------------------------------------
### Modified from:
### ----------------------------------------------------------------------
### Arguments: x is a glm object (mabye can be more generic)
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 8 Aug 2003, 11:21
dev.x <- deviance(x)
df.x <- x$df.residual
round(dev.x/df.x, r)
}
, comment = "08/08/2003")
"dmodes" <-
structure(function(t.ob = NULL, pos = 1, max = NULL, date.sort = TRUE)
{
### Returns a Data frame showing modes of objects in the stated directory.
### Can be sorted by Mode, Length, or Date in addition to default Object
### A variation on this one is modes(). Made for lazy typists to sort by object name
### max is the maximum number of objects returned; useful if ls() is long and
### you're only interested in the first few
###browser()
if(is.null(t.ob)) ls.o <- objects(pos = pos) else ls.o <- objects(pos = pos,
pattern = t.ob)
if(!is.null(max) && max > length(ls.o))
print(paste("Only", length(ls.o), "match", t.ob, "in", search()[pos]))
if(is.null(max) || (max > length(ls.o))) max <- length(ls.o) #
### Rearrange date information
date.chr <- NULL
for(i in seq(ls.o)) {
if(FALSE){
cat(paste("Got down to", ls.o[i], ":\n"))
##browser()
}
date.i <- comment(get(ls.o[i], pos = pos))
if(is.null(date.i))
###
###date.i <- try(comment(get(ls.o[i], pos = pos)))
###if(class(date.i) == "try-error")
date.i <- "NA/NA/NA"
date.chr[i] <- date.i
}
break.at <- occurrence(date.chr, "/", 1:2)
if(is.null(dim(break.at)))
break.at <- matrix(break.at, nrow=1)
first <- break.at[, 1] - 1
second <- break.at[, 2] - 1
day <- as.numeric(substring(date.chr, 1, first))
month <- as.numeric(substring(date.chr, first + 2, second))
year <- as.numeric(substring(date.chr, second + 2))
ob.df <- unfactor(data.frame(Object = ls.o, Day = day, Mon = month, Year =
year))
ob.df <- df.sort(ob.df, "Object", backwards = TRUE)
if(all(is.na(ob.df$Day)))
date.sort <- FALSE ## no point trying to sort
if(date.sort) {
xx <- df.sort(ob.df, c("Day","Mon"))
xx <- df.sort(xx, "Year", backwards = TRUE)[seq(max), ]
}
else xx <- df.sort(ob.df, "Object")
xx$Date <- date.chr[as.numeric(row.names(xx))]
options(warn = -1)
on.exit(options(warn = 0))
out <- list()
dummy <- matrix(nc = 3, nr = 3) # a non-null last list element
out$Object <- c(xx$Object, "dummy")
objs.mat <- cbind(out$Object)
out$Mode <- apply(objs.mat, 1, function(X)
mode(get(X))) #
### Dimension information of dataframes and matrices
dimensions <- apply(objs.mat, 1, function(X)
dim(get(X))[1:2])
y <- dimensions
if(!is.null(y)) {
if(is.list(y)) {
y[unlist(lapply(y, is.null))] <- list(c(NA, NA))
z <- t(matrix(unlist(y), byrow = F, nc = length(y)))
z[is.na(z)] <- "--"
out$Rows <- z[, 1]
out$Cols <- z[, 2]
}
else {
### otherwise it's a matrix which will work differently (not checked in R version)
out$Rows <- y[1, ]
out$Cols <- y[2, ]
}
}
else out$Cols <- out$Rows <- rep("--", length(out$Object)) #
out$Len <- apply(objs.mat, 1, function(X)
length(get(X))) # browser()
out$Date <- xx$Date
dfobj <- apply(objs.mat, 1, function(X)
is.data.frame(get(X))) #
out$Mode[dfobj] <- "dataframe"
out.df <- unfactor(as.data.frame(lapply(out, function(X)
X[seq(max)])), c("Object", "Mode"))
df.to.file(out.df, row.nos = T)
}
, comment = "25/07/2006")
"dmodes.f" <-
structure(function(t.ob = NULL, pos = 1, max = NULL, date.sort = TRUE)
{
### Purpose: lists only the functions
### ----------------------------------------------------------------------
### Modified from: dmodes
### ----------------------------------------------------------------------
### Arguments:
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 10 Sep 2002, 09:55
### Returns a Data frame showing modes of objects in the stated directory.
### Can be sorted by Mode, Length, or Date in addition to default Object
### A variation on this one is modes(). Made for lazy typists to sort by object name
### max is the maximum number of objects returned; useful if ls() is long and
### you're only interested in the first few
if(is.null(t.ob)) ls.o <- objects(pos = pos, all.names =FALSE)
else ls.o <- objects(pos = pos, pattern = t.ob, all.names =FALSE)
if(!is.null(max) && max > length(ls.o))
print(paste("Only", length(ls.o), "match", t.ob, "in", search()[pos]))
if(is.null(max) || (max > length(ls.o))) max <- length(ls.o) #
### Rearrange date information
date.chr <- NULL
for(i in seq(ls.o)) {
date.i <- comment(get(ls.o[i], pos = pos))
if(is.null(date.i))
date.i <- "NA/NA/NA"
date.chr[i] <- date.i
}
break.at <- occurrence(date.chr, "/", 1:2)
if(is.null(dim(break.at)))
break.at_matrix(break.at,nrow=1)
first <- break.at[, 1] - 1
second <- break.at[, 2] - 1
day <- as.numeric(substring(date.chr, 1, first))
month <- as.numeric(substring(date.chr, first + 2, second))
year <- as.numeric(substring(date.chr, second + 2))
ob.df <- unfactor(data.frame(Object = ls.o, Day = day, Mon = month, Year =
year))
ob.df <- df.sort(ob.df, "Object", backwards = TRUE)
if(all(is.na(ob.df$Day)))
date.sort <- FALSE ## no point trying to sort
if(date.sort) {
xx <- df.sort(ob.df, c("Day","Mon"))
xx <- df.sort(xx, "Year", backwards = TRUE)[seq(max), ]
}
else xx <- df.sort(ob.df, "Object")
xx$Date <- date.chr[as.numeric(row.names(xx))]
options(warn = -1)
on.exit(options(warn = 0))
out <- list()
dummy <- matrix(nc = 3, nr = 3) # a non-null last list element
out$Object <- c(xx$Object, "dummy")
objs.mat <- cbind(out$Object)
out$Mode <- apply(objs.mat, 1, function(X)
mode(get(X))) #
### Dimension information of dataframes and matrices
dimensions <- apply(objs.mat, 1, function(X)
dim(get(X))[1:2])
y <- dimensions
if(!is.null(y)) {
if(is.list(y)) {
y[unlist(lapply(y, is.null))] <- list(c(NA, NA))
z <- t(matrix(unlist(y), byrow = F, nc = length(y)))
z[is.na(z)] <- "--"
out$Rows <- z[, 1]
out$Cols <- z[, 2]
}
else {
### otherwise it's a matrix which will work differently (not checked in R version)
out$Rows <- y[1, ]
out$Cols <- y[2, ]
}
}
else out$Cols <- out$Rows <- rep("--", length(out$Object)) #
out$Len <- apply(objs.mat, 1, function(X)
length(get(X))) # browser()
out$Date <- xx$Date
dfobj <- apply(objs.mat, 1, function(X)
is.data.frame(get(X))) #
out$Mode[dfobj] <- "dataframe"
out.df <- unfactor(as.data.frame(lapply(out, function(X)
X[seq(max)])), c("Object", "Mode"))
out.df <- out.df[out.df$Mode == "function", ]
df.to.file(out.df, row.nos = T)
}
, comment = "10/09/2002")
"errorbars" <-
structure(function(x, y, se, yl = y - se, yu = y + se, eps, ...)
{
if(any(is.na(yl)) | any(is.na(yu)))
text(x[is.na(yl) | is.na(yu)], y[is.na(yl) | is.na(yu)], "NA", adj = 1)
segments(x, yl, x, yu, ...)
segments(x - eps, yl, x + eps, yl, ...)
segments(x - eps, yu, x + eps, yu, ...)
}
, comment = "31/06/1996")
"ext.range" <-
structure(function(x, a = 5, ...)
{
## Purpose: Extend the xlim or ylim range of a vector
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments: a is % to extend the range of vector x
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 13 Jan 2006, 09:11
ra <- range(x, ...)
interval <- diff(ra)
adjustment <- interval * a/100
ra + adjustment * c(-1, 1)
}
, comment = "13/01/2006")
"factorize" <-
structure(function(n)
{
p <- n/(z <- 1:ceiling(sqrt(n)))
z <- z[trunc(p) == p]
unique(c(z, rev(n/z)))
}
, comment = "14/12/1998")
"fieller" <-
structure(function(phat, a, b, v11, v12, v22, df.t = Inf, link = "logit")
{
# fiducial interval calculation from Finney p.78
#
# a,b estimated intercept and slope
# v11, v12, v22 variance and covariance of estimates of a and b
# t t statistic
#
if(df.t == Inf) tt <- 1.96 else tt <- qt(0.975, df.t)
g.link <- link.function(link)
a <- g.link(phat) - a
v12 <- v12* -1
m <- a/b
g <- (gg <- (tt/b)^2) * v22
if(g > 0.99 || g < 0) {
return(list(xhat = m, ci = NULL, lower = 1, upper = -1, g = g))
}
xhat0 <- m + g/(1 - g) * m - gg/(1 - g) * v12
vv <- v11 - 2 * m * v12 + m^2 * v22 - g * v11 + gg * v12^2
Ix <- tt/b/(1 - g) * sqrt(vv)
Ix <- abs(Ix)
list(xhat = m, var = vv, lower = xhat0 - Ix, upper = xhat0 + Ix, g = g)
}
, comment = "14/09/2004")
"fill.down" <-
structure(function(x)
{
# Replaces NAs in a numeric vector or blanks in a character one
# with the value preceding it.
# Too lazy to use fill-down in Excel
#
# Warning: Replaces an initial NA with a zero. Might not always be appropriate.
if(is.na(x[1])) x[1] <- 0
stillsome <- T
while(stillsome) {
nas.at <- iff(is.character(x), x == "", is.na(x))
na.indx <- seq(along = x)[nas.at]
x[na.indx] <- x[na.indx - 1]
stillsome <- length(na.indx) > 0
}
x
}
, comment = "16/08/1999")
"find.first" <-
structure(function(x, char = " ")
{
### Finds position of first occurence of a particular character in a character string
if(!is.character(x)) stop(
"find.first() applicable only for character strings")
beg <- NULL
for(k in 1:length(x)) {
leg.vec <- paste(substring(x[k], 1:nchar(x[k]), 1:nchar(x[k])), sep = ""
)
beg[k] <- match(char, leg.vec)
}
beg
}
, comment = "01/06/1996")
"find.pos" <-
structure(function(y, inn)
{
#
# Finds positions of "y" in a vector "inn"
x <- match(inn, y)
seq(along = x)[!is.na(x)]
}
, comment = "30/04/2001")
"find.pos.mat" <-
structure(function(mat = b, test = (b > 1))
{
# Finds positions in matrix of elements which meet test
want <- mat[test]
if(length(want) == 0)
cat("No matches\n")
else pass <- cbind(row(mat)[test], col(mat)[test])
dimnames(pass) <- list(paste("Soln:", 1:dim(pass)[1], sep = ""), c("Row",
"Col"))
pass
}
, comment = "13/05/1997")
"fitconf" <-
structure(function(link = "cloglog", mins, dead, tot, cutoff = NULL, offset, j, leg,
interval = F, pc = c(line = 50, line = 99), pc.spline = NULL, pc.monotone
= c(line = 50, line = 99), span.mono = 1, plot.spline = F, plot.monotone
= NULL, cutoff.mono = 0, xfun = function(x)
x, cm = NULL, numcm = NULL, p.min = NULL, printdat = F, save.resid = F,
cm.strategy = "adjust.later", cm.code = 0, cm.allcodes = NULL, xfun.inv =
NULL, full.robust = F)
{
# link is link function for binomial model
# mins is time of sampling
# dead is number of dead
# total is total number
# cutoff - data at less than this time are not included in the fitting
# offset is the amount added to time to improve linearity
# j is the dataset index number
# log time is used in calculations
# calculations are done using mean zero data
#
# nmid is the number of points suitable for fitting
if(is.null(xfun.inv)) xfun.inv <- function(x)
x
if(all(xfun(exp(1:5)) == (1:5))) {
takelog <- T
xfun.inv <- exp
}
else takelog <- F
#browser()
if(all(xfun(1:8) == (1:8)))
xfun.inv <- function(x)
x
if(all(xfun(1:8) == (1:8)^2)) xfun.inv <- sqrt
# link choices are "identity", "log", "logsurv", "logit", "sqrt",
# "inverse", "probit", "cloglog", "loglog"
inv.fun <- link.inverse(link) # get inverse of link function
link.fun <- link.function(link) # get link function
if(any(dead > tot))
stop(paste("\none or more totals is less than no of dead\n", "Dead:",
paste(dead, collapse = " "), "\nTotal:", paste(tot, collapse = " ")))
mins <- as.character(mins)
cmrows <- mins == cm.code[1]
if(sum(cmrows) == 0 & length(cm.code) > 1) {
cmrows <- mins == cm.code[2]
second.cm <- T
}
else second.cm <- F
x.all <- array(, length(mins))
x.all <- array(, length(mins))
x.all[cmrows] <- - Inf
cutit <- !cmrows & tot > 0
if(!is.null(cm.allcodes)) {
other.rows <- as.logical(match(mins, cm.allcodes, nomatch = 0))
if(any(other.rows))
cutit[other.rows] <- F
x.all[!cutit & !cmrows] <- NA
x.cm <- mins[other.rows]
dead.cm <- dead[other.rows]
tot.cm <- tot[other.rows]
}
else x.cm <- NULL
interp.rows <- cutit | mins == "0"
if(is.null(cutoff.mono))
cutoff.mono <- min(c(0, x.all[cutit]))
mono.rows <- cutit | mins == "0"
x.all[mono.rows] <- as.numeric(mins[mono.rows])
if(any(is.na(x.all[cutit])))
stop(paste("Illegal code", mins[cutit][is.na(x.all[cutit])],
"has appeared."))
mono.rows[mono.rows][x.all[mono.rows] < cutoff.mono] <- F
if(!is.null(pc.monotone) | !is.null(plot.monotone))
monospeak <- paste("Monotone fit: Include x >= ", cutoff.mono, ":", sep
= "")
else monospeak <- ""
if(!is.null(cutoff)) {
cutit[cutit] <- x.all[cutit] >= cutoff
cutspeak <- paste("Line: Include x >= ", cutoff, ":", sep = "")
}
else cutspeak <- ""
if(!is.null(p.min)) {
p <- dead/tot
omitspeak <- paste("Omit until p > ", format(round(p.min, 2)), ":", sep
= "")
testseq <- seq(along = x.all)[p < p.min]
if(length(testseq) > 0)
omit <- 1:max(testseq)
else omit <- 0
}
else {
omit <- 0
omitspeak <- ""
}
cat("\n", j, " ", leg, fill = T)
cutit[omit] <- F
if(is.null(numcm)) numcm <- 0 # numcm may be reset below
if(!is.null(cm)) {
if(numcm == 0)
cmspeak <- paste(" cm (taken as fixed) =")
else cmspeak <- " cm (supplied) ="
cmspeak <- paste(cmspeak, format(round(cm, 3)))
}
else {
if(sum(cmrows) > 0) {
if(any(x.all[cmrows] != - Inf & x.all[cmrows] != 0))
cat("\n*** Warning: Fault in specification of cm rows ***", "\n")
dead0 <- sum(dead[cmrows])
tot0 <- sum(tot[cmrows])
cm <- dead0/tot0
numcm <- tot0
n.cm <- sum(cmrows)
cmspeak <- paste(" cm(obs) =", format(round(cm, 3)), "(from", n.cm,
paste("point", switch((n.cm > 1) + 1,
"",
"s",
), ")", sep = ""))
if(second.cm)
cmspeak <- paste(cmspeak, " cm code (2nd choice) =", cm.code[2])
}
else cmspeak <- "No information is available on cm"
}
if(!is.null(pc.spline)) {
pfrac.interp <- pc.spline/100
yhat.pred <- link.fun(pfrac.interp)
dead.interp <- dead[interp.rows]
tot.interp <- tot[interp.rows]
p.obs <- (dead.interp/tot.interp)
time.interp <- as.numeric(mins[interp.rows])
time1 <- max(time.interp[p.obs == 0])
if(is.na(time1))
time1 <- min(time.interp)
time2 <- min(time.interp[p.obs == 1])
if(is.na(time2))
time2 <- max(time.interp)
time.here <- (time.interp >= time1) & (time.interp <= time2)
dead.interp <- dead.interp[time.here]
tot.interp <- tot.interp[time.here]
p.interp <- (dead.interp + 1/6)/(tot.interp + 1/3)
x.interp <- time.interp[time.here]
}
##
if(!is.null(pc.monotone)) {
time.here <- as.numeric(mins[mono.rows])
ord.x <- order(time.here)
subs.here <- (1:length(mono.rows))[mono.rows][ord.x]
dead.mono <- dead[subs.here]
tot.mono <- tot[subs.here]
p.obs <- (dead.mono/tot.mono)
time.mono <- time.here[ord.x]
}
##
ld.dp <- 2 - floor(log(max(x.all[cutit], na.rm = T))/log(10))
ld.dp <- max(0, ld.dp)
if(is.na(ld.dp))
ld.dp <- 1
if(nchar(monospeak) > 0)
cat(monospeak, "\n")
if(nchar(cutspeak) + nchar(omitspeak) > 0)
cat(cutspeak, " ", omitspeak)
cat(sum(cutit), "points remain_")
cat(cmspeak, "\n")
if(takelog) {
x.trt <- log(x.all[cutit] + offset)
}
else {
x.trt <- xfun(x.all[cutit] + offset)
}
xbar <- mean(x.trt)
x.explan <- x.trt - xbar
dead2 <- dead[cutit]
tot2 <- tot[cutit]
nmid <- length(x.trt[dead2 < tot2])
lt <- rep(-999, length(pc))
b0 <- -999
b1 <- -999
selog50 <- -999
selog99 <- -999
cept <- NA
dev <- NA
df <- NA
dev.robust <- NA
df.robust <- NA
disp <- NA
resid <- NULL
resid.dev <- NULL
lt.se <- rep(NA, length(pc))
cm.est <- cm
cm.n <- c(cm, numcm)
if(!is.null(pc.monotone) && (length(time.mono) >= 2)) {
if(link == "loglog")
link <- "cloglog"
lt.monotone <- pool.adj(time.mono, dead.mono, tot.mono, phat = pc.monotone,
cm = cm, cm.strategy = cm.strategy, cm.allcodes =
cm.allcodes, plotit = plot.monotone, link = link,
xfun = xfun, legend = leg, x.cm = x.cm,
dead.cm = dead.cm, tot.cm = tot.cm, span =
span.mono)
note <- c("( loess)", "( monotone line)")[names(lt.monotone)]
cat(paste(" LT", pc.monotone, ":", format(round(lt.monotone, ld.dp)),
sep = ""), note, "\n")
}
else lt.monotone <- NULL
if(!is.null(pc.spline)) {
lt.spline <- interp.curve(x.interp, p = p.interp, pval = pfrac.interp,
plot.interp = plot.spline, fun = link.fun)
cat(paste(" LT", pc.spline, ":", , sep = ""),
format(round(lt.spline, ld.dp)), "(spline) ", "\n")
}
else lt.spline <- NULL
if(nmid == 2)
full.robust <- F
### browser()
if(nmid >= 2) {
u <- likelihood.ci(cm.n = cm.n, x.explan, xbar, dead2, tot2, link = link,
xfun.inv = xfun.inv, interval = interval, pc = pc,
save.resid = save.resid, cm.strategy = cm.strategy,
full.robust = full.robust)
lt <- u$ld
ld.dp <- 2 - floor(log(lt[length(lt)])/log(10))
if(is.na(ld.dp))
ld.dp <- 1
ld.dp <- max(ld.dp, 0)
b0 <- u$coef[1]
b1 <- u$coef[2]
cept <- b0 - b1 * xbar
lt.se <- u$selog
resid <- u$resid
resid.dev <- u$resid.dev
dev <- u$dev
df <- u$df
dev.robust <- u$dev.robust
df.robust <- u$df.robust
disp <- u$dispersion
linpred <- b0 + b1 * x.explan # comprob <- inv.fun(linpred)
### fitr <- (1 - comprob) * tot2
### fitc <- comprob * tot2
if(printdat) print(cbind(x.all[cutit], dead2, tot2))
### print(cbind(x.all[cutit], x.explan, fitr, dead2, tot2))
### if(min(fitr) < 1) {
### cat("Warning. Some expected numbers 0f dead are < 1",
### fill = T)
### cat("Exp. No. Dead:", format(round(fitr, 3)), fill = T)
###
### }
### if(min(fitc) < 1) {
### cat(
### "Warning. Some expected numbers of survivors are < 1",
### fill = T)
### cat("Exp. Survivors :", format(round(fitc, 3)), fill
### = T)
### }
}
else {
### Estimates of LT50 and LT90 when too few data points are available
###
if(printdat) print(cbind(x.all, dead, tot, cutit))
cat(nmid, " data points give mortalities < 1.0: ")
lt <- NA
cat("No estimate is available", "\n")
}
if(!save.resid) tot2 <- NULL
### Keep totals only if required along with residual information
list(lt = lt, se = lt.se, b0 = cept, b1 = b1, lt.spline = lt.spline,
lt.monotone = lt.monotone, cm = cm.est, dev = dev, df = df, dev.robust
= dev.robust, df.robust = df.robust, dispersion = disp, total = tot2,
resid = resid, resid.dev = resid.dev, times = x.trt)
}
, comment = "14/09/2004")
"flex.contour" <-
structure(function(df = cont.test.df[1:273, (1:3)], transpose = F, label.pos = (13:18)/20,
levs.pc = c(50, 75, 90, 95, 99), lev.cex = 0.6, lev.font = 1, mark = F,
alone = F, ztransform = "logit")
{
# Flexibility in positioning the contour-level labels is the name of the game.
#
# df is dataframe such as produced by gee.
# Requires also is.not.na() in Gems & make.line.loess()
# Not generic yet. Uses only logit transformed % data
xxs <- unique(df[, 1])
yys <- unique(df[, 2])
zzs <- df[, 3]
zzs.mat <- matrix(zzs, nr = length(is.not.na(xxs)))
x.use <- iff(transpose, yys, xxs)
y.use <- iff(transpose, xxs, yys)
if(transpose)
zzs.mat <- t(zzs.mat)
ztransform.fn <- get(ztransform)
lev.pos <- ztransform.fn(levs.pc/100) # (levs.pc are percentages)
lev.txt <- paste(levs.pc)
cont.fit <- contour(x.use, y.use, zzs.mat, save = T, levels = lev.pos)
if(alone)
plot(range(x.use), range(y.use), type = "n")
if(length(label.pos < length(levs.pc)))
label.pos <- rep(label.pos, length(levs.pc))
for(i in 1:length(cont.fit)) {
make.line.loess(cont.fit[[i]], lab = lev.txt[i], mark = mark, label.pos
= label.pos[i], cex = lev.cex, lev.font = lev.font)
}
}
, comment = "03/06/1998")
"flyplot" <-
structure(function(mn = 1:12, ri = 3, cj = 4, line = T, cutoff.mono = NULL, data = NULL,
datafun = NULL, choice = 1, plimits = c(0.1, 0.9995), link = "cloglog",
range.strategy = "page", fit.trend = "line", pc = c(line = 50, line = 99),
span.mono = 1, n.strategy = "mn", xlabels = NULL, xfun = function(x)
x, takelog = F, before = NA, new.page = T, mkh = 0.025, cex = NULL, maint.line
= 1, title.line = 0.25, title.space = 0, xtitle = NULL, ytitle = NULL, ext
= NULL, lt.ld = "LT", datestamp = T, byrow = T, mtext.line = 2)
{
### Title for data subset j is maint
### Prints graphs 4 to a page, starting with graph for data subset j=i
### dset holds all the data
### values of idset (1, 2, ...) identify subsets of dset
### range.strategy may be "individual", or "page", or "all"
### fit.trend may be line, or mono(tone), or mono.line
### mono.line gives monotone fit, then line
### Set pc to decide which LTs will appear on the margin of the plot.
##
options(width = 200)
on.exit(options(width = 80))
trend.4 <- substring(fit.trend, 1, 4)
do.abers <- trend.4 == "mono" | trend.4 == "loes"
oldpar <- par()
on.exit(par(oldpar), T) #
### Set up different par parameters for postscript and motif devices.
if(names(dev.cur()) == "postscript") {
omi <- c(10, 12 + title.space/3, 12 + title.space, 12 + title.space/3)/25.4
mgp <- c(1.75, 0.7, 0)
mai <- c(8, 10, 5, 1)/25.4
}
else {
omi <- c(5, 5 + title.space/3, 7 + title.space, 10 + title.space/3)/25.4
mgp <- c(2.25, 1, 0)
mai <- c(15, 20, 10, 0)/25.4
}
if(new.page & byrow)
par(mai = mai, mgp = mgp, mfrow = c(ri, cj), omi = omi)
else if(new.page & !byrow)
par(mai = mai, mgp = mgp, mfcol = c(ri, cj), omi = omi)
if(is.null(cex) & (ri * cj == 6))
cex <- 0.65
if(!is.null(cex))
par(cex = cex)
final.choice <- choice[length(choice)]
if(!is.null(data)) {
ab <- data[[final.choice]]
leg.ab <- ab$leg.brief
if(is.null(leg.ab))
leg.ab <- ab$legend
choice.call <- ab$choice
if(missing(takelog) & missing(xfun)) {
takelog <- ab$takelog
if(takelog)
xfun <- log
}
}
else {
ab <- NULL
leg.ab <- NULL
choice.call <- choice
}
if(missing(datafun) & !is.null(ab))
datafun <- get(ab$datafun)
if(!is.null(ext))
u <- datafun(choice.call, ext = ext)
else u <- datafun(choice.call)
if(is.null(leg.ab))
leg.ab <- u$leg.brief
if(!is.null(pc)) {
lt.txts <- paste(pc)
lt.prefix <- paste(lt.ld, pc, "=", sep = "")
lt.strategy <- names(pc)
}
else {
lt.strategy <- NULL
lt.txts <- NULL
}
if(missing(ext))
ext <- NULL
if(!is.null(ab) & !is.null(lt.strategy)) {
first.one <- switch(lt.strategy[1],
loess = ,
monotone = ,
mono = "lt.monotone",
line = "lt")
ltmat.leg <- matrix(, dim(ab[[first.one]])[1], length(pc))
for(k in seq(along = lt.strategy)) {
which.lt <- switch(lt.strategy[k],
loess = ,
monotone = ,
mono = "lt.monotone",
line = "lt")
ltmat.leg[, k] <- ab[[which.lt]][, lt.txts[k]]
}
dp <- max(c(3 - ceiling(log(quantile(ab[[first.one]], probs = 0.75,
na.rm = T))/log(10)), 0))
}
else ltmat.leg <- NULL
if(!is.null(ltmat.leg)) {
if(is.null(leg.ab))
leg.ab <- rep("", dim(ltmat.leg)[1])
leg.lt <- apply(ltmat.leg, 1, function(x, z, u, dp)
paste(z, round(x, dp), collapse = "; ", sep = ""),
z = lt.prefix, dp = dp)
legend <- paste(leg.ab, " [", leg.lt, "]", sep = "")
}
else if(!is.null(leg.ab))
legend <- leg.ab
else if(length(lt.txts) > 0)
cat(paste(paste("LT", lt.txts, sep = ""), collapse = ""),
"values are not available, perhaps because pc was not given correctly\n"
)
cutx <- u$cutx
offset <- u$offset
cm <- NULL
if(!is.null(ab)) {
cm.code <- ab$cm.code
cm <- ab$cm
cm.strategy <- ab$cm.strategy
cmtot <- u$cmtot
}
else {
cm.code <- u$cm.code
cm <- u$cm
cmtot <- u$cmtot
cm.strategy <- u$cm.strategy
}
if(is.null(cm.strategy))
cm.strategy <- "adjust.later"
if(is.null(cm.code))
cm.code <- 0
cm.allcodes <- u$cm.allcodes
if(is.null(cm.allcodes))
cm.allcodes <- cm.code
times.all <- as.character(u$time)
tot.all <- u$tot
other.rows <- as.logical(match(times.all, cm.allcodes, nomatch = 0))
keep <- rep(T, length(other.rows))
if(any(is.na(as.numeric(times.all[!other.rows])))) {
cat("\n*** NA(s) appear among non-control times ***\n")
numit <- times.all[!other.rows]
look <- is.na(as.numeric(numit))
print(data.frame(times = numit[look], totals = tot.all[!other.rows][look
], dead = u$dead[!other.rows][look]))
keep[!other.rows][look] <- F
}
if(any(is.na(tot.all) | tot.all == 0)) {
look <- is.na(tot.all) | tot.all == 0
cat("\nNAs or zeros appear among totals\n")
print(data.frame(times = times.all[look], total = tot.all[look], dead =
u$dead[look]))
keep[look] <- F
}
time.txt <- times.all[keep]
times <- array( - Inf, length(time.txt))
use.trt <- !other.rows[keep] | time.txt == "0"
times[use.trt] <- as.numeric(time.txt[use.trt])
xmin <- min(times[use.trt])
idset <- u$id[keep]
tot <- u$total[keep]
dead <- u$dead[keep]
if(missing(xtitle))
xaxtitle <- u$xaxtitle
else xaxtitle <- xtitle
treat <- u$treat
maint <- u$maint
if(is.null(xlabels))
xxlabels <- u$xlabels
else xxlabels <- xlabels
nu <- unique(idset)
nn <- match(nu, idset) # bigstage <- dset[, stagecol]
n <- ri * cj
if(max(mn) > length(nu)) {
maxj <- length(nu)
cat("You have specified plots up to plot no.", max(mn),
"\nwhereas data is available for", maxj, "plots only.\n")
mn <- mn[mn <= maxj]
}
graphnum <- mn
maxj <- max(mn)
if(range.strategy != "individual")
common.range <- switch(range.strategy,
page = mn,
all = 1:length(nu))
todo <- length(nu) - maxj
if(todo > 0)
cat("Datasets", maxj + 1, "to", length(nu), "will remain after this.",
fill = T)
xmin <- min(times)
if(takelog)
xmin <- min(times[times + offset > 0])
xrange <- switch(cm.strategy,
adjust.later = range(times),
abbott = range(times))
ij <- 0
if(range.strategy != "individual") {
xrm <- matrix(, 2, length(common.range))
for(k in common.range) {
ij <- ij + 1
timtim <- times[idset == nu[k]]
xrm[, ij] <- switch(cm.strategy,
adjust.later = range(timtim),
abbott = range(timtim))
}
xrange <- range(xrm)
}
if(!missing(before))
xrange[2] <- before
ij <- 0
cmk <- NULL
for(k in graphnum) {
leg.k <- legend[k]
ij <- ij + 1
j <- nu[k]
take <- idset == j
timesj <- times[take]
deadj <- dead[take]
totj <- tot[take]
time.txj <- time.txt[take]
cmrows <- time.txj == cm.code[1]
second.cm <- F
if(sum(cmrows) == 0 & length(cm.code) > 1) {
cmrows <- time.txj == cm.code[2]
second.cm <- T
}
else second.cm <- F
if(!is.null(cmtot))
numcm <- cmtot[k]
else numcm <- NULL
if(is.null(numcm)) numcm <- 0 # numcm may be reset below
if(is.null(cm)) {
if(sum(cmrows) > 0) {
if(any(timesj[cmrows] != - Inf & timesj[cmrows] != 0))
cat("\n*** Warning: Fault in specification of cm rows ***",
"\n")
dead0 <- sum(deadj[cmrows])
tot0 <- sum(totj[cmrows])
cmk <- dead0/tot0
numcm <- tot0
n.cm <- sum(cmrows)
cmspeak <- paste(" cm(obs) =", format(round(cm, 3)), "(from",
n.cm, paste("point", switch((n.cm > 1) + 1,
"",
"s",
), ")", sep = ""))
if(second.cm)
cmspeak <- paste(cmspeak, " cm code (2nd choice) =", cm.code[2
])
}
else cmspeak <- "No information is available on cm"
}
else {
# !is.null(cm)
cmk <- cm[k]
if(cm.strategy == "abbott") {
cat("Take cm =", format(round(cmk, 3)), "\n")
take <- take & !other.rows
}
}
if(do.abers) {
use.mono <- use.trt[take]
if(!is.null(cutoff.mono))
use.mono[timesj[use.mono] < cutoff.mono] <- F
time.here <- timesj[use.mono]
ord.x <- order(time.here)
subs.here <- (1:length(use.mono))[use.mono][ord.x]
dead.mono <- deadj[subs.here]
tot.mono <- totj[subs.here]
p.mono <- (dead.mono/tot.mono)
time.mono <- time.here[ord.x]
other.take <- other.rows[take]
x.cm <- time.txt[take][other.take]
dead.cm <- deadj[other.take]
tot.cm <- totj[other.take]
if(length(time.mono) < 2)
do.abers <- F
}
if(do.abers) {
if(link == "loglog")
link <- "cloglog"
lt.mono <- pool.adj(time.mono, dead.mono, tot.mono, phat = NULL, cm
= cmk, cm.strategy = cm.strategy, cm.allcodes = cm.allcodes,
plimits = plimits, xtit = xaxtitle, plotit = fit.trend, link =
link, xfun = xfun, legend = leg.k, x.cm = x.cm, dead.cm = dead.cm,
tot.cm = tot.cm, span = span.mono)
}
if(range.strategy == "individual")
xrange <- switch(cm.strategy,
adjust.later = ,
abbott = range(times[take]))
if(diff(range(xxlabels)) > 0)
xlabels <- xxlabels
else xlabels <- pretty(range(times[take]))
if(line) {
ltab <- ab
b <- ltab[["slope"]][k]
a <- ltab[["intercept"]][k]
}
else {
a <- -999
b <- -999
}
startx <- max(xmin, cutx[k])
jk <- switch(n.strategy,
id = j,
mn = k)
if(!do.abers) {
#
simplot(time.txt[take], deadj, totj, link, cm = cmk, cm.strategy =
cm.strategy, xtit = xaxtitle, ytit = ytitle,
main = paste(jk, ": ", "", leg.k, sep = ""), xfun = xfun,
takelog = takelog, line = line,
ab = c(a, b), clip = c(startx, max(timesj)), xlimits = xrange,
xlabels = xlabels, offset = offset, plimits = plimits, mkh = mkh,
title.line = title.line, id = NULL)
}
### if(names(dev.cur()) == "postscript" & nchar(maint) > 0)
### mixed.mtext(texts = maint, side = 3, outer = T, cex = 1.25, line =
### maint.line, adj = 0.5)
### else
mtext(maint, 3, outer = T, cex = 1.4, font = 2, line = maint.line)
}
if(cm.strategy == "abbott")
mtext("Data is plotted following Abott's adjustment", 3, line =
maint.line - 1, outer = T, cex = 0.5)
lt.leg <- paste(paste("lt", pc, ":", names(pc), sep = ""), collapse = "; ")
par(cex = 0.4)
if(datestamp) {
r.txt <- paste(lt.leg, " ", system("date +%d/%m/%y' '%H:%M", TRUE))
l.txt <- system("pwd", TRUE) ## get call to work sometime, deparse(match.call())
mtext(r.txt, side = 1, adj = 1, outer = T, line = mtext.line, cex = .6)
mtext(l.txt, side = 1, adj = 0, outer = T, line = mtext.line, cex = .6)
}
}
, comment = "03/04/2006")
"formula.variables" <-
structure(function(myformula, unwanted = NULL)
{
ops <- c("~", "+", "-", "*", "|", "I", "(", ")", "sqrt", "log", "exp",
"logit", "cloglog", "ang", unwanted)
out <- unlist(myformula)
out[!(out %in% ops)]
}
, comment = "10/10/1997")
"gem.df" <-
structure(list(Object = structure(as.integer(c(26, 65, 67, 89,
80, 60, 25, 94, 10, 9, 99, 73, 53, 45, 39, 77, 41, 61, 19, 47,
76, 75, 74, 31, 5, 63, 43, 21, 42, 37, 17, 64, 55, 88, 98, 28,
97, 6, 86, 2, 4, 93, 72, 14, 59, 36, 24, 85, 70, 69, 29, 57,
90, 83, 20, 82, 16, 38, 12, 11, 84, 18, 46, 1, 44, 23, 34, 13,
7, 95, 91, 52, 58, 56, 87, 71, 51, 96, 8, 78, 40, 22, 54, 62,
27, 32, 68, 15, 33, 3, 48, 35, 66, 81, 49, 92, 50, 79, 30)), .Label = c("add.text",
"agg", "ang.bt", "append.cols", "average.down", "bar.legs", "blank.plot",
"boxplot.fn", "catalogue", "catalogue.gems", "check.fact", "check.na",
"clg", "clipline", "cloglog", "cloglog.bt", "cols.for.genstat",
"comb", "complot", "count.nas", "data.spot", "df.mat", "df.sort",
"df.summary", "df.to.file", "dmodes", "errorbars", "factorize",
"FF", "fieller", "fill.down", "find.first", "find.pos", "find.pos.mat",
"fitconf", "flex.contour", "flyplot", "formula.variables", "gems.df",
"gen.write", "get.col.names", "get.formula", "get.mort.tables",
"glean.blank", "glm", "junk.df", "levels.for.genstat", "likelihood.ci",
"link.function", "link.inverse", "list.mat", "list.to.data.frame",
"logit", "logit.bt", "lookup", "make.df", "make.factors", "make.id",
"make.line.loess", "MASS", "mean.lt", "mean.lt.short", "month.length",
"occurrence", "pid", "plot.lt", "PLS", "pool.adj", "ppaste",
"pseudo.f", "p.sunflowers", "qik.sum", "read.table", "relevel",
"relevel.default", "relevel.factor", "reorder", "rms", "robust.deviance",
"save.new", "sem", "set.screens", "show.nas", "simplot", "s.tapply",
"stdev", "strip.trail", "summarize", "SURV", "tab.file", "table.all",
"tab.mat", "text.bp", "text.lot", "time.to.mort", "tuk.calc",
"unfactor", "Updown", "write.function"), class = "factor"), Mode = structure(as.integer(c(2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)), .Label = c("dataframe",
"function", "numeric"), class = "factor"), Rows = structure(as.integer(c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Label = c("--",
"137", "15"), class = "factor"), Cols = structure(as.integer(c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Label = c("--",
"13", "8"), class = "factor"), Len = c(1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 12, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 13, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1), Date = structure(as.integer(c(85,
66, 18, 76, 76, 53, 48, 10, 10, 6, 81, 78, 78, 78, 75, 30, 33,
61, 60, 50, 8, 8, 8, 41, 22, 82, 32, 47, 44, 28, 28, 71, 71,
65, 46, 36, 64, 9, 59, 37, 74, 49, 34, 77, 7, 7, 4, 26, 43, 25,
16, 11, 70, 5, 5, 80, 67, 23, 58, 58, 55, 63, 29, 27, 19, 45,
31, 17, 12, 69, 52, 52, 14, 2, 62, 35, 1, 68, 13, 83, 38, 51,
57, 54, 86, 3, 84, 40, 24, 39, 79, 79, 15, 21, 20, 73, 42, 72,
56)), .Label = c("01/02/1997", "01/04/1997", "01/06/1996", "01/06/1998",
"01/12/1997", "02/04/2001", "03/06/1998", "03/09/1999", "03/10/1998",
"04/04/2001", "06/01/1998", "06/05/1997", "06/12/1996", "07/04/1997",
"08/10/1995", "09/01/1998", "09/05/1997", "09/05/2001", "09/07/1997",
"09/08/1995", "09/09/1995", "10/08/1999", "10/10/1997", "11/01/1996",
"11/02/1998", "11/04/1998", "11/07/1997", "12/04/1999", "12/07/1997",
"12/07/2000", "13/05/1997", "13/05/1999", "13/06/2000", "13/08/1998",
"14/02/1997", "14/12/1998", "15/09/1998", "15/11/1996", "15/12/1995",
"16/04/1996", "16/08/1999", "16/10/1994", "17/02/1998", "17/04/1999",
"17/05/1997", "17/12/1998", "18/04/1999", "18/04/2001", "18/08/1998",
"18/09/1999", "18/10/1996", "19/04/1997", "20/04/2001", "20/08/1996",
"20/09/1997", "21/03/1994", "21/08/1996", "21/09/1997", "21/09/1998",
"22/04/2000", "22/05/2000", "23/03/1997", "23/07/1997", "23/11/1998",
"23/12/1998", "24/05/2001", "24/11/1997", "24/12/1996", "25/04/1997",
"25/12/1997", "26/01/1999", "26/06/1994", "26/06/1995", "26/08/1998",
"27/03/2001", "27/04/2001", "27/06/1998", "28/03/2001", "29/10/1995",
"29/11/1997", "30/03/2001", "30/07/1999", "30/11/1996", "31/05/1996",
"31/05/2001", "31/06/1996"), class = "factor")), .Names = c("Object",
"Mode", "Rows", "Cols", "Len", "Date"), class = "data.frame", row.names = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
"25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",
"36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46",
"47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57",
"58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68",
"69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79",
"80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90",
"91", "92", "93", "94", "95", "96", "97", "98", "99"), comment = "31/05/2001")
"gems" <-
structure(function(ww = "plot", most = 10)
{
## Purpose: searches Gems for function using dmodes
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 17 Jun 2004, 13:56
gem.pos <- grep("Gems", search())
dmodes(ww, gem.pos, most)
}
, comment = "17/06/2004")
"gen.write" <-
structure(function(df, file)
{
#To create a file readable by genstat from a dataframe
#
write.table(df, file, sep = " ", na = "*", row.names = F, col.names = F, quote = F)
}
, comment = "26/03/2003")
"get.col.names" <-
structure(function(x, patt = ".D")
{
# Gets the names of x that fit a particular pattern patt
col.names <- names(x)
y <- col.names[grep(patt, col.names)]
}
, comment = "13/06/2000")
"get.formula" <-
structure(function(ab.list = ab.second.comb, bit = "cs", link = "cloglog", rnd = 4)
{
# Gets the formula for lines from ab list
ab <- ab.list[[paste(bit)]]
legends <- ab$legend
glean <- get(ab$datafun)
cat(glean()$maint)
cat("\n\n")
out.list <- list()
for(i in 1:len(legends))
out.list[[ab$legend[i]]] <- ppaste(link, "(p) = ", round(ab$intercept[i],
rnd), " + ", round(ab$slope[i], rnd), " * t")
out.list
}
, comment = "17/04/1999")
"get.mort.tables" <-
structure(function(ab.list = ab.second.comb, bit = "cs", intervals = NULL, rnd = 4, file
= "")
{
# Does table of expected mortality for lines from ab list
ab <- ab.list[[paste(bit)]]
legends <- ab$legend
glean <- get(ab$datafun)
cat(glean(choice = bit)$maint) # browser()
cat("\n\n")
if(is.null(intervals))
intervals <- seq(0, round(max(glean(choice = bit)$time)))
out.list <- list(times = intervals)
for(i in 1:len(legends))
out.list[[ab$legend[i]]] <- round(cloglog.bt(ab$intercept[i] + ab$slope[
i] * intervals) * 100, 1)
out.df <- as.data.frame(out.list)
if(file != "") {
write(glean(choice = bit)$maint, file = file)
df.to.file(out.df, file = file, append = T)
}
else out.df
}
, comment = "13/05/1999")
"gimp" <-
structure(function(file = NULL)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 15 Apr 2002, 11:26
if(is.null(file))
system("gimp &")
else system(paste("gimp",file,"&"))
}
, comment = "15/04/2002")
"glean.blank" <-
structure(function(choice = 1)
{
blank.df <- df.sort(blank.df, rev(c("Stage", "Position", "Rep")))
attach(blank.df)
on.exit(detach("blank.df"))
idset <- make.id(Time)
cutx <- NULL
leg.brief <- unique(paste(Stage, Position, " Rep", Rep, sep = ""))
maint <- paste(degree(40), "heat treatment of some bug or other")
xlabels <- c(0, 0)
xaxtitle <- "Time (minutes)"
list(id = idset, times = Time, total = Dead + Live, dead = Dead, cutx =
cutx, offset = 0, xaxtitle = xaxtitle, maint = maint, legend = leg.brief,
xlabels = xlabels, takelog = F)
}
, comment = "09/07/1997")
"gmt" <-
structure(function(x, format = "%Y-%m-%d")
{
### Purpose: Converts a text string of yyyy-mm-dd format to POSIXct at tz = GMT
### ----------------------------------------------------------------------
### Modified from: (useful to use with ddif() to get days difference
### ----------------------------------------------------------------------
### Arguments: modify format for zillions of different date formats
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 20 Feb 2004, 09:41
x1 <- strptime(x, format = format)
as.POSIXct(x1, tz = "GMT")
}
, comment = "05/03/2004")
"gnum" <-
structure(function(file = NULL)
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 15 Apr 2002, 11:26
if(is.null(file))
system("gnumeric &")
else system(paste("gnumeric",file,"&"))
}
, comment = "28/06/2002")
"gv" <-
structure(function(file = NULL)
{
## Purpose: Starts Ghostview in the current directory
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 21 Mar 2002, 13:18
if(is.null(file))
system("gv &")
else system(paste("gv",file,"&"))
}
, comment = "28/06/2002")
"iff" <-
structure(function(test, true.result, false.result)
{
if(test)
true.result
else false.result
}
, comment = "15/07/1997")
"ipsofactor" <-
structure(function(x)
{
## Purpose: Makes a factor with levels in the current order
## ----------------------------------------------------------------------
## Arguments:
## x: character vector
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 6 Mar 2002, 10:47
y <- factor(x, levels = unique(x))
y
}
, comment = "28/06/2002")
"junk.df" <-
structure(list(Fruit = c("apple", "pear", "plum", "cherry", "apricot",
"apple", "pear", "plum", "cherry", "apricot", "apple", "pear",
"plum", "cherry", "apricot"), January = c(4, 15, 10, 11, 14,
5, 8, 18, 11, 15, 18, 14, 11, 5, 14), February = c(8, 13, 6,
5, 13, 12, 14, 9, 7, 13, 11, 9, 11, 14, 5), March = c(11, 15,
11, 13, 12, 13, 14, 15, 16, 15, 10, 16, 14, 10, 16), April = c(22,
15, 9, 13, 16, 21, 8, 18, 14, 22, 13, 18, 23, 10, 14), May = c(25,
11, 15, 18, 24, 15, 9, 11, 24, 17, 25, 23, 13, 20, 20), Week = c("A",
"A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C",
"C")), .Names = c("Fruit", "January", "February", "March", "April",
"May", "Week"), row.names = c("1", "2", "3", "4", "5", "6", "7",
"8", "9", "10", "11", "12", "13", "14", "15"), comment = "31/05/2001", class = "data.frame")
"kullas" <-
structure(function()
{
### Purpose: Basic Trellis colours for lattice plots
### ----------------------------------------------------------------------
### Modified from: col.spray in /home/hrapgc/Rstuff/garry/industry
### ----------------------------------------------------------------------
### Arguments:
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 16 Jan 2004, 13:52
lab.cex <- 1.2
list(background = list(col="transparent"),
bar.fill = list(col="#c8ffc8"),
box.rectangle = list(col = "darkgreen"),
box.umbrella = list(col = "darkseagreen", lty = 1),
box.dot = list(col = "darkorchid", cex = 1),
dot.line = list(col="#e8e8e8"),
dot.symbol = list(col="orchid", cex = .8, pch = 1),
plot.line = list(col="darkgreen"),
plot.symbol = list(col="orchid", cex = .8, pch = 1),
plot.polygon = list(alpha = 1, col="darkseagreen", border = "black",
lty = 1, lwd = 1),
par.xlab.text = list(cex = lab.cex),
par.ylab.text = list(cex = lab.cex),
axis.components = list(
left = list(tck = .6, pad1 = .5),
top = list(tck = .6, pad1 = .65),
right = list(tck = .6, pad1 = .5),
bottom = list(tck = .6, pad1 = .65)),
regions = list(col = heat.colors(100)),
reference.line = list(col="turquoise", lty = 3),
superpose.line = list(col = c("darkgreen", "red", "royalblue",
"orange", "brown", "turquoise", "orchid"),
lty = rep(1,7), lwd = rep(1, 7)),
superpose.symbol = list(pch = c(1:7), cex = rep(.9, 7),
col = c("darkgreen", "red", "royalblue", "orange", "brown",
"turquoise", "orchid"))
)
}
, comment = "03/08/2006")
"lat.set" <-
structure(function(setting = par.xlab.text, subsetting = cex, to = 1.3)
{
## Purpose: Change trellis settings in one move
## ----------------------------------------------------------------------
## Arguments:
## setting -- which setting is to be changed
## (a name of a list from trellis.settings)
## subsetting -- which part of that list
## to -- what are we changing it to
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Date: 19 Apr 2002, 10:06
y <- as.character(substitute(setting))
z <- as.character(substitute(subsetting))
x <- trellis.par.get(y)
x[[z]] <- to
trellis.par.set(y, x)
}
, comment = "07/06/2002")
"lattice" <-
structure(function()
{
### Installs grid and lattice library
library(grid, lib.loc = "/home/hrapgc/Rstuff/library/")
library(lattice, lib.loc = "/home/hrapgc/Rstuff/library/")
}
, comment = "02/10/2001")
"lerrorbars" <-
structure(function(x, y, se, yl = y - se, yu = y + se, eps, ...)
{
## Purpose: Lattice errorbars
## ----------------------------------------------------------------------
## Modified from: errorbar
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Marcus Davy , Creation date: 31 Aug 2005, 10:58
if(any(is.na(yl)) | any(is.na(yu)))
ltext(x[is.na(yl) | is.na(yu)], y[is.na(yl) | is.na(yu)], "NA",
adj = 1)
lsegments(x, yl, x, yu, ...)
lsegments(x - eps, yl, x + eps, yl, ...)
lsegments(x - eps, yu, x + eps, yu, ...)
}
, comment = "31/08/2005")
"levels.for.genstat" <-
structure(function(x, no.quotes = T)
{
# Returns a string of the factor levels to use in genstat
# paste from between " marks
# change the opening and closing " to ' if quotes are wanted
paste(unique(x), collapse = ifelse(no.quotes, ", ", "', '"))
}
, comment = "27/03/2003")
"likelihood.ci" <-
structure(function(cm.n = c(0, 0), heatval, xbar, dead, tot, offset = 0, link = "cloglog",
xfun.inv = exp, interval = T, pc = c(50, 99), dp = 1, save.resid = F,
cm.strategy = "adjust.later", full.robust = T)
{
# Uses robust glm
# heatval is the x-variable.
# xbar is a (usually centred) origin of measurement for heatval
cm <- cm.n[1]
numcm <- cm.n[2]
pfrac <- pc/100
if(cm.strategy == "adjust.later")
pval <- cm + (1 - cm) * pfrac
else pval <- pfrac
n <- length(heatval)
xvec <- heatval
link.fun <- link.function(link) # inv.fun <- link.inverse(link)
resid <- NULL
resid.dev <- NULL
if(cm.strategy == "adjust.later") {
p <- dead/tot # ymat <- cbind(dead, tot - dead)
reflect.logsurv <- 1
if(link == "logsurv") {
ymat <- cbind(tot - dead, dead)
reflect.logsurv <- -1
}
if(link == "log" | link == "logsurv")
fam <- quasi(link = log, variance = "mu(1-mu)")
else fam <- call("quasibinomial", link)
pnx <- data.frame(p.obs = p, x = heatval)
u.call <- call("glm", p.obs ~ x, data = pnx, family = fam, weight = tot)
u <- eval(u.call)
if(full.robust) {
r1 <- residuals(u, type = "deviance")
here.r <- abs(r1) != max(abs(r1))
u2.call <- call("glm", p.obs ~ x, family = fam, data = pnx[here.r, ],
weight = tot[here.r])
u2 <- eval(u2.call)
assign("startr", predict(u2, newdata = pnx), frame = 1)
u3.call <- call("glm", p.obs ~ x, family = fam, data = pnx, start =
startr, weight = tot)
u <- eval(u3.call)
}
if(any(abs(u$weights) == Inf)) {
u$weights[abs(u$weights) == Inf] <- 0
# Fudge to ensure calculations continue
}
uu <- summary.glm(u, correl = F)
if(save.resid) {
resid.dev <- residuals(u, type = "deviance") * reflect.logsurv
resid <- residuals(u, type = "response") * reflect.logsurv
}
b0 <- uu$coef[1, 1] * reflect.logsurv
b1 <- uu$coef[2, 1] * reflect.logsurv
cm.est <- NULL # phat <- inv.fun(b0 + b1 * xvec)
fit.p <- fitted(u)
if(link == "logsurv")
fit.p <- 1 - fit.p
rob.info <- robust.deviance(fit.p, tot, dead)
dev.robust <- rob.info[[1]]
df.robust <- rob.info[[2]]
dev0 <- uu$deviance
disp0 <- uu$dispersion
df0 <- uu$df[2]
}
else if(cm.strategy == "abbott") {
uu <- abbott(cmobs = cm, numcm = numcm, dose = heatval, dead = dead,
total = tot, link = link, max.iter = 12, save.resid = save.resid)
b0 <- uu$coef[1]
b1 <- uu$coef[2]
cm.est <- uu$cm
dev0 <- uu$dev
df0 <- length(heatval) - 2
if(df0 > 0)
disp0 <- dev0/df0
else disp0 <- NA
if(save.resid) {
resid.dev <- uu$resid.dev
resid <- uu$resid
}
dev.robust <- NULL
df.robust <- NULL
}
else stop(paste("Unrecognized cm strategy:", cm.strategy))
cat("Dev. =", format(round(dev0, 2)), "(df", df0, ")", " [Disp =", format(
round(disp0, 3)), "]", " [Dev =", format(round(dev.robust, 2)), "(df",
df.robust, ")]", fill = T)
if(!is.na(disp0)) {
if(disp0 < 1)
het <- 1
else het <- disp0
}
else het <- NA
if(!is.na(disp0)) {
if(disp0 > 1)
df.t <- df0
else df.t <- Inf
}
else df.t <- NA
v11 <- (uu$cov.un[1, 1]) * het
v22 <- (uu$cov.un[2, 2]) * het
v12 <- (uu$cov.un[1, 2]) * het
ld <- xfun.inv((link.fun(pval) - b0)/b1 + xbar) - offset
dp <- 2 - floor(log(ld[length(ld)])/log(10))
if(is.na(dp))
dp <- 1
dp <- max(dp, 0)
selog <- array(, length(pc))
j <- 0
if(!is.na(v22))
for(percent in pval) {
j <- j + 1
if(b1/sqrt(v22) > 0.5) {
ci <- fieller(percent, b0, b1, v11, v12, v22, df.t = df.t, link =
link)
if(is.null(ci$var))
ci$var <- NA
ld[j] <- xfun.inv(ci$xhat + xbar) - offset
selog[j] <- sqrt(ci$var)
cat(paste("LT", pc[j], ":", sep = ""), format(round(ld[j], dp)),
" ", fill = F)
if(interval)
if(ci$upper > ci$lower) {
cat("95% C.I. is", format(round(xfun.inv(ci$lower + xbar) -
offset, dp)), " to", format(round(xfun.inv(ci$upper +
xbar) - offset, dp)), fill = T)
}
else cat(" * CI could not be calculated (g =", paste(round(ci$
g, 5), ") *", sep = ""), fill = T)
}
else cat(paste("LT", pc[j], ":", sep = ""), format(round(ld[j], dp
)), " ", fill = T)
}
else {
for(i in seq(along = pc))
cat(paste(" LT", pc[i], ":", sep = ""), format(round(ld[i], dp)))
cat("\n")
}
cept <- b0 - b1 * xbar
se0 <- sqrt(v11)
se1 <- sqrt(v22)
cept <- b0 - b1 * xbar
var.cept <- v11 + v22 * xbar^2 - 2 * xbar * v12
rho.cept.b <- (v12 - xbar * v22)/sqrt(var.cept)/se1
se <- c(sqrt(var.cept), se1)
rho <- v12/(se0 * se1)
cat(" Coeffs:", format(round(c(cept, b1), 4)), "[SE", format(round(se, 3)
), " r=", format(round(rho.cept.b, 5)), "]", "\n")
list(coef = c(b0, b1), selog = selog, vcov = c(v11, v12, v22), df = df0, cm
= cm.est, dev = dev0, dev.robust = dev.robust, df.robust = df.robust,
ld = ld, resid = resid, resid.dev = resid.dev)
}
, comment = "14/12/2004")
"link.function" <-
structure(function(link)
{
link.options <- c("identity", "log", "logsurv", "logit", "sqrt", "inverse",
"probit", "loglog", "cloglog", "help")
if(link == "help")
cat("\nOptions are:", paste(link.options, sep = " "), "\n")
link.int <- charmatch(link, link.options)
if(is.na(link.int))
stop("Invalid link type")
else if(link.int == 0)
stop("Ambiguous link type")
g <- switch(link,
identity = function(p)
p,
log = function(p)
log(p),
logsurv = function(p)
- log(1 - p),
logit = function(p)
{
p <- ifelse(p > 1, NA, ifelse(p < 0, NA, p))
p1 <- ifelse(p > 1 - 9.9999999999999998e-17, 1 -
9.9999999999999998e-17, ifelse(p < 1.0000000000000001e-16,
1.0000000000000001e-16, p))
log(p1/(1 - p1))
}
,
sqrt = function(p)
sqrt(p),
inverse = function(p)
1/p,
probit = qnorm,
loglog = ,
cloglog = function(p)
{
p <- ifelse(p > 1, NA, ifelse(p < 0, NA, p))
p1 <- ifelse(p > 1 - 9.9999999999999998e-8, 1 -
9.9999999999999998e-8, ifelse(p < 1.0000000000000001e-8,
1.0000000000000001e-8, p))
log((log(1 - p1)*-1))
}
)
g
}
, comment = "16/09/2002")
"link.inverse" <-
structure(function(link = "help")
{
link.options <- c("identity", "log", "logsurv", "logit", "sqrt", "inverse",
"probit", "loglog", "cloglog", "help")
if(link == "help")
cat("\nOptions are:", paste(link.options, sep = " "), "\n")
link.int <- charmatch(link, link.options)
if(is.na(link.int))
stop(paste("Invalid link type: \nOptions are:", paste(link.options,
collapse = " ")))
else if(link.int == 0)
stop("Ambiguous link type")
f <- switch(link,
identity = function(x)
x,
log = function(x)
exp(x),
logsurv = function(x)
1 - exp( - x),
logit = function(x)
{
z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, x)))
z/(z + 1)
}
,
sqrt = function(x)
x^2,
inverse = function(x)
1/x,
probit = pnorm,
loglog = ,
cloglog = function(x)
{
z <- exp(ifelse(x > 80, 80, ifelse(x < -80, -80, x)))
1 - exp(.Uminus(z))
}
)
f
}
, comment = "16/10/1994")
"list.mat" <-
structure(function(x)
{
# To convert a list of vectors into a matrix
#
# x must be a named list of "named vectors"
rownames <- names(x)
colnames <- names(x[[1]])
mat <- matrix(nr = length(rownames), nc = length(colnames), dimnames = list(
rownames, colnames))
for(i in rownames)
mat[i, ] <- unlist(x[[i]])
mat
}
, comment = "01/02/1997")
"list.to.data.frame" <-
structure(function(x)
{
# 1) Test list element lengths are the same
lengths <- unlist(lapply(x, length))
len <- mean(lengths)
if(all(lengths == len)) {
attr(x, "class") <- "data.frame"
attr(x, "row.names") <- 1:len
}
else {
stop("list elements are not of same length")
}
x
}
, comment = "19/04/1997")
"logit" <-
structure(function(p)
{
pres <- (p > 0) & (p < 1)
p[pres] <- log(p[pres]/(1 - p[pres]))
n <- length(p[!pres])
p[!pres] <- rep(NA, n)
p
}
, comment = "28/03/2001")
"logit.bt" <-
structure(function(x)
{
p <- exp(x)/(1 + exp(x))
p
}
, comment = "04/10/2001")
"lookup" <-
structure(function(x, y, z = NULL)
{
# Looks up a legend to find the corresponding elements of the input vector
#
# y has all possible values x could take;
# z has all corresponding unabbreviated values of interest
# This function returns those bits of z which correspond to the input values x
#
# if z is NULL, y must be a matrix of two columns
if(is.null(z) && dim(y)[2] != 2) stop("z must be specified OR y a matrix\n"
)
if(!is.null(dim(y))) {
z <- y[, 2]
y <- y[, 1]
}
names(z) <- y
z[x]
}
, comment = "26/01/1999")
"make.df" <-
structure(function(text.file, ...)
{
# Making a dataframe from a text file and using tab separators
read.table(text.file, T, as.is = T, sep = "\t", row.names = NULL, ...)
}
, comment = "01/04/1997")
"make.factors" <-
structure(function(df, fcols = 1, ordered = F, ...)
{
### To change specified columns into factors. By default levels are reordered
### in alphabetical order. Keeping that order makes it an ordered factor (ordered = T)
### which does some strange things in glms and the like, but it's useful for ordering
### factors for trellis plotting.
###
### Without named columns, will not work unless fcols are sequential (col nos. are used:
### somehow that makes a difference)
###
### Vectors can also be used, in which case setting sorted to T is the same
### as using factor()
###
onecol <- is.null(dim(df)) #
if(onecol) {
vec <- unfactor(df)
if(!ordered)
df <- factor(vec, ...)
else df <- ordered(vec, labels = unique(vec, ...))
}
else {
if(is.character(fcols))
fcols <- find.pos(fcols, names(df))
for(i in fcols) {
if(!ordered) {
df[[i]] <- factor(df[[i]], ...)
}
else df[[i]] <- ordered(df[[i]], levels = unique(df[[i]]), ...)
}
}
df
}
, comment = "06/01/1998")
"make.id" <-
structure(function(x)
{
# Creates one id value for each row. Value increases by 1 for each
# element where the value of x is less than for the previous element.
dx <- c(0, diff(x))
u <- rep(0, length(x))
u[dx < 0] <- 1
cumsum(u) + 1
}
, comment = "07/04/1997")
"make.line.loess" <-
structure(function(xy, lab, label.pos = 0.9, cex = par()$cex, lev.font = 1, lty = 1, lwd
= par()$lwd, inclusion = "ALL", mark = T, mark.extra = "")
{
overlap <- function(x, y, xl, yl)
{
over <- F
i <- 1
while(i <= length(xl)) {
first.pair.match <- (1:length(xl[[i]]))[(xl[[i]] == x[1]) & yl[[i]] ==
y[1]]
if(is.na(diff(x)))
browser()
if((diff(x) == 0) & (diff(y) == 0)) {
if(length(first.pair.match) > 0)
over <- T
}
else {
second.pair.match <- (1:length(xl[[i]]))[(xl[[i]] == x[2]) & yl[[i
]] == y[2]]
names(first.pair.match) <- rep("f", length(first.pair.match))
names(second.pair.match) <- rep("s", length(second.pair.match))
pair.match <- sort(c(first.pair.match, second.pair.match))
adj.pairs <- (1:(length(pair.match) - 1))[diff(pair.match) == 1]
adj.pairs <- adj.pairs[!is.na(adj.pairs)]
if(length(adj.pairs) > 0)
if((names(pair.match)[adj.pairs[1]] == "f") & (names(pair.match
)[adj.pairs[1] + 1] == "s"))
over <- T
}
i <- i + 1
}
over
}
assign("overlap", overlap, 0)
pos <- function(posn, x, y)
{
if(posn == 1)
length(x)
else {
dx <- c(0, diff(x))
dy <- c(0, diff(y))
lengths <- sqrt(dx^2 + dy^2)
cum.lengths <- cumsum(lengths)
full.length <- cum.lengths[length(cum.lengths)]
main.index <- ((1:length(cum.lengths))[cum.lengths > posn *
full.length])[1] - 1
left.over <- posn * full.length - cum.lengths[main.index]
main.index + left.over/lengths[main.index + 1]
}
}
assign("pos", pos, 0)
label.line <- function(label.pos, x, y, lab, cex, i, font)
{
if(label.pos >= 0) {
j <- pos(label.pos, x, y)
index <- trunc(j)
if(!is.na(diff(x[index + 0:1])) & !is.na(diff(y[index + 0:1])))
if(!overlap(x[index + 0:1], y[index + 0:1], x.keep, y.keep)) {
j <- j - index
tx <- x[index] + j * (x[index + 1] - x[index])
ty <- y[index] + j * (y[index + 1] - y[index])
text(tx, ty, lab, cex = cex, font = font)
x.keep[[i - 1]] <- x
y.keep[[i - 1]] <- y
assign("x.keep", x.keep, 1)
assign("y.keep", y.keep, 1)
}
}
}
plot.line <- function(x, y, lty, lwd, inclusion)
{
find.point <- function(inclusion, x, y, start)
{
p <- pos(inclusion, x, y)
index <- trunc(p)
p <- p - index
tx <- x[index] + p * (x[index + 1] - x[index])
ty <- y[index] + p * (y[index + 1] - y[index])
list(tx = tx, ty = ty, index = index + start)
}
if(is.character(inclusion)) {
lines(x, y, lty = lty, lwd = lwd)
T
}
else if(is.null(inclusion))
F
else if(inclusion[1] == inclusion[2])
F
else {
if(inclusion[1] > inclusion[2])
inclusion <- c(inclusion[1], 1, 0, inclusion[2])
x.save <- x
y.save <- y
for(i in 1:(length(inclusion) %/% 2)) {
x <- x.save
y <- y.save
start <- find.point(inclusion[1 + (i - 1) * 2], x, y, T)
finish <- find.point(inclusion[2 + (i - 1) * 2], x, y, F)
x <- c(start$tx, x[start$index:finish$index], finish$tx)
y <- c(start$ty, y[start$index:finish$index], finish$ty)
lines(x, y, lty = lty, lwd = lwd)
}
T
}
}
## Main ##
xy$x <- c(NA, xy$x, NA)
na <- (1:length(xy$x))[is.na(xy$x)]
xy$x <- xy$x[ - (na[diff(na) == 1])]
xy$y <- c(NA, xy$y, NA)
xy$y <- xy$y[ - (na[diff(na) == 1])]
na <- (1:length(xy$x))[is.na(xy$x)]
x.keep <- list()
y.keep <- list()
assign("x.keep", x.keep, 1)
assign("y.keep", y.keep, 1)
if(length(label.pos) < length(na) - 1)
label.pos <- c(label.pos, rep(label.pos[1], length(na) - length(
label.pos)))
cat("Lev ", lab, mark.extra, " Lines ", max(length(na) - 1, 0), "\n", sep
= "")
if(length(na) > 0)
for(i in 2:length(na)) {
x <- xy$x[(na[i - 1] + 1):(na[i] - 1)]
y <- xy$y[(na[i - 1] + 1):(na[i] - 1)]
old.x <- x
old.y <- y
x <- x[!is.na(old.x) & !is.na(old.y)]
y <- y[!is.na(old.x) & !is.na(old.y)]
if(length(x) > 0) {
if(plot.line(x, y, lty, lwd, inclusion[[i - 1]]))
label.line(label.pos[i - 1], x, y, paste(lab, ifelse(mark,
paste(mark.extra, "-", i - 1, sep = ""), ""), sep = ""), cex,
i, font = lev.font)
}
}
}
, comment = "03/06/1998")
"mean.lt" <-
structure(function(start.list = ab.soluble, choose = 1, lt = 99, intv = 95, new.order = F,
leg.beg = NULL, leg.end = 4, insect = "interesting critters", fit = NULL,
two.tables = F, borders = F, omit = NULL)
{
## Requires LT data in an ab. type list with sublist names such as "total"
## referring to LTs at those locations
##
## LT50s calculated by mono fit: LT99s calculated by line fit
## Has to group species to take mean, and then restore to the original
## order of species (otherwise, it's very simple)
## The name "species" really refers to location and life stage in
## some cases: it is what appears on the individual mortality plots
##
## leg.end will usually be used to remove "Rep_" from hames. If there are more than
## 9 reps in the legend name, it won't work.
##
## leg.beg is the number of characters to omit from the legend text in determining what
## groups similar trials
## leg.beg by default is the postion of the first blank in the legend name: That is
## to allow the practice of having trial numbers in that legend: if there's no space
## in the legend, it won't work.
ablist <- as.character(substitute(start.list))
leng.name <- nchar(ablist)
interval <- intv/100
if(is.character(choose)) {
choose.i <- (1:length(names(start.list)))[names(start.list) ==
as.character(choose)]
}
else (choose.i <- choose)
ab <- start.list[[choose.i]]
legs <- as.character(ab$legend)
### Seems to be necessary for correct class
if(is.null(leg.beg)) {
### Find the first space in the legend names:
leg.beg <- NULL
for(k in 1:length(legs)) {
leg.vec <- paste(substring(legs[k], 1:nchar(legs[k]), 1:nchar(legs[k]
)), sep = "")
leg.beg[k] <- match(" ", leg.vec)
}
}
spec <- substring(legs, leg.beg + 1, nchar(legs) - leg.end)
include <- rep(TRUE, length = (length(spec))) #
### Are there any we wish to exclude? (22/1/2001)
if(!is.null(omit))
include[grep(omit, spec)] <- FALSE
spec <- spec[include]
uniq.spec <- unique(spec)
sortspec <- sort(spec)
if(is.null(fit))
fit <- ifelse(lt == 50, "monotone", "line")
if(fit == "line")
lt.vector <- ab$lt[, paste(lt)]
else lt.vector <- ab$lt.monotone[, paste(lt)]
type <- paste(insect, " [", choose, "] (", fit, " fit)", sep = "")
lt.vector[lt.vector < 0 | lt.vector == Inf] <- NA # (22/1/2001)
lt.vector <- lt.vector[include] # (22/1/2001)
lts <- round(as.numeric(lt.vector), 5)
lt.1 <- matrix(round(lts, 1), ncol = 1, dimnames = list(spec, NULL))
lt.2 <- NULL #
### Sort to arrange species together
for(i in uniq.spec) {
lt.2 <- c(lt.2, lt.1[dimnames(lt.1)[[1]] == i, ])
}
### Change vector into single column matrix
lt.3 <- matrix(lt.2, nc = 1, dimnames = list(names(lt.2), NULL)) #
### Have vector of grouped lts. Continue finding means, etc.
av.log.lt <- tapply(log(lts), spec, mean, na.rm = T)
my.var <- function(x) var(x, na.rm = T)
var.log.lt.i <- tapply(log(lts), spec, my.var) #
### Replace any missing variances with zero
var.log.lt <- ifelse(is.na(var.log.lt.i), 0, var.log.lt.i)
sem.log.lt <- tapply(log(lts), spec, sem, na.rm = T)
reps <- tapply(lts, spec, function(x)
length(x[!is.na(x)]))
deg.pool <- sum(reps[reps > 0] - 1)
ind.var.free <- var.log.lt * (reps - 1)
pool.var.free <- sum(ind.var.free)/deg.pool
delta <- qt(1 - (1 - interval)/2, deg.pool) * sqrt(pool.var.free/reps)
log.up <- av.log.lt + delta
log.low <- av.log.lt - delta
lt.mean <- round(exp(av.log.lt), 1)
upper <- round(exp(log.up), 1)
lower <- round(exp(log.low), 1)
sem <- round(sem.log.lt * lt.mean, 3) #
### Make matrix of means, etc
y <- cbind(reps, lt.mean, lower, upper, sem)
x <- y[uniq.spec, ] #
### Put individual LTs into matrix
rn <- dimnames(lt.3)[[1]]
urn <- unique(rn)
cn <- max(table(rn))
mm <- matrix(rep("---", length(urn) * cn), nc = cn, dimnames = list(urn,
paste("Rep", 1:cn, sep = "")))
for(i in urn) {
### Not all rows have same number of reps
kk <- lt.2[names(lt.2) == i]
for(k in 1:length(kk))
mm[i, k] <- kk[k]
}
### If original order of plots is not wanted, sort into new order
if(new.order) {
mm <- mm[order(dimnames(mm)[[1]]), ]
x <- x[order(dimnames(x)[[1]]), ]
}
### Output to screen
if(two.tables) {
cat(c(paste("Separate LT", lt, "s for ", type, ":", sep = ""), "\n\n"))
tab.mat(mm)
cat("\n\n")
cat(c(paste("Mean LT", lt, "s and ", intv, "% c.i.:", sep = ""), "\n\n")
)
tab.mat(x)
cat("\n\n")
}
else {
one.table <- cbind(mm, x)
cat(c(paste("Separate LT", lt, "s with mean and ", intv, "% c.i. for ",
type, ":", sep = ""), "\n\n"))
tab.mat(one.table)
}
if(borders)
print.char.matrix(unfactor(as.data.frame(one.table)), col.names = T)
}
, comment = "14/12/2004")
"mean.lt.short" <-
structure(function(lt.vector, lt = 99, intv = 95)
{
# Shortened version of mean.lt()
interval <- intv/100
lts <- round(as.numeric(lt.vector), 5)
lts <- lts[!is.na(lts)]
reps <- length(lts) # Change vector into single column matrix
lt.3 <- matrix(lts, nc = 1) #
# Have vector of grouped lts. Continue finding means, etc.
av.log.lt <- mean(log(lts), na.rm = T)
var.log.lt <- my.var(log(lts)) #
# Replace any missing variances with zero
sem.log.lt <- sem(log(lts), na.rm = T)
deg <- length(lts) - 1
delta <- qt(1 - (1 - interval)/2, deg) * sqrt(var.log.lt/reps)
log.up <- av.log.lt + delta
log.low <- av.log.lt - delta
lt.mean <- round(exp(av.log.lt), 1)
upper <- round(exp(log.up), 1)
lower <- round(exp(log.low), 1)
sem <- round(sem.log.lt * lt.mean, 3) #
# Make vector of means, etc
y <- matrix(c(reps, lt.mean, lower, upper, sem), nr = 1)
dimnames(y) <- list(NULL, c("reps", "lt.mean", "lower", "upper", " sem"))
# Put individual LTs into vector
mm <- matrix(round(lts, 1), nr = 1)
dimnames(mm) <- list(NULL, paste("Rep", 1:length(mm), sep = ""))
# Output to screen
cat(c(paste("Separate LT", lt, "s ", sep = ""), "\n\n"))
tab.mat(mm)
cat("\n\n")
cat(c(paste("Mean LT", lt, "s and ", intv, "% c.i.:", sep = ""), "\n\n"))
#if(dim(y)[1]==1)
#x_matrix(x,nc=5,dimnames=dimnames(y))
tab.mat(y)
cat("\n\n")
}
, comment = "20/08/1996")
"occurrence" <-
structure(function(x, char = " ", positions = 1)
{
### Finds positions of specified occurrences of a character (char) in a
### character string (x), or vector thereof
#
### Default is to find the first occurrence. A vector of occurrences can be given
if(!is.character(x)) x <- as.character(x)
z <- as.list(x)
find.gaps <- function(a, b)
{
size <- 1:nchar(a)
size[substring(a, 1:nchar(a), 1:nchar(a)) == b]
}
y <- lapply(z, function(w)
find.gaps(w, char))
longest <- max(sapply(y, length))
names(y) <- paste(1:length(y)) #
# if the number of occurrences is not the same for all elenents...
if(any(diff((sapply(y, length)))) > 0) {
# Fix up odd ones
short.y <- names(y[sapply(y, length) < longest])
for(i in short.y) {
short.y.len <- length(unlist(y[[i]]))
y[[i]] <- c(unlist(y[[i]]), rep(NA, (longest - short.y.len)))
}
}
if(max(positions) > longest) {
positions <- positions[!(positions > longest)]
cat(paste("No more than", longest, "of such character in any element\n")
)
}
t(matrix(unlist(y), byrow = F, nc = length(y)))[, positions]
}
, comment = "26/01/1999")
"p.sunflowers" <-
structure(function(x, y, number, size = 0.125, add = FALSE, rotate = F, pch = 16, ...)
{
## Purpose: Produce a 'sunflower'-Plot
## -------------------------------------------------------------------------
## Arguments: x,y: coordinates;
## number[i] = number of times for (x[i],y[i]) [may be 0]
## size: in inches add : Should add to a previous plot ?
## rotate: randomly rotate flowers? further args: as for plot(..)
## -------------------------------------------------------------------------
## Authors: Andreas Ruckstuhl, Werner Stahel, Martin Maechler, Tim Hesterberg
## Date : Aug 89 / Jan 93, March 92, Jan, Dec 93, Jan 93
## Examples: p.sunflowers(x=sort(2*round(rnorm(100))), y=round(rnorm(100),0))
## 16:12, 16 Nov 2006 (NZDT)[[User:Sven|Sven]] p.sunflowers(rnorm(100),rnorm(100),number=rpois(n=100,lambda=2),
## rotate=T, main="Sunflower plot")
n <- length(x)
if(length(y) != n)
stop("x & y must have same length !")
if(missing(number)) {
orderxy <- order(x, y)
x <- x[orderxy]
y <- y[orderxy]
first <- c(T, (x[-1] != x[ - n]) | (y[-1] != y[ - n]))
x <- x[first]
y <- y[first]
number <- diff(c((1:n)[first], n + 1))
}
else {
if(length(number) != n)
stop("number must have same length as x & y\n!")
x <- x[number > 0]
y <- y[number > 0]
number <- number[number > 0]
}
n <- length(x)
if(!add) {
axislabels <- match(c("xlab", "ylab"), names(list(...)))
if(!is.na(axislabels[1]))
xlab <- list(...)[[axislabels[1]]]
else xlab <- deparse(substitute(x))
if(!is.na(axislabels[2]))
ylab <- list(...)[[axislabels[2]]]
else ylab <- deparse(substitute(y))
plot(x, y, xlab = xlab, ylab = ylab, type = "n", ...)
}
nequ1 <- number == 1
if(any(nequ1))
points(x[nequ1], y[nequ1], pch = pch, csi = size * 1.25)
if(any(!nequ1))
points(x[!nequ1], y[!nequ1], pch = pch, csi = size * 0.80000000000000004
)
i.multi <- (1:n)[number > 1]
if(length(i.multi)) {
ppin <- par()$pin
pusr <- par()$usr
xr <- (size * abs(pusr[2] - pusr[1]))/ppin[1]
yr <- (size * abs(pusr[4] - pusr[3]))/ppin[2]
i.rep <- rep(i.multi, number[number > 1])
z <- NULL
for(i in i.multi)
z <- c(z, 1:number[i] + if(rotate) runif(1) else 0)
deg <- (2 * pi * z)/number[i.rep]
segments(x[i.rep], y[i.rep], x[i.rep] + xr * sin(deg), y[i.rep] + yr *
cos(deg))
}
}
, comment = "14/02/1997")
"pdfn" <-
structure(function(x, ps2pdf = TRUE, cex = .5, adj = 1, side = 1, line = -1, ...)
{
### Purpose: Print file name on a plot
### ----------------------------------------------------------------------
### Modified from:
### ----------------------------------------------------------------------
### Arguments: ps2pdf -- are we making a postscript name into a PDF name?
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 24 Feb 2004, 14:55
par(new = T)
par(omi = rep(0,4), mai = rep(0,4), ...)
if(ps2pdf)
x <- ppaste(substring(x, 1,nchar(x)-1), "df")
blank.plot()
mtext(x, side = side, line = line, cex = cex, adj = adj)
}
, comment = "28/07/2004")
"pid" <-
structure(function()
### current PID
system("echo $PPID",intern=T)
, comment = "24/05/2001")
"plot.lt" <-
structure(function(lt = 99, datafun = get.leaf, choose = "lbam", fun = function(x)
x, poly = T, poly.limits = list(NULL), poly.keep.x = list(NULL),
poly.plot.limits = list(NULL), deg = 2, splin = F, wide.legend = 1,
high.legend = 1, expand.y = 0, print.statistics = T, x.limits = NULL,
y.limits = NULL, predict.x = NULL, confidence = 0.95, compare = NULL,
fit.robust = F, plot.ci = F, plotit = T, plot.means = T, plot.points = F,
plot.single.sem = F, plot.groups.sem = T, legend.sem = T, legend.linespace
= 1.25, errors = "aov", include = NULL, parallel = NULL, mkh = 0.08,
datestamp = T, prefix = "", prefix.y = "", legend.header = "", legend.note
= "Limits shown are one SE\neither side of the mean", bars.pretext =
"2 SE", na.means = T, wts = NULL, style = list(lwd = 1, plotch = c(0, 1, 2,
5, 3, 4, 16, 18, 6:15), cex.labels = 1, cex.title = 1.15, mai.xtra = c(0.25,
0.5, 0.25, 0.5)), mtext.line = -1)
{
### Means, for each temperature in each list item, are calculated and plotted
### An overall SD is calculated for each set of list items,
### as identified in groups.
### Set na.means to F to omit plotting of means when data for a mean relies
### partly on NAs
### Examples
### plot.lt(lt=50,datafun=get.air.CA.lt,poly=F, poly.limits=list(c(32.5,40)),
### splin=T,wide.legend=0.67)
### plot.lt(lt=50,datafun=get.co2.o2.lt,poly=c(T,F),deg=1, splin=F,wide.legend=0.67)
### plot.lt(lt=99,datafun=get.co2.o2.lt,poly=c(T,F),deg=2, splin=F,wide.legend=0.67)
### plot.lt(lt=50,datafun=get.pretreat.lt,poly=T,deg=2, splin=F,wide.legend=0.67)
#
### aov.errors = T causes the within x mean square to be used for SEs etc.
#
### lt= : gives lt % to be printed on the y-axis label
#
### datafun: This function collects up the data, labelling information etc.,
### and delivers it in a standard form to plot.lt.
#
### fun: Transformation function for the y-axis. Usually identity or log.
#
### poly: T if a line or polynomial is required.
#
### poly.limits: used in connection with splin=T when a specific polynomial
### (usually a quadratic) is required over the part of the range specified
### by poly.limits, with a spline for the remainder.
###
### deg: degree for polynomial(s)
#
### Note: Parameters poly, poly.limits and deg are expanded, if necessary,
### into a vector or (for poly.limits) a list.
### wide.legend, high.legend: point at which to start printing first legend item
### Specify as a fraction of the x- or y-distance.
### expand.y: e.g. 25% will increase the upper y-limit to give 25% extra space
### for legend etc.
#
### print.statistics: whether to print out coeffs, ses etc of fitted polynomials
### plot.means, plot.points: means and/or points
### plot.single.sem, plot.groups.sem: strategy for calculating sem
### errors: strategy for calculation of error mean square.
### "aov" calculates the mean square from the "within x" sum of squares.
### Alternatives are lack.of.fit (deviations from the fitted polynomial)
### and total (aov and lack.of.fit combined).
### include: This is a vector of booleans; T where a list element is included.
#
#
par(new = F)
options(width = 180)
var.new <- function(x)
{
pres <- !is.na(x)
var(x[pres])
}
date.txt <- unix("date '+DATE: %d/%m/%y %H:%Mh'")
docs <- paste(unix("pwd"), as.character(as.name(match.call())), date.txt,
collapse = "")
print(docs)
if(!plotit) {
if(missing(plot.means))
plot.means <- F
if(missing(plot.points))
plot.points <- F
}
plotch <- style$plotch
lwd <- style$lwd
cex.title <- style$cex.title
cex.labels <- style$cex.labels
mai.xtra <- style$mai.xtra
if(!is.null(confidence) & !is.null(compare)) {
cat("\nValues set both for ci's & for comparisons.",
"\nCode for both is not implemented at this point.",
"\nWill do comparisons only.\n")
confidence <- NULL
}
if(plot.points | plot.means | plotit)
plot.some <- T
else plot.some <- F
xtra <- numeric(0)
linetype.len <- 0
oldpar <- par()
on.exit(par(oldpar))
par(mai = par()$mai + mai.xtra)
u <- datafun(lt = lt, choice = choose)
leg <- u$leg
ltraw.list <- u$lt
xval.list <- u$xval
xlab <- u$xlab
xaxlab <- u$xaxlab
yaxlab <- u$yaxlab
ylab <- u$ylab
groups <- u$groups
maint <- u$maint
if(is.list(maint))
maint <- paste(maint[["prefix"]], " LT~v~c.8~.", lt, "~h0~c1~.", , " ",
maint[["main"]], sep = "")
else maint <- paste(prefix, maint)
mean.list <- lapply(xval.list, unique)
### Set up structure that has the right shape
num.list <- mean.list
uniqx.list <- mean.list
nlist <- length(xval.list)
if(is.numeric(include)) {
tmp <- rep(F, nlist)
tmp[include] <- T
include <- tmp
}
else if(is.null(include))
include <- rep(T, nlist)
else if(length(include) < nlist)
include <- c(include, rep(F, nlist - length(include)))
include.parallel <- rep(F, nlist)
if(!is.null(parallel)) {
include.parallel[parallel] <- T
xval.parallel <- list()
pres.parallel <- list()
lt.parallel <- list()
n.parallel <- sum(include.parallel)
num.parallel <- rep(0, n.parallel)
wts.parallel <- rep(0, n.parallel)
df.parallel <- rep(0, n.parallel)
leg.parallel <- leg[include.parallel]
}
if(is.null(groups))
groups <- as.list((1:nlist)[include])
ngroups <- length(groups)
if(nlist > 1 & length(deg) == 1)
deg.pol <- rep(deg, nlist)
else deg.pol <- deg
if(nlist > 1 & length(poly) == 1)
poly <- rep(poly, nlist)
if(nlist > 1 & length(splin) == 1)
splin <- rep(splin, nlist)
if(!is.list(poly.limits))
poly.limits <- list(poly.limits)
if(!is.list(poly.keep.x))
poly.keep.x <- list(poly.keep.x)
if(nlist > 1 & length(poly.limits) == 1)
poly.limits <- rep(poly.limits, nlist)
if(nlist > 1 & length(poly.keep.x) == 1)
poly.keep <- rep(poly.keep.x, nlist)
if(!is.list(poly.plot.limits))
poly.plot.limits <- list(poly.plot.limits)
if(nlist > 1 & length(poly.plot.limits) == 1)
poly.plot.limits <- rep(poly.plot.limits, nlist)
if(!is.null(wts))
if(nlist > 1 & length(wts) == 1)
wts <- rep(wts, nlist)
xlab <- u$xlab
if(sum(include) > length(plotch))
plotch <- c(plotch, letters)
vm.aov <- array(, nlist)
vm.other <- array(, nlist)
df.other <- array(, nlist)
numdf <- array(, nlist)
nummin <- array(, nlist)
nummax <- array(, nlist)
lt.list <- ltraw.list
for(i in 1:nlist)
lt.list[[i]] <- fun(ltraw.list[[i]])
mat.range <- matrix(unlist(lapply(ltraw.list, range, na.rm = T)), nrow = 2)
unlist.raw <- unlist(ltraw.list)
yrange <- range(mat.range[, include], na.rm = T)
log.y <- F
if(all(fun(exp(1:4)) == 1:4)) {
log.y <- T
mat.range <- matrix(unlist(lapply(ltraw.list, function(x, na.rm = T)
range(x[x > 0]), na.rm = T)), nrow = 2)
yrange[1] <- min(mat.range, na.rm = T)
yrange[1] <- min(unlist.raw[unlist.raw > 0], na.rm = T)
}
if(is.null(y.limits)) {
funrange <- fun(yrange)
y.limits <- yrange
}
else funrange <- fun(y.limits)
funrange[2] <- funrange[2] + diff(funrange) * expand.y
### Often set expand.y = 0.05 or 0.1, to allow room for legend
if(is.null(yaxlab)) {
ylabels <- pretty(y.limits)
if(all(fun(exp(1:4)) == 1:4)) {
splity <- sqrt(y.limits[1] * y.limits[2])
ylabels2 <- pretty(c(splity, max(y.limits)), 4)
ylabels1 <- pretty(c(min(y.limits), splity), 4)
ylabels <- sort(unique(c(ylabels1, ylabels2)))
if(ylabels[1] == 0)
ylabels[1] <- mean(ylabels[1:2])
}
}
else ylabels <- yaxlab
ncharwid <- max(nchar(paste(ylabels)))
ylabadd <- (ncharwid - 1.5) * 0.25
ytitladd <- (1.25 * ncharwid - 2) * 0.125
if(plot.some) {
if(is.null(x.limits))
x.limits <- range(unlist(xval.list), na.rm = T)
plot(x.limits, funrange, xlim = x.limits, type = "n", xlab = "", xaxt =
"n", ylab = "", yaxt = "n")
if(is.null(xaxlab))
xaxlab <- pretty(x.limits)
axis(1, at = xaxlab)
axis(2, at = fun(ylabels), labels = paste(ylabels), mgp = c(3 + ytitladd,
1 + ylabadd, 0))
chh <- par()$cxy[2]
chw <- par()$cxy[1]
uxy <- par()$usr
}
outvals <- vector("list", sum(include))
names(outvals) <- names(lt.list)[include]
i.plot <- 0
j.parallel <- 0
for(i in 1:nlist)
if(include[i]) {
u.x <- NULL
bt.y <- NULL
i.plot <- i.plot + 1
deg.i <- deg.pol[i]
pres <- !is.na(lt.list[[i]])
if(all(!pres))
next
temps <- xval.list[[i]][pres]
ltval <- lt.list[[i]][pres]
if(plot.points)
points(temps, ltval, pch = plotch[i.plot], mkh = 0.75 * mkh, lwd
= lwd)
vv <- tapply(ltval, list(temps), var)
my <- tapply(ltval, list(temps), mean)
tab.temps <- table(temps)
utemps <- as.numeric(names(tab.temps))
num <- as.numeric(tab.temps)
num.list[[i]] <- num
uniqx.list[[i]] <- utemps
if(!na.means) {
my.test.na <- tapply(lt.list[[i]], list(xval.list[[i]]), function(
x)
any(is.na(x)))
utemps.all <- unique(xval.list[[i]])
utemps.na <- utemps.all[my.test.na]
if(length(utemps.na) > 0) {
u.match <- match(utemps.na, utemps, nomatch = 0)
if(any(u.match > 0))
my[u.match] <- NA
}
mean.list[[i]] <- my
}
if(plot.means)
points(utemps, my, pch = plotch[i.plot], mkh = mkh, lwd = 2)
nummin[i] <- min(num[!is.na(my)])
nummax[i] <- max(num[!is.na(my)])
sd <- sqrt(vv)
sem <- sd/sqrt(num)
df.aov.rss <- sum(num - 1)
numdf[i] <- df.aov.rss
ms.aov <- sum(vv * (num - 1), na.rm = T)/df.aov.rss
vm.aov[i] <- ms.aov
sd.aov <- sqrt(ms.aov)
cat("\n", leg[i], " ", date(), fill = T)
av.se <- data.frame(rbind(format(round(my, 2)), format(round(sem, 2)),
num))
row.names(av.se) <- c("Means are:", "sem's are:", "Nos.")
names(av.se) <- paste(utemps)
print(av.se)
cat("sds are", format(round(sd, 2)))
if(df.aov.rss > 0)
cat(" pooled:", format(round(sd.aov, 3)), " df =", df.aov.rss,
fill = T)
if(log.y) {
btmeans <- exp(my)
btsem <- btmeans * sem
bt.se <- data.frame(rbind(format(round(btmeans, 1)), format(round(
btsem, 2))))
row.names(bt.se) <- c("BT means:", "Appr sem:")
names(bt.se) <- paste(utemps)
print(bt.se)
}
cat("\n", fill = T)
rtemps <- range(temps)
utemps <- uniqx.list[[i]]
poly.i <- poly[i]
if(poly.i)
linetype.len <- 3.5
splin.i <- splin[i]
if(poly.i & length(utemps) <= deg.i) {
if(length(utemps) > 1)
cat("*** Insufficient points for polynomial of degree", deg.i,
fill = T)
poly.i <- F
}
if(splin.i & length(utemps) < 4) {
if(length(utemps) > 1)
cat(
"*** Insufficient points for spline -- at least 4 are required",
fill = T)
splin.i <- F
}
poly.lim <- poly.limits[[i]]
poly.keep <- poly.keep.x[[i]]
if(is.null(poly.lim))
poly.lim <- rtemps
else if(any(is.na(match(poly.lim, utemps))))
cat(
"One or both endpoints for the polynomial fit do not match actual data values",
fill = T)
if(length(poly.lim) == 1 & is.null(poly.keep)) {
cat("Only one limit has been given; cannot fit quadratic", fill =
T)
poly.lim <- rep(poly.lim, 2)
}
if(is.null(poly.keep)) {
keep <- (utemps >= poly.lim[1] & utemps <= poly.lim[2])
poly.keep <- utemps[keep]
}
else keep <- as.logical(match(utemps, poly.keep, nomatch = F))
nu <- (1:length(utemps))[keep]
y.means <- my[nu]
x.means <- utemps[nu]
yfit <- my
use.quad <- as.logical(match(temps, poly.keep, nomatch = F))
x.quad <- temps[use.quad]
y.quad <- ltval[use.quad]
num.quad <- num[nu]
df.quad <- sum(num.quad - 1)
if(df.quad == 0)
ms.quad <- 0
else ms.quad <- sum(vv[nu] * (num.quad - 1), na.rm = T)/df.quad
ems <- ms.quad
df.t <- df.quad
if(poly.i) {
poly.legend <- "Polynomial"
if(fit.robust)
poly.legend <- "Polynomial (robust fit)"
if(!is.null(poly.keep.x[[i]]))
cat(poly.legend, "of degree", deg.i, "through points",
poly.keep, fill = T)
else cat(poly.legend, "of degree", deg.i, "from points", poly.lim[
1], "to", poly.lim[2], fill = T)
uu <- poly.fit(x.quad, y.quad, num.quad, deg.i, ms.quad, df.quad,
errors, print.statistics = print.statistics, predict.x =
predict.x, fit.robust = fit.robust, wts = wts)
ypred.se <- uu$pred$se.fit
df.t <- uu$df.ms
xtra <- uu$xtra
ems <- uu$ems
vm.other[i] <- ems
df.other[i] <- df.t
df.t0 <- uu$pred$df
ems0 <- uu$pred$ms
if(!is.null(df.t0))
if(df.t != df.t0)
cat("\n ******* df.t =", df.t, " df.t0 =", df.t0, "\n")
if(!is.null(ems0)) if(ems != ems0) {
cat("****** Warning ******\n")
cat("ems =", ems, " ems0 =", ems0, "\n\n")
}
### vanilla ms; depends on errors
b <- uu$b
if(include.parallel[i]) {
j.parallel <- j.parallel + 1
xval.parallel[[j.parallel]] <- x.quad
lt.parallel[[j.parallel]] <- y.quad
num.parallel[j.parallel] <- sum(num[nu])
wts.parallel[j.parallel] <- 1/ems
df.parallel[j.parallel] <- df.t
}
if(!is.null(predict.x)) {
ypred <- uu$pred$fit
u.x <- predict.x
if(!is.null(confidence)) {
alpha <- 1 - confidence
c.txt <- paste(round(confidence * 100))
c.txt1 <- NULL
t.stat <- qt(1 - alpha/2, df.t)
}
else if(!is.null(compare)) {
alpha <- 1 - compare
t.stat <- qt(1 - alpha/2, df.t)/sqrt(2)
c.txt <- format(round((2 * pt(t.stat, df.t) - 1) * 100, 1))
c.txt1 <- paste(round((1 - compare) * 100))
}
ci.low <- ypred - ypred.se * t.stat
ci.high <- ypred + ypred.se * t.stat
if(log.y) {
bt.y <- exp(ypred)
bt.sem <- bt.y * ypred.se
xx <- data.frame(rbind(format(round(bt.y, 1)), format(signif(
exp(ci.low), 3)), format(signif(exp(ci.high), 3)), format(
signif(bt.sem, 2))))
dimnames(xx) <- list(c("BT preds:", paste(c.txt, "pc low:",
sep = ""), paste(c.txt, "pc high:", sep = ""),
"Approx sem:"), paste(u.x))
print(xx)
if(!is.null(c.txt1))
cat("\nNon-overlap of ci's is equivalent to a test for",
"\nno difference at the", paste(c.txt1, "pc", sep = ""
), "level.\n")
}
else {
cat("\nTemps :", format(round(u.x, 1)), fill = T)
cat("Pred vals:", format(round(ypred, 1)), fill = T)
cat("Appr. sem:", format(round(ypred.se, 2)), fill = T)
bt.y <- ypred
}
}
sem[nu] <- sqrt(ems/num.quad)
if(length(nu) < length(utemps) & errors == "aov") {
ms.other <- sum(vv[ - nu] * (num[ - nu] - 1), na.rm = T)/sum(
num[ - nu] - 1)
sem[ - nu] <- sqrt(ms.other/num[ - nu])
}
}
if(plot.single.sem)
errorbars(utemps, my, sem, 0.25 * chw)
if(poly.i) {
if(is.null(poly.plot.limits[[i]]))
xpoints <- pretty(utemps[nu], 50)
else xpoints <- pretty(poly.plot.limits[[i]], 50)
yfit[nu] <- rep(b[1], length(nu))
yhat <- rep(b[1], length(xpoints))
}
if(poly.i)
if(deg.i > 0)
for(k in 2:(deg.i + 1)) {
yfit[nu] <- yfit[nu] + b[k] * utemps[nu]^(k - 1)
yhat <- yhat + b[k] * xpoints^(k - 1)
}
if(plot.some) {
if(poly.i & !splin.i)
lines(xpoints, yhat, lty = i.plot, lwd = lwd)
if(plot.ci) {
lines(spline(u.x, ci.low), lty = i.plot)
lines(spline(u.x, ci.high), lty = i.plot)
}
}
if(splin.i) {
### spline; perhaps polynomial through points between poly.limits
if(diff(rtemps) > diff(range(utemps[nu]))) cat(
"Extend polynomial to the whole range of x-values as a spline curve.",
fill = T)
if(length(utemps) >= 3) {
uus <- spline(utemps, yfit)
xpoints <- uus$x
yhat <- uus$y
if(length(xtra) > 3) {
xtrafit <- approx(uus$x, uus$y, xout = xtra)$y
if(log.y)
xtrafit <- exp(xtrafit)
names(xtrafit) <- paste(xtra)
cat("\n Predictions from spline fit:\n")
print(round(xtrafit, 2))
}
if(plot.some)
lines(xpoints, yhat, lty = i.plot)
}
}
if(!is.null(u.x) & !is.null(bt.y))
outvals[[i.plot]] <- cbind(u.x, bt.y)
}
### Parallel curve analysis, if requested
### At present, lines are the only option (no plots)
if(!is.null(parallel)) {
cat("\nParallel line or curve analysis", fill = T)
cat("Factors are:", leg.parallel, fill = T)
deg.here <- max(deg.pol[include.parallel])
x.all <- unlist(xval.parallel)
y.all <- unlist(lt.parallel)
wts.all <- rep(wts.parallel, num.parallel)
id.parallel <- rep(1:n.parallel, num.parallel)
f.parallel <- factor(id.parallel, labels = leg.parallel)
u.parallel <- lm(lt ~ C(stage, treatment) + temp, data = list(lt = y.all,
temp = x.all, stage = f.parallel, wts = wts.all, deg.here = deg.here),
weights = wts)
u.one <- lm(lt ~ temp, data = list(lt = y.all, temp = x.all, stage =
f.parallel, wts = wts.all, deg.here = deg.here), weights = wts)
print(summary(u.parallel, correl = F))
print(anova(u.one, u.parallel), test = "F")
}
semin <- array(, ngroups)
semax <- array(, ngroups)
groupsd <- array(, nlist)
for(k in 1:ngroups) {
sublist <- groups[[k]] ### N.B. May have >1 list item per group
leg.k <- paste(leg[sublist], collapse = "; ")
if(!any(include[sublist])) {
cat("** List nos.", sublist, "have been specified",
"\nOne or more of these is not in the include list", fill = T)
break
}
nmin <- min(nummin[sublist])
nmax <- max(nummax[sublist])
sd.aov <- sqrt(sum(vm.aov[sublist] * numdf[sublist], na.rm = T)/sum(
numdf[sublist]))
sd.other <- sqrt(sum(vm.other[sublist] * df.other[sublist], na.rm = T)/
sum(df.other[sublist]))
if(is.na(sd.other))
next
if(errors == "aov" & sum(numdf[sublist]) == 0 & poly.i)
errors.here <- "lack.of.fit"
else errors.here <- errors
sd.points <- switch(errors.here,
aov = sd.aov,
lack.of.fit = ,
total = sd.other)
groupsd[k] <- sd.points
if(length(sublist) > 1) {
cat("\nGroup together lists", sublist, fill = T)
cat(" Overall pooled sd =", format(round(sd.points, 3)), fill = T)
}
semax[k] <- sd.points * sqrt(1/nmin)
semin[k] <- sd.points * sqrt(1/nmax)
cat("\n", leg.k, ":")
if(nmax > nmin)
cat("sem = Max:", format(round(semax[k], 3)), " Min:", format(
round(semin[k], 3)), fill = T)
else cat("sem =", format(round(semax[k], 3)), fill = T)
}
if(plot.some) {
if(plot.groups.sem)
if(ngroups > 1) {
se.tmp <- semax
se.tmp[is.na(se.tmp)] <- 0.55 * chh
se.last <- se.tmp[ngroups]
se.two <- se.tmp + c(0, se.tmp[ - ngroups]) + 0.4 * chh
if(length(legend.linespace) == 1) legend.linespace <- se.two/chh
### Can get spacing from parameter legend.linespace by making this a vector
}
x.legend <- uxy[1] + diff(uxy[1:2]) * wide.legend + chw
if(wide.legend == 1)
x.legend <- uxy[2] - linetype.len * chw - max(nchar(leg)) * chw
y.legend <- uxy[3] + diff(uxy[3:4]) * high.legend - chh
if(high.legend == 1) {
y.legend <- uxy[4] - 0.5 * chh
}
if(high.legend == 0)
y.legend <- uxy[1] + se.last + sum(se.two)
xpos <- rep(x.legend, length(leg))
ypos <- y.legend - cumsum(legend.linespace) * chh
}
j0 <- 0
if(plot.groups.sem)
for(k in 1:ngroups) {
sd.here <- groupsd[k]
sublist <- groups[[k]]
if(legend.sem) {
j0 <- j0 + 1
sem.max <- semax[k]
sem.min <- semin[k]
errorbars(xpos[j0] - 1.25 * chw, ypos[j0], sem.max, eps = chw/2)
pretext.pos <- xpos[j0] - 2.5 * chw
if(!is.na(sem.min) & sem.min < 0.75 * sem.max) {
errorbars(xpos[j0] - 3 * chw, ypos[j0], sem.min, eps = chw/2)
pretext.pos <- xpos[j0] - 4.25 * chw
}
if(j0 == 1)
text(pretext.pos, ypos[j0], bars.pretext, adj = 1)
}
for(jk in sublist) {
mv <- mean.list[[jk]]
num <- num.list[[jk]]
sem <- sd.here/sqrt(num)
uniqx <- uniqx.list[[jk]]
if(!legend.sem)
errorbars(uniqx, mv, sem, eps = chw/2)
}
}
if(!plot.some)
return()
note.pos <- ypos[j0] - 1 * chh
if(plot.groups.sem)
note.pos <- note.pos - semax[length(semax)]
if(nchar(legend.note) > 0)
text(xpos[1], note.pos, legend.note, adj = 0, cex = 0.75)
oldmgp <- par(mgp = c(2.75 + ytitladd, 0.5, 0))$oldmgp
if(names(dev.cur()) == "postscript")
mixed.mtext(texts = paste(prefix.y, "LT~v~c.8~.", lt, "~h0~c1~. ", ylab,
sep = ""), side = 2, line = 2.5 + ytitladd, adj = 0.5, cex =
cex.labels)
else mtext(text = paste(prefix.y, "LT", lt, " ", ylab, sep = ""), side = 2,
line = 2.5 + ytitladd, adj = 0.5, cex = cex.labels)
par()$oldmgp
if(length(leg) < nlist) {
cat("Insufficient legends", fill = T)
leg <- c(leg, rep("", nlist - length(leg)))
}
else if(length(leg) > nlist)
cat("NB: There are superfluous legends", fill = T)
if(length(leg) < nlist)
cat("Insufficient legends", fill = T)
j0 <- 0
if(sum(include) > 1) {
for(i in 1:nlist)
if(include[i]) {
j0 <- j0 + 1
points(xpos[j0], ypos[j0], pch = plotch[j0], cex = 0.9, mkh = mkh)
if(poly[i])
lines(c(xpos[j0] + 1.5 * chw, xpos[j0] + (linetype.len + 1.5) *
chw), c(ypos[j0], ypos[j0]), lty = j0)
if(names(dev.cur()) == "postscript")
mixed.text(xpos[j0] + (linetype.len + 2) * chw, ypos[j0], leg[i
], adj = 0, cex = 0.8)
else text(xpos[j0] + (linetype.len + 2) * chw, ypos[j0], leg[i],
adj = 0)
}
}
if(is.null(xlab)) {
xlab <- ""
cat("No setting was provided for xlab\n")
}
if(names(dev.cur()) == "postscript") {
mixed.mtext(texts = paste(xlab, sep = ""), side = 1, line = 2.25, adj =
0.5, cex = cex.labels)
if(maint != "")
mixed.mtext(texts = maint, side = 3, line = 1.75, adj = 0.5, cex =
cex.title)
}
else {
mtext(paste(xlab, sep = ""), 1, line = 2.25, cex = cex.labels)
mtext(text = maint, side = 3, line = 1.75, adj = 0.5, cex = cex.title)
}
oldcex <- par(cex = 0.4)
r.txt <- unix("date +%d/%m/%y' '%H:%M")
l.txt <- paste(unix("pwd"), as.character(as.name(match.call())), collapse
= "")
### docs <- paste(unix("pwd"), as.character(as.name(match.call())),
### date.txt, collapse = "")
if(datestamp) {
mtext(r.txt, side = 1, adj = 1, outer = T, line = mtext.line)
mtext(l.txt, side = 1, adj = 0, outer = T, line = mtext.line)
}
options(width = 80) # docs
invisible()
}
, comment = "08/10/1995")
"plot.resid" <-
structure(function(xx)
{
## Purpose: Plots fitted values against the residuals for aovs, lms, etc
## ----------------------------------------------------------------------
## Modified from:
## ----------------------------------------------------------------------
## Arguments: xx is an aov, glm, lm type object
## ----------------------------------------------------------------------
## Author: Patrick Connolly, Creation date: 26 Aug 2005, 16:14
with(xx, plot(fitted.values, residuals))
}
, comment = "26/08/2005")
"pool.adj" <-
structure(function(x, r, n, phat, cm, cm.strategy = "adjust.later", cm.code = NULL,
cm.allcodes = NULL, plimits = c(0.0025, 0.9995), plotit = "loess",
link = "loglog", xfun = function(x) x, legend = "", x.cm = NULL,
dead.cm = NULL, tot.cm = NULL, xtit = "Time", ytit = NULL,
mkh = 0.025, extrap.low = 0.25, extrap.high = 0.25, span = 1,
deg = 1, lo.lty = 2, lwd = 1, n.lo = 50, title.line = 0.25)
{
### pool.adj is called both by fitconf and by flyplot
###
### extrap.low (extrap.high) is fraction
### of range allowed for extrapolation below (above)
### These are halved if there are just two points
if(!is.null(plotit)) {
plotit[plotit == "loess"] <- "mono.loess"
plotit[plotit == "mono"] <- "mono.loess"
show.curve <- as.logical(match(c("mono.line", "mono.loess"), plotit,
nomatch = 0))
}
else show.curve <- rep(F, 2)
if(!is.null(phat)) {
phat.strategy <- names(phat)
new.lt <- T
lt <- array(NA, length(phat))
nam.lt <- names(phat)
ab <- vector("list", length(phat))
}
else {
lt <- NULL
ab <- NULL
new.lt <- F
phat.strategy <- ""
}
plot.some <- any(show.curve)
if(all(xfun(exp(1:8)) == 1:8)) {
xfun.inv <- exp
takelog <- T
}
else {
xfun.inv <- function(x)
x
takelog <- F
}
if(any(diff(x)) < 0)
stop("** Monotone fitting requires that data is ordered by x-values **")
x0 <- x
r0 <- r
n0 <- n
p0 <- r0/n0
if(any(diff(x) == 0)) {
r <- tapply(r, x, sum)
n <- tapply(n, x, sum)
x <- as.numeric(names(r))
}
if(takelog)
if(any(x <= 0)) {
x.0 <- x > 0
x <- x[x.0]
r <- r[x.0]
n <- n[x.0]
}
g <- link.function(link)
nn <- length(n)
p <- r/n
if(is.null(cm)) {
cat("\nNo control mortality was supplied.\n")
cat("\nFor monotone estimation set cm = 0\n")
}
if(cm.strategy == "adjust.later" & !is.null(phat))
ptarg <- cm + ((1 - cm) * phat)/100
id <- 1:nn
repeat {
i <- 2
while(i <= nn && p[id[i]] >= p[id[i - 1]]) i <- i + 1
if(i <= nn) {
x[id[i - 1]] <- x[id[i]] <- (n[id[i - 1]] * x[id[i - 1]] + n[id[i]] *
x[id[i]])/(n[id[i - 1]] + n[id[i]])
r[id[i - 1]] <- r[id[i]] <- r[id[i]] + r[id[i - 1]]
n[id[i - 1]] <- n[id[i]] <- n[id[i]] + n[id[i - 1]]
p[id[i - 1]] <- p[id[i]] <- r[id[i]]/n[id[i]]
idi <- id[i]
while(i <= nn && id[i] == idi) {
id[i] <- id[i - 1]
i <- i + 1
}
}
else break
}
n1 <- max((1:nn)[p == 0])
if(!is.na(n1))
id <- id[id >= n1]
use <- unique(id)
ux <- x[use]
ux.tran <- xfun(ux)
if(!is.null(phat))
ptarg.interp <- ptarg
r.use <- r[use]
n.use <- n[use]
p.use <- r.use/n.use
if(cm.strategy == "abbott") {
here.use <- p.use >= cm
use <- use[here.use]
ux <- ux[here.use]
ux.tran <- xfun(ux)
p.use <- (p.use[here.use] - cm)/(1 - cm)
}
xx <- cbind(r[use], (n - r)[use])
p.max <- max(p.use)
p.min <- min(p.use)
done.loess <- F
show.line <- show.curve[1]
show.loess <- show.curve[2]
nx <- sum(p.use < 1)
n01 <- sum(p.use > 0 & p.use < 1)
phat.strategy[phat.strategy == "mono"] <- "loess"
phat.strategy[phat.strategy == "mono.loess"] <- "loess"
calc.loess <- any(as.logical(match("loess", phat.strategy, nomatch = 0)))
if(length(ux) >= 4 & nx >= 3 & (show.loess | calc.loess) & n01 >= 2) {
gp <- g(p.use)
df.lo <- data.frame(y = p.use, x = ux.tran, n = n.use)
browser()
u <- switch(link,
identity = gam(y ~ lo(x, span = span, degree = deg), family =
binomial("identity"), data = df.lo, weight = n), logsurv = ,
log = gam(y ~ lo(x), family = binomial("log"), data = df.lo,
weight = n, span = span, degree = deg),
logit = gam(y ~ lo(x, span = span, degree = deg),
family = binomial("logit"), data = df.lo, weight = n),
sqrt = gam(y ~ lo(x, span = span, degree = deg),
family = binomial("sqrt"), data = df.lo, weight = n),
inverse = gam(y ~ lo(x, span = span, degree = deg),
family = binomial("inverse"), data = df.lo, weight = n),
probit = gam(y ~ lo(x, span = span, degree = deg),
family = binomial("probit"), data = df.lo, weight = n),
cloglog = ,
loglog = gam(y ~ lo(x, span = span, degree = deg),
family = binomial("cloglog"), data = df.lo, weight = n))
df.new <- data.frame(x = pretty(ux.tran, n.lo))
hat.u <- predict(u)
u.lo <- loess(y ~ x, data = list(y = hat.u, x = ux.tran), span = span,
deg = deg, span = span) # If no. of points is 4 and span<1,
### the fitted curve has a tendency to oscillate up and down
hat.loess <- predict(u.lo, newdata = df.new)
range.loess <- range(hat.loess, na.rm = T)
done.loess <- T
}
else done.loess <- F
for(i in seq(along = phat)) {
if(n01 < 2)
break
pc <- ptarg[i]
yhat <- g(pc)
get.lt <- as.logical(match(c("line", "loess"), phat.strategy[i], nomatch
= 0))
lt.line <- get.lt[1]
lt.loess <- get.lt[2]
pc.interp <- ptarg.interp[i]
yhat.interp <- g(pc.interp)
alt.line <- F
if(!done.loess || yhat.interp < range.loess[1] || yhat.interp >
range.loess[2])
alt.line <- T
mn.use <- 1:length(use)
if(any(r.use < n.use))
n.ones <- max(mn.use[r.use < n.use]) + 1
else n.ones <- mn.use[1]
take.above <- (r.use > n.use * pc) & (mn.use <= n.ones)
take.low <- r.use <= n.use * pc
if(lt.line | (alt.line & lt.loess))
if(sum(take.above) >= 1) {
mn.above <- mn.use[take.above]
n.above <- min(c(max(mn.above), min(mn.above) + 2))
n.low <- max(c(1, n.above - 3))
if(n.low == 1)
n.above <- min(c(max(mn.above), min(mn.above) + 3))
}
else if(sum(take.low) >= 1) {
mn.low <- mn.use[take.low]
n.above <- max(mn.low)
n.low <- max(c(mn.low[1], n.above - 3))
}
else {
lt.line <- F
alt.line <- F
}
if(lt.line | alt.line) {
p.here <- p.use[n.low:n.above]
not.one <- sum(p.here < 1)
g.here <- g(p.here)
g.interp <- g(pc.interp)
n.span <- n.above - n.low + 1
eps.low <- extrap.low * diff(range(g.here))
eps.high <- extrap.high * diff(range(g.here))
if(n.span == 2) {
eps.low <- eps.low/2
eps.high <- eps.high/2
}
if(n.span >= 2 & not.one > 1)
show.line <- ((max(g.here) > pc.interp - eps.high & min(g.here) <
pc.interp) | (min(g.here) < pc.interp + eps.low &
max(g.here) > pc.interp))
else {
show.line <- F
lt.line <- F
alt.line <- F
}
}
if(lt.line | (alt.line & lt.loess)) {
rn <- xx[n.low:n.above, ]
xv <- xfun(ux[n.low:n.above])
u <- switch(link,
identity = glm(y ~ x, family = binomial("identity"),
data = list(y = rn, x = xv)), logsurv = ,
log = glm(y ~ x, family = binomial("log"),
data = list(y = rn, x = xv)),
logit = glm(y ~ x, family = binomial("logit"),
data = list(y = rn, x = xv)),
sqrt = glm(y ~ x, family = binomial("sqrt"),
data = list(y = rn, x = xv)),
inverse = glm(y ~ x, family = binomial("inverse"),
data = list(y = rn, x = xv)),
probit = glm(y ~ x, family = binomial("probit"),
data = list(y = rn, x = xv)),
cloglog = ,
loglog = glm(y ~ x, family = binomial("cloglog"),
data = list(y = rn, x = xv)))
if(!is.na(u$coef[2]) & (lt.line | (alt.line & lt.loess)) & u$coef[2] >
0) {
lt[i] <- xfun.inv((yhat - u$coef[1])/u$coef[2])
ab[[i]] <- u$coef[1:2]
nam.lt[i] <- "line"
}
else {
lt[i] <- NA
ab <- NULL
}
}
if(done.loess && lt.loess && min(hat.loess, na.rm = T) <= yhat && max(
hat.loess, na.rm = T) >= yhat) {
here.hat <- !is.na(hat.loess)
lt[i] <- xfun.inv(approx(x = hat.loess[here.hat],
y = df.new$x[ here.hat], xout = yhat.interp)$y)
nam.lt[i] <- "mono"
}
}
if(new.lt & !is.null(phat)) {
p.txt <- paste("LT", phat, "=", format(round(lt, 1)), sep = "")
leg <- paste(legend, " [", paste(p.txt, collapse = "; "), "]", sep = "")
names(lt) <- nam.lt
}
else leg <- legend
if(plot.some) {
if(!is.null(x.cm)) {
x0 <- c(paste(x0), x.cm)
r0 <- c(r0, dead.cm)
n0 <- c(n0, tot.cm)
}
### NB x0, r0, n0 do not appear below
simplot(x = x0, resp = r0, tot = n0, fun = link, cm.strategy =
cm.strategy, cm = cm, cm.code = cm.code,
cm.allcodes = cm.allcodes, xtit = xtit, ytit = ytit, mkh = mkh,
xfun = xfun, main = leg, title.line = title.line,
plimits = plimits, points.lwd = lwd)
yp <- g(plimits)
for(j in seq(along = ab))
abline(ab[[j]], lty = i)
if(done.loess) {
here <- hat.loess >= yp[1] & hat.loess <= yp[2]
lines(df.new$x[here], hat.loess[here], lty = lo.lty, lwd = lwd)
}
}
lt
}
, comment = "01/10/2004")
"ppaste" <-
structure(function(...)
{
paste(..., sep = "")
}
, comment = "11/02/1998")
"print.char.matrix" <-
structure(function (x, file = "", col.name.align = "cen", col.txt.align = "right",
cell.align = "cen", hsep = "|", vsep = "-", csep = "+", row.names = TRUE,
col.names = FALSE, append = FALSE, top.border = TRUE, left.border = TRUE)
{
### To print a data frame or matrix to a text file or screen
### and having names line up with stacked cells
###
### First, add row names as first column (might be removed later)
rownames <- dimnames(x)[[1]]
x <- cbind(rownames, x)
dimnames(x)[[1]] <- seq(nrow(x))
### Set up some padding functions:
###
pad.left <- function(z, pads) {
### Pads spaces to left of text
padding <- paste(rep(" ", pads), collapse = "")
paste(padding, z, sep = "")
}
pad.mid <- function(z, pads) {
### Centres text in available space
padding.right <- paste(rep(" ", pads%/%2), collapse = "")
padding.left <- paste(rep(" ", pads - pads%/%2), collapse = "")
paste(padding.left, z, padding.right, sep = "")
}
pad.right <- function(z, pads) {
### Pads spaces to right of text
padding <- paste(rep(" ", pads), collapse = "")
paste(z, padding, sep = "")
}
### (Padding happens on the opposite side to alignment)
pad.types <- c("left", "mid", "right")
names(pad.types) <- c("right", "cen", "left")
pad.name <- pad.types[col.name.align]
pad.txt <- pad.types[col.txt.align]
pad.cell <- pad.types[cell.align]
### Padding character columns
### Need columns with uniform number of characters
pad.char.col.right <- function(y) {
### For aligning text to LHS of column
col.width <- nchar(y)
biggest <- max(col.width)
smallest <- min(col.width)
padding <- biggest - col.width
out <- NULL
for (i in seq(y)) out[i] <- pad.right(y[i], pads = padding[i])
out
}
pad.char.col.left <- function(y) {
### For aligning text to RHS of column
col.width <- nchar(y)
biggest <- max(col.width)
smallest <- min(col.width)
padding <- biggest - col.width
out <- NULL
for (i in seq(y)) out[i] <- pad.left(y[i], pads = padding[i])
out
}
pad.char.col.mid <- function(y) {
### For aligning text to centre of column
col.width <- nchar(y)
biggest <- max(col.width)
smallest <- min(col.width)
padding <- biggest - col.width
out <- NULL
for (i in seq(y)) out[i] <- pad.mid(y[i], pads = padding[i])
out
}
### which functions to use this time.
pad.name.fn <- get(paste("pad.", pad.name, sep = ""))
pad.txt.fn <- get(paste("pad.char.col.", pad.txt, sep = ""))
pad.cell.fn <- get(paste("pad.", pad.cell, sep = ""))
###
### Remove troublesome factors
x <- as.data.frame(x)
fac.col <- names(x)[sapply(x, is.factor)]
for (i in fac.col) x[, i] <- I(as.character(x[, i]))
### ARE ANY LINE BREAKS IN ANY COLUMNS?
break.list <- list()
for (i in seq(nrow(x))) {
x.i <- unlist(x[i, ])
rows.i <- sapply(strsplit(unlist(x[i, ]), "\n"), length)
rows.i[rows.i < 1] <- 1
break.list[[i]] <- rows.i
}
break.row <- sapply(break.list, function(x) any(x > 1))
names(break.row) <- seq(nrow(x))
xx <- x
if (any(break.row)) {
### add in extra row/s
xx <- NULL
reprow <- lapply(break.list, unique)
for (k in seq(nrow(x))) {
x.k <- unlist(x[k, ])
x.k[x.k == ""] <- " "
if (break.row[k]) {
l.k <- strsplit(x.k, "\n")
add.blanks <- max(break.list[[k]]) - break.list[[k]]
names(l.k) <- names(add.blanks) <- seq(length(l.k))
if (any(add.blanks > 0)) {
for (kk in names(add.blanks[add.blanks > 0]))
l.k[[kk]] <- c(l.k[[kk]], rep(" ", add.blanks[kk]))
}
l.k.df <- as.data.frame(l.k)
names(l.k.df) <- names(x)
xx <- rbind(xx, as.matrix(l.k.df))
}
else xx <- rbind(xx, x.k)
}
row.names(xx) <- paste(rep(row.names(x), sapply(reprow,
max)), unlist(reprow), sep = ".")
### Make an index for the rows to be printed
rn <- row.names(xx)
rnb <- strsplit(rn, "\\.")
rpref <- codes(factor(sapply(rnb, function(z) z[1])))
}
else rpref <- seq(nrow(x))
x <- as.data.frame(xx)
### Character columns need different treatment from numeric columns
char.cols <- sapply(x, is.character)
if (any(char.cols))
x[char.cols] <- sapply(x[char.cols], pad.txt.fn)
### Change numeric columns into character
if (any(!char.cols))
x[!char.cols] <- sapply(x[!char.cols], format)
### now all character columns each of which is uniform element width
###
### Lining up names with their columns
### Sometimes the names of columns are wider than the columns they name,
### sometimes vice versa.
###
names.width <- nchar(names(x))
if (!col.names)
names.width <- rep(0, length(names.width))
cell.width <- sapply(x, function(y) max(nchar(as.character(y))))
### (the width of the characters in the cells as distinct
### from their names)
name.pads <- cell.width - names.width
cell.pads <- -name.pads
name.pads[name.pads < 0] <- 0
cell.pads[cell.pads < 0] <- 0
pad.names <- name.pads > 0
pad.cells <- cell.pads > 0
### Pad out the column names if necessary:
if (any(pad.names)) {
stretch.names <- names(x)[pad.names]
for (i in stretch.names) {
names(x)[names(x) == i] <- pad.name.fn(i, name.pads[i])
}
}
### likewise for the cells and columns
if (any(pad.cells)) {
stretch.cells <- names(x)[pad.cells]
for (j in stretch.cells) x[, j] <- pad.cell.fn(x[, j],
cell.pads[j])
}
### Remove row names if not required
if (!row.names)
x <- x[-1]
### Put the column names on top of matrix
if (col.names) {
mat2 <- rbind(names(x), as.matrix(x))
if(row.names) # Don't remove first column name if there are no rownames
mat2[1, 1] <- paste(rep(" ", cell.width[1]), collapse = "") # get rid of "rownames" from first column
}
else mat2 <- as.matrix(x)
mat.names.width <- nchar(mat2[1, ])
### character string to separate rows
space.h <- ""
for (k in seq(mat.names.width)) {
space.h <- c(space.h, rep(vsep, mat.names.width[k]),
csep)
}
line.sep <- paste(c(ifelse(left.border, csep, ""), space.h),
collapse = "")
if (col.names)
rpref <- c(0, rpref, 0)
else rpref <- c(rpref, 0)
### print to screen or file
if (top.border) {
write(line.sep, file = file, append = append)
append <- TRUE
}
for (i in 1:nrow(mat2)) {
if (left.border)
write(paste(paste(c("", mat2[i, ]), collapse = hsep),
hsep, sep = ""), file = file, append = append)
else write(paste(paste(mat2[i, ], collapse = hsep), hsep,
sep = ""), file = file, append = append)
append <- TRUE
### print separator if row prefix is not same as next one
if (rpref[i] != rpref[i + 1])
write(line.sep, file = file, append = TRUE)
}
}
, comment = "19/10/2006")
"pseudo.f" <-
structure(function(y, show.aov = F)
{
# Calculates pseudo F-statistics from an anova of a glm object y
x <- summary(y)
bits <- unlist(dimnames(x)[[1]])
resid.df <- x[rev(bits)[1], "Resid. Df"]
out.df <- NULL
for(i in bits[-1]) {
f.stat <- (x[i, "Deviance"]/x[i, "Df"])/(x[rev(bits)[1], "Resid. Dev"]/
resid.df)
prf <- 1 - round(pf(f.stat, x[i, "Df"], resid.df), 3)
out.df <- rbind(out.df, data.frame(Factor = i, "F value" = round(f.stat,
4), Df1 = x[i, "Df"], Df2 = resid.df, Prob = prf))
}
if(show.aov) {
print(x)
cat("\n")
}
out.df
}
, comment = "30/04/2001")
"qik.sum" <-
structure(function(df, cols = 1:4, ordered = F)
{
### Does a quick summary, first coercing specified cols to factors and making a matrix
### into a dataframe.
#
### Requires make.factors() also
df <- as.data.frame(df)
summary(make.factors(df, cols, ordered = ordered))
}
, comment = "13/08/1998")
"relevel.default" <-
structure(function(x, ref, ...)
stop("relevel only for factors")
, comment = "03/09/1999")
"relevel.factor" <-
structure(function(x, ref, ...)
{
lev <- levels(x)
nlev <- length(lev)
if(is.character(ref))
ref <- match(ref, lev)
if(is.na(ref))
stop("ref must be an existing level")
if(ref < 1 || ref > nlev)
stop(paste("ref =", ref, "must be in 1 :", nlev))
factor(x, levels = lev[c(ref, seq(along = lev)[ - ref])], labels = lev[c(
ref, seq(along = lev)[ - ref])])
}
, comment = "03/09/1999")
"remove.shading" <-
structure(function()
{
# To get rid of shading on shingles above plots
strip.background <- trellis.par.get("strip.background")
strip.background$col <- 0
trellis.par.set("strip.background", strip.background)
strip.shingle <- trellis.par.get("strip.shingle")
strip.shingle$col <- 0
trellis.par.set("strip.shingle", strip.shingle)
}
, comment = "09/07/2001")
"rms" <-
structure(function(x, na.rm = T, zero.rm = F)
{
if(na.rm)
x <- x[!is.na(x)]
if(zero.rm)
x <- x[abs(x) > 0]
sqrt(mean(x^2))
}
, comment = "30/11/1996")
"robust.deviance" <-
structure(function(phat, n, r, min.expected = 1)
{
### phat : vector of fitted probabilities
### n : binomial totals
### r : binomial successes
### min.expected : groups combined so that minimum number expected
eps <- 1.0000000000000001e-15
if(any(n < 0))
stop("can't have 0 or negative binomial sample size")
rhat <- phat * n
ord <- order(phat)
ind <- match(1:length(n), ord)
cum.rhat <- cumsum(rhat[ord])
g1 <- cum.rhat < min.expected
g1 <- c(T, g1)[1:length(n)]
revcum.rhat <- rev(cumsum(rev((n - rhat)[ord])))
g3 <- revcum.rhat < min.expected
g3 <- c(g3, T)[-1]
g2 <- !(g1 | g3)
df <- sum(g2) ### Really, sum(g2)+2-2; the endpoints are always F
grp1 <- rep(NA, length(n))
grp1[g1] <- 1
grp1[g2] <- 2:(sum(g2) + 1)
grp1[g3] <- sum(g2) + 2
grp <- grp1[ind]
grp
r.grp <- tapply(r, grp, sum)
n.grp <- tapply(n, grp, sum)
phat.grp <- tapply(phat * n, grp, sum)/tapply(n, grp, sum)
qhat.grp <- 1 - phat.grp
pobs <- r.grp/n.grp
qobs <- 1 - pobs
D1 <- ifelse(pobs < eps, log((pobs/phat.grp)^r.grp), r.grp * log(pobs/
phat.grp))
D2 <- ifelse(qobs < eps, log((qobs/qhat.grp)^(n.grp - r.grp)), (n.grp -
r.grp) * log(qobs/qhat.grp))
dev <- 2 * sum(D1 + D2)
list(dev = dev, df = df)
}
, comment = "26/06/1994")
"s.tapply" <-
structure(function(x, y, z, ...)
{
# Counteracts tapply's propensity to return in sorted index order
# instead of the order in which they were
tapply(x, y, z, ...)[(unique(y))]
}
, comment = "30/04/2001")
"save.new" <-
structure(function(neuobj, gem.pos = NULL)
{
# For saving a new function to the Gems database
if(is.null(gem.pos)) gem.pos <- grep("Gems", search())
fun.name <- deparse(substitute(neuobj))
assign(fun.name, neuobj, pos = gem.pos)
save(list = ls(all = TRUE, pos = gem.pos), file = substring(search()[
gem.pos], 6))
}
, comment = "27/04/2001")
"sem" <-
structure(function(x, na.rm = F)
{
if(na.rm)
x <- x[!is.na(x)]
if(length(x) == 0)
NA
else sqrt(var(x)/length(x))
}
, comment = "09/09/1995")
"set.screens" <-
structure(function(across = 2, down = 3, left.margin = 10, top.margin = 20, height = 80,
width = 90, ygap = -10, xgap = -1, overwrite = T, centre.x = T, centre.y =
T, byrow = T)
{
# Uses split.screen for setting up screens: alternative to par(mfrow/mfcol)
uu <- ps.region/72 * 25.4 # Change points to mm
x.all <- uu[3] - uu[1]
y.all <- uu[4] - uu[2] #
# Check if there is room for what has been set
spare.x <- x.all - across * width - xgap * (across - 1)
if(spare.x < 0)
stop("\nInsufficient space:\nCheck left.margin, width and xgap")
spare.y <- y.all - down * height - ygap * (down - 1)
if(spare.y < 0) stop(
"\nInsufficient space:\nCheck top.margin, height and ygap") #
# Set up coordinates for the screens
top <- top.margin + ifelse(centre.y, spare.y/(down + 1), 0)
left <- left.margin + ifelse(centre.x, spare.x/(across + 1), 0)
N <- across * down
first.x.corners <- c(left - uu[1], left - uu[1] + width)/x.all
first.y.corners <- c(uu[4] - top - height, uu[4] - top)/y.all #browser()
if(N == 1)
pos.mat <- matrix(c(first.x.corner, first.y.corner), nr = 1)
else {
pos.mat <- NULL
if(byrow) {
for(j in 1:down) {
y.corners.j <- first.y.corners - (ygap + height)/y.all * (j - 1)
for(i in 1:across) {
x.corners.i <- first.x.corners + ((width + xgap) * (i - 1))/
x.all
pos.mat <- rbind(pos.mat, c(x.corners.i, y.corners.j))
}
}
}
else {
for(i in 1:across) {
x.corners.i <- first.x.corners + ((width + xgap) * (i - 1))/x.all
for(j in 1:down) {
y.corners.j <- first.y.corners - (ygap + height)/y.all * (j - 1
)
pos.mat <- rbind(pos.mat, c(x.corners.i, y.corners.j))
}
}
}
split.screen(pos.mat, erase = !overwrite)
}
}
, comment = "29/11/1997")
"show.nas" <-
structure(function(df)
{
### identifies which columns of a dataframe have nas.
#
sapply(df, function(x)
any(is.na(x)))
}
, comment = "01/12/1997")
"simplot" <-
structure(function(x, resp, tot, fun, cm.strategy = "adjust.later", cm = 0, cm.code =
NULL, cm.allcodes = NULL, xtit = "Time in Min", main = "", mkh = 0.025,
xfun = function(x) x, takelog = F, line = F, ab = c(NA, NA),
clip = NULL, xlimits = c(0, 0),
xlabels = c(0, 0), offset = 0, title.line = 1, plimits = c(0.0025, 0.9995),
ytit = NULL, new = F, lty = 1, plot.char = c(0, 1, 2, 5, 6:16),
id = NULL, points.lwd = 1)
{
if(missing(mkh))
mkh <- par()$cin[2]/3
if(new)
par(new = T)
g <- link.function(fun)
m <- length(x)
if(is.null(id))
id <- 1
if(length(id) == 1)
id <- rep(id, length(x))
if(takelog)
xfun <- log
if(all(xfun(exp(1:8)) == (1:8))) takelog <- T
### *****************************************************************
if(is.null(cm.code))
cm.code <- 0
if(is.null(cm.allcodes))
cm.allcodes <- cm.code
other.rows <- match(x, cm.allcodes, nomatch = 0)
cm.here.codes <- unique(cm.allcodes[other.rows])
cmrows <- as.logical(other.rows)
if(cm.strategy == "adjust.later" & !takelog)
here.trt <- !cmrows | x == "0"
else here.trt <- !cmrows
x.trt <- as.numeric(x[here.trt])
resp.trt <- resp[here.trt]
tot.trt <- tot[here.trt]
id.trt <- id[here.trt]
if(missing(xlimits))
xlimits <- range(x.trt)
if(takelog) {
if(any(x.trt + offset < 0)) {
cat(
"*** Warning: The following x-values have been omitted because\n x"
)
if(offset != 0)
cat("+", offset)
cat(" < 0:", fill = T)
cat(format(x.trt[x.trt + offset < 0]), fill = T)
}
if(any(x.trt + offset <= 0)) {
x.logable <- x.trt + offset > 0
resp.trt <- resp.trt[x.logable]
tot.trt <- tot.trt[x.logable]
x.trt <- x.trt[x.logable]
id.trt <- id.trt[x.logable]
}
if(xlimits[1] + offset <= 0)
xlimits[1] <- min(x.trt[x.trt + offset > 0], na.rm = T)
}
if(diff(range(xlabels)) == 0)
xlabels <- pretty(xlimits)
if(takelog)
xlabels <- xlabels[xlabels > 0]
cat(xlabels, " ") # xlabels specifies numeric labels
p <- resp.trt/tot.trt
prange <- range(p)
if(prange[1] < 0 | prange[2] > 1) {
cat("Error: Smallest p is", prange[1], fill = T)
cat(" Largest p is", prange[2], fill = T)
break
}
cm0 <- switch(cm.strategy,
adjust.later = 0,
abbott = cm)
here <- p > cm0 & p < 1
p.adj <- (p - cm0)/(1 - cm0)
ylow <- g(plimits[1])
yhigh <- g(plimits[2])
eps <- (yhigh - ylow) * 0.01
ylow <- ylow - eps
yhigh <- yhigh + eps
labp <- c(1, 5, 25, 50, 75, 95, 99, 99.9)
labp <- labp[labp/100 > plimits[1] & labp/100 < plimits[2]]
aty <- g(labp/100)
dist <- yhigh - ylow
if(is.null(ytit)) {
fun.text <- fun
ytit <- paste("% Mortality,", fun.text, "scale")
}
par(mkh = mkh)
u.id <- unique(id.trt)
xlim.here <- xfun(xlimits + offset)
if(takelog)
xlim.here <- xlim.here - c(0.115 * diff(xlim.here), 0)
x.cpos <- xlim.here[1]
if(sum(here) > 0) {
plot(xfun(x.trt + offset)[here], g(p.adj[here]), xlim = xlim.here, ylim
= c(ylow, yhigh), xlab = "", ylab = ytit, axes = F, type = "n")
}
else plot(xfun(xlimits + offset), c(ylow, yhigh), xlim = xlim.here, xlab =
"", ylab = ytit, axes = F, type = "n")
midhigh <- 0.5 * (yhigh + par()$usr[4])
midlow <- 0.5 * (ylow + par()$usr[3])
chw <- par()$cxy[1]
chh <- par()$cxy[2]
j <- 0
for(i in u.id) {
j <- j + 1
here.i <- here & i == id.trt
g.i <- g(p.adj[here.i])
if(any(here.i) & any(g.i < ylow | g.i > yhigh)) {
cat("\n*** One or more points is outside limits set by plimits.***\n"
)
cat("\n*** plimits =", format(round(plimits, 4)), "***", "\n")
cat("\n*** Points are:\n")
if(any(g.i < ylow)) {
cat(format(round(p.adj[g.i < ylow], 4)), "\n")
g.i[g.i < ylow] <- ylow
}
if(any(g.i > yhigh)) {
cat(format(round(p.adj[g.i > yhigh], 4)), "\n")
g.i[g.i > yhigh] <- yhigh
}
}
if(any(here.i))
points(xfun(x.trt + offset)[here.i], g(p.adj[here.i]), pch =
plot.char[j], mkh = mkh)
x100 <- x.trt[p == 1 & i == id.trt]
x0 <- x.trt[p == 0 & i == id.trt]
x0.100 <- c(x0, x100)
y0.100 <- c(rep(midlow, length(x0)), rep(midhigh, length(x100)))
if(length(x0.100) > 0) {
points(xfun(x0.100 + offset), y0.100, pch = ".", mkh = mkh)
points(xfun(x0.100 + offset), y0.100, pch = plot.char[j], mkh = mkh)
}
}
if((length(cm.here.codes) > 0) & cm.strategy == "adjust.later") {
u.cm <- unique(cm.here.codes)
for(cm.i in u.cm) {
### here.cm <- x == cm.i ## Try character 8/9/04
here.cm <- as.character(cm.i) == x
p.cm <- resp[here.cm]/tot[here.cm]
here.cm <- here.cm[p.cm > 0]
p.cm <- p.cm[p.cm >= 0]## add in equality
#browser()
if(sum(here) > 0)
text(rep(x.cpos, length(p.cm)), g(p.cm), rep(cm.i, length(p.cm)))
}
}
if(line) {
b <- ab[2]
a <- ab[1]
if(is.na(b))
b <- -999
if(is.null(clip))
clip <- xlimits
if(b > 0)
clipline(xfun(clip + offset), c(ylow, yhigh - chh/3), a, b)
}
xmax <- max(x.trt)
epsx <- 0.0125 * diff(par()$usr[1:2])
epsy <- 0.0125 * diff(par()$usr[3:4])
if(takelog) {
x.cut1 <- x.cpos + epsx
x.cut2 <- x.cpos + 2.5 * epsx
here.pos <- xfun(xlabels + offset) > x.cut2 & xlabels <= xmax
}
else here.pos <- xlabels <= xmax
xpos <- xlabels[here.pos]
if(length(xpos) <= 2)
xpos <- xlabels[xlabels <= xmax]
### if(FALSE) {
### if(names(dev.cur()) == "postscript" & nchar(main) > 0)
### mixed.mtext(side = 3, line = title.line, texts = main, cex = par()$cex *
### 1.2, adj = 0.5)
### else
### }
mtext(side = 3, line = title.line, text = main, cex = par()$cex * 1.2)
title(xlab = xtit)
axis(1, at = xfun(xpos + offset), labels = F)
xusr1 <- par()$usr[1]
axis(1, at = xfun(xpos + offset), labels = paste(xpos), mgp = c(2, 0.5, 0))
if(takelog) {
axis(1, at = c(par()$usr[1], x.cut1), labels = F, tck = 0, mgp = c(2,
0.5, 0))
axis(1, at = c(x.cut2, par()$usr[2]), labels = F, tck = 0, mgp = c(2,
0.5, 0))
lines(x.cut1 + c( - epsx, epsx), par()$usr[3] + 2 * c( - epsy, epsy),
xpd = T)
lines(x.cut2 + c( - epsx, epsx), par()$usr[3] + 2 * c( - epsy, epsy),
xpd = T)
}
else axis(1, at = par()$usr[1:2], labels = F, tck = 0, mgp = c(2, 0.5, 0))
pr <- c(1, 5, 20, 50, 80, 99)
extreme.pr <- c(0, 100)
npr <- length(pr)
if(all(g(pr/100) == pr/100))
identity <- T
else identity <- F
if(identity)
pr <- c(0, 25, 50, 75, 100)
else pr <- pr[g(pr/100) > ylow & g(pr/100) < yhigh]
axis(2, at = g(pr/100), labels = paste(pr))
if(identity) {
midhigh <- g(1)
axis(2, at = (0:20)/20, labels = F)
}
if(!identity) {
axis(2, at = c(ylow + 2 * epsy, yhigh - 2 * epsy), labels = F, tck = 0)
for(k in 1:2) {
y.end <- c(ylow, yhigh)[k]
axis(2, at = c(y.end, par()$usr[2 + k]), labels = F, tck = 0)
}
kk <- (1:2)[c(length(x0) > 0, length(x100) > 0)]
for(k in kk) {
y.mid <- c(midlow, midhigh)[k]
axis(2, at = y.mid, tck = 0, labels = paste(extreme.pr[k]))
y.end <- c(ylow, yhigh)[k]
k <- k + 1
abline(h = y.end, lty = 3)
lines(xusr1 + c( - epsx, epsx), c(y.end - 1.5 * epsy, y.end + 1.5 *
epsy), xpd = T)
## 1/5/06 fixed up long running error with this: ifelse(k == 2, 4, 0) * epsy +
lines(xusr1 + c( - epsx, epsx), ifelse(k == 2, 4, 0) * epsy +
c(y.end - 3.5 * epsy, y.end - 0.5 * epsy), lty = lty, xpd = T)
}
}
else axis(2, at = par()$usr[3:4], labels = F, tck = 0)
box(bty = "7")
}
, comment = "01/05/2006")
"stdev" <-
structure(function(x, na.rm = F)
{
if(na.rm)
x <- x[!is.na(x)]
if(length(x) == 0)
NA
else sqrt(var(x))
}
, comment = "21/09/1998")
"strip.trail" <-
structure(function(x, trail = " ")
{
### This function strips away trailing characters of specified value
### from character vectors: everything from the first occurrence of the trail
### character is removed
begin.spaces <- find.first(x, trail)
need.strip <- !is.na(begin.spaces)
x[need.strip] <- substring(x[need.strip], 1, begin.spaces[need.strip] - 1)
x
}
, comment = "23/03/1997")
"subsets" <-
structure(function(n, r, v = 1:n) # works in S or R
{
## Purpose: Gets permutations of vectors
## ----------------------------------------------------------------------
## Modified from: Rnews (Vol 2 I think)
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Bill Venables , Creation date: 6 Jun 2000
if(r <= 0) vector(mode(v), 0) else
if(r >= n) v[1:n] else {
rbind(cbind(v[1], Recall(n-1, r-1, v[-1])),
Recall(n-1, r, v[-1]))
}
}
, comment = "06/06/2000")
"summarize" <-
structure(function(data, y.cols = NULL, split.cols = NULL, rep.col = "Rep", file = "",
append = FALSE, help = FALSE, bp.width = 30, between = FALSE, dec = 3,
select = c("num", "mean", "sem", "min", "max", "boxplots"),
df.out = FALSE, pad = "mid", show.type = FALSE)
{
### Keep Stephen's beginning:
### Requires text.bp(), df.summary(), and df.to.file()
summary.names <- c("N", "Missing", "Mean", "Median", "Trimmed Mean",
"Std. Dev.", "SEM", "Min", "Max", "1st Quartile",
"3rd Quartile", "Sum", "Boxplots")
summary.types <- c("num", "missing", "mean", "median", "trm", "stdev",
"sem", "min", "max", "q1", "q3", "sum", "boxplots")
help.df <- (data.frame(Select.this = summary.types, to.get.this =
summary.names))
if(help) {
cat("\n In the argument \"select\"\n")
print(help.df)
return()
}
### Make character vectors of y.cols and split.cols if necessary:
if(is.numeric(y.cols))
y.cols <- names(data)[y.cols]
if(is.numeric(split.cols)) split.cols <- names(data)[split.cols] #
### If between reps are needed, first find means for each rep
if(between & show.type) {
if(file == "")
cat("\n\t16:12, 16 Nov 2006 (NZDT)16:12, 16 Nov 2006 (NZDT)~ BETWEEN REPLICATES 16:12, 16 Nov 2006 (NZDT)16:12, 16 Nov 2006 (NZDT)~\n\n")
else write("\n\t16:12, 16 Nov 2006 (NZDT)16:12, 16 Nov 2006 (NZDT)~ BETWEEN REPLICATES 16:12, 16 Nov 2006 (NZDT)16:12, 16 Nov 2006 (NZDT)~\n\n", file =
file, append = append)
append <- TRUE
data <- df.summary(data = data, y.cols = y.cols,
split.cols = c(split.cols, rep.col), decimals = 6)
}
else {
if(show.type){
if(file == "")
cat("\n\t********* WITHIN REPLICATES *********\n\n")
else write("\n\t********* WITHIN REPLICATES *********\n\n", file = file,
append = append)
append <- TRUE
}
}
s.tapply <- function(x, y, z, ...)
tapply(x, y, z, ...)[(unique(y))]
q1 <- function(x)
quantile(x, na.rm = TRUE)[2]
q3 <- function(x)
quantile(x, na.rm = TRUE)[4]
no.missing <- function(x)
length(x[is.na(x)])
len.rm <- function(x)
length(x[!is.na(x)])
round.stdev <- function(x, y = dec)
round(stdev(x, na.rm = TRUE), y)
round.mean <- function(x, y = dec)
round(mean(x, na.rm = TRUE), y)
round.sem <- function(x, y = dec)
round(sem(x, na.rm = TRUE), y)
round.max <- function(x, y = dec)
round(max(x, na.rm = TRUE), y)
round.min <- function(x, y = dec)
round(min(x, na.rm = TRUE), y)
round.median <- function(x, y = dec)
round(median(x, na.rm = TRUE), y)
add.boxplot <- function(x, y.range = range.j, width = bp.width)
text.bp(x, y.range, width) #
### Functions are now set up:
#
### Make named vector of functions to be used:
use.functions <- c("len.rm", "no.missing", "round.mean", "round.median",
"trm", "round.stdev", "round.sem", "round.min", "round.max",
"q1", "q3", "sum", "add.boxplot")
names(use.functions) <- summary.types
names(summary.names) <- summary.types #
### Use dataframes instead of matrices:
if(is.matrix(data))
data <- unfactor(as.data.frame(data))
else data <- unfactor(data) ### avoid pesky factors
data <- df.sort(data, rev(split.cols)) #
### Create index and adjust for single split columns if necessary:
indx.df <- data.frame(data[, split.cols])
names(indx.df) <- split.cols ### Necessary for single split.cols
attach(indx.df)
indx <- NULL
for(i in split.cols)
indx <- paste(indx, get(i), sep = ":")
detach("indx.df")
out.df <- data.frame(indx.df[match(unique(indx), indx), split.cols])
names(out.df) <- split.cols ### Necessary for single split.cols
dimnames(out.df) <- list(paste(1:nrow(out.df)), names(out.df))
for(j in y.cols) {
range.j <- range(data[, j], na.rm = TRUE)
for(k in select)
out.df[[summary.names[k]]] <- as.vector(s.tapply(data[, j], indx,
get(use.functions[k])))
if(df.out) {
names(out.df) <- c(split.cols, paste(j, select, sep = "."))
return(out.df)
}
else {
### Get rid of weird row number names:
dimnames(out.df)[[1]] <- seq(nrow(out.df))
out.df <- unfactor(out.df)
numeric.cols <- !sapply(out.df, function(x)
any(is.na(as.numeric(x))))
## add in as.vector to cope with R-1.8.1 26/11/03 (and above)
out.df[numeric.cols] <- as.vector(lapply(out.df[numeric.cols], as.numeric))
underlines <- paste(rep("=", nchar(j)), collapse = "")
if(file == "") {
cat(paste("\n ", j, "\n", sep = ""))
cat(paste(" ", underlines, "\n", sep = ""))
print(out.df)
}
else {
write(paste("\n", j), file = file, append = append)
append <- TRUE # hereafter appending will always be wanted
write(paste(" ", underlines, sep = ""), file = file, append = TRUE)
df.to.file(out.df, file = file, append = append, pad = pad)
}
}
}
invisible()
}
, comment = "17/02/2004")
"tab.file" <-
structure(function(x, file, rounding = 3, ...)
{
### To write a dataframe or matrix to a tab delimited file
if(!is.null(rounding)) {
if(is.data.frame(x))
x <- as.matrix(x)
x <- round(x, rounding)
}
write.table(x, file, sep = "\t", ...)
}
, comment = "25/12/1997")
"tab.mat" <-
structure(function(mat)
{
rownames <- dimnames(mat)[[1]]
colnames <- dimnames(mat)[[2]]
cat(c("\t", paste(colnames, collapse = "\t")), "\n")
for(i in 1:dim(mat)[1]) {
cat(paste(c(rownames[i], mat[i, ]), collapse = "\t"), "\n")
}
}
, comment = "26/06/1995")
"table.all" <-
structure(function(x, own.seq = seq(0, 3, by = 0.5), minx = 0, maxx = NULL)
{
### returns a table that indicates a zero if the value is absent
### 1a) if else when specifying own sequence
if(!is.null(own.seq)) {
same.intervals <- mean(diff(own.seq))
if(all(diff(own.seq) == same.intervals)) {
table(c(x, own.seq)) - 1
}
else {
stop("Sequence steps not all the same")
}
}
else {
### 1b) when own sequence is NULL
### 2) if else when specifying a maximum for x
if(!is.null(maxx)) {
table(c(x, minx:maxx)) - 1
}
else {
table(c(x, minx:max(x, na.rm = T))) - 1
}
}
}
, comment = "19/04/1997")
"ternary" <-
structure(function(X, pch = par("pch"), lcex = 1, add = FALSE,
ord = 1:3, space = 1, ...)
{
### Purpose: Draws ternary plots
### ----------------------------------------------------------------------
### Modified from: MASS Skye dataset
### ----------------------------------------------------------------------
### Arguments: space: how much space needed to allow for long apex names?
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 7 Apr 2004, 10:01
X <- as.matrix(X)
if(any(X) < 0) stop("X must be non-negative")
s <- drop(X %*% rep(1, ncol(X)))
if(any(s <= 0)) stop("each row of X must have a positive sum")
if(max(abs(s - 1)) > 1e-6) {
warning("row(s) of X will be rescaled")
X <- X / s
}
X <- X[, ord]
s3 <- sqrt(1/3)
if(!add)
{
oldpty <- par("pty")
on.exit(par(pty = oldpty))
par(pty = "s")
plot(c(-s3, s3) * space, c(0.5 - s3, 0.5 + s3) * space, type = "n",
axes = FALSE, xlab = "", ylab = "")
polygon(c(0, -s3, s3), c(1, 0, 0), density = 0)
lab <- NULL
if(!is.null(dn <- dimnames(X))) lab <- dn[[2]]
if(length(lab) < 3) lab <- as.character(1:3)
eps <- 0.05 * lcex
text(c(0, s3 + eps * 0.7, -s3 - eps * 0.7),
c(1 + eps, -0.1 * eps, - 0.1 * eps), lab, cex = lcex)
}
points((X[, 2] - X[, 3]) * s3, X[, 1], pch = pch, ...)
}
, comment = "15/04/2004")
"text.bp" <-
structure(function(x, y.range, width = 45, no.data.string = "NA")
{
### Does a "text boxplot" to use in summarize()
x <- x[!is.na(x)]
if(length(x) < 1)
return(no.data.string)
x.quan <- (boxplot(x, plot = F)$stats[, 1])
x.outliers <- boxplot(x, plot = F)$out
#browser()
if(length(x.outliers)>0)
out.pos <- round(approx(seq(y.range[1], y.range[2], length = width), 1:
width, xout = x.outliers)$y)
qps <- round(approx(seq(y.range[1], y.range[2], length = width), 1:width,
xout = x.quan)$y) #
### Have quantile positions...
### Make allowances for overlapping positions of quantiles
#
qps.start <- qps
if(prod(diff(qps[-3])) == 0) {
if(qps[2] - qps[1] < 1)
qps[2] <- qps[1] + 1
if(qps[5] - qps[4] < 1)
qps[4] <- qps[4] - 1
if(qps[4] - qps[2] < 1) {
if(qps[5] < width)
qps[4:5] <- qps[4:5] + 1
else qps[1:2] <- qps[1:2] - 1
}
}
bar.vec <- rep(" ", width)
bar.vec[qps[1]] <- "("
bar.vec[qps[2]] <- "["
bar.vec[qps[4]] <- "]"
bar.vec[qps[5]] <- ")"
if(qps[2] - qps[1] > 1)
bar.vec[(qps[1] + 1):(qps[2] - 1)] <- "."
if(qps[4] - qps[2] > 1)
bar.vec[(qps[2] + 1):(qps[4] - 1)] <- "-"
if(qps[5] - qps[4] > 1) bar.vec[(qps[4] + 1):(qps[5] - 1)] <- "." #
### Add in outliers if present
if(length(x.outliers)>0)
bar.vec[out.pos] <- "*"
bar.tex <- paste(bar.vec, collapse = "")
full.bar.tex <- paste(":", bar.tex, ":", sep = "")
full.bar.tex
}
, comment = "18/06/2001")
"text.lot" <-
structure(function(outfile = "write.fn.src")
{
# creates source file to list all functions
src.comm <- ppaste("source('", outfile, "')") # to use src file
oblist <- ls(pos = 1)
new.function <- 0 # no of functions to write
app <- FALSE
today <- system("date +%d/%m/%Y", TRUE)
for(i in oblist) {
date.fn <- comment(get(i))
if(is.function(get(i)) & (date.fn == today | date.fn == "") ) {
source.line <- ppaste("write.function(", i, ")")
write(source.line, file = outfile, append = app)
new.function <- new.function + 1
app <- TRUE
}
}
cat(ppaste("Now you can write ", new.function, " function",
ifelse(new.function > 1, "s", "")," using:\n", src.comm,
"\n"))
}
, comment = "04/02/2003")
"tidy" <-
structure(function()
{
### Purpose: Source file that adds dates to undated functions and keeps code.
### ----------------------------------------------------------------------
### Modified from: incorporates catalogue and few others
### ----------------------------------------------------------------------
### Arguments:
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 17 Oct 2002, 11:06
source("/home/hrapgc/Rstuff/Gems/tidyup")
}
, comment = "17/10/2002")
"time.to.mort" <-
structure(function(ab, link = "cloglog", pc = 0.99996830000000003)
{
### For calculating time to (by default) Probit9 mortality, or any other %
#
### Note: if the x-axis has used a transformation before fitting the line,
### the reverse transformation must be made to this output.
#
f <- link.function(link)
targ <- f(ab$cm + (1 - ab$cm) * pc)
lt <- (targ - ab$intercept)/ab$slope
lt
}
, comment = "25/04/1997")
"tuk.calc" <-
structure(function(comp = 3, error.df = 50, se = 0.2029, interval = 0.94999999999999996,
sed.logical = T)
{
if(sed.logical) {
qtukey(interval, comp, error.df) * (se/sqrt(2))
}
else {
qtukey(interval, comp, error.df) * se
}
}
, comment = "24/12/1996")
"tuk.groups" <-
structure(function(xx, y, sing.top.ok = TRUE, print.mat = FALSE, conf = 0.95,
full.names = FALSE, print.grid = FALSE, print.sim.df = FALSE)
{
### Purpose: Get Tukey HSD groupings from an aov
### ----------------------------------------------------------------------
### Modified from: Various -- terminology is from palatability of plants
### ----------------------------------------------------------------------
### Arguments: sing.top.ok Sometimes need to fudge last value (Why??)
### not entirely trustable -- use with caution
### can't handle factor levels containing a '-' character
### xx: an anova object
### y: which factor are we interested in?
### sing.top.ok: Is it OK to have a singular top value? (fudge?)
### full.names: long names will maker colnames too long
### print.mat: Do we want to see the matrix with astrixes?
### print.grid: Do we want the table printed with grids?
### print.sim.df: Print the dataframe with levels that
### do not have an HSD with some other?
### ----------------------------------------------------------------------
### Author: Patrick Connolly, Creation date: 20 Aug 2003, 14:07
xx.tuk <- TukeyHSD(xx, ordered = T, conf.level = conf)
tuk.df <- as.data.frame(xx.tuk[[y]])
tuk.df$HSD[tuk.df$lwr < 0 ] <- "*"
## Get order of most palatable host
palatable <- names(sort(model.tables(xx, type = "means")$tables[[y]]))
plants <- rownames(tuk.df)
plants.list <- strsplit(plants, "-")
tuk.df$Plant1 <- unlist(lapply(plants.list, "[", 2))
tuk.df$Plant2 <- unlist(lapply(plants.list, "[", 1))
same.groups <- list()
tuk.same <- tuk.df[tuk.df$lwr < 0 , ]
tuk.diff <- tuk.df[tuk.df$lwr >= 0 , ]
if(print.sim.df)
print(tuk.same)
if(print.mat){ ## Print similarity matrix
diff.mat <- matrix("", nrow = length(palatable), ncol = length(palatable))
dimnames(diff.mat) <- list(palatable, palatable)
for(i in palatable){
for(j in palatable){
if(i != j ){
first <- match(i, tuk.df$Plant1)
second <- match(j, tuk.df$Plant2)
whichrow <- paste(j, i, sep = "-")
diff.mat[i, j] <- tuk.df[whichrow, ]$HSD
}
}
}
diag(diff.mat) <- "∑"
diff.mat[is.na(diff.mat)] <- ""
diff.df <- as.data.frame(t(diff.mat))
if(!full.names)
names(diff.df) <- LETTERS[seq(palatable)]
if(print.grid)
print.char.matrix(as.data.frame(diff.df), col.names = T)
else print(as.data.frame(diff.df))
cat("\n")
}
for(i in palatable){
same.i1 <- tuk.same[tuk.same$Plant1 == i, "Plant2"]
same.i2 <- tuk.same[tuk.same$Plant2 == i, "Plant1"]
same.groups[[i]]$max <- same.groups[[i]]$min <- match(i, palatable)
if(length(same.i1) > 0)
same.groups[[i]]$max <- match(rev(same.i1)[1], palatable)
if(length(same.i2) > 0)
same.groups[[i]]$min <- match(same.i2, palatable)
}
### Find corners in similarity matrix
mins <- sapply(same.groups, function(x) x$min[1])
maxs <- sapply(same.groups, function(x) x$max[1])
low <- rep(seq(unique(mins)), table(mins))
#
if(!sing.top.ok){# might be inaccurate to stop singular top value
if(rev(low)[1] != rev(low)[2])
low[length(low)] <- low[length(low) - 1]
}
high <- rep(seq(unique(maxs)), table(maxs)) # same length as palatable
if(length(low) != length(high))
return()
group <- list()
for(i in seq(palatable))
group[[palatable[i]]] <- paste(sort(unique(letters[seq(low[i], high[i])])),
collapse = "")
### return a named vector
groups <- sapply(group, function(x) x)
groups
}
, comment = "18/03/2005")
"unfactor" <-
structure(function(x, cols = "every.factor.column", notify = F)
{
#Converts factors back to numerics or character vectors
#
# A dataframe (or vector) of the same dimensions and names is returned
# Default is to change all factor columns
# Alternatively, specific columns can be specified
#
# notify: Notify that no factors need changing?
# Vectors have to be handled differently:
onecol <- is.null(dim(x))
if(is.null(dim(x))) {
x <- as.data.frame(x)
cols <- 1
}
# Handling data frames:
if(is.numeric(cols))
change.cols <- cols
else {
if(cols == "every.factor.column") {
fact.cols <- check.fact(x)
change.cols <- find.pos(T, fact.cols) #
# (a numerical vector of column numbers that are to be changed)
}
else change.cols <- find.pos(cols, names(x))
}
if(length(change.cols) == 0 & notify)
cat(paste("No factors to change in", as.character(substitute(x)), "\n"))
else {
options(warn = -1)
on.exit(options(warn = 0))
for(i in change.cols) {
levels.i <- levels(x[[i]])
level.nos <- as.numeric(levels.i)
numbers <- ifelse(any(is.na(level.nos)), F, T)
if(numbers)
x[, i] <- as.numeric(I(as.character(x[, i])))
else x[, i] <- I(as.character(x[, i]))
}
}
if(onecol)
x <- as.vector(unlist(x))
x
}
, comment = "30/04/2001")
"wipe" <-
structure(function (x)
{
## Purpose: Removes the file in /tmp/ and sets source attribe to NULL
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Zaphod Beeblebrox, Date: 2 Mar 2002, 17:09
ob.name <- deparse(substitute(x))
attr.com <- ppaste("attr(", ob.name, ", \"source\") <- NULL\n")
tmp.com <- ppaste("rm /tmp/hrapgc.", ob.name, ".R\n")
try(system(tmp.com))
write.com <- ppaste("write.function(", ob.name, ")")
cat(" If you want to keep the code, now run these commands:\n")
cat(paste(write.com, attr.com, sep = " ;"))
}
, comment = "28/06/2002")
"write.function" <-
structure(function(x)
{
### writes a function to file in the ./scr/ directory
fun.name <- as.character(substitute(x))
filedate <- comment(get(fun.name))
location <- system("pwd", TRUE)
file.name <- ppaste(location,"/.src/R.", fun.name)
write(paste(" Listing of:", fun.name), ncolumns = 1, file = file.name)
write(paste(" Located in:", location), file = file.name, ncolumns = 1, append = TRUE)
write(paste("Last updated:", filedate,
"\n**************************************\n"),
file = file.name, append = TRUE)
dump(fun.name, file = file.name, append = TRUE)
}
, comment = "22/08/2002")
