1
$\begingroup$

For purposes of validating a forecast model, I'd like to compare a GLARMA model that I developed to a null model that includes the same autocorrelation effects but lacks the environmental data. GLARMA are like ARMA, but they allow for count data as the independent variable.

Here's a reference for this model type.

Dunsmuir, William T. M. and Scott, David J. (2015) The glarma Package for Observation-Driven Time Series Regression of Counts. Journal of Statistical Software, 67(7), 1–36. http://dx.doi.org/10.18637/jss.v067.i07

library(glarma)
library(tidyverse)
data(Polio)
y <- Polio[, 2]
X <- as.matrix(Polio[, 3:8])
glarmamod <- glarma(y, X, thetaLags = c(1,2,5), type = "Poi", method = "FS",
                    residuals = "Pearson", maxit = 100, grad = 1e-6)
glarmamod

Call: glarma(y = y, X = X, type = "Poi", method = "FS", residuals = "Pearson", 
    thetaLags = c(1, 2, 5), maxit = 100, grad = 1e-06)

GLARMA Coefficients:
  theta_1    theta_2    theta_5  
0.2184597  0.1272311  0.0872861  

Linear Model Coefficients:
       Intcpt          Trend      CosAnnual      SinAnnual  CosSemiAnnual  SinSemiAnnual  
    0.1299754     -3.9283714     -0.0991262     -0.5308445      0.2111276     -0.3932302  

Degrees of Freedom: 167 Total (i.e. Null);  159 Residual
Null Deviance: 343.0004 
Residual Deviance: 250.6179 
AIC: 536.7052

X_null <- X[,1, drop = FALSE]
glarmamod_null <- glarma(y, X_null, thetaLags = c(1,2,5), type = "Poi", method = "FS",
                    residuals = "Pearson", maxit = 100, grad = 1e-6)
glarmamod_null

Call: glarma(y = y, X = X_null, type = "Poi", method = "FS", residuals = "Pearson", 
    thetaLags = c(1, 2, 5), maxit = 100, grad = 1e-06)

GLARMA Coefficients:
  theta_1    theta_2    theta_5  
0.2700613  0.1629822  0.1112386  

Linear Model Coefficients:
   Intcpt  
0.1579249  

Degrees of Freedom: 167 Total (i.e. Null);  164 Residual
Null Deviance: 343.0004 
Residual Deviance: 280.207 
AIC: 555.5405

However, I'm struggling to do this with my actual data and model.

events <- c(2, 0, 1, 1, 2, 2, 0, 2, 0, 0, 0, 2, 2, 2, 1, 0, 0, 1, 1, 1, 
1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 
0, 0, 0, 1, 0, 3, 2, 0, 1, 2, 3, 1, 2, 1, 0, 0, 2, 2, 2, 2, 0, 
2, 2, 1, 0, 2, 1, 2, 0, 2, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 17, 9, 5, 3, 2, 3, 
6, 3, 3, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 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, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 
0, 0, 2, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 
1, 0, 0, 0, 0, 0, 2, 0, 1, 0, 0, 0)

concentration <- c(-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -1.2, -1.8, -0.8, 
-3.9, -3.9, -0.7, -3.9, -2.1, -1.3, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -2.3, -3.9, 
-2.8, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -2.7, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-2.8, -3.9, -2.4, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -1.6, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -2.5, -1.3, 0.1, 0.8, 2, 
-2.8, -2.4, 0, -0.3, -0.7, 1.4, -1.5, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -1.6, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, 
-3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9, -3.9)

events_ts <- ts(events)
concentration_mtx <- cbind(intercept = 1, concentration) %>% as.matrix()

My model with data works fine.

glarma(
  y = events_ts,
  X = concentration_mtx,
  phiLags = NULL,
  thetaLags = c(1,6)
)

Call: glarma(y = events_ts, X = concentration_mtx, phiLags = NULL, 
    thetaLags = c(1, 6))

GLARMA Coefficients:
  theta_1    theta_6  
0.1301675  0.1773473  

Linear Model Coefficients:
 intercept      concentration  
-0.9438270   0.5214789  

Degrees of Freedom: 261 Total (i.e. Null);  258 Residual
Null Deviance: 452.4841 
Residual Deviance: 224.4038 
AIC: 435.2015

However, the null equivalent fails with an error.

concentration_mtx_null <- concentration_mtx[,1, drop = FALSE]

glarma(
  y = events_ts,
  X = concentration_mtx_null,
  phiLags = NULL,
  thetaLags = c(1, 6)
)

Error in glarma(y = events_ts, X = concentration_mtx_null, phiLags = NULL, : Fisher Scoring fails to converge from the initial estimates.

Why does my model not converge in the absence of environmental data while the built in polio model does?

Is there any way I can get a null glarma model for my data that is the same as my working model but lacks environmental data?

Thanks for any advice.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
0
$\begingroup$

The help file for this function posts a similar scenario.

## Score Type (GAS)  Residuals, Newton Raphson
## Note: Newton Raphson fails to converge from GLM initial estimates.
## Setting up the initial estimates by ourselves
init.delta <- glarmamod$delta
beta <- init.delta[1:6]
thetaInit <- init.delta[7:9]

glarmamod <- glarma(y, X, beta = beta, thetaLags = c(1, 2, 5),
                    thetaInit = thetaInit, type ="Poi", method = "NR",
                    residuals = "Score", maxit = 100, grad = 1e-6)

Thus, they advocate setting up initial estimates for theta. In the example above, one can copy thetaInit from the theta estimates initial non-null function.

Thus running:

glarma(
  y = events_ts,
  X = concentration_mtx_null,
  phiLags = NULL,
  thetaLags = c(1, 6),
  thetaInit = c(0.07373976 , 0.15219858  )
)

Works as expected.

Call: glarma(y = events_ts, X = concentration_mtx_null, phiLags = NULL, 
    thetaLags = c(1, 6), thetaInit = c(0.07373976, 0.15219858))

GLARMA Coefficients:
  theta_1    theta_6  
0.1627414  0.1072906  

Linear Model Coefficients:
 intercept  
-0.6770703  

Degrees of Freedom: 261 Total (i.e. Null);  259 Residual
Null Deviance: 452.4841 
Residual Deviance: 684.3455 
AIC: 561.5263
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.