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

[irenerlover]

I had posted the same question in Stackoverflow.
https://stackoverflow.com/questions/77269879/how-to-extract-the-result-from-gnm-object-by-using-broomtidy-statsconfint-wi/77270019?noredirect=1#comment136222463_77270019

jared_mamrot https://stackoverflow.com/users/12957340/jared-mamrot and I had attempted different methods, including
for loop,lapply(), Map() and purrr::map(), and still could not extract the result from broom::tidy.

I know that broom::tidy.gnm just directly uses the profile.gnm to obtain the estimate and confidence intervals from the gnm object. jared_mamrot suspects that the error related to how the function handles the args.

Here is the sample dataset and the code I used:

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?