helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

How to extract the overdispersion parameters in regression models assuming negative binomial distribution #74

Closed ghost closed 1 year ago

ghost commented 1 year ago

I'm working on regression models on count data and am assuming negative binomial distribution. It seems that the fitting process ran successfully, but how can I extract the overdispersion(often called phi?) from the fitSSM object?

helske commented 1 year ago

See my reply in this issue: https://github.com/helske/KFAS/issues/20, the overdispersion parameters are part of model$u.

ghost commented 1 year ago

I think I've misunderstood model$u as the offset term ( I found a non-english website falsely describing model$u as an offset term.) Then, is there a way to add an offset term to the model?

helske commented 1 year ago

Due to poor design choices model$u corresponds to different things depending on the respective distribution, e.g. it's an offset in Poisson, a number of trials in binomial, and a dispersion parameter in negative binomial. If you want to use offset in a negative binomial case you can add log(offset) as an explanatory variable to the model, and set the corresponding initial state distribution to constant N(1, 0). So something like this:

model <- SSModel(
  y ~ -1 + SSMregression(~ x + log_offset, a1 = c(0, 1), P1inf = diag(c(1, 0))), 
  distribution = "negative binomial")
ghost commented 1 year ago

Thank you very much. I'm having difficulty reconstructing my model following your instruction

In my model, I want time-varying intercepts as well as time-varying coefficients for up to 2 variables I used models as below: Intercept term: "-1" in your example to "SSMtrend(1, Q=NA)" Regression term: ~x1(+x2)+log_offset, a1=c(0, (0,) 1), P1inf=diag(1, (1,) 0), Q=diag(c(NA, (NA,) 0) update function: updatefn <- function(pars, model){ model$Q[1,1,] = exp(pars[1]) model$Q[2,2,] = exp(pars[2]) (model$Q[4,4,] = exp(pars[3])) model$u[,] = exp(pars[4]) return(model) } when I run 1 or 2-variable models, I get model$u with like 10^12 order, which seems not reasonable for a dispersion parameter.

Please let me know if I'm doing something wrong here.

helske commented 1 year ago

Hard to know without the actual reproducible example, but I would double-check that the elements of Q and P1inf match with the order of the states (e.g. by checking the state names in model$a1). It could also be due to bad initial values for the optimization which causes the dispersion parameter to wander toward infinity.

ghost commented 1 year ago

Sorry that I cannot give you details about my data. At least I confirmed that the Q and P1 inf match, and that the coefficients for the offset term are fixed to 1 (as expected). At this point, seems like I have to find another way to do this analysis. (At least I've successfully obtained regression results of time-varying regression when using log-transformed "rates" instead of handling data as count data)

Thank you very much for your help anyway.

helske commented 1 year ago

Assuming there is no errors in your code and no bugs in KFAS, then what your results mean is that essentially there is no overdispersion and you could just use Poisson distribution, as such a large dispersion parameter means that NB tends towards Poisson with mean=var. Running a Poisson version of the model should then lead to very similar results.

You could also test the bssm package to see whether the its NB model gives similar results (it also directly supports offsets).