[6541] in s-news-athena
Summary: contours in triangles
daemon@ATHENA.MIT.EDU (MARCEL.BAUMGARTNER@01.CHLSNR.nestr)
Wed Mar 1 02:25:50 1995
Date: Wed, 1 Mar 1995 08:10:53 +0100
From: MARCEL.BAUMGARTNER@01.CHLSNR.nestrd.ch
To: s-news@utstat.toronto.edu
------------------------------ Start of body part 1
------------------------------ Start of body part 2
Sorry if this is a repost!
Hi,
some days ago I posted a query regarding contour plots in
triangles, especially useful when displaying results of ternary
mixture designs.
Thanks to everybody who replied. I think there is a need for some
routines for doing mixture analysis. I was told that not even
S+DOX includes such functions.
It turned out that writing what I need was not that difficult,
and in fact it is a very nice exercice in geometry and linear algebra.
In what follows, I include the I wrote. The functions are
very basic, a lot more can be done. Please feel free to contact me
for comments and improvements.
Marcel Baumgartner
Marcel.Baumgartner@chlsnr.nestrd.ch
*********************************************************************
*********** Plotting routines for analyzing 3-way *****************
*********** mixture designs *****************
*********************************************************************
Set par(pty="s") to get real equilateral triangles (and print
in portrait mode).
To plot an equilateral triangle:
================================
vertices <- matrix(c(0,0,.5,sqrt(3)/2,1,0,0,0),ncol=2,byrow)
plot(vertices,type="l",axes=F,xlab="",ylab="")
Top corner: component 1
Left corner: component 2
Right corner: component 3
To add one or several design points:
====================================
points(x%*%triangle.trans,pch=0)
where x is a design matrix (ncol=3, sum of rows = constant,
mixture values between 0 and 1) and
triangle.trans <- matrix(c(0.5,sqrt(3)/2,0,0,1,0),ncol=2,byrow=T)
To plot contours of the fitted polynomial (using, for example, lm()):
=====================================================================
Fill in the estimated coefficients in the following template
(this is not at all user-friendly, but I see no other way yet):
z.template <- function(x, y)
{
aux <- z.trans(x, y)
x1 <- aux$x1
x2 <- aux$x2
x3 <- aux$x3
#Linear model
b1 * x1 + b2 * x2 + b3 * x3
#Quadratic model
# b1 * x1 + b2 * x2 + b3 * x3 +
b12 * x1 * x2 + b13 * x1 * x3 + b23 * x2 * x3
#Full Cubic model
# b1 * x1 + b2 * x2 + b3 * x3 +
# b12 * x1 * x2 + b13 * x1 * x3 + b23 * x2 * x3 +
# d12 * x1 * x2 * (x1 - x2) + d13 * x1 * x3 * (x1 -x3) +
# d23 * x2 * x3 * (x2 - x3) +
# b123 * x1 * x2 * x3
#Special Cubic model
# b1 * x1 + b2 * x2 + b3 * x3 +
# b12 * x1 * x2 + b13 * x1 * x3 + b23 * x2 * x3 +
# b123 * x1 * x2 * x3
}
where
z.trans <- function(x, y, constant = 1)
{
x1 <- constant * ((2 * y)/sqrt(3))
x2 <- constant * ( - x - y/sqrt(3) + 1)
x3 <- constant * (x - y/sqrt(3))
list(x1 = x1, x2 = x2, x3 = x3)
}
The "constant" option is useful when your components do not sum up to
1, but to a fixed constant less than 1.
The following function does the rest:
function(z, npoints = 40, ..., levels, manual = F, labels = paste("X", 1:3, sep
= ""))
{
#Written by Marcel Baumgartner, Marcel.Baumgartner@chlsnr.nestrd.ch
#z: a function in the format of z.template
#npoints: number of points for contours
#levels: desired levels for contour
#manual: if TRUE, then levels can be put manually using locator()
#labels: labels for mixture components
x.con <- seq(0, 1, length = npoints)
y.con <- seq(0, sqrt(3)/2, length = npoints)
if(missing(levels))
con.out <- contour(x.con, y.con, outer(x.con, y.con, z), save
= T, plot = F)
else con.out <- contour(x.con, y.con, outer(x.con, y.con, z), save = T,
plot = F, levels = levels)
n <- length(con.out)
par(pty = "s")
plot(c(0.5, 0.5), type = "n", xlim = c(0, 1), ylim = c(0, 1), xlab = "",
ylab = "", bty = "n", axes = F)
for(i in 1:n) {
x <- con.out[[i]]$x
y <- con.out[[i]]$y
test <- rep(T, length(x))
test[y > sqrt(3)/2 | y/x > sqrt(3) | y + sqrt(3) * x > sqrt(3)] <-
F
lines(x[test], y[test])
if(!manual) {
place <- length(x[test][!is.na(x[test])])
place <- round(median(1:place),0)
text(x[place], y[place], names(con.out)[i])
}
}
vertices <- matrix(c(0, 0, 0.5, sqrt(3)/2, 1, 0, 0, 0), ncol = 2, byrow
= T)
lines(vertices)
text(matrix(c(0.5, sqrt(3)/2 * 1.05, -0.05, -0.05, 1.05, -0.05), ncol
= 2, byrow = T), labels)
if(manual) {
cat("Put manually contour labels: \n")
print(names(con.out))
for(i in 1:n)
text(locator(1), names(con.out)[i])
}
invisible(con.out)
}
An example (see page 50 of J. A. Cornell, "Experiments with Mixtures:
Designs, Models, and the Analysis of Mixture Data", John Wiley,
1981):
z.p _ function(x, y)
{
aux <- z.trans(x, y)
x1 <- aux$x1
x2 <- aux$x2
x3 <- aux$x3 #Linear model
# b1 * x1 + b2 * x2 + b3 * x3 #Quadratic model
# b1 * x1 + b2 * x2 + b3 * x3 +
# b12 * x1 * x2 + b13 * x1 * x3 + b23 * x2 * x3 #Full Cubic model
# b1 * x1 + b2 * x2 + b3 * x3 +
# b12 * x1 * x2 + b13 * x1 * x3 + b23 * x2 * x3 +
# d12 * x1 * x2 * (x1 - x2) + d13 * x1 * x3 * (x1 -x3) +
# d23 * x2 * x3 * (x2 - x3) +
# b123 * x1 * x2 * x3
#Special Cubic model
1.8 * x1 + 25.4 * x2 + 28.6 * x3 - 34.8 * x1 * x2 - 48.4 * x1 * x3 -
94.4 * x2 * x3 + 624.6 * x1 * x2 * x3
}
motif()
plot.mixture(z.p,levels=c(5,10,15,20))
Improvements to be made:
o Localisation of level labels
o Allow z to use for example coef(model) instead of manually entering
the model coefficients
------------------------------ End of body part 2
