Linearity of dominance hierarchies

Here, I’m posting two sets of codes:
(1) A code to produce a coefficient K, based on work by Kendall, as reported in Appleby (1983).
(2) A code to produce a linearity index, h, and a modification of this index by deVries (1995), often referred to as h’.

For a thorough explanation of what these values mean, please refer to the citations below.

Note that I have

This is a set of codes to implement a test of linearity as described by Appleby (1983). It uses an index called K, based on work on paired comparisons by Kendall.

library(igraph) dat=read.csv(file.choose(),header=TRUE,row.names=1,check.names=FALSE) # read .csv file m=as.matrix(dat) g=graph.adjacency(m,mode="directed",weighted=TRUE,diag=FALSE) E(g)$width=E(g)$weight totint=m+t(m) # This set of codes reproduces Appleby's (1983) test of linearity  # First step is to change dyads that have unknown relationship =0.5 mod=m for (i in 1:nrow(m)){     for (j in 1:nrow(m)){         if (m[i,j]>m[j,i]) mod[i,j]=1             else if (m[i,j]==m[j,i])              mod[i,j]=mod[j,i]=0.5             else              mod[i,j]=0         if (totint[i,j]==0) mod[i,j]=0.5                 }         } diag(mod)=0 N=nrow(m) Si=rowSums(mod) d=(N*(N-1)*(2*N-1)/12)-(0.5*sum(Si^2)) df=(N*(N-1)*(N-2)/((N-4)^2)) chi=(8/(N-4))*((N*(N-1)*(N-2))/24-(d+0.5))+df pAppleby=1-pchisq(chi,df=df) maxd=ifelse(N%%2==1, (N^3-N)/24, (N^3-4*N)/24) K= 1-d/maxd K pAppleby

This is a set of codes to implement the test of linearity of a dominance matrix, as outlined in de Vries (1995). Detailed notes to come later.

library(igraph) dat=read.csv(file.choose(),header=TRUE,row.names=1,check.names=FALSE) # read adjacency matrix from .csv file m=as.matrix(dat) g=graph.adjacency(m,mode="directed",weighted=TRUE,diag=FALSE) E(g)$width=E(g)$weight plot.igraph(g,vertex.label=V(g)$name,layout=layout.fruchterman.reingold, vertex.color="white",edge.color="black",vertex.label.color="black") totint=m+t(m) N=nrow(m) V0=degree(g,mode="out") rawh=(12/((N^3)-N))*sum((V0-((N-1)/2))^2) #This calculates the original Landau's h value ## This set is the modified test of linearity a la de Vries (1995). There are three major steps:  #Step 1) randomly fill in the null relationships such that all individuals either wins (=1), loses (=0) to each individual. Known ties are denoted as 0.5 for both individuals. You then calculate the h-value for this tournament -- this is the h0 value. The "modified Landau's h" as denoted by de Vries (1995) is the mean value of these h0 values in 10,000 randomizations.  #Step 2) Create a completely random tournament in which all individuals either win or lose to each individual. Calculate the h-value for this, which is the hr value.  #Step 3) Compare hr and h0: p-value is the number of times hr is bigger or equal to h0 in 10,000 simulations.  h0=vector(length=10000) hr=vector(length=10000) t=0 for (k in 1:10000){ newmat=m for (i in 1:N){ for (j in 1:N){ if (totint[i,j]>0) if (m[i,j]>m[j,i]) newmat[i,j]=1 else if (m[i,j]==m[j,i])  newmat[i,j]=newmat[j,i]=0.5 else  newmat[i,j]=0 else if (j>i){ newmat[i,j]=sample(c(0,1),1) newmat[j,i]=abs(newmat[i,j]-1)}}} diag(newmat)=0 V=rowSums(newmat) h0[k]=(12/((N^3)-N))*sum((V-((N-1)/2))^2) nm=matrix(nrow=N,ncol=N) for (i in 1:nrow(m)){ for (j in 1:nrow(m)){ if (j>i){ nm[i,j]=sample(c(0,1),1) nm[j,i]=abs(nm[i,j]-1)}}} diag(nm)=0 Vr=rowSums(nm) hr[k]=(12/((N^3)-N))*sum((Vr-((N-1)/2))^2) if (hr[k]>=h0[k]) t=t+1} hmod=mean(h0) p=t/10000 cat(" Landau's h= ",rawh,"\n","modified Landau's h= ",hmod,"\n","p-value from simulations= ",p) hist(hr,xlim=c(0,1),xlab="Landau h values from simulation") abline(v=hmod,lty=3,lwd=1.5)

References:

de Vries, Han. 1995. An improved test of linearity in dominance hierarchies containing unknown or tied relationships. Animal Behavior. 50: 1375-1389.

Share this: