レキシスモデルから粗発生率テーブルを作成する

Question

研究のフォローアップデータをモデル化するために、Rの「Epi」パッケージを使用しています。 私はLexisモデルを宣言したり、Poissonと（生存パッケージと組み合わせて）Cox回帰を実行することに問題はありません。

最初のデータレビューの一部として、Rのレキシスモデル（任意のポアソン/コックスモデルにあらかじめ適合する）のデータから粗調整されていない発生率/イベントレートのテーブルを作成する簡単な方法を見つける必要があります。 - のない出力

#Generic Syntax Example 
total <-cbind(tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum),tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum)) 
    #Add up the number of events within the stratifying variable 
    #Add up the amount of follow-up time within the stratifying the variable 
    rates <- tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum)/tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum)*10^3 
    #Given rates per 1,000 person years 
    ratetable <- (cbind(totals,rates)) 

#Specific Example based on the dataset 
totals <-cbind(tapply(lexis_model$lex.Xst,lexis_model$grade,sum),tapply(lexis_model$lex.dur,lexis_model$grade,sum)) 
rates <- tapply(lexis_model$lex.Xst,lexis_model$grade,sum)/tapply(lexis_model$lex.dur,lexis_model$grade,sum)*10^3 
ratetable <- (cbind(totals,rates)) 
ratetable 

       rates 
1 90 20338.234 4.4251630 
2 64 7265.065 8.8092811 
#Shows number of events, years follow-up, number of events per 1000 years follow-up, stratified by the stratifying variable 

注この粗製未調整/絶対速度：

私は私はこれを実行し、探索データ解析の一部として変数により層別化することを可能にする符号化されたアプローチを発見しましたポアソンモデル。上記のコードが実際に望ましい出力を生成することは分かりましたが、私は人々がレキシスデータセットを取ってこれを出力できるコマンドを認識しているかどうかを見たいと思っていました。私はEpiとepitoolsパッケージで利用可能なコマンドを見てきました。何も見逃しているかもしれませんが、これを行うための明白な方法は見当たりませんでした。

これは非常に一般的なことですが、単純にレキシスデータセットと階層化変数を指定することで、誰かがパッケージ/関数を認識できるかどうか疑問に思っていました（実際には、一回で上に）。

理想的には出力は以下のようなものになります（私はRの賛成から離れて移動しようとしていますSTATAから取られている！）：

最初の20行のコピーまたはのように実際のデータは（データが既にエピパッケージを使用して語彙モデルにして置かれているので、すべての関連語彙の変数がある）はこちら： https://www.dropbox.com/s/yjyz1kzusysz941/rate_table_data.xlsx?dl=0

Answer 1

私は単にのようなtidyverse Rパッケージを使用してこれを行うだろう：

library(tidyverse) 
lexis_model %>% 
    group_by(grade) %>% 
    summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>% 
    mutate(rate = sum_Xst/sum_dur*10^3) -> rateable 
rateable 

# A tibble: 2 x 4 
# grade sum_Xst sum_dur  rate 
# <dbl> <int>  <dbl> <dbl> 
# 1  1  2 375.24709 5.329821 
# 2  2  0 92.44079 0.000000 

そして、あなたは、関数自身にこれをラップすることができます：

rateFunc <- function(data, strat_var) 
{ 
    lexis_model %>% 
    group_by_(strat_var) %>% 
    summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>% 
    mutate(rate = sum_Xst/sum_dur*10^3) 
} 

あなたがして呼び出します：tidyversesummariseとmutateの組み合わせを使用すると、それは非常に簡単ですので、

rateFunc(lexis_model, "grade") 

これは便利ですテーブルに要約統計量を追加してください。

EDIT： 質問に明確化した後、これはrateコマンドを使用してpopEpiパッケージ使用して行うことができます。

popEpi::rate(lexis_model, obs = lex.Xst, pyrs = lex.dur, print = grade) 

# Crude rates and 95% confidence intervals: 
# grade lex.Xst lex.dur  rate  SE.rate  rate.lo rate.hi 
# 1:  1  2 375.2472 0.00532982 0.003768752 0.001332942 0.0213115 
# 2:  2  0 92.4408 0.00000000 0.000000000 0.000000000  NaN