search

jomulder / BFpack

BFpack can be used for testing statistical hypotheses using the Bayes factor, a Bayesian criterion originally developed by sir Harold Jeffreys.
https://bfpack.info/
14 stars 4 forks source link

Bf.default #6

Closed Jaeoc closed 5 years ago

Jaeoc commented 5 years ago

Hi Joris,

Here's the function, let me know if this was what you had in mind.

Best, Anton

cjvanlissa commented 5 years ago

Hi Anton! Could you please clarify what cases this function is intended to cover?

Jaeoc commented 5 years ago

Hi Caspar,

Joris asked me to write a function similar to bain.default (https://github.com/cjvanlissa/bain/blob/master/R/bain_methods.R) but with:

This is the bare bones version of what I understood this entailed.

cjvanlissa commented 5 years ago

Dear Anton, thank you for clarifying. I meant: What type of input is the BF.default function expecting? In bain, this is a named vector, named covariance matrix, sample size, number of group parameters and number of joint parameters. It is an interface that higher-order functions like bain.lm call after extracting the relevant information from the lm object.

Jaeoc commented 5 years ago

Ah, I see, thanks for clarifying. I see what you mean now that I look at the bain.lm.

Joris and I were originally discussing to write a function more like that that could be applied to all models that would use gaussian measures (correlations, glm, glmer), but this is not currently how at least the correlation function is written and Joris didn't want to rewrite it right now (I don't think the other functions exist yet).

The 'default' S3 methods are to my understanding in general intended to be applied to any object of a class for which there is no other relevant methods. The BF.default method thus takes the same basic input as any BF method should [according to the draft paper Joris send around]: a fitted model object (x), a character vector (hypothesis) to feed to parse_hypothesis(), and a vector of priors. It then attempts to extract data using lm extractor functions (e.g., coef()) and applies the gaussian approach to evaluate the hypotheses.

I have tested it so far only with glm objects, so the method could easily be called just BF.glm. I am not sure what other classes Joris expects it to work with.

jomulder commented 5 years ago

Hi Anton, thanks for this. Before I accept this pull request, could you check whether the function is also able to get an object of the class "lavaan" (using the cfa function or sem function in the lavaan package), "coxph" (using the coxph function from the survival package), and an object of class "glmerMod" using the glmer function in the lme4 package?

cjvanlissa commented 5 years ago

I think this might be going wrong then, considering there is already a BF.lm, and glm is an lm object. Thus, BF.default will never get called for a glm object. Same for a lavaan object; that will just trigger BF.lavaan.

jomulder commented 5 years ago

Ai indeed a glm also has a lm as second class. I did not yet check this. The point is that BF.lm can only handle a "regular" lm object, i.e., a non-glm object. Any ideas how to deal with this?

On Wed, Jun 5, 2019 at 4:43 PM C. J. van Lissa notifications@github.com wrote:

I think this might be going wrong then, considering there is already a BF.lm, and glm is an lm object. Thus, BF.default will never get called for a glm object. Same for a lavaan object; that will just trigger BF.lavaan.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/jomulder/BFpack/pull/6?email_source=notifications&email_token=AFZBLCA2RDCVEN7CK5J6IBDPY7GHVA5CNFSM4HSIQVA2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODW753KY#issuecomment-499113387, or mute the thread https://github.com/notifications/unsubscribe-auth/AFZBLCHPE55QYPTTCMG5O6DPY7GHVANCNFSM4HSIQVAQ .

jomulder commented 5 years ago

Btw there is no BF.lavaan function (yet), at least as far as I know.

On Wed, Jun 5, 2019 at 4:47 PM Joris Mulder jomulder@gmail.com wrote:

Ai indeed a glm also has a lm as second class. I did not yet check this. The point is that BF.lm can only handle a "regular" lm object, i.e., a non-glm object. Any ideas how to deal with this?

On Wed, Jun 5, 2019 at 4:43 PM C. J. van Lissa notifications@github.com wrote:

I think this might be going wrong then, considering there is already a BF.lm, and glm is an lm object. Thus, BF.default will never get called for a glm object. Same for a lavaan object; that will just trigger BF.lavaan.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/jomulder/BFpack/pull/6?email_source=notifications&email_token=AFZBLCA2RDCVEN7CK5J6IBDPY7GHVA5CNFSM4HSIQVA2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODW753KY#issuecomment-499113387, or mute the thread https://github.com/notifications/unsubscribe-auth/AFZBLCHPE55QYPTTCMG5O6DPY7GHVANCNFSM4HSIQVAQ .

cjvanlissa commented 5 years ago

The point is that BF.lm can only handle a "regular" lm object, i.e., a non-glm object. Any ideas how to deal with this?

Yes - methods are called for classes in order. So an object of class c("glm", "lm") will trigger BF.glm. If that doesn't exist, it will trigger BF.lm. If that doesn't exist, it will trigger BF.default.

So the solution is to make a function for the higher-order class (glm).

Jaeoc commented 5 years ago

I think this might be going wrong then, considering there is already a BF.lm, and glm is an lm object. Thus, BF.default will never get called for a glm object.

Good point Caspar, missed that.

Yes, a simple solution would be to just convert this function to an internal function instead (basically rename it to something like gaussian_estimator()) and then pass this to separate higher order methods for the classes glm, lavaan, coxph, glmerMOD.

So essentially something like:

BF.glm <- function(x, hypothesis, prior, ...){

gaussian_estimator(x, hypothesis, prior, ...)

}

Possibly with different model extraction between class methods if the same functions don't work (e.g., coef() does not seem to quite work the same way on glmerMOD objects as on glm).

I can probably do that tomorrow. Maybe it would be preferable (safer) to not have a default method, so that the BF() function simply throws an error when it receives an unrecognized class?

jomulder commented 5 years ago

Sounds like a good idea that should be fine at least for now. So let's go for that! Note that gaussian_estimator() does the same thing as the bain function except that gaussian_estimator() is more stable as it just calls dmvnorm and pmvnorm while the bain function sometimes gives slightly different outcomes because a sampling algorithm in Fortran resulting in some Monte Carlo errors. For the arguments of gaussian_estimator(), we could do the same as for bain.default (Anton, this can be found on Caspar's GIthub page) which requires a named mean vector, posterior covariance matrix, and sample size. These elements will then extracted using different functions depending on the class of the model object. I will talk next week with Herbert about merging bain and BFpack. Cheers guys, Joris

On Wed, Jun 5, 2019 at 8:14 PM Jaeoc notifications@github.com wrote:

Good point Caspar, missed that.

Yes, a simple solution would be to just convert this function to an internal function instead (basically rename it to something like gaussian_estimator()) and then pass this to separate higher order methods for the classes glm, lavaan, coxph, glmerMOD.

So essentially something like:

BF.glm <- function(x, hypothesis, prior, ...){

gaussian_estimator(x, hypothesis, prior, ...)

}

Possibly with different model extraction between class methods if the same functions don't work (e.g., coef() does not seem to quite work the same way on glmerMOD objects as on glm).

I can probably do that tomorrow. Maybe it would be preferable (safer) to not have a default method, so that the BF() function simply throws an error when it receives an unrecognized class?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/jomulder/BFpack/pull/6?email_source=notifications&email_token=AFZBLCBH5OK7PMGRAKBQUHLPY766XA5CNFSM4HSIQVA2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODXASARQ#issuecomment-499195974, or mute the thread https://github.com/notifications/unsubscribe-auth/AFZBLCAZ5OJTUIBRQMCWO23PY766XANCNFSM4HSIQVAQ .

cjvanlissa commented 5 years ago

That sounds good! I think we're now on the same page :)

Jaeoc commented 5 years ago

I have now added the methods for the above classes. Seems to work pretty well overall according to my quick testing, but some things to note:

All methods basically take a fitted model object (x), extract a named mean vector, the covariance matrix and the number of observations. These are then passed together with any hypotheses and priors to the internal Gaussian_estimator() function.

cjvanlissa commented 5 years ago

Can you update your version of bain to my development version using install_github("cjvanlissa/bain")? I believe it should parse hypotheses for lavaan parameters.

Jaeoc commented 5 years ago

Just tried, but still getting parsing errors:


Error: Could not parse the names of the 'estimates' supplied to parse_hypothesis(). Estimate names must start with a letter or period (.), and can be a combination of letters, digits, period and underscore (_).
The estimates violating these rules were originally named: 'visual=~x2', 'visual=~x3', 'textual=~x5', 'textual=~x6', 'speed=~x8', 'speed=~x9', 'x1~~x1', 'x2~~x2', 'x3~~x3', 'x4~~x4', 'x5~~x5', 'x6~~x6', 'x7~~x7', 'x8~~x8', 'x9~~x9', 'visual~~visual', 'textual~~textual', 'speed~~speed', 'visual~~textual', 'visual~~speed', 'textual~~speed'.
After parsing, these parameters are named: 'visual~~visual', 'textual~~textual'
cjvanlissa commented 5 years ago

Weird. I fixed it somewhere; after the 13th I'll find the code for you. In the meanwhile, you can maybe build your functions and just gsub the names for the time being?

jomulder commented 5 years ago

Hi both, sounds good to have something that at least works for now. Let me now when I can/should merge the pull request. Cheers, Joris