bayes_fosr with gibbs returns "fosr" string instead of a model

[kastman]

Hi all,

For some reason, using bayes_fosr with est.method = 'VB' is returning a model just fine, but with est.method = 'Gibbs' I only get the string "fosr" as the return value, e.g. mod <- bayes_fosr(func ~ preds, est.method = 'Gibbs'); mod == 'fosr'. Additionally, Gibbs was working previously on this same system and I can't think of anything that has changed on it; I reinstalled refund but have no luck.

Any suggestions on troubleshooting? Any logfiles I can add? Sorry, I'm not familiar with the guts of the package so I'm not sure how to go deeper. Thanks,

[jeff-goldsmith]

Thanks for noting the bug, although I haven't been able to recreate it (and don't see anything obvious in the code that would lead to this kind of issue). I tried running a modified model from the help file and it worked as expected:

library(refund)

data(DTI)
Gibbs = bayes_fosr(cca ~ pasat, data = DTI, Kt = 10,
  est.method = "Gibbs", N.iter = 500, N.burn = 200)

If you can check whether the above example works for you and / or post an example that has the behavior you describe, that would be really helpful.

[kastman]

I just checked, and the above example returns a model correctly while my own model still fails, so it must be something to do with my model specification. Using identical parameters to the test, I get:

Error in gibbs_mult_fpca(formula = formula, data = data, ...) : 
  unused argument (kt = 10)

This suggested that the actual estimation function being used successfully wasn't gibbs_mult_fpca. Looking at the switch statement in bayes_fosr, it appears that my code including the random effect is using gibbs_mult_fpca while the example is using gibbs_cs_fpca. I looked at the DTI example, but couldn't find an appropriate combination to try to make the DTI estimation break with a random effect glancing at the example quickly (I got various indexing errors, shape collapsing errors, etc.)... So I suspect that the problem is somewhere in gibbs_mult_fpca, though I'm not sure where specifically.

Thanks again for your help with this. Is it worthwhile to gist the data? Here is my environment, running in RStudio:

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin16.1.0 (64-bit)
Running under: macOS  10.13.3

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] fda.usc_1.3.0  rpart_4.1-10   mgcv_1.8-15    nlme_3.1-128   MASS_7.3-45    fda_2.4.7      Matrix_1.2-7.1
[8] refund_0.1-16 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.13     magic_1.5-6      munsell_0.4.3    colorspace_1.3-2 lattice_0.20-34  grpreg_3.1-2    
 [7] rlang_0.1.2      minqa_1.2.4      plyr_1.8.4       tools_3.3.2      parallel_3.3.2   grid_3.3.2      
[13] gtable_0.2.0     gamm4_0.2-5      lme4_1.1-13      lazyeval_0.2.0   tibble_1.3.4     nloptr_1.0.4    
[19] ggplot2_2.2.1    RLRsim_3.1-3     pbs_1.1          scales_0.5.0     boot_1.3-18     

[julia-wrobel]

@kastman just from looking at the error you're getting, is it possible in your model specification you are using the argument 'kt = 10' rather than 'Kt = 10'?

[jeff-goldsmith]

Adding to this, although the DTI dataset doesn't have a random effects structure, you could try adding a subject ID at least for the purposes of troubleshooting:

data(DTI)
DTI = subset(DTI, select = c(cca, case, pasat))
DTI = DTI[complete.cases(DTI),]
DTI$subj = sort(sample(1:20, size = 334, replace = TRUE))

fit = bayes_fosr(cca ~ pasat + re(subj), data = DTI, Kt = 10)

The code above seems to run and produce reasonable results, at least.

If you're willing to gist the data / make the full example public that could be helpful.

[kastman]

@julia-wrobel Good catch - you're right that I was using the Kt arg incorrectly. Thanks!

@jeff-goldsmith Unfortunately, the minimal example w/ DTI does work for me, so suspect this must be about my data. The sampler does complete until the end (i.e., a higher value of N.iter does increase the time spent sampling, so I get the feeling that the chain is being sampled correctly and something just isn't being passed off right or returning early).

My use-case is a functional analysis of hand-gripper data, with subjects gripping to different conditions. Here's a gist with a basic breaking model. I do hope to release the whole entirety of the dataset as soon as it's published, but this is a small sample of real grip-courses with the conditions and subjects shuffled. https://gist.github.com/kastman/305542adf248a763ffcbb42756fb3755 Thanks in advance for looking at what I could possibly be doing wrong!

[jeff-goldsmith]

Thanks for creating a gist, @kastman ! I took a look, and ran into exactly the same issue using the data you posted. Oddly, running the function after saving the version in the developer branch to the local environment worked just fine -- still not exactly sure why, since I don't think there are major differences between the devel and master branch versions.

In any case, can you try installing the developer version using the code below and see if that fixes your problem?

# install.packages("devtools")
library(devtools)
devtools::install_github("refunders/refund", ref="devel")

[kastman]

Thanks for getting back to me so fast, @jeff-goldsmith. It looks like that worked, whatever it is! Thanks again for the prompt response.

I'm closing this issue for now - just for the future that was somewhere between the release version 0.1-16 and commit #33a6b4c.

[kastman]

Hi @jeff-goldsmith and @julia-wrobel ,

Apologies for dragging up the issue again, for something that's only tangential; it just doesn't quite deserve its own issue.

I'm trying to run my model again after re-registration, which had been working successfully with no problems. For some reason any time I include an interaction term my design matrix is now singular: Error in solve.default(as.numeric((A + IJ * D/2)/(b.q.sigma.me)) * sig.X + : system is computationally singular: reciprocal condition number = 1E-20, although none of the columns I'm including are highly correlated. Do you have any advice on how to look into why the interaction is breaking the design matrix? I'm totally flumoxed as to what has changed. This is the case both with installing HEAD from github with devtools and also 0.16. Thanks in advance for any advice,