Hi R Gurus & Lurkers... Thanks in advance to anyone who is willing to
tackle this! Bryan
I have been implementing the graphical manova method described in "An
Introduction to Ggobi" (from the Ggobi web site). A stand alone working
code is appended below. The code is almost the same as described in the
"Introduction document", with one bug fix. A quick summary of the code is
that it takes one's data, and fits an ellipsoid to it at a requested
confidence level, then color codes everything for display. If you don't
have ggobi installed, remove the ggobi from the last line and just use
>graphic_manova(response,class)
You will probably have to comment out the last 4 lines of the graphic_manova
function as well to avoid trivial errors.
Here's the R question: If the variable "class" has more than two levels
(factors) in it, the code executes but runs into an error because
cis = lapply(sub.groups, combined, cl=cl)
creates "cis" with a bunch of NA's, which then cause havoc when one tries to
do any matrix operations on it (not surprisingly). The NA's follow an
interesting pattern: the ellipsoid points are generated for the first two
dimensions (pc1 and pc2), but NA's are generated for the third dimension
(pc3). So cis contains the 3 original data dimensions, 1000 added ellipsoid
points to go with pc1, and 1000 added ellipsoid points to go with pc2, and
1000 NA's to go with pc3.... I don't see why the third set of data is any
different than the first two, and the first two execute correctly.
********
# generate sample data
pc1=rnorm(20, sd=1)
pc2=rnorm(20, mean = 10, sd=2)
pc3=rnorm(20, sd=3)
class=factor(c("group 1","group 1","group 1","group 2","group 2","group
2","group 2","group 2","group 2","group 2","group 1","group 1","group
1","group 1","group 2","group 2","group 2","group 2","group 1","group 2"),
ordered=TRUE)
response=cbind(pc1, pc2, pc3)
# Now generate confidence ellipsoids using the method described
# in "An Introduction to RGGOBI" with minor modifications
# Define 3 functions to do the heavy lifting
# First: a function that generates a random set of points on a sphere
# centered on the mean of the passed data, skewed to match the variance
# of the passed data (which turns the sphere into an ellipsoid),
# and adjusted in size to match the desired confidence level.
ellipse = function(data, npoints=1000, cl=0.95, mean=colMeans(data),
cov=var(data),n=nrow(data))
{
norm.vec = function(x) x/sqrt(sum(x^2))
p = length(mean)
ev = eigen(cov)
sphere = matrix(rnorm(npoints*p), ncol=p)
cntr = t(apply(sphere, 1, norm.vec)) # normalized sphere
cntr = cntr %*% diag(sqrt(ev$values)) %*% t(ev$vectors) # ellipsoid of
correct shape
Fcrit = qf(cl, p, n-p)
scalefactor = sqrt((p*(n-1))/(n*(n-p)))*Fcrit
cntr = cntr*scalefactor # ellipsoid of correct size
if (!missing(data)) # only relevant when no data passed
colnames(cntr) = colnames(data)
cntr+rep(mean, each=npoints)
}
# Next a function that combines the original data with the generated
ellipsoid
combined = function(data, cl=0.95)
{
dm = data.matrix(data)
ellipse = as.data.frame(ellipse(dm, npoints=1000, cl=cl))
both = rbind(data, ellipse)
both$SIM = factor(rep(c(FALSE,TRUE),c(nrow(data),1000)))
both
}
# Now a function to separate the dataset into categories
graphic_manova = function(data, catvar, cl=0.68)
{
sub.groups = data.frame(cbind(data,catvar))
sub.groups = split(sub.groups,catvar)
cis = lapply(sub.groups, combined, cl=cl)
df = as.data.frame(do.call(rbind, cis))
df$var = factor(rep(names(cis), sapply(cis, nrow)))
g = ggobi(df)
glyph_type(g[1]) = c(6,1)[df$SIM] # makes dots of ellipsoids tiny
glyph_color(g[1]) = df$var # properly colors the two groups
invisible(g)
}
# Now actually do the computations & plot the data!
# ggobi(combined(response)) # This is a debugging check point
ggobi(graphic_manova(response,class))
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
