これは、stackoverflow での最初の投稿です。ガイドラインに従って質問してみました。
私は現在限界構造モデルに取り組んでいます。このトピックに関する Hernan の記事のほとんどを読み、彼の著書What ifと提供された R コードに取り組みました。ただし、私が見つけたすべての例はかなり単純です。たとえば、患者は追跡調査の開始時に治療 A を受けるか受けないかのいずれかであり、その間ずっとその治療を続けます。学んだことを自分のデータセットに適用するのは難しいと思います。
私のデータセットは、イベントまでの時間に関するデータです。フォローアップ期間中、患者は治療 A を受けることができます。一部の患者はフォローアップの開始時にすでに治療 A を受けていますが、他の患者はフォローアップの後半に開始するか、まったく受けていません。さらに、治療 A を受けている患者は、経過観察中に治療を中止し、治療を再開することができます。つまり、時間変動治療です。
ここでは、イベントまでの時間データの大まかな MSM のコードを含むデータセットの例を示します。Cox PH 回帰の代わりに、オッズ比が Cox モデルのハザード比に類似しているプールされたロジスティック回帰を使用します。簡単にするために、打ち切りの重みは計算しませんでした。
set.seed(1948)
idvar<-1:100 #100 subjects
v1<-sample(c(0,1), replace=TRUE, size=length(idvar)) #One baseline variable
time<-sample(1:100, replace=TRUE,size=length(idvar)) #Follow-up time
event<-sample(c(0,1), replace=TRUE, size=length(idvar)) #End of follow-up event. If 0 patients either reached administrative censoring or a competing event
df<-data.frame(idvar,v1,time,event) #Create dataframe
library(splitstackshape)
df.expand<-expandRows(df,"time", drop = F) #Create dataframe with 1 row per 1 unit of time per subject
df.expand$time_1 <- sequence(rle(df.expand$idvar)$lengths)-1 #Calculates times in expanded dataframe
df.expand$event_true <- as.factor(ifelse(df.expand$time_1==df.expand$time-1 & #Only event at last follow-up
df.expand$event==1, 1, 0))
df.expand$treatment<-sample(c(0,1), replace=TRUE,size=nrow(df.expand)) #Time-varying treatment
df.expand$v2<-sample(c(0,1), replace=TRUE,size=nrow(df.expand)) #Time-varying covariate
#Calculate weights. Need help here
pred_dnom = predict(glm(treatment~factor(v1) + factor(v2), #Time-varying covariates only in denominator
data = df.expand,
family = binomial(link = "logit")), type = "response")
pred_num = predict(glm(treatment~factor(v1),
data = df.expand,
family = binomial(link = "logit")), type = "response")
df.expand$sw <- ifelse(df.expand$treatment == 0, ((1 - pred_num) / (1 - pred_dnom)),
(pred_num / pred_dnom))
##Normally I would also calculate censoring weights and multiply those with the treatment weights##
fit1<-glm(event_true==1~treatment+factor(time_1),
data = df.expand, weights = sw,
family = binomial(link = "logit"))
ただし、上記の例では、重みは以前の治療履歴と時間によって変化する交絡因子の治療履歴を考慮していません。だから私の質問は、これらの歴史を時変治療にどのように組み込むかということです. Fewellらによるこの記事。患者が治療を中止して再開することができない場合に、時変治療のためにそうする方法を説明しています。私は引用します：
各被験者の各月までの完全な治療履歴の確率を推定するために ( 3の分母)、各月に観察された治療の推定確率を時間の経過とともに累積的に乗算します。各被験者の最初の推定確率はそのまま残されます。それ以外の場合は、現時点での推定確率に前の時点での推定確率を掛けます。[Fewell Z, Hernán MA, Wolfe F, Tilling K, Choi H, Sterne JAC. 限界構造モデルを使用した時間依存交絡の制御。スタタジャーナル。2004;4(4):402-420。ドイ:10.1177/1536867X0400400403]
彼らの記事のコードに見られるように、患者が治療を開始すると、重みは 1 に設定されます。多くの人が重みを計算するために使用するもう 1 つの方法は、ipw パッケージです。ただし、このパッケージは、患者が治療を開始したときにも重みを 1 に設定し、その後も重みは 1 のままです。
ipw パッケージの作成者の 1 人も、Grafféo などの記事に寄稿しました。これは、ipw パッケージを拡張して、時間によって変化する重みを許可します。彼らは、治療を受けていないときに治療を開始するための重みと、治療中に治療を中止するための重みを計算することによってこれを行います。(記事のサポート情報として利用可能なコード) しかし、彼らは治療の歴史/時変交絡因子を考慮していませんでした.
より正確には、交絡治療曝露の全履歴を考慮したのではなく、時間 t での対象の曝露は、時間 t で行われた交絡治療のみに依存していると考えただけでした。次に、関心のある露出が変更された最初の時点までの重みを計算する代わりに、すべての時点、つまり、関心ステータスの処理のすべての変更とすべてのイベント時間について重みが計算されます。[Grafféo, N, Latouche, A, Geskus, RB, Chevret, S. 治療加重の逆確率を使用した時変曝露のモデル化。生体認証ジャーナル。2018; 60: 323–332。https://doi.org/10.1002/bimj.201600223]
彼らは次のことを提案します
交絡因子の歴史は考慮しませんでした。これは、関心のある曝露と交絡因子との関係におけるスイッチごとに異なるモデルを使用して対処できます。このモデル化オプションを研究するには、さらなる作業が必要です。最も簡単な方法は、以前の暴露のカウンターを共変量として治療モデルに含めることです。[Grafféo, N, Latouche, A, Geskus, RB, Chevret, S. 治療加重の逆確率を使用した時変曝露のモデル化。生体認証ジャーナル。2018; 60: 323–332。https://doi.org/10.1002/bimj.201600223]
そのため、Fewell による方法には、治療と時変変数の履歴が組み込まれていますが、治療開始後に重みを変更することはできません。Grafféo による方法では、重みを経時的に変化させることができますが、治療や時間によって変化する変数の履歴は組み込まれていません。したがって、私が実際に望んでいるのは、これらの方法を組み合わせることですが、それがどのように可能になるかは絶対にわかりません.
私の質問が明確で（そして正しく尋ねられて）、誰かが提案してくれることを願っています。ありがとう！