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

```{r}
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

```{r}
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.

```{r}
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.
```{r}
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.
```{r}
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.