0

我有一个名为的数据集bjmd，看起来像这样（简化）：

      rte   year   y  obs
22037 46001  1     0   1
22042 46001  2     4   3
22047 46001  3     5   3
22202 46002  1    11   1
22207 46002  2    14   1
22212 46002  3     6   1
22140 46003  1     5   6
22141 46003  2     2   6
22142 46003  3     6   6

我想运行一个循环来glm对每个不同的rte（46001、46002、46003）进行分析。在每个rte中，有多个years，它们都需要包含在glm分析中。从每条路线的glm测试中，我正在使用坡度并创建另一个表格，其中路线和坡度作为列。这就是我想要的样子：

rte    slope
46001   x
46002   y
46003   z

这是我想出的for循环代码：

route<-with(bjmd,unique(rte))
slope<-with(bjmd,numeric(length(unique(rte))))
table<-data.frame(route,slope)
for (i in unique(as.factor(bjmd$rte))) {
  data<-subset(bjmd, rte=='i')
  slope[i] <- coef(summary(glm(y ~  year+obs,
                               family = poisson(link=log),data=data)))[2,1]
  table[i,2] <-paste(slope[i])
})
table

这段代码有问题，因为我的斜率不断得到 0 值：

  route slope
1 46001     0
2 46002     0
3 46003     0

有人可以帮忙指出我在哪里搞砸了吗？

4

1 回答 1

1

不需要循环；只需split根据 . 将您的数据集分成组rte。然后用 为每个组拟合一个模型lapply。

lapply(split(bjmd, bjmd$rte), function(dat) glm(y ~ year + obs, data=dat))

您还可以使用交互项一次性对所有内容进行建模。预测值将相同，但残差偏差 df 以及因此的 P 值将不同。哪种方法更适合您的需求取决于您的项目。

glm(y ~ (year + obs) * factor(rte), data=bjmd)
于 2013-07-06T15:48:39.760 回答