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
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?
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?