• file: matrix_tau.r
• document page: 31

code:

#----read data------------------------------------------------------------------
Data <- read.table("D:/Heidi/Master/R_Daten/SummaryStatistics/tau.txt", skip=5, fill=TRUE)

Columns <- ncol(Data)+1
Row <- nrow(Data)

x <- 6
n <- 2

#----read data line by line-----------------------------------------------------
tauMatrix <- as.matrix(scan("D:/Heidi/Master/R_Daten/SummaryStatistics/tau.txt",
                             what=double(0), skip=x, nlines=1, nmax=n))
tauMatrix <- cbind(t(tauMatrix), matrix(NA, ncol=(Columns-n), nrow=1))
tauMatrix <- tauMatrix[,2:Columns]

n <- n + 1
x <- x + 1

while(n<(Row+1)){
  nextrow <- as.matrix(scan("D:/Heidi/Master/R_Daten/SummaryStatistics/tau.txt",
                       what=double(0), skip=x, nlines=1, nmax=n))
  nextrow <- cbind(t(nextrow), matrix(NA, ncol=(Columns-n), nrow=1))
  nextrow <- nextrow[,2:Columns]

  tauMatrix <- rbind(tauMatrix, nextrow)

  n <- n + 1
  x <- x + 1
}

a <- ncol(tauMatrix)
b <- nrow(tauMatrix)

x <- c(1:a)
y <- c(1:b)

#----draw plot------------------------------------------------------------------
library(fields)
ColorRamp <- rgb( seq(1,0,length=256),  # Red
                  seq(1,0,length=256),  # Green
                  seq(1,1,length=256))  # Blue

#----Mirror matrix (left-right)----
mirror.matrix <- function(x) {
  xx <- as.data.frame(x);
  xx <- rev(xx);
  xx <- as.matrix(xx);
  xx;
}

#----Rotate matrix 270 clockworks----
rotate270.matrix <- function(x) {
  mirror.matrix(t(x))
}

tauMatrix <- rotate270.matrix(tauMatrix)

#----draw matrix plot----
image.plot(x,y, tauMatrix, col=ColorRamp, main="Divergence times allowing for unequal
              population sizes (tau)", xlab="Population", ylab="Population",
              legend.args=list( text="tau", cex=1.0, side=4, line=2.2), axes = FALSE)
       contour(tauMatrix, add = TRUE)
       axis(1, at = c(1:a))
       axis(2, at = c(1:b), labels=c(b:1))
       box()