Converting state matrices to adjacency matrices to thresholded matrices

Question

I am currently taking part in research utilizing Naked Mole Rats. I am utilizing data analysis techniques to answer the overall question, "Can we detect certain types of naked mole rat behavior utilizing statistical techniques such as hamming distance and similarity graphs?".

Certain statistical techniques, such as hamming distance have been utilized to detect sleeping patterns, however, that is a pretty blatant type of behavior in my opinion. Whereas, analyzing particular behavior after the birth of a new pup is much harder to detect.

My main programming language is Java and I am new to R. I have written two scripts that are particularly important to the research. The first one, which I will post here, reads in a state matrix with 59 columns and an undetermined amount of rows per file.

Each column signifies a rat, each row signifies one time window, and each data point represents the location the animal is found at

One row of the data looks like this written on one line (an example of a file can be found here:

2   6   3   1   5   2   4   5   7   1   1   1   1   1   2   1   2   1   2   3   

1   2   2   2   1   1   2   1   1   3   2   2   2   1   2   2   1   1   2   3   

7   2   2   2   2   2   1   2   1   2   2   1   1   2   2   4   1   3   4

The program below takes in this matrix, creates an adjacency matrix for the specific file (ie a matrix 59 x 59) that counts how many time windows rat 1 has been with rat 2, and how many time windows rat 1 has been with rat 3, etc. Additionally, the program then creates a thresholded matrix to see which animals spend more then, for example, 50% percent of their time together.

#This version of code works with multiple data files in a particular       folder, 
#performs some matrix calculations on them and outputs the new matrices with their corresponding titles 
#into that same folder - currently there are 165 stateFiles within this particular folder which 
#means 330 additional files were created in that same folder in about a minutes worth of time. 

stateFiles <- list.files("/home/zackymo/Desktop/birth_data_zach_calcs_75/", pattern="*state.txt", full.names=TRUE)
thres <- .75

for(file in stateFiles){
  rawData = read.table(file)
  df <- data.frame(rawData)
  incidence <- matrix(rep(0, ncol(df)*ncol(df)), nrow=ncol(df))
  thresholded <- matrix(rep(0, ncol(df)*ncol(df)), nrow=ncol(df)) 

  #Set the diagonal = to constants due to the fact that both are symmetric matrices 
  diag(incidence) <- nrow(df)
  diag(thresholded) <- 1

  #These loops turn the state matrix into a boolean matrix, and then counts up the true and falses, adds them to the incidence matrix, and then creates the thresholded matrix. 
  for (i in 1:(ncol(df)-1)) {
    for (j in (i+1):ncol(df)) {
      incidence[i,j] = incidence[j,i] = sum(df[,i] == df[,j])
      if(incidence[i,j]/nrow(df) >= thres){
        thresholded[i,j] = thresholded[j,i] = 1
      }#should these conditions be reversed for effiency? How can I check that? 
      else{
        thresholded[i,j] = thresholded[j,i] = 0
      }
    }
  }
  #Data preparation for output file - transposing matrices 
  tincidence <- t(incidence) 
  tthresholded <- t(thresholded)

  #Manipulating the fileName to be able to create a new file w/ a meaningful file name
  #I may be able to split by / then edit the paste below f[1] + f[2] +   f[3] + f[4] + f[5]
  f <- strsplit(file, split='.', fixed=TRUE)
  f <- unlist(f)

  #This will not work -> because you end up 
  #Setting the current directory for the output files to be placed, and never 
  #setwd("/home/zackymo/Desktop/birth_data1_out_files")

  #Creating the various output files 
 write(tincidence, file = paste(f[1],"toInc", sep = ""),
        ncolumns = ncol(df),
        sep = " ") 

  write(tthresholded, file = paste(f[1],"toThresholded", sep = ""),
        ncolumns = ncol(df),
        sep = " ") 
}

I was wondering if someone would be able to critique this code in terms of following R programming standards, along with the general logic and efficiency behind the program.

Answer 1

I'll give it a shot, as I read through from the top. I'll skip the comments but keep them, this is an excellent habit.

You are doing well identifying your globals and leaving them at the top like you did:

stateFiles <- list.files("/home/zackymo/Desktop/birth_data_zach_calcs_75/", pattern="*state.txt", full.names=TRUE)
thres <- .75

I would only have put the input directory in its own variable workingDir so you have more direct access to it if you decide to change it later. Also, you could have chosen to setwd(workingDir) and used full.names = FALSE within list.files.

rawData = read.table(file)
df <- data.frame(rawData)

You are using both the <- and = operators for assignments. I'd recommend you pick the one you prefer and stick to it for consistency (for info, the community standard is to use <-).

read.table returns a data.frame so df <- data.frame(rawData) was not needed. Your data is only numeric (even integer) so you could have chosen to store it into a matrix rather than a data.frame: data.matrix(read.table(...)). Matrices can make many operations faster, though not the ones you are using here.

incidence <- matrix(rep(0, ncol(df)*ncol(df)), nrow=ncol(df))
thresholded <- matrix(rep(0, ncol(df)*ncol(df)), nrow=ncol(df)) 

You have written ncol(df) six times in two lines... Store it in a well named variable (maybe numRats?). Same goes for nrow(df). Also, rep() was not needed, you could just use R's recycling rule:

 incidence <- matrix(0, nrow = numRats, ncol = numRats)

Now you are doing:

 diag(incidence) <- nrow(df)
 diag(thresholded) <- 1

Have a look at ?diag, you could have initialized you two matrices directly as follows:

incidence <- diag(nrow(df), ncol(df))
threshold <- diag(1, ncol(df))

Now looking at:

  if(incidence[i,j]/nrow(df) >= thres){
    thresholded[i,j] = thresholded[j,i] = 1
  }#should these conditions be reversed for effiency? How can I check that? 
  else{
    thresholded[i,j] = thresholded[j,i] = 0
  }

Notice how you are repeating code by writing thresholded[i,j] = thresholded[j,i] = a couple times? You could have done:

thresholded[i,j] = thresholded[j,i] = if (incidence[i,j]/nrow(df) >= thres) 1 else 0

You could even have let R do the conversion from logical to integer:

thresholded[i,j] = thresholded[j,i] = incidence[i,j]/nrow(df) >= thres

Regardless, let's take a step back. The construction of the thresholded matrix is so simple once you have incidence that you should just do the following all the way at the end:

thresholded <- as.integer(incidence >= thres * nrow(df))

(I would even encourage you to drop the as.integer and keep a matrix of booleans). Also, since this is so easy to compute on the fly for any given threshold value, I would discourage you from storing the output in a file. You would only need to save the incidence matrix. This way, you will save disk space and read/write computation times.

Next,

tincidence <- t(incidence) 
tthresholded <- t(thresholded)

seems completely unnecessary since the two matrices are symmetric by construction.

Next, your string manipulation for the output filenames seems a bit complicated. If I get it correctly, you just want to change the extension from txt to toInc? If so, I think you should use sub() as follows:

outFile <- sub("txt$", "toInc", file)

Last, stay away from write for writing your output to a file. It is a really old function; the fact you have to provide ncol as an input tells me it was designed with on-screen printing in mind. Instead, you could use write.table to be consistent with your using of read.table at the beginning of your script:

write.table(incidence, file = outFile, row.names = FALSE, col.names = FALSE)

Last, I'd like to take an extra step back and offer alternatives to the way you computed the incidence matrix. Your code is rather optimal, despite a double for loop. I think you can make it simpler and maybe more robust, at only small costs to computation times.

First alternative: instead of treating the diagonal separately, include it to the for loops by changing the indices to:

for (i in 1:ncol(df)) {
    for (j in i:ncol(df)) {

This also makes your code robust to the special case where df would have a single column (I let you see how your code would break in that case...)

The second alternative would be to replace the double for loop with a call to outer. It is a little slower but see how short and elegant your code becomes:

dat <- data.matrix(read.table(file))
n <- ncol(dat)
incidence <- outer(1:n, 1:n, FUN = function(i, j) colSums(dat[, i] == dat[, j]))
thresholded <- incidence >= thres * nrow(dat)

and done! I hope it helps. Good luck with your project.