Independently standardize and write

[joannarifkin]

Hi there!

This is a follow up on #36 and on #40. I can use https://github.com/kdm9's extra parallelized code to generate a collection of .Robj results for different runs, but I'm having trouble figuring out how they be loaded to run standardize.xvals and write.xvals.

I've tried iterating through loading the files and adding each to a list and then run standardize.xvals and write.xvals on that, but I don't think it's doing the trick. It sounds like it should be possible, from the text below, but I'm not doing it right.

Any thoughts?

Thanks!

Joanna

          Ok so I _think_ this is salvageable.  Inside the `x.validation` call there are 4 main things that happen:

1. it makes the data partitions
2. it performs all the cross-validation replicates (parallelized, if that was specified)
3. it standardizes the cross-validation results
4. it writes out the output

Your run got truncated part of the way through Step 2 above, but I think you should be able to finish the remaining runs by loading the data.partitions object created in your initial run and then calling the internal function x.validation.rep (using the ::: syntax to call an unexported function; e.g., conStruct:::x.validation.rep). Once you've completed all the runs, I think you should be able finish out Steps 3 and 4 above. Inside x.validation, those steps are accomplished via:

names(x.val) <- paste0("rep_", 1:n.reps)
x.val <- lapply(x.val, standardize.xvals)
save(x.val,file=paste0(prefix,".xval.results.Robj"))
write.xvals(x.val,prefix)

I think you should be able to take all the same steps, again using the ::: syntax to call standardize.xvals and write.xvals, which are unexported functions in the conStruct package.

Let me know if you hit snags along the way - I'm happy to help troubleshoot!

Originally posted by @gbradburd in #40 (comment)

[gbradburd]

hmm yeah this should be doable. can you show the snippet of code where you're attempting this, and, if possible, attach a compressed .Robj file of one of the xval runs from @kdm9's extra-parallelized code?

[joannarifkin]

Totally!

Here are some .Robj files and code.

"Archive.zip" contains spatial and non-spatial results for K=1 and K=2 from the demo data, generated using @kdm9's code. I also have stuff from my own data, but was using the demo data to test all these data wrangling steps because it can be compared to the original.

 `
files <- list.files(pattern="example_reg_.*conStruct.results.Robj", full.names=TRUE, recursive=FALSE) #List files
results_from_files=list() #initiate empty lists to store files and index to files
results_index<-list()
for (f in files){
#print(f)
name<-as.character(f)
#print(strsplit(name, "/|_|rep|K"))
namesplit<-strsplit(name, "/|_|rep|K")
f_rep<-as.numeric(namesplit[[1]][6])
#print(f_rep)
f_K<-as.numeric(namesplit[[1]][7])
#print(f_K)
f_mdl<-namesplit[[1]][4]
#print(f_mdl)
index_list<-list(f_rep, f_mdl, f_K)
names(index_list)<-c("rep", "mdl", "K")
load(f)
print(index_list)
conStruct.results<-standardize.xvals(conStruct.results)
results_index<-append(results_index, list(index_list))
names(conStruct.results)<-paste(f_rep, f_mdl, f_K, sep="_")
results_from_files <- append(results_from_files, conStruct.results)
}`

I've tried a few different versions of "standardize.xvals" and "write.xvals," or their components, on the resulting lists of files, using .Robj files from both the csh and regular version of x.validation, but don't seem to be parsing the input in the right way to average the correct things.

Thanks and let me know what else I can do to help!

Archive.zip

Below: csh.xvals command to generate the example files.

`
 data(conStruct.data)
 my.csh.xvals <- csh.x.validation(train.prop = 0.9,
                          n.reps = 4,
                          K = 1:3,
                          freqs = conStruct.data$allele.frequencies,
                          data.partitions = NULL,
                          geoDist = conStruct.data$geoDist,
                          coords = conStruct.data$coords,
                          prefix = "csh_version",
                          n.iter = 1e2,
                          make.figs = FALSE,
                          save.files = TRUE,
                          parallel = FALSE,
                          n.nodes = NULL)`

[gbradburd]

It looks like the .Robjs in the Archive.zip file were the data.block and conStruct.results files that are saved as a by-product of an x.validation run. The next step in the cross-validation algorithm is to calculate the likelihood of the "testing" data partition given the parameters estimated in each of those model runs. Those likelihoods (one for each model run in that cross-validation replicate, i.e. on a particular testing data partition) need to be standardized for comparison against the same from a different replicate. It looks like @kdm9's csh.x.validation function returns a data.frame with the mean likelihoods (averaged over MCMC iterations) for each model run on each replicate (the mean.fit column).

So, the way I would do this is:

First, define a helper function to turn the tibble output by @kdm9's csh.x.validation into a list in the same format as that output by conStruct::x.validation:

tibble2list <- function(K,n.reps,csh.xvals){
	xvals.list <- vector(mode="list",length=n.reps)
	names(xvals.list) <- paste0("rep_", 1:n.reps)
	for(i in 1:n.reps){
			xvals.list[[i]][["sp"]] <- vector(mode="list",length=K)
				names(xvals.list[[i]][["sp"]]) <- 1:K
			xvals.list[[i]][["nsp"]] <- vector(mode="list",length=K)
				names(xvals.list[[i]][["nsp"]]) <- 1:K
		for(k in 1:K){
			xvals.list[[i]][["sp"]][[k]] <- csh.xvals$mean.fit[csh.xvals$rep==i & csh.xvals$mdl=="sp" & csh.xvals$K==k]
			xvals.list[[i]][["nsp"]][[k]] <- csh.xvals$mean.fit[csh.xvals$rep==i & csh.xvals$mdl=="nsp" & csh.xvals$K==k]
		}
	}
	return(xvals.list)
}

Then make another function for standardizing the output of csh.x.validation:

stdize.csh.xvals <- function(K,n.reps,csh.xvals,prefix){
	xvals.list <- tibble2list(K,n.reps,csh.xvals)
	std.xvals <- lapply(xvals.list,conStruct:::standardize.xvals)
	conStruct:::write.xvals(std.xvals,prefix)
	return(std.xvals)
}

I tried that locally and it worked for me:

data(conStruct.data)
my.csh.xvals <- constructhelpers::csh.x.validation(train.prop = 0.9,
						                          n.reps = 4,
						                          K = 1:3,
						                          freqs = conStruct.data$allele.frequencies,
						                          data.partitions = NULL,
						                          geoDist = conStruct.data$geoDist,
						                          coords = conStruct.data$coords,
						                          prefix = "csh_version",
						                          n.iter = 1e2,
						                          make.figs = FALSE,
						                          save.files = TRUE,
						                          parallel = FALSE,
						                          n.nodes = NULL)
stdize.csh.xvals(K=3,n.reps=4,csh.xvals=my.csh.xvals,prefix="TESTING")

Let me know if it works for you!

[joannarifkin]

Hi Gideon,

Thanks! It works when I run the constructhelpers x.validation and save it, but I'm having trouble figuring out how to load the existing files into a tibble that has the same structure.

I'm guessing the smart money here is just to rerun the parallel versions and save the "my.csh.xvals" object for downstream processing, right? Rather than figure out how to reload it from the .Robj files?

Thanks again!

Joanna

[gbradburd]

Hm yeah if you only have access to the conStruct.results and data.block objects that are output by csh.x.validation (with the save.files argument set to TRUE), then I don't think you can reconstruct the cross-validation scores. So, yeah smart money is probably to re-run the cross-validation and save the my.csh.xvals for downstream steps. Sorry!

[joannarifkin]

Great, thanks! Does having the modelfit objects help at all? I'm thinking saving the object is probably where it's at. Thanks again! Joanna Rifkin PhD they/them Computational biologist

…

On Wed, Feb 14, 2024 at 12:54 PM Gideon ***@***.***> wrote: Hm yeah if you only have access to the conStruct.results and data.block objects that are output by csh.x.validation (with the save.files argument set to TRUE), then I don't think you can reconstruct the cross-validation scores. So, yeah smart money is probably to re-run the cross-validation and save the my.csh.xvals for downstream steps. Sorry! — Reply to this email directly, view it on GitHub <#68 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFPL6CTA7O6XI7DZL2VKUG3YTT24JAVCNFSM6AAAAABC4XKGO2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNBUGMZDEMZXHA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[gbradburd]

If I understand what you're asking, I think you can safely set the save.files argument to FALSE because you don't need the model.fit or conStruct.results or data.block objects that are output as intermediates. Just saving the xvals for standardization should be enough. Is that what you're asking?

[joannarifkin]

Ah, no, I was asking if having the model.fit objects in addition to conStruct.results and data.block made it possible to reconstruct the cross-validation scores, but I think I've answered my own question because the training data are still gone. Got it! Also, if it's helpful I can share my SGE and .r scripts for running this on an HPC. I'm guessing downstream users might want that? Cheers, Joanna

…

On Wed, Feb 14, 2024 at 4:14 PM Gideon ***@***.***> wrote: If I understand what you're asking, I think you can safely set the save.files argument to FALSE because you don't need the model.fit or conStruct.results or data.block objects that are output as intermediates. Just saving the xvals for standardization should be enough. Is that what you're asking? — Reply to this email directly, view it on GitHub <#68 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFPL6CQVEV7RLHB5AEFGPFDYTUSMJAVCNFSM6AAAAABC4XKGO2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNBUGYYTEOBRGI> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Joanna Rifkin PhD they/them Computational biologist

[gbradburd]

yeah sharing the scripts would be great! Thanks Joanna!

[joannarifkin]

Sorry for the delay - here are HPC scripts and the .yml for a conda environment.

This is a wrapper shell script. It's written in SGE but could easily be converted for SLURM

#!/bin/bash
#$ -cwd #launch from current working directory
#$ -V
#$ -M name@email.com
#$ -m bea #email when complete or crash
#$ -q your_queue.q
#$ -pe smp 64 #allocate resources

Job Starts Here

source activate ~/.conda/envs/conStruct

R CMD BATCH --no-save conStruct_parallel.r

This is the R script itself

library(conStruct)
library(parallel)
library(foreach)
library(doParallel)
library(constructhelpers)
setwd("/your/favorite/directory")

freqs<-readRDS("~data_matrix.RData") #a data matrix saved as an RData file
pop_geospheredist<-as.matrix(read.csv("geosphere_dist.csv"), header=F) #distance matrix generated with geosphere
pop_coords<-as.matrix(read.csv("pop_coords.csv"))

cl <- makeCluster(64,type="FORK")
registerDoParallel(cl)

my.xvals <- csh.x.validation(train.prop = 0.8,
n.reps = 8,
K = 1:10,
freqs = freqs,
data.partitions = NULL,
geoDist = pop_geospheredist,
coords = pop_coords,
prefix = "name/outfiles",
n.iter = 1e4,
make.figs = TRUE,
save.files = TRUE,
parallel = TRUE,
n.nodes = 64)

save(my.xvals, file = "xvals.RData")
print("DONE")

stopCluster(cl)

Then, one can simply load xvals.RData (possibly still on an HPC), depending how big it is, run the standardizations above, and optionally also inspect the layer contributions.

Hope that's helpful!

conStruct.yml.txt