Re: [MORPHMET] pairwise matrix of vector angles in R
Hi David,
What you touched on is the art of statistical computing. You used a logical
function to only calculate angles between non-identical vectors and avoid NaN
values; our geomorph function turns off warnings and waits until the end and
replaces what should be computational 0 values with actual 0 values. Either
way accomplishes the same goal. In our case, the method is faster, which is
important when performing the same operation over thousands of permutations, as
in a hypothesis test that generated distributions of null angles. Another
benefit is the tcrossprod function maintains row names so that there is no
confusion about the order of vectors, especially with large pairwise matrices.
Although your functions resemble more the conceptual goal one wants to
accomplish, they are painfully slow over many permutations. I have
experimented with a function like yours (uses loops to fill in values of a
pairwise matrix) and found that it is actually better to create a distance
matrix, populate it, and then coerce the distance matrix into a symmetric
matrix. Doing this has two advantages. First, you are not duplicating
operations, like finding the angle between vectors 1 and 3 and then again
between vectors 3 and 1. Second, when coercing a distance matrix into a
symmetric matrix, the diagonal is automatically 0. Here is your code augmented
to do this:
# matrix where each vector is a row
vec.mat <- rbind(...)
# empty distance matrix-like matrix for pairwise angles
angle.mat <- as.dist(matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat)))
# angle between vectors function, in degrees
vec.angle <- function(v1,v2){
acos((t(v1)%*%v2)/(sqrt(sum(v1^2))*sqrt(sum(v2^2*180/pi}
# create a matrix of unique pairwise comaprisons
pw.mat <- combn(nrow(vec.mat), 2)
# fill the all pairwise angles matrix
for(i in 1:ncol(pw.mat)){
angle.mat[i] <- vec.angle(vec.mat[pw.mat[1,i],], vec.mat[pw.mat[2,i],])
}
# turn distance-like matrix into symmetric matrix
angle.mat <- as.matrix(angle.mat)
This approach is computationally MUCH more efficient without having to trap
potential warnings or errors. What we use in geomorph is even more efficient
than this but has less of a conceptual connection to the goal (one has to
recognize that we are first calculating inverse vector length as weights from
the diagonal of a cross-product matrix calculations, multiplying vectors by
their weights to make them unit length, then obtain the outer-product matrix,
having inner-products between unit length vectors as elements, which is vector
correlation by definition). If you were to rip open our geomorph code - and
that is possible, thanks to R - you will find it is replete with many
computational efficiencies that might seem decoupled from the purest
definitions of the statistics we calculate. Statistical computing is both
cumbersome and elegant this way.
Happy computing!
Mike
> On Nov 5, 2017, at 4:48 PM, David Katz wrote:
>
> Hi everyone,
>
> Apologies to the group. In retrospect, despite the disclaimer in my first
> email, better to have error checked before sending (as Mike noted, I wrote
> acos not cos, and I borrowed from old code that assumed the vectors were
> already unit length; I also made mistakes in how the loop was set up;
> overall, a mess).
>
> I might as well provide the corrected, checked code. Afterward, I explain why
> the function to compute the angle between vectors (vec.angle) requires a test
> for whether the vectors being compared are duplicates of one another
> ("any(duplicated(.))")...
>
> # matrix where each vector is a row
> vec.mat <- rbind(...)
>
> # empty square matrix for pairwise angles
> angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
>
> # angle between vectors function, in degrees
> vec.angle <- function(v1,v2){
> ifelse(any(duplicated(rbind(v1, v2))),0,
> acos((t(v1)%*%v2)/(sqrt(sum(v1^2))*sqrt(sum(v2^2*180/pi)}
>
> # fill the all pairwise angles matrix
> for(i in 1:nrow(vec.mat)){for(j in 1:nrow(vec.mat)){
> angle.mat[i, j] <- vec.angle(v1=vec.mat[i,], v2=vec.mat[j,])}}
>
> The cosine of the angle between two identical vectors (call it CTI) is 1. But
> R will occasionally compute a value greater than one. This is not immediately
> apparent: If CTI is called, the value displayed in the console will be 1;
> however, use of the logical operator CTI==1 returns FALSE, and the logical
> CTI > 1 returns TRUE. I think it's a machine precision issue. Because cosine
> ranges from -1 to 1, there's no such thing as an arccosine for a value
> greater than 1. R will return Nan for acos(CTI). Incorporating
> any(duplicated(.)) into the vec.angle function assures that identical vectors
> will have an angle between them of 0. This will be so whether a vector is
> being compared to itself (i = j, hence deposited along the diagonal of the
> matrix of angles), or to an identical vector in another row (i != j, hence
Re: [MORPHMET] pairwise matrix of vector angles in R
Hi everyone,
Apologies to the group. In retrospect, despite the disclaimer in my first
email, better to have error checked before sending (as Mike noted, I wrote
acos not cos, and I borrowed from old code that assumed the vectors were
already unit length; I also made mistakes in how the loop was set up;
overall, a mess).
I might as well provide the corrected, checked code. Afterward, I explain
why the function to compute the angle between vectors (vec.angle) requires
a test for whether the vectors being compared are duplicates of one another
("any(duplicated(.))")...
# matrix where each vector is a row
vec.mat <- rbind(...)
# empty square matrix for pairwise angles
angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
# angle between vectors function, in degrees
vec.angle <- function(v1,v2){
ifelse(any(duplicated(rbind(v1, v2))),0,
acos((t(v1)%*%v2)/(sqrt(sum(v1^2))*sqrt(sum(v2^2*180/pi)}
# fill the all pairwise angles matrix
for(i in 1:nrow(vec.mat)){for(j in 1:nrow(vec.mat)){
angle.mat[i, j] <- vec.angle(v1=vec.mat[i,], v2=vec.mat[j,])}}
The cosine of the angle between two *identical* vectors (call it CTI) is 1.
But R will *occasionally* compute a value greater than one. This is not
immediately apparent: If CTI is called, the value displayed in the
console will be 1; however, use of the logical operator CTI*==*1 returns
FALSE, and the logical CTI* > *1 returns TRUE. I think it's a machine
precision issue. Because cosine ranges from -1 to 1, there's no such thing
as an arccosine for a value greater than 1. R will return Nan for acos(CTI)
. Incorporating any(duplicated(.)) into the vec.angle function assures that
identical vectors will have an angle between them of 0. This will be so
whether a vector is being compared to itself (*i = j, *hence deposited
along the diagonal of the matrix of angles), or to an identical vector in
another row (*i != j, *hence deposited in an off-diagonal cell).
I think the line diag(vc) = 0 in Mike's vec.angle.matrix code is meant to
do the same thing for his diagonal values.
Again, sorry for any confusion I caused. Thanks again to Mike.
Best, David
On Fri, Nov 3, 2017 at 10:46 AM, Mike Collyer wrote:
> Dear David, and others,
>
> Be careful with the code you just introduced here. There are a couple of
> mistakes. First, vectors need to be unit length. You code does not
> transform the vectors to unit length. Second, it’s the arccosine, the
> cosine of the vector inner product that finds the angle. Third, though
> less of an error but more of a precision issue, there is no need to round
> 180/pi, as you have done. This is what R will return if you ask it to
> divide 180 by pi, but the four decimal places are for your visual benefit.
> By using 180/pi instead of 57.2958, the results will be more precise. This
> last part is not much issue but the first two are big problems.
>
> Here is a demo using two x,y vectors - 0,1 and 1,1 - which we know has an
> angle of 45 degrees between them in a plane. I calculate the angles using
> your method and one with unit length vectors and arccosine of their inner
> products, which is what we use in geomorph.
>
> > v1 <- c(1,0)
> > v2 <- c(1,1)
> > V <- rbind(v1, v2)
> >
> > # Katz method
> > vec.mat <- V
> > angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
> > vec.angle <- function(v1,v2){cos(t(v1)%*%v2)*57.2958}
> > for(i in 1:nrow(vec.mat)) {
> + for(j in 1:nrow(vec.mat)) {
> + angle.mat[i, j] <- vec.angle(vec.mat[i,], vec.mat[j,])
> + }
> + }
> >
> > angle.mat
> [,1] [,2]
> [1,] 30.95705 30.95705
> [2,] 30.95705 -23.84347
> >
> > geomorph:::vec.ang.matrix(V, type = "deg")
>v1 v2
> v1 0 45
> v2 45 0
>
> As you can see, using the cosine and not transforming the vectors to unit
> length finds results that suggest there is some non-0 degree angle between
> a vector and itself, in addition to the incorrect angle between vectors.
>
> You indicated that you were patching together code, so it might be an
> oversight, but it is an important distinction for others who might use the
> code. Here are the “guts” of the geomorph functions in case somebody wants
> to reconcile the points I just made with the set-up in your code, which is
> similar.
>
> > geomorph:::vec.cor.matrix
> function(M) {
> M <- as.matrix(M)
> w <- 1/sqrt(diag(tcrossprod(M))) # weights used to make vectors unit
> length
> vc = tcrossprod(M*w)# outer-product matrix finds inner-products between
> vectors for all elements
> options(warn = -1) # turn off warnings for diagonal elements
> vc # vector correlations returned
> }
>
> > geomorph:::vec.ang.matrix
> function(M, type = c("rad", "deg", "r")){
> M <- as.matrix(M)
> type <- match.arg(type)
> if(type == "r") {
> vc <- vec.cor.matrix(M) # as above
> } else {
> vc <- vec.cor.matrix(M)
> vc <- acos(vc) # finds angles for vector correlations
> diag(vc) = 0 # Make sure computational 0s are true 0s
> }
> if(t
Re: [MORPHMET] pairwise matrix of vector angles in R
Dear David, and others,
Be careful with the code you just introduced here. There are a couple of
mistakes. First, vectors need to be unit length. You code does not transform
the vectors to unit length. Second, it’s the arccosine, the cosine of the
vector inner product that finds the angle. Third, though less of an error but
more of a precision issue, there is no need to round 180/pi, as you have done.
This is what R will return if you ask it to divide 180 by pi, but the four
decimal places are for your visual benefit. By using 180/pi instead of
57.2958, the results will be more precise. This last part is not much issue
but the first two are big problems.
Here is a demo using two x,y vectors - 0,1 and 1,1 - which we know has an angle
of 45 degrees between them in a plane. I calculate the angles using your
method and one with unit length vectors and arccosine of their inner products,
which is what we use in geomorph.
> v1 <- c(1,0)
> v2 <- c(1,1)
> V <- rbind(v1, v2)
>
> # Katz method
> vec.mat <- V
> angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
> vec.angle <- function(v1,v2){cos(t(v1)%*%v2)*57.2958}
> for(i in 1:nrow(vec.mat)) {
+ for(j in 1:nrow(vec.mat)) {
+ angle.mat[i, j] <- vec.angle(vec.mat[i,], vec.mat[j,])
+ }
+ }
>
> angle.mat
[,1] [,2]
[1,] 30.95705 30.95705
[2,] 30.95705 -23.84347
>
> geomorph:::vec.ang.matrix(V, type = "deg")
v1 v2
v1 0 45
v2 45 0
As you can see, using the cosine and not transforming the vectors to unit
length finds results that suggest there is some non-0 degree angle between a
vector and itself, in addition to the incorrect angle between vectors.
You indicated that you were patching together code, so it might be an
oversight, but it is an important distinction for others who might use the
code. Here are the “guts” of the geomorph functions in case somebody wants to
reconcile the points I just made with the set-up in your code, which is similar.
> geomorph:::vec.cor.matrix
function(M) {
M <- as.matrix(M)
w <- 1/sqrt(diag(tcrossprod(M))) # weights used to make vectors unit length
vc = tcrossprod(M*w)# outer-product matrix finds inner-products between
vectors for all elements
options(warn = -1) # turn off warnings for diagonal elements
vc # vector correlations returned
}
> geomorph:::vec.ang.matrix
function(M, type = c("rad", "deg", "r")){
M <- as.matrix(M)
type <- match.arg(type)
if(type == "r") {
vc <- vec.cor.matrix(M) # as above
} else {
vc <- vec.cor.matrix(M)
vc <- acos(vc) # finds angles for vector correlations
diag(vc) = 0 # Make sure computational 0s are true 0s
}
if(type == "deg") vc <- vc*180/pi # turns radians into degrees
vc
}
These functions have some extras that would not pertain to the general solution
but are meant to trap warnings and round 0s for users, but they should not get
in the way of understanding.
Cheers!
Mike
> On Nov 3, 2017, at 11:33 AM, David Katz wrote:
>
> I think this does it (but please check; I quickly stuck together two
> different pieces of code)...
>
> # matrix where each vector is a row
> vec.mat <- ...
> #Compute group*loci matrix of mean microsatellite lengths
> angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
> # angle function (radians converted to degrees)
> vec.angle <- function(v1,v2){cos(t(v1)%*%v2)*57.2958}
> # angles-to-matrix loop
> for(i in 1:nrow(vec.mat))
> {
> for(j in 1:nrow(vec.mat))
> {
> # angle bw vec1 and vec2
> angle.mat[i, j] <- vec.angle(vec.mat[i,], vec.mat[j,])}
> }
> return(angle.mat)
> }
>
> On Fri, Nov 3, 2017 at 2:38 AM, andrea cardini > wrote:
> Dear All,
> please, does anyone know if there's an R package that, using a matrix with
> several vectors (e.g., coefficients for allometric regressions in different
> taxa), will compute the pairwise (all possible pairs of taxa) matrix of
> vector angles?
>
> Thanks in advance for any suggestion.
> Cheers
>
> Andrea
>
>
> --
>
> Dr. Andrea Cardini
> Researcher, Dipartimento di Scienze Chimiche e Geologiche, Università di
> Modena e Reggio Emilia, Via Campi, 103 - 41125 Modena - Italy
>
> tel. 0039 059 2058472
>
> Adjunct Associate Professor, School of Anatomy, Physiology and Human Biology,
> The University of Western Australia, 35 Stirling Highway, Crawley WA 6009,
> Australia
>
> E-mail address: alcard...@gmail.com ,
> andrea.card...@unimore.it
> WEBPAGE: https://sites.google.com/site/alcardini/home/main
>
>
> FREE Yellow BOOK on Geometric Morphometrics:
> http://www.italian-journal-of-mammalogy.it/public/journals/3/issue_241_complete_100.pdf
>
>
>
> ESTIMATE YOUR GL
Re: [MORPHMET] pairwise matrix of vector angles in R
I think this does it (but please check; I quickly stuck together two
different pieces of code)...
# matrix where each vector is a row
vec.mat <- ...
#Compute group*loci matrix of mean microsatellite lengths
angle.mat <- matrix(NA,nrow=nrow(vec.mat),ncol=nrow(vec.mat))
# angle function (radians converted to degrees)
vec.angle <- function(v1,v2){cos(t(v1)%*%v2)*57.2958}
# angles-to-matrix loop
for(i in 1:nrow(vec.mat))
{
for(j in 1:nrow(vec.mat))
{
# angle bw vec1 and vec2
angle.mat[i, j] <- vec.angle(vec.mat[i,], vec.mat[j,])}
}
return(angle.mat)
}
On Fri, Nov 3, 2017 at 2:38 AM, andrea cardini wrote:
> Dear All,
> please, does anyone know if there's an R package that, using a matrix with
> several vectors (e.g., coefficients for allometric regressions in different
> taxa), will compute the pairwise (all possible pairs of taxa) matrix of
> vector angles?
>
> Thanks in advance for any suggestion.
> Cheers
>
> Andrea
>
>
> --
>
> Dr. Andrea Cardini
> Researcher, Dipartimento di Scienze Chimiche e Geologiche, Università di
> Modena e Reggio Emilia, Via Campi, 103 - 41125 Modena - Italy
>
> tel. 0039 059 2058472
>
> Adjunct Associate Professor, School of Anatomy, Physiology and Human
> Biology, The University of Western Australia, 35 Stirling Highway, Crawley
> WA 6009, Australia
>
> E-mail address: alcard...@gmail.com, andrea.card...@unimore.it
> WEBPAGE: https://sites.google.com/site/alcardini/home/main
>
> FREE Yellow BOOK on Geometric Morphometrics:
> http://www.italian-journal-of-mammalogy.it/public/journals/3
> /issue_241_complete_100.pdf
>
> ESTIMATE YOUR GLOBAL FOOTPRINT: http://www.footprintnetwork.or
> g/en/index.php/GFN/page/calculators/
>
> --
> MORPHMET may be accessed via its webpage at http://www.morphometrics.org
> --- You received this message because you are subscribed to the Google
> Groups "MORPHMET" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to morphmet+unsubscr...@morphometrics.org.
>
>
--
David C. Katz, Ph.D.
Postdoctoral Fellow
Benedikt Hallgrimsson Lab
University of Calgary
Research Associate
Department of Anthropology
University of California, Davis
ResearchGate profile
Personal webpage
Re: [MORPHMET] pairwise matrix of vector angles in R
Andrea, If you already have the matrix of coefficients, geomorph has an internal function that can do what you seek to do. You can try this geomorph:::vec.ang.matrix(myMatrix, type = “r”) # for vector correlations geomorph:::vec.ang.matrix(myMatrix, type = “rad”) # for vector angles in radians geomorph:::vec.ang.matrix(myMatrix, type = “deg”) # for vector angles in degrees This function is an internal function for the advanced.procD.lm function, where one could test the angular differences among the allometric vectors of taxa. It assumes the matrix is arranges such that the rows are taxa vectors. Please note that if one uses advanced.procD.lm to perform this test, the matrix of pairwise angles is available as part of the results. We have an example for a homogeneity of slopes test in the help file that should be analogous to your example. There are two ways to obtain the results summary(HOS) # HOS is the test name; the pairwise angle matrix is part of the summary HOS$obs.slopes.angles Cheers! Mike > On Nov 3, 2017, at 4:38 AM, andrea cardini wrote: > > Dear All, > please, does anyone know if there's an R package that, using a matrix with > several vectors (e.g., coefficients for allometric regressions in different > taxa), will compute the pairwise (all possible pairs of taxa) matrix of > vector angles? > > Thanks in advance for any suggestion. > Cheers > > Andrea > > > -- > > Dr. Andrea Cardini > Researcher, Dipartimento di Scienze Chimiche e Geologiche, Università di > Modena e Reggio Emilia, Via Campi, 103 - 41125 Modena - Italy > tel. 0039 059 2058472 > > Adjunct Associate Professor, School of Anatomy, Physiology and Human Biology, > The University of Western Australia, 35 Stirling Highway, Crawley WA 6009, > Australia > > E-mail address: alcard...@gmail.com, andrea.card...@unimore.it > WEBPAGE: https://sites.google.com/site/alcardini/home/main > > FREE Yellow BOOK on Geometric Morphometrics: > http://www.italian-journal-of-mammalogy.it/public/journals/3/issue_241_complete_100.pdf > > ESTIMATE YOUR GLOBAL FOOTPRINT: > http://www.footprintnetwork.org/en/index.php/GFN/page/calculators/ > > -- > MORPHMET may be accessed via its webpage at http://www.morphometrics.org > --- You received this message because you are subscribed to the Google Groups > "MORPHMET" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to morphmet+unsubscr...@morphometrics.org. > -- MORPHMET may be accessed via its webpage at http://www.morphometrics.org --- You received this message because you are subscribed to the Google Groups "MORPHMET" group. To unsubscribe from this group and stop receiving emails from it, send an email to morphmet+unsubscr...@morphometrics.org.
