如何在forloop中使用broom：：tidy/stats：：confint从gnm对象中提取结果？

nc1teljy 于 2023-10-13 发布在其他

关注(0)|答案(1)|浏览(123)

我尝试使用以下代码从gnm对象中提取系数：

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

数据集和部分代码来自本文：
https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/1471-2288-14-122#Sec13
数据集：https://static-content.springer.com/esm/art%3A10.1186%2F1471-2288-14-122/MediaObjects/12874_2014_1140_MOESM2_ESM.zip
产品编号：https://static-content.springer.com/esm/art%3A10.1186%2F1471-2288-14-122/MediaObjects/12874_2014_1140_MOESM1_ESM.docx
运行forloop后，弹出以下错误：

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

其中i = 1。
然而，当我一个接一个地运行for循环时：

> 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

我能准确无误地得出结果.有什么我错过的吗？

r

来源：https://stackoverflow.com/questions/77269879/how-to-extract-the-result-from-gnm-object-by-using-broomtidy-statsconfint-wi

1条答案

按热度按时间

uelo1irk1#

伟大的问题;感谢您提供相关细节。我开始研究gnm()的源代码，但无法找出问题的根本原因。我还尝试将for循环转换为函数，并使用lapply()，Map()和purrr::map()，得到了相同的结果。
作为一个潜在的解决方法，也许你可以在循环中运行for (i in 2:12)，然后“手动”添加i = 1的有问题的结果，例如。

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

创建于2023-10-11使用reprex v2.0.2
这能解决你的问题吗？

赞(0）回复(0）举报 2023-10-13

我来回答

如何在forloop中使用broom：：tidy/stats：：confint从gnm对象中提取结果？

1条答案

相关问题

热门标签

最新问答