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の中には複数のがあり、それらすべてを分析yearに含める必要があります。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 に答える