How to extract the result from gnm object by using broom::tidy/stats::confint within a forloop?

Question

I tried to use the following code to extract the coefficients from the gnm object:

library(tidyverse)
library(gnm)
library(broom)
library(haven)

data = read_dta('londondataset2002_2006.dta')

data$ozone10 <- data$ozone/10

# GENERATE MONTH AND YEAR
data$month  <- as.factor(months(data$date))
data$year   <- as.factor(format(data$date, format="%Y") )
data$dow    <- as.factor(weekdays(data$date))
data$stratum <- as.factor(data$year:data$month:data$dow)

data <- data[order(data$date),]

# FIT A CONDITIONAL POISSON MODEL WITH A YEAR X MONTH X DOW STRATA
modelcpr = vector(mode = 'list',length = 12)

for(i in 1:12){
  
  
  modelcpr1 <- gnm(numdeaths ~ ozone10 + 
                               lag(temperature,i), data=data, family=poisson,
                               eliminate=factor(stratum))
  
  
  modelcpr[[i]] = broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
    select(estimate, conf.low, conf.high) %>% as.matrix %>% unname
}

#vs
#only for i = 1

modelcpr1 <- gnm(numdeaths ~ ozone10 + 
                             lag(temperature,1), data=data, family=poisson,
                             eliminate=factor(stratum))

#broom::tidy
broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
            select(estimate, conf.low, conf.high) %>% as.matrix %>% unname

The dataset and part of the code are from this paper:

https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/1471-2288-14-122#Sec13

dataset: https://static-content.springer.com/esm/art%3A10.1186%2F1471-2288-14-122/MediaObjects/12874_2014_1140_MOESM2_ESM.zip

code: https://static-content.springer.com/esm/art%3A10.1186%2F1471-2288-14-122/MediaObjects/12874_2014_1140_MOESM1_ESM.docx

---

After running the forloop, the following error popped up:

modelcpr = vector(mode = 'list',length = 12)

for(i in 1:12){
  
  
  modelcpr1 <- gnm(numdeaths ~ ozone10 + 
                               lag(temperature,i), data=data, family=poisson,
                               eliminate=factor(stratum))
  
  
  modelcpr[[i]] = broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
    select(estimate, conf.low, conf.high) %>% as.matrix %>% unname
}

Error in profile.gnm(object, which = parm, alpha = 1 - level, trace = trace) : 
  profiling has found a better solution, so original fit had not converged
In addition: Warning message:
In sqrt((deviance(updated) - fittedDev)/disp) : NaNs produced

where i = 1.

However, when I run the for loop one by one:

> modelcpr1 <- gnm(numdeaths ~ ozone10 + 
+                              lag(temperature,1), data=data, family=poisson,
+                              eliminate=factor(stratum))
> 
> #broom::tidy
> broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
+             select(estimate, conf.low, conf.high) %>% as.matrix %>% unname
         [,1]      [,2]     [,3]
[1,] 1.000446 0.9988817 1.002013

I can obtain the result without any error. Is there any thing that I missed?

Answer 1

Great question; thanks for including the relevant details. I started working through the source code of gnm() but couldn't figure out the root cause of the problem. I also tried converting the for-loop into a function and using lapply(), Map() and purrr::map(), and got the same result.

As a potential workaround, perhaps you could run for (i in 2:12) in the loop and then add the problematic result for i = 1 'manually', e.g.

library(tidyverse)
# install.packages("gnm")
library(gnm)
library(broom)
library(haven)

data = read_dta('/Users/jared/Desktop/londondataset2002_2006.dta')

data$ozone10 <- data$ozone/10

# GENERATE MONTH AND YEAR
data$month  <- as.factor(months(data$date))
data$year   <- as.factor(format(data$date, format="%Y") )
data$dow    <- as.factor(weekdays(data$date))
data$stratum <- as.factor(data$year:data$month:data$dow)

data <- data[order(data$date),]

# FIT A CONDITIONAL POISSON MODEL WITH A YEAR X MONTH X DOW STRATA
modelcpr = list()

for(i in 2:12){
  modelcpr1 <- gnm(numdeaths ~ ozone10 + lag(temperature, i),
                   data=data, family=poisson,
                   eliminate=factor(stratum))
  modelcpr[[i]] <- broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
    select(estimate, conf.low, conf.high) %>% as.matrix %>% unname
}
#> Warning: The `tidy()` method for objects of class `gnm` is not maintained by the broom team, and is only supported through the `glm` tidier method. Please be cautious in interpreting and reporting broom output.
#> 
#> This warning is displayed once per session.

modelcpr1 <- gnm(numdeaths ~ ozone10 + 
                   lag(temperature,1), data=data, family=poisson,
                 eliminate=factor(stratum))

modelcpr[[1]] <- broom::tidy(modelcpr1,conf.int = T, exponentiate = T) %>% slice(n()) %>% 
  select(estimate, conf.low, conf.high) %>% as.matrix %>% unname

modelcpr
#> [[1]]
#>          [,1]      [,2]     [,3]
#> [1,] 1.000446 0.9988817 1.002013
#> 
#> [[2]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9977508 0.9962128 0.9992911
#> 
#> [[3]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9959252 0.9959004 0.9959574
#> 
#> [[4]]
#>           [,1]    [,2]      [,3]
#> [1,] 0.9964251 0.99639 0.9964698
#> 
#> [[5]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9967111 0.9966801 0.9967492
#> 
#> [[6]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9960485 0.9960293 0.9960723
#> 
#> [[7]]
#>           [,1]     [,2]      [,3]
#> [1,] 0.9958503 0.995833 0.9958719
#> 
#> [[8]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9955268 0.9955123 0.9955448
#> 
#> [[9]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9958603 0.9958431 0.9958819
#> 
#> [[10]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9961307 0.9961111 0.9961547
#> 
#> [[11]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9956007 0.9955877 0.9956168
#> 
#> [[12]]
#>           [,1]      [,2]      [,3]
#> [1,] 0.9947129 0.9947042 0.9947233

^{Created on 2023-10-11 with reprex v2.0.2}

Would that solve your problem?