[R] superimpose histogram on biplot

[Kevin Wright]

I know that this is a late follow-up, but I just found this posting
today while searching for something else.

Two of my colleagues and myself have explored the use of Mosaic plots
for displaying the amount of variation captured by each axis.  See the
following paper:

@Article{Laffont2007,
  author =       {Jean-Louis Laffont and Mohamed Hanafi and Kevin Wright},
  title =        {Numerical and Graphical Measures to Facilitate the
                  Interpretation of GGE Biplots},
  journal =      {Crop Science},
  year =         2007,
  volume =       47,
  pages =        {990-996},
  annote =       {File: Biplots}
}

Here is the example code at the end of the paper.

Y = data.frame(E1 = c(50, 55, 65, 50, 60, 65, 75),
E2 = c(67, 71, 76, 80, 82, 89, 95),
E3 = c(90, 93, 95, 102, 97, 106, 117),
E4 = c(98, 102, 105, 130, 135, 137, 133),
E5 = c(120, 129, 134, 138, 151, 153, 155))
rownames(Y) = c("G1", "G2", "G3", "G4", "G5",
"G6", "G7")
Y = scale(Y, center = TRUE, scale = FALSE)
YG = matrix(rowMeans(Y)) %*% rep(1,ncol(Y))
YGE = Y - YG
Ysvd = svd(Y) # Singular value decomposition
G = U = Ysvd$u
H = Ysvd$v %*% diag(Ysvd$d)
# Formula: SSGk = diag(U´ YG YG´ U)
SSGk = diag(crossprod(crossprod(YG,U)))
# Formula: SSGEk = diag(U´ YGE YGE´ U)
SSGEk = diag(crossprod(crossprod(YGE,U)))
# Proportion of SS for each axis
TSS = sum(Y^2)
AxisProp = round(100*(SSGk + SSGEk)/TSS,0)
# First, the environment-focused GGE biplot
biplot(G[,1:2], H[,1:2], xlabs = rownames(Y),ylabs =
colnames(Y),
xlab = paste("Axis 1: ",AxisProp[1],"%"),
ylab = paste("Axis 2: ",AxisProp[2],"%"))
# Then the mosaic plot
mosaicdata = data.frame(G = SSGk, GE = SSGEk)
mosaicplot(mosaicdata, main = "", off = c(5,1),
cex.axis = 1, col = c("dimgray", "snow2"))

Kevin Wright


On Tue, Apr 1, 2008 at 11:37 AM, Stéphane Dray
<dray at biomserv.univ-lyon1.fr> wrote:
> Hi Jennie,
> if you have questions about ade4, you can use the adelist
> http://listes.univ-lyon1.fr/wws/info/adelist
>
> ade4 has the add.scatter.eig function which is exactly what you need. Try
>
> add.scatter.eig(nipmat$eig, xax=1, yax=2)
>
> Cheers,
>
>
>
> Jennie Lavine wrote:
>> Hi all,
>>
>> I've been trying to figure out how to superimpose a histogram on a
>> biplot that shows the relative contribution of each axis.  I have
>> been using the NIPALS function (http://biomserv.univ-lyon1.fr/~dray/
>> files/softwares/nipals.R) to run principal component analyses.  Here
>> is a toy example.
>>
>> source("http://biomserv.univ-lyon1.fr/~dray/files/softwares/nipals.R")
>> mat=matrix(runif(100,0,1), ncol=10, nrow=10)
>> nipmat=nipals(mat, nf=9)
>> scatter(nipmat)
>>
>> In the plot generated by the above "scatter" command, there is a
>> histogram in the upper left corner.  I want to know how to
>> superimpose that histogram on a similar plot, such as the following:
>>
>> groups = as.factor(c(rep(1,2), rep(2,4), rep(3,4)))
>> s.arrow(dfxy=nipmat$co[,1:2]*8, sub="Day 10", possub="bottomleft",
>> csub=3)
>> s.class(dfxy=nipmat$li[,1:2], fac=groups, cellipse=2, axesell=F,
>> cstar=0 , col=c(2:3), add.plot=T)
>>
>> I can create the histogram using:
>>
>> plot(nipmat$eig/sum(nipmat$eig), type='h')
>>
>> but I don't know how to superimpose it on the above graph.
>>
>> Thanks for any  help!
>>
>> Best,
>> Jennie
>>
>>
>> Jennie Lavine
>> PhD Candidate
>> Center for Infectious Disease Dynamics
>> Penn State University
>> 505 ASI Building
>> (814)863-1815
>>
>>
>>       [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at 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.
>>
>>
>>
>
> --
> Stéphane DRAY (dray at biomserv.univ-lyon1.fr )
> Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - Lyon I
> 43, Bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France
> Tel: 33 4 72 43 27 57       Fax: 33 4 72 43 13 88
> http://biomserv.univ-lyon1.fr/~dray/
>
> ______________________________________________
> R-help at 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.
>