# different types of plots for data mining
prototypes <- (x,y=,fun=) {
# basically a cleaned-up version of "aggregate"
# that is specialized for one "by" factor
(!(y) && (y) == "function") {
fun <- y
y <-
}
((y)) {
resp <- (x)
y <- x[resp]]
(!(y)) ("no class variable")
x <- not(x,resp)
} {
resp <- "Group.1"
}
y = (y)
((fun,)) {
# fast implementation of special case
(indicators.factor(y) %*% (x))
} ((fun,)) {
# fast implementation of special case
f <- indicators.factor(y)
f <- scale.rows(f, 1/(f))
(f %*% (x))
} (T) {
# faster than below
xs <- ((x),y)
xp <- (xs,(xl) (xl,2,fun))
((xp))
} {
yl <- (y)
(yl) <- resp
x <- (x,yl,fun)
(x) <- (y)
x[setdiff((x),resp)]
}
}
hclust.prototypes <- (x,) {
resp <- (x)
px <- not(prototypes(x),resp)
((px)^2,members=(x[resp]]),)
}
reorder.svd.data.frame <- (x,dims=(1,2)) {
("using svd to order\n")
i <- (x,)
nx <- x[i]
nx = na.dummy(nx,0)
s <- (nx - rep.row((nx),(nx)))
v1 = s$u,1]
s <- (nx - rep.col((nx),(nx)))
v2 = s$v,1]
# remove flip ambiguity by alphabetic sort
((x)[which.min(v1)] > (x)[which.max(v1)]) {
v1 <- -v1
}
((x)[which.min(v2)] > (x)[which.max(v2)]) {
v2 <- -v2
}
ord = ((v1),(v2))
(d 1:((x))) {
(d %in% dims) {
} {
ord[d]] = 1:(x)[d]
}
}
((!i)) {
j = 1:(i)
j[i] = (i)[ord[2]]]
ord[2]] = j
}
ord
}
reorder.hc.data.frame <- (x,dims=(1,2)) {
("using hclust to order\n")
ord = ()
(1 %in% dims) {
# rows
hc = ((x)^2,method="ward")
ord[1]] = hc$
} {
ord[1]] = 1:(x)
}
(2 %in% dims) {
# columns
hc = (((x))^2,method="ward")
ord[2]] = hc$
} {
ord[2]] = 1:(x)
}
ord
}
apply_order <- (x, ord) {
x <- x[ord[1]],]
x <- x,ord[2]]]
x
}
#' Data Image
#'
#' Each value in a data matrix is represented by a colored pixel.
#'
#' @param x a matrix or data frame.
#' @param reorder If \code{TRUE}, the rows and columns are reordered so that
#' the data along each row and each column follows a linear trend. Can also be
#' a vector of two logicals, indicating separately whether to reorder rows
#' and/or columns.
#' @param scale If \code{TRUE}, each column is scaled to have minimum 0 and
#' maximum 1.
#' @param col a vector of colors to use for representing values.
#' @param ... extra arguments to \code{\link{image.default}}.
#' @author Tom Minka
#' @export
#' @seealso \code{\link{YR.colors}}
#' @references M. C. Minnotte and R. W. West. "The Data Image: A Tool For
#' Exploring High Dimensional Data Sets." Proceedings of the ASA Section on
#' Statistical Graphics, 1998. \url{http://citeseer.nj.nec.com/72085.html}
#' \url{http://math.usu.edu/~minnotte/research/pubs.html}
#' @examples
#'
#' data(USJudgeRatings)
#' data_image(USJudgeRatings,col=RC.colors(32))
#'
<- (x,las=1,=T,=T,) {
#mar=c(2.5,8.5,0,0.1)
x = (x)
x <- x[sapply(x,)]
() x <- apply.df(x,scale.range)
( %% 2 == 0) {
x = apply.df(x,)
}
( > 2) {
x = (x)
s = linear.profiles(x)$s
s = (s)
x = scale.cols(x,s)
}
(() == 1) = (,2)
(()) {
ord <- reorder.hc.data.frame(x,())
x <- apply_order(x, ord)
}
image.table((x),las=las,)
(ord)
}
#' Reordered and scaled star plot
#'
#' Automatically reorders and scales the variables to make a
#' readable star plot.
#'
#' @param x a matrix or data frame. Each row makes one star.
#' @param proto If \code{TRUE} and the response variable is a factor,
#' then only the mean value for each response is plotted.
#' See \code{\link{prototypes}}.
#' @param reorder If \code{TRUE}, the rows and columns are reordered so
#' that the data along each row and each column follows a linear trend.
#' Can also be a vector of two logicals, indicating separately
#' whether to reorder rows and/or columns.
#' @param scale If \code{TRUE}, each column is scaled to have minimum
#' 0 and maximum 1.
#' @param... additional arguments for \code{\link{stars}}.
#' @author{Tom Minka}
#' @seealso \code{\link{stars}}
#' @examples
#' data(iris)
#' star_plot(iris)
#' @export
<- (x,draw.segments=T,proto=T,=T,=T,) {
# reorder = c(T,F) will only reorder cases
x = (x)
resp <- (x)
((x[resp]]) && proto) {
("showing median of each",resp,"\n")
f <- x[resp]
x <- x[sapply(x,)]
x <- (x,f,)
(x) <- (f[1]])
}
x <- x[sapply(x,)]
() {
x <- apply.df(x,scale.range)
x = apply.df(x,)
}
(() == 1) = (,2)
(()) {
ord <- reorder.svd.data.frame(x,())
x <- apply_order(x, ord)
}
#opar <- par(mar=c(0.05,0,0,0.1))
#on.exit(par(opar))
stars(x, draw.segments=draw.segments, =F, )
}
parallel.cases <- (x,yscale=("linear","range","none"),
xscale=("equal","linear","none"),
yaxt=(yscale=="none") "s" "n",) {
x <- (x)
y <- x
# assume no factors
y <- (y)
x <- 1:(y)
yscale <- (yscale)
xscale <- (xscale)
(yscale=="linear") {
# do this to subtract means
# scaling y here has no effect
y <- (y)
(xscale != "none") {
# want yij*sj = aj + bj*xi
# aj subtracts mean of yij
# then xi is top eigenvector
s = (y)
x = s$u,1]
s = (1/s$v,1])
(xscale == "equal") {
x <- rank.stable(x)
# now act as if x was fixed and solve for scaling again
xscale <- "none"
}
xscale <- "none"
}
(xscale=="none") {
# solve for sj when xj is known
xt <- x - (x)
(T) {
# scale to standardize residuals
r = (y*y) - (xt %*% y)^2/(xt*xt)
r = r/(y)
s = 1/(r)
} {
# scale to match svd solution
r = (xt %*% y)/(xt*xt)
s = 1/r
}
}
# never any reversals
y <- y * rep.row(s,(y))
} {
(yscale == "range") {
# map to [0,1]
y <- y - rep.row((y,2,),(y))
y <- y / rep.row((y,2,),(y))
}
(xscale != "none") {
# want yij = aj + bj*xi
# subtract mean in each col to remove aj,
# then xi is top eigenvector
s <- (y)
x <- s$u,1]
(xscale == "equal") x <- rank.stable(x)
}
}
# remove x flip ambiguity by alphabetic sort
((y)[which.min(x)] > (y)[which.max(x)]) {
x <- -x
}
i <- (x)
y <- reorder.rows(y,i)
x <- x[i]
(x,y,xlab="",ylab="",yaxt=yaxt,)
i
}
#' Parallel-coordinate plot
#'
#' @description Plots a data matrix using parallel coordinates.
#' @param x a data matrix or frame. Only numeric columns will be
#' plotted.
#' @param yscale describes how the columns should be scaled.
#' @param xscale describes how the columns should be placed.
#' @param proto If \code{TRUE} and the response variable is a factor,
#' then only the mean value for each response is plotted.
#' See \code{\link{prototypes}}.
#' @param ...Extra arguments to \code{\link{labeled.curves}}.
#' @details
#' Let the rows of the matrix be \emph{cases} and the columns be
#' \emph{dimensions}.
#' Each dimension is laid out as a vertical axis, and the value for each
#' case is plotted as a point on that axis. Points corresponding to the
#' same case are connected by lines. Thus each case is represented by a
#' curve or \emph{profile}, moving horizontally from dimension to dimension.
#'
#' To enhance readability, the dimensions should usually be shifted to have
#' a common center and scaled to have
#' comparable units. This is controlled by \code{yscale}.
#' If \code{yscale="linear"}, the dimensions are automatically shifted
#' and scaled
#' by the linear profiles method. That is, the profiles are
#' made to be as straight as possible.
#' This may result in negative scaling for a dimension,
#' i.e. the axis is reversed so that moving up means a lower value.
#' When this happens, the dimension name is prepended with a minus sign.
#' If \code{yscale="range"}, the dimensions are shifted and scaled (by a
#' positive number) so that the data ranges from 0 to 1 on each axis.
#' If \code{yscale="none"}, the dimensions are plotted as-is.
#'
#' Another important choice is the ordering and spacing of the axes
#' on the plot.
#' If \code{xscale="linear"}, the axes are placed according to the
#' linear profiles method.
#' If \code{xscale="equal"}, the axes are placed similarly to
#' \code{"linear"} but with the constraint that they must be equally
#' spaced.
#' If \code{xscale="none"}, the axes are placed in the order
#' that they appear in the matrix.
#'
#' If the data frame has a categorical response, the profiles are colored
#' according to the response. The color scheme can be adjusted using the
#' arguments to \code{\link{labeled.curves}}.
#' @author Tom Minka
#' @references A. Inselberg and B. Dimsdale.
#' "Parallel coordinates: A tool for visualizing multidimensional geometry."
#' Proc. of Visualization 90, p. 361-78, 1990.
#'
#' E. J. Wegman.
#' "Hyperdimensional Data Analysis Using Parallel Coordinates."
#' JASA 85:664-675, 1990.
#' @seealso \code{\link{parallel.cases}},\code{\link{star.plot}}
#' @examples
#' data(iris)
#' parallel_plot(iris)
#' parallel_plot(iris,proto=F,labels=NULL)
#' @export
<- (x,yscale=("linear","range","none"),flipy=F,
xscale=("equal","linear","none"),proto=T,flipx=F,
yaxt=(yscale=="none") "s" "n",=1,) {
x <- (x)
resp <- (x)
((x[resp]])) {
(proto) {
("showing median of each",resp,"\n")
f <- x[resp]
x <- x[sapply(x,)]
x <- (x,f,)
(x) <- ((f[1]]))
<- 1:(x)
} {
# plot cases, but color by class
<- ((x[resp]]))
}
x <- not(x,resp)
} {
<- 1:(x)
}
y <- x
# assume no factors
y <- (y)
x <- 1:(y)
yscale <- (yscale)
xscale <- (xscale)
(yscale=="linear") {
# do this to subtract means
# scaling y here has no effect
y <- (y)
# bugfix for scale: NaNs are created when a column has zero variance
y[is.nan(y)] = 0
(xscale != "none") {
r = linear.profiles(y,=); x = r$x; s = r$s
(xscale == "equal") {
x <- rank.stable(x)
# now act as if x was fixed and solve for scaling again
xscale <- "none"
}
}
(xscale=="none") {
s = linear.profiles(y,x=x,=)$s
}
# minimize number of reversals
((s<0) > (s>0)) s <- -s
(flipy) s = -s
y <- y * rep.row(s,(y))
reversed <- (s<0)
} {
(yscale == "range") {
# map to [0,1]
y <- y - rep.row((y,2,),(y))
y <- y / rep.row((y,2,)+eps,(y))
}
(xscale != "none") {
x = linear.profiles(y,s=1,=)$x
(((x))) ("internal: x is nan")
(xscale == "equal") x <- rank.stable(x)
}
reversed = F
}
# remove x flip ambiguity by alphabetic sort
((y)[which.min(x)] > (y)[which.max(x)]) {
x <- -x
}
(flipx) x = -x
# must do this after the sort
((reversed)) {
("axis reversed for",(y)[reversed],"\n")
(y)[reversed] <- ((y)[reversed], (s) ("-",s,sep=""))
}
i <- (x)
y <- reorder.cols(y,i)
x <- x[i]
(x,(y),=,xlab="",ylab="",yaxt=yaxt,)
}
improve.s <- (A,s,n) {
d = (s)
((A) != d || (A) != d) error("A is wrong size")
(j 1:d) {
b = (A[j,-j] %*% s-j])/A[j,j]
s1 = 0.5*(-b + (b^2 + 4*n/A[j,j]))
s2 = 0.5*(-b - (b^2 + 4*n/A[j,j]))
j1 = A[j,j]*(s1 + b)^2 - n*(s1^2)
j2 = A[j,j]*(s2 + b)^2 - n*(s2^2)
((j1,j2))
s[j] = (j1 < j2) s1 s2
}
s
}
linear.profiles2 <- (y,x=,s=) {
y = (y)
n = (y)
d <- (y)
r = linear.profiles(y,x,s); x = r$x; s = r$s
s = (d)
#print(s)
#print(x)
<- (d) - (1/d,(d,d))
(iter 1:10) {
old.x = x
(iter2 1:100) {
old.s = s
xt <- x - (x)
P = - (xt,xt)/(xt*xt)
# Hadamard product
A = P * ((y) %*% y)
# function value
J = (s) %*% A %*% s - n*((s^2))
(J[1])
s = improve.s(A,s,(y))
activity = ((s - old.s))
#print(activity)
(activity < 1e-2)
}
((x)) {
# improve x
x = (scale.cols(y,s) %*% )$v,1]
#print(x)
}
activity = ((x - old.x))
(activity < 1e-4)
}
s = s * ((1/s^2))
(iter,"iters\n")
(x=x,s=s)
}
linear.profiles <- (y,x=,s=,=1) {
y = (y)
#y = scale(y)
((x) && (s)) {
# Hartigan's linear profile algorithm
# want yij*sj = ai + bi*xj
# yij = ai*(1/sj) + bi*(xj/sj)
# so 1/sj = v[j,1], xj/sj = v[j,2]
# the constraint is sum_j 1/sj^2 = 1
# scaling by singular values is irrelevant
sv <- (y)
x <- sv$v,2]/sv$v,1]
s <- 1/sv$v,1]
} (==1 && (s)) {
# solve for sj when xj is known
sv = (y)
s = 1/sv$v,1]
} ((s)) {
# solve for sj when xj is known
# this comes from solving for (ai,bi) analytically and substituting
# to get eigenvector problem for sj
# the constraint is sum_j sj^2 = 1
d <- (y)
xt <- x - (x)
a <- (d) - (1/d,(d,d)) - (xt,xt)/(xt*xt)
# Hadamard product
a <- a * ((y) %*% y)
s <- (a)
# want the smallest eigenvector of a
s <- s$u,d]
} ((x)) {
# want yij = ai + bi*xj
# subtract mean in each row to remove ai,
# then xj is top eigenvector
sv <- (y - rep.col((y),(y)))
x <- sv$v,1]
}
(x=x,s=s)
}
# abbreviate strings so that max width is w, and strings remain unique
abbreviate.width <- (s,w,cex=("cex")) {
old.max <- w
{
wid <- (s,cex=cex)
((wid)<=w)
i <- wid > w
j <- (wid)
# this may cause a collision
s[i] <- (s[i],(s[j])-1)
}
s
}
# minka: use auto.mar and auto.layout
stars <-
(x, full = , = , radius = ,
= (x)[1]], locations = ,
= , len = 1,
key.loc = , key.labels = (x)[2]], key.xpd = ,
= , = , flip.labels = ,
draw.segments = , col.segments = 1:n.seg,
col.stars = ,
axes = , frame.plot = axes,
main = , = , xlab = "", ylab = "",
cex = 0.8, lwd = 0.25, lty = ("lty"), xpd = ,
mar = (("mar"),auto.mar(main,xlab,ylab,axes)),
)
{
# technically, strheight doesn't work until we've made axes
#char.height <- max(strheight(labels,cex=cex))
char.height <- cex*("cxy")[2]
cell.size <- () {
((flip.labels)) flip.labels = T
ff <- (2.1,2.1)
# minka: adjust loc for labels
(!()) {
ff[2] <- ff[2] + char.height*(1+flip.labels)
}
ff
}
if ((x))
x <- (x)
if (!(x))
("x must be a matrix or a data frame")
if (!(x))
("data in x must be numeric")
n.loc <- (x)
n.seg <- (x)
if ((locations)) { ## Default (x,y) locations matrix
# minka: use auto.layout
(()) {
ff = cell.size()
= ff[2]/ff[1]
= auto.layout(n.loc,)
}
(!() && (flip.labels))
flip.labels <- [2] * (()) > 30
ff = cell.size()
locations <- expand.grid(ff[1] * 1:[2],
ff[2] * [1]:1)[1:n.loc, ]
}
{
if ((locations) && (locations) == 2) {
## all stars around the same origin
locations <- ((locations[1],n.loc),
(locations[2],n.loc))
(!() && n.loc > 1)
("labels don't make sense for a single location")
<-
}
{
if ((locations))
locations <- (locations)
if (!(locations) || (locations) != 2)
("locations must be a 2-column matrix.")
if (n.loc != (locations))
("number of rows of locations and x must be equal.")
}
((flip.labels))
flip.labels <- # have no grid
}
xloc <- locations,1]
yloc <- locations,2]
## Angles start at zero and pace around the circle counter
## clock-wise in equal increments.
angles <-
(full)
(0, 2*, =n.seg+1)-(n.seg+1)]
if (draw.segments)
(0, , =n.seg+1)-(n.seg+1)]
(0, , =n.seg)
if ((angles) != n.seg)
("length(angles) must be the same as ncol(x)")
## Missing values are treated as 0
if () {
x <- (x, 2, (x)
(x - (x, na.rm = ))/((x, na.rm = )))
# minka: sqrt because we are depicting area
x = (x)
} {
# just scale so max is 1, as assumed by rest of code
x = x/(x)
}
## Missing values are treated as 0
x[is.na(x)] <- 0
mx <- (x <- x * len)
# minka: place key automatically at top
# do this now for computing ylim
(!(key.loc) && (key.loc)) {
key.loc <- ((xloc), (yloc) + 2.1 + char.height)
}
(()) {
<- (xloc) + (-mx, mx)
# minka: should make room for labels too - use range.text
(!(key.loc)) {
# minka: make room for the key
w <- ((key.labels,cex=cex,="inches"))/("pin")[1]
xlim.key = extend.inches((key.loc[1]-len,key.loc[1]+len),w,w)
[1] <- ([1],xlim.key[1])
[2] <- ([2],xlim.key[2])
}
}
(()) {
<- (yloc) + (-mx, mx)
(!()) {
# minka: make room for the bottom row labels
[1] <- (yloc - len - char.height)
}
(!(key.loc)) {
# minka: make room for the key
[1] <- ([1],key.loc[2]-len-char.height)
[2] <- ([2],key.loc[2]+len+char.height)
}
}
deg <- / 180
## The asp argument keeps everything (the symbols!) square
op <- (mar = mar, xpd = xpd) ; ((op))
(0, ="n", , =, =,
main = main, = , xlab = xlab, ylab=ylab, = 1, axes = axes)
s.x <- xloc + x * ((angles), (n.loc,n.seg))
s.y <- yloc + x * ((angles), (n.loc,n.seg))
if ( draw.segments ) {
aangl <- (angles, (full)2* )
for (i 1:n.loc) { ## for each location, draw a segment diagram
px <- py <- ()
for (j 1:n.seg) {
k <- (from = aangl[j], to = aangl[j+1], = 1*deg)
px <- (px, xloc[i], s.x[i,j], x[i,j]*(k) + xloc[i], )
py <- (py, yloc[i], s.y[i,j], x[i,j]*(k) + yloc[i], )
}
(px, py, = col.segments, lwd=lwd, lty=lty)
}
} # Segment diagrams
{ # Draw stars instead
for (i 1:n.loc) {
(s.x[i,], s.y[i,], lwd=lwd, lty=lty, = col.stars[i])
if (radius)
((xloc[i],n.seg),
(yloc[i],n.seg),
s.x[i,], s.y[i,], lwd=lwd, lty=lty)
}
}
(!()) {
## vertical text offset from center
y.off <- mx * ((full) 1 0.1)
(flip.labels)
# minka: bug fix
y.off <- y.off + cex*("cxy")[2] * (((n.loc)-1)%%2)
##DBG cat("mx=",format(mx),"y.off:"); str(y.off)
(xloc, yloc - y.off, , cex=cex, adj=(0.5, 1))
}
if ( !(key.loc) ) { ## Draw unit key
## usually allow drawing outside plot region:
(xpd = key.xpd) # had `xpd' already above
key.x <- len * (angles) + key.loc[1]
key.y <- len * (angles) + key.loc[2]
if (draw.segments) {
px <- py <- ()
for (j 1:n.seg) {
k <- (from = aangl[j], to = aangl[j+1], = 1*deg)
px <- (px, key.loc[1], key.x[j], len * (k) + key.loc[1], )
py <- (py, key.loc[2], key.y[j], len * (k) + key.loc[2], )
}
(px, py, = col.segments, lwd=lwd, lty=lty)
}
{ # draw unit star
(key.x, key.y, lwd=lwd, lty=lty)
if (radius)
((key.loc[1],n.seg), (key.loc[2],n.seg),
key.x, key.y, lwd=lwd, lty=lty)
}
## Radial Labeling -- should this be a standalone function ?
lab.angl <- angles +
(draw.segments) (angles[2] - angles[1]) / 2 0
label.x <- 1.1 * len * (lab.angl) + key.loc[1]
label.y <- 1.1 * len * (lab.angl) + key.loc[2]
## Maybe do the following without loop {need not use adj but ..)!
for (k 1:n.seg) {
text.adj <-
(## horizontal
if (lab.angl[k] < 90*deg || lab.angl[k] > 270*deg) 0
if (lab.angl[k] > 90*deg && lab.angl[k] < 270*deg) 1
0.5,
## vertical
if (lab.angl[k] <= 90*deg) (1 - lab.angl[k] / (90*deg)) /2
if (lab.angl[k] <= 270*deg)
(lab.angl[k] - 90*deg) / (180*deg)
## lab.angl[k] > 270*deg
1 - (lab.angl[k] - 270*deg) / (180*deg)
)
(label.x[k], label.y[k],
= key.labels[k], cex = cex, adj = text.adj)
}
} # Unit key is drawn and labelled
if (frame.plot) ()
(locations)
}
# modified by Tom Minka
# abbreviation, better spacing, better handling of key
my.stars <-
(x, full = , = , radius = ,
= (x)[1]],
locations = , = , = , len = 1,
= , abbrev=T,
key.loc = , key.labels = , key.xpd = ,
draw.segments = , axes = ,
cex = 0.8, lwd = 0.25, )
{
if ((x))
x <- (x)
if (!(x))
("x must be a matrix or a data frame")
if (!(x))
("data in x must be numeric")
n.loc <- (x)
n.seg <- (x)
deg <- / 180 # segments only
seg.colors <- (!()) 1:n.seg
if ((locations)) { # make loc matrix
md <- ((n.loc)) # => md^2 >= n.loc
loc.x <- 2.1*((1:md)-1)
loc.y <- (loc.x)
# minka: adjust loc for labels
(!()) {
# technically, strheight doesn't work until we've made axes
spc <- ((,cex=cex))
loc.y <- ((2.1+spc)*((1:md)-1))
}
loc <- expand.grid(loc.x, loc.y)[1:n.loc, ]
}
{
if ((locations) && (locations) == 2) {
## all stars around the same origin
loc <- ((locations[1],n.loc),
(locations[2],n.loc))
(!() && n.loc > 1)
("labels don't make sense for a single location")
<-
}
{
if ((locations))
locations <- (locations)
if (!(locations) || (locations) != 2)
("locations must be a 2-column matrix.")
loc <- (locations)
if (n.loc != (loc))
("number of rows of locations and x must be equal.")
}
}
## Angles start at zero and pace around the circle counter
## clock-wise in equal increments.
angles <-
(full)
(0, 2*, =n.seg+1)-(n.seg+1)]
if (draw.segments)
(0, , =n.seg+1)-(n.seg+1)]
(0, , =n.seg)
if ((angles) != n.seg)
("length(angles) must be the same as ncol(x)")
## Missing values are treated as 0
x[is.na(x)] <- 0
if () {
# minka: use full range of [0.1,1]
x <- (x,2,(x,2,),FUN="-")
x <- (x,2,(x,2,)/0.9,FUN="/")
x <- (x,2,0.1,FUN="+")
## Columns of 0s will put NAs in x, next line gets rid of them
x[is.na(x)] <- 0
# minka: sqrt because we are depicting area
x = (x)
}
x <- x * len
# minka: place key automatically at top
(!(key.loc) && (key.loc)) {
key.loc <- ((loc[1]]), (loc[2]]) + 2.5)
}
((key.labels)) key.labels <- (x)[2]]
(()) {
<- (loc,1] + len, loc,1] - len)
(!(key.loc)) {
w <- ((key.labels,cex=cex))
[1] <- ([1],key.loc[1]-1-w)
[2] <- ([2],key.loc[1]+1+w)
}
}
(()) {
<- (loc,2] + len, loc,2] - len)
(!()) {
spc <- ((,cex=cex))
[1] <- (loc,2] - 1 - spc)
}
(!(key.loc)) {
h <- ((key.labels,cex=cex))
[1] <- ([1],key.loc[2]-1-h)
[2] <- ([2],key.loc[2]+1+h)
}
}
## The asp argument keeps everything square
(0, ="n", , =, =,
xlab="", ylab="", = 1, axes = axes)
# minka: abbreviate labels if necessary
(!() && abbrev) {
<- abbreviate.width(,2,cex=cex)
}
if ( draw.segments ) {
## for each location, draw a segment diagram
for ( i 1:n.loc ) {
poly.x <-
poly.y <-
start.x <- x[i,] * ( angles ) + loc[i,1]
start.y <- x[i,] * ( angles ) + loc[i,2]
### FIXME : we can do without the following inner loop !
for (j 1:n.seg) {
poly.x <- (poly.x,loc[i,1],start.x[j])
poly.y <- (poly.y,loc[i,2],start.y[j])
next.angle <-
(j < n.seg)
angles[j+1]
else ((full) 360 180) * deg
k <- (from = angles[j], to = next.angle, = deg)
poly.x <- (poly.x, x[i,j] * ( k ) + loc[i,1], )
poly.y <- (poly.y, x[i,j] * ( k ) + loc[i,2], )
}
(poly.x, poly.y, = seg.colors, lwd=lwd)
if (!())
(loc[i,1], loc[i,2] - (full)(x) 0.1 * (x),
[i], cex=cex, adj=(0.5,1))
}
} # Segment diagrams
{ # Draw stars instead
for ( i 1:n.loc ) {
temp.x <- x[i,] * ( angles ) + loc[i,1]
temp.y <- x[i,] * ( angles ) + loc[i,2]
if (radius)
((loc[i,1],n.seg),
(loc[i,2],n.seg),
temp.x, temp.y, lwd=lwd)
((temp.x, temp.x[1]),
(temp.y, temp.y[1]), lwd=lwd)
if (!())
(loc[i,1], loc[i,2] - (full)(x) 0.1 * (x),
[i], cex=cex, adj=(0.5,1))
}
}
if ( !(key.loc) ) { ## Draw unit key
## allow drawing outside plot region (inside figure region):
op <- (xpd = key.xpd); ((op))
key.r <- len * 0.8
key.x.coord <- ( angles ) * key.r + key.loc[1]
key.y.coord <- ( angles ) * key.r + key.loc[2]
if ( draw.segments ) {
key.x <-
key.y <-
for (j 1:n.seg){
key.x <- (key.x,key.loc[1],key.x.coord[j])
key.y <- (key.y,key.loc[2],key.y.coord[j])
k <- angles[j] + deg
next.angle <-
if (j < n.seg) angles[j+1]
else ((full) 360 180) * deg
while (k < next.angle) {
key.x <- (key.x, key.r * ( k ) + key.loc[1])
key.y <- (key.y, key.r * ( k ) + key.loc[2])
k <- k + deg
}
key.x <- (key.x, key.r * ( next.angle ) + key.loc[1], )
key.y <- (key.y, key.r * ( next.angle ) + key.loc[2], )
}
(key.x, key.y, = seg.colors, lwd=lwd)
}
{ # draw unit star
if ( radius )
((key.loc[1],n.seg), (key.loc[2],n.seg),
key.x.coord, key.y.coord, lwd=lwd)
((key.x.coord, key.x.coord[1]),
(key.y.coord, key.y.coord[1]), lwd=lwd)
}
lab.angl <- angles +
(draw.segments) (angles[2] - angles[1]) / 2 0
label.r <- len*0.9
label.x <- ( lab.angl ) * label.r + key.loc[1]
label.y <- ( lab.angl ) * label.r + key.loc[2]
##-- FIXME : Do the following witout loop !
for (k 1:n.seg) {
text.adj <-
(## horizontal
if (lab.angl[k] < 90*deg || lab.angl[k] > 270*deg) 0
if (lab.angl[k] > 90*deg && lab.angl[k] < 270*deg) 1
0.5,
## vertical
if (lab.angl[k] <= 90*deg) (1 - lab.angl[k] / (90*deg)) /2
if (lab.angl[k] <= 270*deg)
(lab.angl[k] - 90*deg) / (180*deg)
## lab.angl[k] > 270*deg
1 - (lab.angl[k] - 270*deg) / (180*deg)
)
(label.x[k], label.y[k],
= key.labels[k], cex = cex, adj = text.adj)
}
} # Unit key is drawn and labelled
()
}
qp <- (A,b,tol=1e-10) {
# returns the vector x which minimizes x'x subject to Ax=b.
#library(MASS)
# standardize the problem
s = ((A * A))
A = A/rep.col(s,(A))
b = b/s
d = (A)
x = (0,(d,1))
ws = ()
(iter 1:100) {
# find most violated constraint
s = (A %*% x) - b
# constraints in the working set must be satisfied
# this is to avoid problems with numerical inaccuracy
s[ws] = 0
i = (s)
(s[i] >= -tol)
(i %in% ws) {
(s[i])
("double violated constraint")
}
ws = (ws,i)
(T) {
# solve equality-constrained problem
Aw = A[ws,,=F]
bw = b[ws,=F]
(T) {
y = (Aw %*% (Aw), bw)
} {
y = ginv(Aw %*% (Aw)) %*% bw
}
x = (Aw) %*% y
((y<=eps)) {
# remove y<0 from working set
ws = ws[y>eps]
}
}
#cat("working set:",ws,"\n",fill=60)
}
x
}
# not used
collision = (x,y,w,h) {
# returns T iff the first box collides with any of other boxes
n = (x)
(i 2:n) {
(((x[1]-x[i]) < (w[1]+w[i])/2) &&
((y[1]-y[i]) < (h[1]+h[i])/2)) (T)
}
F
}
find.gap <- (=0,h0,y,h,symmetric=T) {
# returns the smallest position of a box of size h0 which will not collide
# with boxes of size h centered on y.
# sort y
ord = (y)
y = y[ord]
h = h[ord]
n = (y)
# compute all gaps
dy = (y)-(h[1:(n-1)]+h[2:n])/2
# dy must be positive since there are no collisions among y's
i = (dy > h0)
# upper end
y.new = (y[n]+(h[n]+h0)/2,)
((i) > 0) {
# this isn't quite optimal
# lower bound
y.new = (y.new, y[i]+(h[i]+h0)/2)
i.next = (i+1,n)
# upper bound
y.new = (y.new, y[i.next]-(h[i.next]+h0)/2)
}
(symmetric) {
# lower end
y.new = (y.new, y[1]-(h[1]+h0)/2)
}
(!symmetric) {
y.new = y.new[y.new >= ]
((y.new) == 0) ("empty")
}
y.new[which.min((y.new - ))]
}
move.collisions.vertical <- (x,w,h,=0,symmetric=T) {
# find the minimal vertical shifts so that boxes of width w and height h
# centered on x do not overlap.
n = (x)
((w) == 1) w = (w,n)
((h) == 1) h = (h,n)
(() == 1) = (,n)
y = (n)
= (h)/4
# loop in reverse order
y[n] = [n]
(i in (n-1):1) {
# find nearby boxes (those for which overlap is possible)
= (i+1):n
j = [abs(x[i]-x[is]) < (w[i]+w[is])/2]
# now only need to avoid vertical overlap
((j) == 0) y[i] = [i]
y[i] = find.gap([i],h[i],y[j],h[j],symmetric)
(y[i] < [i] && !symmetric) {
(y[i])
([i])
(h[i])
(y[j])
(h[j])
("")
}
}
y
}
move.collisions <- (x,w) {
# move x's the minimal amount so that boxes of width w centered on x
# do not overlap.
# units of x and w must match (preferably inches).
((x) == 1) (x)
(F) {
# jitter x
x = x + 2*eps*((x))
x = x + ((x))
}
ord = (x)
x = x[ord]
w = w[ord]
# construct constraint matrix
# constraints are (x[j]+dx[j])-(x[i]+dx[i]) > (w[j]+w[i])/2
# or dx[j]-dx[i] > (w[j]+w[i])/2 -(x[j]-x[i])
# for all j>i
n = (x)
nc = n*(n-1)/2
A = (0,(nc,n))
b = (0,(nc,1))
index = 1
(i 1:(n-1)) {
(j in (i+1):n) {
A[index,j] = 1
A[index,i] = -1
b[index] = (w[j]+w[i])/2 - (x[j]-x[i])
index = index + 1
}
}
(A) = (x)
x = x + (qp(A,b))
# restore original order
x[ord] = x
x
}
absdist <- (x,y=x) {
(rep.row(y,(x)) - rep.col(x,(y)))
}
clip <- (x) {
(x > 0)*x
}
line.of.sight <- (bbox,bbox.seg,dx,dy,) {
n = (dx)
e = (0,n)
(i 1:n) {
bbox.i = bbox[setdiff(1:n,i),]
s = bbox.i - rep.row(bbox.seg[i,],(bbox.i))
s,(1,3)] = -s,(1,3)]
m1 = (s,1,)
bbox.i = bbox.i(m1 > 0),,=F]
((bbox.i) > 0) {
#print(nrow(bbox.i))
dm = ((dy[i],0,-dx[i],0, 0,dy[i],-dx[i],0,
dy[i],0,0,-dx[i], 0,dy[i],0,-dx[i]),(4,4))
s = (bbox.i %*% dm) - [i]
m1 = (s,1,)
m2 = (s,1,)
e[i] = ((m1*m2 < 0)*(-m1,m2)^2)
}
}
(e)
}
move.collisions2 <- (x,y,w,h,cost=(1,(x)),los=F) {
# move (x,y) the minimal amount so that boxes of width w and height h
# centered on (x,y) do not overlap.
# units of (x,y) and (w,h) must match (preferably inches).
# cost is a vector specifying how immovable each box is.
# cost[i] = 0 means box i will not be moved.
# returns the new (x,y) values, as a matrix.
ok = !(x) & !(y) & !(w) & !(h)
x0=x;y0=y;w0=w;h0=h;cost0=cost
x=x[ok];y=y[ok];w=w[ok];h=h[ok];cost=cost[ok]
n = (x)
cost2 = (cost,cost)
n2 = 2*n
i = ((1:n2)[cost2 > 0],(1:n2)[cost2 == 0])
ir = 1:n2
ir[i] = 1:n2
nz = (cost2 == 0)
wij = (rep.row(w,n)+rep.col(w,n))/2
(wij) = 0
hij = (rep.row(h,n)+rep.col(h,n))/2
(hij) = 0
cost.fcn <- (dxy) {
# returns the "badness" of a given c(dx,dy) vector
dxy = (dxy,(0,nz))[ir]
dx = dxy[1:n]
dy = dxy[n+(1:n)]
x = x+dx
y = y+dy
# want total movement to be small
e = dx*dx+dy*dy
# moving beyond the original box is especially bad
e = e/100 + clip((dx) - w/2)^2 + clip((dy) - h/2)^2
e = (cost*e)
# total overlap
eo = clip(wij-absdist(x))*clip(hij-absdist(y))
eo = (eo^2)
= 0
(los) {
# want to have an unblocked "line of sight" to the original position
bbox = (x-w/2,x+w/2,y-h/2,y+h/2)
= x*dy - y*dx
bbox.seg = (x-dx,x,y-dy,y)
bbox.seg = ((bbox.seg,1],bbox.seg,2]),
(bbox.seg,1],bbox.seg,2]),
(bbox.seg,3],bbox.seg,4]),
(bbox.seg,3],bbox.seg,4]))
= line.of.sight(bbox,bbox.seg,dx,dy,)
}
#cat(e," ",eo," ",el,"\n")
e + 100*eo + 100*
}
dxy = (0,n2-nz)
(F) {
res = (dxy,cost.fcn,control=(maxit=10000))
dxy = res$
("optim code: ",res$convergence,"\n")
("iters: ",res$counts,"\n")
("fcn value: ",res$value,"\n")
}
(T) {
res = (cost.fcn,dxy)
dxy = res$estimate
("nlm code: ",res$code,"\n")
("iters: ",res$iterations,"\n")
("fcn value: ",res$minimum,"\n")
}
# see if it is really optimal
#dxy[1] = 0
#print(cost.fcn(dxy))
dxy = (dxy,(0,nz))[ir]
dx = dxy[1:n]
dy = dxy[n+(1:n)]
x = x+dx
y = y+dy
x0[ok]=x;y0[ok]=y;x=x0;y=y0
(x,y)
}
panel.hist <- (x, )
{
# helper function for pairs(), from examples
# usage: pairs(x,diag.panel=panel.hist)
# put histograms on the diagonal
usr <- ("usr"); ((usr))
(usr = (usr[1:2], 0, 1.5) )
h <- (x, = )
breaks <- h$breaks; nB <- (breaks)
y <- h$counts; y <- y/(y)
(breaks-nB], 0, breaks-1], y, )
}
# default color for partial residuals
()$col.res <- 1
# default lowess span
#formals(panel.smooth)$span <- 4/5
#formals(lowess)$f <- 3/4
# returns a model frame where the responses have been replaced by their
# residuals under the fit
residual.frame <- (object,) {
(()) {
x <- my.model.frame(object)
resp <- (x)
((object,"glm")) {
p <- object$fitted.value
x[resp]] <- (object$y - p)/p/(1-p)
}
x[resp]] <- (object)
}
{
x <-
resp <- (object)
((object,"glm") || (object,"tree")) {
x[resp]] = misclassResiduals(object,)
} {
# expand p to include NAs
p <- (object,x)[rownames(x)]
(p) <-
x[resp]] <- x[resp]] - p
#attr(x,"terms") <- terms(object)
}
x = my.model.frame(((resp,"~.")),x)
}
x
}
partial.residuals <- (object,omit,condensed=T) {
resp <- (object)
pred <- predictor.vars(object)
(condensed) {
# remove the predictor and refit the model
x <- ()
(i 1:(omit)) {
fit <- (object,(".~. -",omit[i]))
#cat(pred[i],coef(fit),"\n")
x[i]] <- (fit)
}
(x) <- omit
} {
# do not refit the model
x <- ((object,="partial"),check.names=F)[omit]
r = (object)
(i 1:(x)) {
(x[i]]) = (r)
}
}
x
}
partial.residual.frame <- (object,omit,) {
x = (object)
resp <- (x)
x[resp]] = partial.residuals(object,omit,)[1]]
pred = predictor.vars(((resp,"~",omit)))
x[c(pred,resp)]
}
expand.frame <- (,v) {
# also see expand.model.frame
((v) == 0) ()
vs = (v,":")
(i 1:(v)) {
x = [vs[i]]]
(((x,))) {
as.contrasts = (f) ((f)) (f==(f)[2]) f
x = apply.df(x,as.contrasts)
}
[v[i]]]] <- (x,1,)
}
}
#' Plot predictors versus response
#'
#' Makes a matrix of pairwise scatterplots with lowess-type trend lines.
#'
#' @author{Tom Minka}
#'
#' @aliases predict_plot predict_plot.data.frame predict_plot.formula
#' @param formula a formula specifying the response and predictor variables
#' @param data,x a data frame with at least two columns
#' @param partial a model from which to compute partial residuals (used
#' by \code{\link{predict_plot.lm}}).
#' @param mcol,mlwd If plotting partial residuals of an \code{lm},
#' the color and width of the model predictions.
#' @param layout a vector \code{c(rows,cols)} specifying the desired
#' layout of panels. Otherwise chosen automatically based on the size
#' of the plotting window.
#' @param highlight a logical vector specifying which predictors to highlight.
#' @param se If \code{TRUE}, show standard errors in linecharts.
#' @param scol,slwd color and width of trend lines.
#' @param span,degree,family parameters for the trend line (see \code{loess}).
#' @param rtype how a factor response should be handled when drawing a
#' trend line.
#' @param identify.pred A character vector of predictor names for which
#' to interactively \code{\link{identify}} points. If \code{TRUE},
#' done for all predictors.
#' @param mar margins within each panel
#' @param xaxt,yaxt arguments to \code{par}
#' @param col plotting color for symbols
#' @param asp Aspect ratio for each panel. If \code{"auto"}, the aspect
#' ratio is chosen automatically based on the trend line and
#' \code{\link{auto.aspect}}.
#' @param given,given.lab,nlevels,pretty,key,bg,color.palette,pch.palette used for conditioning plots.
#' @param main,xlab,ylab axis labels.
#' @param ... extra arguments passed to \code{predict_plot.data.frame} or
#' \code{\link{plot}}.
#' @return
#' If the predictor is numeric, makes a scatterplot with loess line on top.
#' If the predictor is a factor, makes a \code{\link{linechart}}.
#'
#' @seealso \code{\link{loess}}, \code{\link{model.plot}}
#' @examples
#' data(Cars)
#' predict_plot(Price ~ ., CarsT)
#' fit <- lm(Price ~ ., CarsT)
#' predict_plot(Price~ ., CarsT, partial=fit)
#' # same thing using predict_plot.lm
#' predict_plot(fit, partial = TRUE)
#' @export
<- () ("predict_plot")
#' @export
predict_plot.default <- (x,y,,xlab,ylab,) {
= (x,y)
(!()) () =
((xlab)) xlab = ((x))
((ylab)) ylab = ((y))
(,xlab=xlab,ylab=ylab,)
}
std.error = (x) (x)/((x))
#' @export
<- (x,=,partial=,condensed=T,
rtype=("prob","logit","probit"),
highlight,se=T,
identify.pred=F,
mcol=0,mlwd=2,
scol="red",slwd=2,
span=2/3,degree=2,="symmetric",
mar=,bg=("bg"),
xaxt=("xaxt"),yaxt=("yaxt"),
=1,=,
given=,given.lab=,=2,
=T,key=!(given),
color.palette=,
pch.palette=default.pch,
main=,xlab,ylab=,) {
rtype <- (rtype)
resp <- (x)
((degree) && (x[resp]])) degree = 0
pred <- predictor.terms(x)
x <- expand.frame(x,(pred,predictor.vars(x)))
(!(partial)) {
((highlight)) highlight <- predictor.terms(partial)
(!(given.lab)) {
(given.lab %in% predictor.vars(partial)) {
partial = (partial,(".~.-",given.lab))
}
}
pres <- partial.residuals(partial,pred,condensed=condensed)
((ylab)) ylab = ("Partial",resp)
}
{
((highlight)) highlight <-
((ylab)) ylab = resp
}
(!(given)) {
# key==1 means separate panel, key==2 means on each panel
(key==2 && (main)) main = given.lab
b =
(!(given)) {
# make given a factor
b <- (given,,=)
given <- rcut(given,b)
} = ((given))
((color.palette)) color.palette = color.palette()
lev = 1:
level.color = color.palette((lev-1) %% (color.palette))+1]
level.pch <- pch.palette(((lev-1) %/% (color.palette)) %% (pch.palette))+1]
}
(()) <- auto.layout((pred) + (key==1))
opar <- (bg=bg)
(( != (1,1))) {
opar = (opar,(mfrow=))
}
(!(mar)) {
((mar)) {
mar = auto.mar(main=main,xlab=((xlab)) "y" xlab,
ylab=((ylab)) "y" ylab,xaxt=xaxt,yaxt=yaxt,key=(key==2))
}
opar <- (opar,(mar=mar))
}
((opar))
(i pred) {
(!(partial)) {
# sometimes pres is longer than x, because of NAs
x[resp]] <- pres[i]][rownames(x)]
((x[resp]])) (("partial of",i,"not found"))
}
((x[i]]) && !(x[i]])) {
x[i]] <- (x[i]],(x[resp]]))
}
# plot points
((given)) {
# no given
(!(scol)) {
# smooth curve
#k <- is.finite(x[,i]) & is.finite(x[[resp]])
# keep warnings about infinite points
k <- !(x,i]) & !(x[resp]])
fit =
(scol != 0) {
((x[resp]])) {
# factor response
(((x[resp]]))==2 && !(x[i]])) {
(rtype=="prob") {
fit = loprob(x[k,i], x[k,resp], degree=degree, span=span)
(F) {
# piecewise constant approx
(tree)
= (x=x[k,i],y=x[k,resp])
tr = tree(y~x,)
partition.tree(tr,add=T,=F)
}
} {
fit <- loprob(x[k,i], x[k,resp])
p <- fit$y-1
p <- (rtype,logit=(p/(1-p)),probit=(p))
fit$y <- p+1.5
(x[i]],2*(x[resp]])-3,xlab=i,ylab=rtype)
}
}
} {
# numeric response
(!(x[i]])) {
(!((x[k,i])))
((i,"has infinite values.\nCheck your transformation."))
(!((x[k,resp])))
((resp,"has infinite values.\nCheck your transformation."))
((k) > 1500) {
# save time on dense plots
fit = (x[k,i], x[k,resp], f=span)
} {
fit = lowess2(x[k,i], x[k,resp], f=span, degree=degree, =)
}
}
}
}
}
real.asp = ((,"auto")) (fit)
((x[resp]])) {
# factor response
(rtype=="prob") {
((x[i]])) {
# should sort if the factor is not ordered
#mosaicplot(table(x[[i]], x[[resp]]), xlab=i, ylab=ylab, ...)
tab = (x[resp]],x[i]])
((tab)) = (resp,i)
(row.probs(tab,se=T), xlab=" ",)
} {
(x[i]], x[resp]], xlab="", ylab=ylab, )
}
}
} ((x[i]])) {
# factor predictor, numeric response
= rtable(((resp,"~",i)),x)
se = rtable(((resp,"~",i)),x,fun=std.error) * 1.64
(,se,=,ylab=ylab,xlab=" ",)
} {
# numeric response and predictor
(x[i]], x[resp]], xlab="", ylab=ylab, main=main,
=real.asp, =, )
(!(partial) && (partial,"lm")) {
a <- (partial)
#a0 = a["(Intercept)"]
#if(is.null(a0)) a0 = 0
a0 = -(x[i]])*a[i]
(mcol != 0 && !(a[i])) (a0,a[i],=mcol,lwd=mlwd)
}
}
((xlab)) {
(i %in% highlight) (xlab=i,font.lab=2)
(xlab=i)
} (xlab=xlab)
(!(fit)) (fit,=scol,lwd=slwd)
# identify
((identify.pred == T) || (i %in% identify.pred)) {
(x[k,i],x[k,resp],=(x)[k])
}
} ((x[i]])) {
# given
fmla = ((resp,"~",i,"+",given.lab))
x[given.lab]] = given
= rtable(fmla,x)
(se) {
rt.se = rtable(fmla,x, fun=std.error) * 1.64
rt.se[is.na(rt.se)] = (x[resp]])
} rt.se =
((),(rt.se),xlab=" ",
xaxt=xaxt,yaxt=yaxt,row.labels=,)
font.lab = (i %in% highlight) 2 1
(xlab=i,font.lab=font.lab)
} {
# condition on the given variable
()
()
(main=main,ylab=ylab,=(key==1) 1)
(i %in% highlight) (xlab=i,font.lab=2)
(xlab=i)
((x[i]])) {
# setup window for factor
<- ((x[i]]))
(=(0.5,+0.5),=(x[resp]]))
(xaxt != "n") (1,1:,=(x[i]]))
(yaxt != "n") (2)
(("jittering",i,"\n"))
} {
# setup window for numeric
(=(x[i]],na.rm=T),=(x[resp]],na.rm=T))
(xaxt != "n") (1)
(yaxt != "n") (2)
}
lev <- (given)
(g 1:(lev)) {
color <- level.color[g]
<- level.pch[g]
val <- lev[g]
k <- !(x,i]) & !(x,resp])
k = k & (given == val)
((x[i]])) {
<- (((x[k,i]))-0.5)/5
((x[k,i])+, x[k,resp], =color,
=, )
} {
(x[k,i], x[k,resp], =color, =, )
}
((x[resp]])) {
fit = loprob(x[k,i], x[k,resp])
} {
fit = lowess2(x[k,i], x[k,resp], f=span, degree=degree)
}
(fit,=color,lwd=slwd)
}
}
}
# show legend for the given variable, as a separate panel
# this is better than putting a color key on each plot
(!(given)) {
(key==1) {
()
(!(given.lab)) {
font.lab = (given.lab %in% highlight) 2 1
(xlab=given.lab,font.lab=font.lab)
}
lev = (given)
y <- ((lev)+0.02)
(i 1:(lev)) {
color <- color.palette[i]
val <- lev[i]
(0.5,0.5+y[i]-y[length(y)]/2,val,=color,adj=0.5)
(((level.pch))>1)
(0.6+(val)/2,0.5+y[i],=level.pch[i])
}
} (key==2) {
i = 1:((given))
((b))
(color.palette[i],level.pch[i],(given))
(color.palette[i],level.pch[i],breaks=b)
}
}
}
#' Plot predictors versus residuals.
#'
#' Makes a matrix of pairwise scatterplots with lowess-type trend lines.
#'
#' Partial residuals are computed by fitting a new model with the
#' predictor removed, which is different from \code{residuals(type="partial")}.
#' @author Tom Minka
#' @param object the output of \code{lm}.
#' @param data a data frame to use instead of \code{model.frame(object)}.
#' @param partial If \code{TRUE}, plot partial residuals instead of residuals.
#' @param ylab axis label.
#' @param ... extra arguments for \code{\link{predict_plot.data.frame}}.
#' @return
#' A plot similar to \code{\link{predict_plot}}, but where the vertical
#' axis is residuals. These plots can be used to judge which predictors
#' should be added to the model.
#' @seealso \code{\link{predict_plot}}
#' @examples
#' # see the examples for predict_plot
#' @export
<- (object,,partial=F,ylab=,) {
(!partial) {
((ylab)) ylab = ("Residual",(object))
(()) {
res <- residual.frame(object)
} {
res <- residual.frame(object,)
}
("plotting residuals\n")
(res,highlight=predictor.terms(object),ylab=ylab,)
} {
(()) <- (object)
= my.model.frame((((object),"~.")),)
("plotting partial residuals\n")
(,partial=object,ylab=ylab,)
}
}
given.vars <- () {
rhs <- [3]]
((rhs) && ((rhs[1]]) == "|")) {
# remove givens from formula
(rhs[3]])
}
}
remove.given <- () {
rhs <- [3]]
[3]] <- rhs[2]]
}
#' @exportS3Method
<- (,=(),) {
# formula has givens?
given = given.vars()
(!(given)) {
= remove.given()
(()) g <- (given,env=)
g <- [given]]
((g))
(("variable \"",given,"\" not found",sep=""))
((,,
given=g,given.lab=given,))
}
(F) {
expr <- (expand = F)
expr$ <-
expr$ <-
expr[1]] <- ("model.frame.default")
x <- (expr, ())
} {
# formula has its own environment
x <- (,,=)
}
(x,)
}
# only used once in lab 10
# expand.cross <- function(object) {
# resp <- response.var(object)
# pred <- predictor.vars(object)
# # quadratic terms only
# pred2 <- c()
# for(i in pred) {
# for(j in pred) {
# if(match(i,pred) < match(j,pred)) pred2 = c(pred2,paste(i,j,sep=":"))
# }
# }
# pred = c(pred,pred2)
# formula(paste(resp,"~",paste(pred,collapse="+")))
# }
# step.up <- function(object,scope=expand.cross(object)) {
# step(object,list(upper=scope,lower=formula(object)))
# }
#' Contour plot matrix
#'
#' @examples
#' fit <- lm(Fertility ~ ., data=swiss)
#' interact.plot(fit)
#' @export
<- (object, ) ("interact.plot")
#' @export
interact.plot.lm <- (object,,partial=F,main=,) {
(!partial) {
(()) {
res <- residual.frame(object)
} {
res <- residual.frame(object,)
}
("plotting residuals\n")
((main)) main = "Residuals"
interact.plot.data.frame(res,highlight=predictor.terms(object),main=main,)
} {
(()) <- (object)
= my.model.frame((((object),"~.")),)
("plotting partial residuals\n")
interact.plot.data.frame(,partial=object,main=main,)
}
}
#' @export
interact.plot.data.frame <- (x,ypred,partial=,highlight,span=0.75,
scol="red",slwd=2,=("*",":"),
xaxt="n",yaxt="n",se=T,
bg=("bg"),nlev=8,main=,) {
# type=":" means partials remove cross term only
# type="*" means partials remove everything
#library(modreg)
= ()
resp <- (x)
pred = predictor.vars(x)
((ypred)) ypred = pred
((bg) & ((x[pred],))) bg = (0.65)
n = (pred)
mar=(0.05,0.05,0,0.05)
((ypred == pred)) = (n,n)
= auto.layout(n*(ypred))
opar = (mfrow=,mar=mar,bg=bg)
((opar))
(!(partial)) {
((highlight)) highlight <- predictor.terms(partial)
((main)) {
main = ("Partial residuals\n",
( == ":") "cross-term" "both predictors","absent",sep="\n")
}
}
((highlight)) highlight <-
(py ypred) {
(px pred) {
(px == py) {
()
()
y = 0.5
((x[py]])) {
auto.legend((x[py]]),((x[py]])))
#y = 0.5
}
(0.5,y,py,font=(py %in% highlight) 2 1)
} (!(x[py]]) && ((py,pred) < (px,pred))) {
()
((py,pred) == 1 && (px,pred) == 2)
(0.5,0.5,main)
} {
# do this to convert the names so formula(y) will work
y = (x[c(px,py,resp)])
(!(partial)) {
predxy = (px,,py,sep="")
fit = (partial,(".~.-",predxy))
y[resp]] = (fit)
}
mar = (0.05,0.05,0,0.05)
((y[px]]) && (y[py]])) {
= rtable((y),y)
xlab = (px %in% ypred) "" px
(se) {
rt.se = rtable((y),y,fun=std.error) * 1.64
rt.se[is.na(rt.se)] = (x[resp]])
} rt.se =
((),(rt.se),xlab=xlab,ylab="",yaxt=yaxt,=,)
} ((y[px]])) {
()
} ((y[py]])) {
fmla = ((resp,"~",px,"|",py))
(fmla,y,xlab="",ylab="",xaxt=xaxt,yaxt=yaxt,key=F)
} ((y[resp]])) {
(y,
xlab="",ylab="",mar=mar,key=F,axes=F,n=nlev,)
} {
(((y),y,span=span),
xlab="",ylab="",zlab="",mar=mar,key=F,axes=F,n=nlev,)
}
predxy = (px,":",py,sep="")
predyx = (py,":",px,sep="")
(((predxy,predyx) %in% highlight)) {
(="red")
} {
()
}
}
}
}
}
(!("my.factor.scope")) old.factor.scope =
my.factor.scope <- (, scope)
{
<- scope$
add <- scope$add
(() && !()) {# have base model
nmdrop <- ()
facs <-
(()) {
nmfac <- ()
where <- (nmdrop, nmfac, 0)
((!where)) ("lower scope is not included in model")
facs <- , -where, = ]
nmdrop <- nmfac-where]
} nmdrop <- ()
# minka: allow non-hierarchy
} nmdrop <- (0)
(!(add)) nmadd <- (0)
{
nmfac <- ()
nmadd <- (add)
(!(nmfac)) {
where <- (nmfac, nmadd, 0)
((!where)) ("upper scope does not include model")
nmadd <- nmadd-where]
add <- add, -where, = ]
}
# minka: allow non-hierarchy
}
( = nmdrop, add = nmadd)
}
#replaceInNamespace("factor.scope",my.factor.scope)
#replaceInNamespace("factor.scope",old.factor.scope)
