same Log Likelihood of different SSModel

[JonBrommer]

Hi

I find KFAS returns the same log-likelihood for different SSModel objects applied to the same data but assuming different observation processes. I think below is a reproducible example that illustrates the issue.

library(KFAS)
library(MASS)
#generate data
n.sites=5 #number of sites
T.steps=28 #time steps
n<-rep(log(20),T.steps) #starting value for state (true) population mean
sigma_p<-.3 #process variance
r.true=.03 #population growth rate on log scale
r<-rnorm((T.steps-1),mean=r.true,sigma_p) #stochastic pgr includes process variance
for (t in 2:T.steps) { # process on log scale
  n[t]<-n[t-1]  + r[t-1]
}
y=NULL
for (t in 1:T.steps) { #observation on datascale
  y[t]<-sum(rnegbin(n.sites,exp(n[t]),3)) #summed up (total) with theta = 3
}
# --- SSM settings
Zt <- matrix(c(1, 0), 1, 2)
Tt <- matrix(c(1, 0, 1, 1), 2, 2)
Rt <- matrix(c(1, 0), 2, 1)
Qt <- matrix(NA)
P1inf <- diag(2)
#poisson
model_poisson <- SSModel(y ~ -1 +
                           SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, P1inf = P1inf),
                         distribution = "poisson", u = n.sites)
#----------- NB version
model_nb <- SSModel(y ~ -1 +
                      SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, P1inf = P1inf),
                    distribution = "negative binomial", u = n.sites)
#log likelihood makes no sense. always returns same value
logLik(model_poisson)
logLik(model_nb)

[helske]

That log-likelihood value is -.Machine$double.xmax^0.75 i.e. basically -infinity. This is bevause you have NA value in Qt:

is.SSModel(model_nb, na.check=TRUE, return.logical=FALSE)
Error in is.SSModel(model_nb, na.check = TRUE, return.logical = FALSE) : 
  System matrices (excluding Z) contain NA or infinite values, covariance matrices contain values larger than 1e+07

[JonBrommer]

Right. For the Gaussian SSM the logL is returned when fitting the model, but not for non-Gaussian. A sort of work around (continuing from the previous example) would be to estimate Qt and re-run the model with that, which produces log likelihood values as Qt != NA . Does that make sense (for this simple case) or a dumb idea?

#use fitted parameters to infer logL
#--poisson
fit_poisson <- fitSSM(model_poisson, inits = c(0, 0), method = "BFGS")
out_poisson <- KFS(fit_poisson$model)
#--nb
fit_nb <- fitSSM(model_nb, inits = c(0, 0), method = "BFGS")
out_nb <- KFS(fit_nb$model)
#re-run model with estimated parameters
model_poisson <- SSModel(y ~ -1 +
                      SSMcustom(Z = Zt, T = Tt, R = Rt, Q = out_poisson$model$Q, P1inf = P1inf),
                    distribution = "negative binomial", u = n.sites)
model_nb <- SSModel(y ~ -1 +
                      SSMcustom(Z = Zt, T = Tt, R = Rt, Q = out_nb$model$Q, P1inf = P1inf),
                    distribution = "negative binomial", u = n.sites)
logLik(model_nb)-logLik(model_poisson)

[helske]

What do you mean that the log-likelihood is not returned for non-Gaussian? From the fitSSM the final log-likelihood value is in -fit_poisson$optim.out$value, or you can call logLik(fit_poisson$model) as you do with logLik(fit_poisson$model). So there's no need to recreate the model with SSModel as it is already in fit_poisson$model.

[JonBrommer]

Thanks for your reply. I got confused because for the gaussian case the $model is not needed to get the logL. Below example that illustrates the different calls to get the logL in gaussian and poisson. Much obliged

data("alcohol")
deaths <- window(alcohol[, 2], end = 2007)
population <- window(alcohol[, 6], end = 2007)
Zt <- matrix(c(1, 0), 1, 2)
Ht <- matrix(NA)
Tt <- matrix(c(1, 0, 1, 1), 2, 2)
Rt <- matrix(c(1, 0), 2, 1)
Qt <- matrix(NA)
a1 <- matrix(0, 2, 1)
P1 <- matrix(0, 2, 2)
P1inf <- diag(2)
model_gaussian <- SSModel(deaths / population ~ -1 +
                                 + SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1,
                                             + P1inf = P1inf), H = Ht)
model_poisson <- SSModel(deaths ~ -1 +
                           SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, P1inf = P1inf),
                         distribution = "poisson", u = population)
fit_gaussian <- fitSSM(model_gaussian, inits = c(0, 0), method = "BFGS")
fit_poisson <- fitSSM(model_poisson, inits = c(0, 0), method = "BFGS")
logLik(fit_gaussian) #returns logL
logLik(fit_poisson) #producs error
logLik(fit_poisson$model) #slightly different syntax, returns logL
logLik(fit_gaussian$model) #works also

[helske]

Hmm, that is strange, if I run your code I get an error from logLik(fit_gaussian) as well as I would expect as fit_gaussian is just a list for which there is no logLik method.

[JonBrommer]

You are right. I started from clean slate and I now understand the correct call to logLik() includes the fitted model $model . Thanks again for clarification