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 <dck...@ucdavis.edu> 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
> 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 <mlcoll...@gmail.com
> <mailto:mlcoll...@gmail.com>> 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(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 <dck...@ucdavis.edu
>> <mailto:dck...@ucdavis.edu>> 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 <alcard...@gmail.com
>> <mailto:alcard...@gmail.com>> 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
>> <https://maps.google.com/?q=Emilia,+Via+Campi,+103+-+41125+Modena+-+Italy&entry=gmail&source=g>
>> 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
>> <https://maps.google.com/?q=35+Stirling+Highway,+Crawley+WA+6009,+Australia&entry=gmail&source=g>
>>
>> E-mail address: alcard...@gmail.com <mailto:alcard...@gmail.com>,
>> andrea.card...@unimore.it <mailto:andrea.card...@unimore.it>
>> WEBPAGE: https://sites.google.com/site/alcardini/home/main
>> <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
>>
>> <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/
>> <http://www.footprintnetwork.org/en/index.php/GFN/page/calculators/>
>>
>> --
>> MORPHMET may be accessed via its webpage at http://www.morphometrics.org
>> <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
>> <mailto:morphmet%2bunsubscr...@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 <https://www.researchgate.net/profile/David_Katz29>
>> Personal webpage
>> <https://davidckatz.wordpress.com/>
>>
>> --
>> MORPHMET may be accessed via its webpage at http://www.morphometrics.org
>> <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
>> <mailto: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 <https://www.researchgate.net/profile/David_Katz29>
> Personal webpage
> <https://davidckatz.wordpress.com/>
--
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.
