3

私の質問は、Stata のワイブル回帰によって推定された係数から導出された遷移確率の標準偏差 (SD) の計算に関するものです。

移行確率は、白血病患者の疾患進行を 90 日間 (約 10 年) の 40 サイクルにわたってモデル化するために使用されています。対応するマルコフ サイクル確率とその SD を使用してパラメーターを近似できるベータ分布を作成するには、確率の SD (マルコフ モデルの実行中に変化する) が必要です。次に、これらの分布を使用して確率的感度分析を行います。つまり、単純な確率 (サイクルごとに 1 つ) の代わりに使用され、そこからランダムに引き出されて、モデルの費用対効果の結果の堅牢性を評価できます。

とにかく、イベント生存データまでの時間を使用して、回帰分析を使用して、遷移確率を生成する方程式に差し込むことができる係数を推定しました。例えば...

. streg, nohr dist(weibull)

        failure _d:  event
   analysis time _t:  time

Fitting constant-only model:

Iteration 0:   log likelihood = -171.82384
Iteration 1:   log likelihood = -158.78902
Iteration 2:   log likelihood = -158.64499
Iteration 3:   log likelihood = -158.64497
Iteration 4:   log likelihood = -158.64497

Fitting full model:
Iteration 0:   log likelihood = -158.64497  

Weibull regression -- log relative-hazard form 

No. of subjects =           93                     Number of obs   =        93
No. of failures =           62
Time at risk    =        60250
                                                   LR chi2(0)      =     -0.00
Log likelihood  =   -158.64497                     Prob > chi2     =         .

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------
-------------+------
       _cons |  -4.307123   .4483219    -9.61   0.000    -5.185818   -3.428429
-------------+----------------------------------------------------------
-------------+------
       /ln_p |  -.4638212   .1020754    -4.54   0.000    -.6638854    -.263757
-------------+----------------------------------------------------------
-------------+------
           p |    .628876   .0641928                      .5148471    .7681602
         1/p |   1.590139   .1623141                      1.301812    1.942324

次に、p と _cons、時間 (つまり、マルコフ サイクル数) の t、サイクルの長さの u (通常は 1 年、私は白血病患者と仕事をしているので、私の場合は 90 日です) を使用する方程式 () で確率を作成します。再発または死亡などのイベントを起こす可能性が非常に高い人）。

したがって、ラムダ = p、ガンマ = (exp(_cons))

gen result = (exp((lambda*((t-u)^ (gamma)))-(lambda*(t^(gamma)))))

gen transitions = 1-result

変動性に目を向けると、まず係数の標準誤差を計算します

. nlcom (exp(_b[_cons])) (exp(_b[/ln_p]))

       _nl_1:  exp(_b[_cons])
       _nl_2:  exp(_b[/ln_p])

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------
-------------+------
       _nl_1 |   .0116539   .0044932     2.59   0.009     .0028474    .0204604
       _nl_2 |   .6153864    .054186    11.36   0.000     .5091838     .721589

しかし、私が本当に求めているのは、遷移値の標準誤差です。たとえば、

nlcom (_b[transitions])

しかし、これは機能せず、私が使用している本には、この追加情報を取得するためのヒントがありません. 近づく方法についてのフィードバックは大歓迎です。

4

2 に答える 2

1

この回答は正しくありませんが、有用なコメントがいくつか生成されたのでそのままにしておきます

sysuse auto, clear gen u = 90 +rnormal()

セット シード 1234 キャプチャ プログラム ドロップ _all

program define myprog , rclass tempvar result reg turn disp /* ここで -streg- ステートメントを置き換えます */ gen result' = _b[disp]*u sumresult' return scalar sd = r(sd) end

ブートストラップ sdr = r(sd): myprog estat ブートストラップ、bc パーセンタイル

注: ブートストラップされたプログラムでは、新しい変数 (結果) を一時的に定義する必要があります。そうしないと、ブートストラップ複製ごとに変数が新たに作成されるため、gen ステートメントでエラーが発生します。

于 2014-03-11T02:23:21.520 に答える
1

2 番目の回答: 2014-03-23

更新: 2014-03-26負の確率を修正しました: Emily のコードの転記でエラーが発生しました。また、Austin Nichols ( http://www.stata.com/statalist/archive/2014-03/msg00002.htmlnlcom )によって Statalist で提案されているように、それを示しています。Austin のコードに 1 つの修正を加えました。

ブートストラップは依然として解決の鍵です。ターゲット量は、 から推定されたパラメーターの非線形の組み合わせに基づく式によって計算された確率stregです。の後に返される行列には推定値が含まれていないため、 は標準誤差を推定しません。e(b)stregnlcomこれは、ブートストラップにとって理想的な状況です。標準的なアプローチが採用さmyprogれています。パラメーターを推定するプログラムを作成します。それからbootstrapそのプログラム。

この例では、t 値の範囲の遷移確率 pt が推定されます。ユーザーは、範囲の最小値と最大値、およびtスカラーを設定する必要がありuます。おそらく興味深いのは、推定されたパラメータの数が可変であるため、foreach内にステートメントが必要なことですmyprog。また、bootstrapによって返される推定値のリストで構成される引数が必要myprogです。このリストもforeachループで構築されます。

/* set u  and minimum and maximum times here */
scalar u = 1
local tmin = 1
local tmax = 3

set linesize 80

capture program drop _all

program define myprog , rclass
syntax anything
streg , nohr  dist(weibull)
scalar lambda = exp(_b[ln_p:_cons])
scalar gamma =exp(_b[_t:_cons])

forvalues t = `1'/`2'{
scalar p`t'= 1 -  ///
(exp((lambda*((`t'-u)^(gamma)))-(lambda*(`t'^(gamma)))))
return scalar p`t' = p`t'
}
end


webuse cancer, clear
stset studytime, fail(died)
set seed 450811

/* set up list of returned probabilities for bootstrap */
forvalues t = `tmin'/`tmax' {
local p`t' = "p" + string(`t')
local rp`t'= "`p`t''" + "=" +  "("+ "r(" + "`p`t''" +"))"
local rlist =  `"`rlist' `rp`t''"'
}

bootstrap `rlist', nodots: myprog `tmin' `tmax'
forvalues t = `tmin'/`tmax' {
qui streg, nohr dist(weibull)
nlcom 1 - ///
(exp((exp(_b[ln_p:_cons])*((`t'-u)^(exp(_b[_t:_cons]))))- ///
(exp(_b[ln_p:_cons])*(`t'^(exp(_b[_t:_cons]))))))
}

結果： 

Bootstrap results                               Number of obs      =        48
                                                Replications       =        50
      command:  myprog 1 3
           p1:  r(p1)
           p2:  r(p2)
           p3:  r(p3)

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          p1 |   .7009447   .0503893    13.91   0.000     .6021834    .7997059
          p2 |   .0187127    .007727     2.42   0.015     .0035681    .0338573
          p3 |   .0111243   .0047095     2.36   0.018     .0018939    .0203548
------------------------------------------------------------------------------

/* results of  nlcom */
------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |   .7009447   .0543671    12.89   0.000      .594387    .8075023
-------------+----------------------------------------------------------------
       _nl_1 |   .0187127   .0082077     2.28   0.023     .0026259    .0347995
-------------+----------------------------------------------------------------
       _nl_1 |   .0111243   .0049765     2.24   0.025     .0013706    .0208781
------------------------------------------------------------------------------
于 2014-03-22T20:52:04.890 に答える