Error in map_interaction_qtl

[ivlachos]

Hi,

Step 1 finished correctly and then used the kernel to map eQTLs.

The process seems to be passing correctly different milestones (I had to add a few utility functions from dox that were missing) but ends with an error:

[a great number of genes and variants examined, all with normal termination] .... .... Optimization terminated normally: Convergence detected: relative gradient magnitude is below tolerance

Error in { : task 1 failed - "task 1 failed - "invalid 'length' argument"" Calls: map_interaction_qtl -> %do% ->

Any ideas?

[davidaknowles]

Can you let me know which utility functions from dox you needed, would be nice to add them here.

Did you set checkpoint_dir? I would recommend doing that so you can a) check whether it running successfully for genes before the one where it crashes b) restart without having to rerun all the successfully run genes.

R isn't giving a very informative error message because the issue is occurring inside a parallel %dopar% loop. "invalid 'length' argument" is what happens when NULL gets passed to numeric(). The only place numeric is called in map_interaction_qtl is on line 169: init$beta=c(init$beta,numeric(ncol(interact))) Here interact is the design matrix for the full model with interaction, which I suspect ends up being NULL here for some reason. Would you be able to make to make a MWE? (I realize this may be challenging if you can't share the data easily).

[ivlachos]

Hi,

I am sourcing utils.R from dox, since the one(s) that were missing were found there.

It's been quite a while and I can't recall which were missing but I can remove the source call and see if it nags about something.

Based on your response, I'll overload map_interaction_qtl with a local version and trace. The code of the wrapper is quite straightforward. I'm not sure how I can access the stanmodel object from the package to pass it to optimize.

The example is feasible but it might be challenging due to privacy constraints. I could try to create a random dataset with same dimensions and see if I can replicate the error.

Thanks

[davidaknowles]

Ok thanks, I guess I had utils.R from dox sourced when testing suez. I'll work out what needs adding.

Stan based packages access the compiled stan models through stanmodels, which is private... but one of the good/bad things about R is you can still access private package functions/members using :::, so suez:::stanmodels$suez_step_2 is what you can pass to optimizing.

[ivlachos]

I believe the origin of the issue is here:

interact=model.matrix(~geno:condition,data=anno) interact=interact[,3:ncol(interact)]

the model matrix (first call) created is an N*3 matrix (intercept, Genotype:ConditionA, Genotype:ConditionB) The second call reduces the matrix into a vector, so a few lines below:

init$beta=c(init$beta,numeric(ncol(interact)))

It breaks since it requests ncol from a vector.

[davidaknowles]

Ah, the old drop=F issue. So does interact=interact[,3:ncol(interact),drop=F] fix it?

I think the key missing function from utils.R is fix_diag, which I've now added to master.

[ivlachos]

I already added it and and I've overloaded the original one. I'm creating a shorter example to see quickly if it's a fix.