I am wanting to plot pca loadings onto sites superimposed on a contour
map. There are a maximum of 20 sites at which loadings might appear -
however due to the nature of my data, missing data has meant that some
stations have not been included in some of the pca.
For example, I am performing pca for 5 different variables at 20 different
sites for a number of different time intervals and want to perform a cross
comparison to see where key trends lie. Where data is not available for
any particular station at any given time interval I am still keen to
analyse the principal components on the others.
I have been using the following code to produce the plots I want:
for (i in 1:4) {
image(east.grid,north.grid,t(map.matrix),col=my.colors,
axes=TRUE,xlab="",ylab="")
text(stations$x.axis,stations$y.axis,labels=as.character(stations$station))
title(paste("Component",i))
pc <- p6.5am.pca$rotation[,i]
pc.pos <- (pc > 0)
if (any(pc.pos)) {
symbols(stations$x.axis[pc.pos],stations$y.axis[pc.pos],
circles=scale*pc[pc.pos],lwd=2,fg="red",inches=FALSE,add=TRUE)
}
if (any(!pc.pos)) {
symbols(stations$x.axis[!pc.pos],stations$y.axis[!pc.pos],
circles=-scale*pc[!pc.pos],lwd=2,lty=2,fg="blue",inches=FALSE,
add=TRUE)
}
box()
}
My problem is this - due to the large number of plots required I would
like to be able to batch process the results, but in the case of applying
the above code when one (or more) station(s) is missing from the pca I get
an error message as the x and y lengths differ. Is there a way I can get the
code
to simply return a zero loading at a site if it hasn't been included in
the calculation of the pca? (or is there another better way around this
problem?)
I hope that this question makes some sense.
Any ideas?
Thanks,
laura
