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> 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> 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>
> 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, 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 <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.
>
>
>
--
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.
