PLN-team / PLNmodels

A collection of Poisson lognormal models for multivariate count data analysis
https://pln-team.github.io/PLNmodels
GNU General Public License v3.0
54 stars 18 forks source link

Questions. #93

Closed brgew closed 1 year ago

brgew commented 1 year ago

Hi,

I have some questions related to PLNmodels v 1.0.0.

o when I run the example given in the vcov.PLNfit() documentation, I get NULL . That is,

> data(trichoptera)
> trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
> myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera)
> vcov(myPLN) ## variance-covariance of B
NULL

    I wonder if this is the expected result. Am I making a mistake?

o it looks like the PLNnetwork_param() config_post options to jackknife and bootstrap exit with errors. I wonder     if I am missing something or making another mistake?

Thank you for the valuable package!

jchiquet commented 1 year ago

Hi @brgew and thank you for your feedback on PLNmodels.

when I run the example given in the vcov.PLNfit() documentation, I get NULL . I wonder if this is the expected result. Am I making a mistake?

Indeed, this is now the expected result : we use to send back the approximated variational variance-covariance estimator, but it highly underestimates the true variance of the parameters and should not be used blindly. However, I agree that this behavior is strange. I will discuss with @mahendra-mariadassou to decide whether we keep the vcov method or purely suppress it.

If you really want your old 'vcov' behavior for the moment, you need,

> data(trichoptera)
> trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
> myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, 
                              control = PLN_param(config_post = list(variational_var = TRUE)))
> vcov(myPLN) # variational variance-covariance of B

To get more trustful estimator of the variance, use jackknife or bootstrap estimation( like you seems to do regarding your second question)

> data(trichoptera)
> trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
> myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, 
                               control = PLN_param(config_post = list(jackknife = TRUE)))
> standard_error(myPLN, "jackknife", "B") #  estimated variance of B

It looks like the PLNnetwork_param() config_post options to jackknife and bootstrap exit with errors. I wonder if I am missing something or making another mistake?

That's a bug, I can reproduce it, thank for pointing it

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLNnetwork(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, 
                                          control = PLNnetwork_param(config_post = list(jackknife = TRUE)))
Error in (function (data, params, config) :
unused arguments (Y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 1, 3, 4, 2, 2, 1, 4, 2, 1, 0, 1, 4, 3, 2, 2, 3, 6, 7, 5, 3, 0, 1, 3, 1, 1, 0, 0, 0, 12, 7, 12, 15, 6, 4, 18, 4, 5, 1, 4, 3, 3, 6, 5, 4, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 8, 32, 176, 69, 14, 4, 29, 8, 2, 3, 3, 40, 20, 500, 142, 44, 31, 30, 71, 20, 5, 8, 2, 10, 9, 14, 2, 4, 1173, 2671, 33, 62, 220, 29, 60, 21, 3, 5, 18, 49, 33, 71, 28, 37, 103, 11, 17, 27, 0, 0, 4, 2, 1, 0, 1, 0, 0, 0, 0, 2, 2, 19, 9, 1, 4, 6, 5, 0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 1, 4, 8, 20, 6, 3, 2, 0, 0, 0, 0, 0, 1, 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,

Will fix this soon on github, sorry about that. We indeed changed many things, and the jackknife/bootstrap thing was not highly tested.

brgew commented 1 year ago

Hi,

I am grateful for the feedback and insights! Thank you again for the valuable package and your efforts.