6

Stan で堅牢なロジスティック回帰 (Robit) を実行したいと考えています。このモデルは、Gelman & Hill の「Data Analysis Using Regression and Multilevel Methods」(2006 年、pp. 124) で提案されていますが、実装方法がわかりません。Stan の Github リポジトリとリファレンス マニュアルを確認しましたが、残念ながらまだ混乱しています。通常のロジスティック回帰をモデル化するために使用したコードを次に示します。たとえば、自由度 7 の分布でエラーが発生するようにするには、何を追加すればよいでしょうか? ひょっとして、マルチレベル解析を行っても同じ手順になるのでしょうか?

library(rstan)

set.seed(1)
x1 <- rnorm(100)  
x2 <- rnorm(100)
z <- 1 + 2*x1 + 3*x2      
pr <- 1/(1+exp(-z))       
y <- rbinom(100,1,pr)  

df <- list(N=100, y=y,x1=x1,x2=x2)

# Stan code
model1 <- '
data {                          
  int<lower=0> N;          
  int<lower=0,upper=1> y[N];  
  vector[N] x1;         
  vector[N] x2;
}
parameters {
  real beta_0;     
  real beta_1;        
  real beta_2; 
}
model {
  y ~ bernoulli_logit(beta_0 + beta_1 * x1 + beta_2 * x2);
}
'
# Run the model
fit <- stan(model_code = model1, data = df, iter = 1000, chains = 4)
print(fit)

ありがとうございました！

4

4 に答える 4

4

Luc Coffeng がStan のメーリング リストでこの回答を送ってくれたので、ここに追加する必要があると思いました。彼は言った：

「Robit 回帰のベースとして GLM を採用します。標準誤差項を に置き換えるだけで、e ~ student_t(7, 0, sigma_e)スケールが適切sigma_e ~ cauchy(0, 2)だと思います ((-5,5) の逆ロジットがほとんどをカバーするため、とにかく 5 を超えることはありません)。 [0,1] 間隔の. t-エラーのスケールに加えて, パラメータとして t-エラーの df を指定することもできます. 推奨されるコードについては以下を参照してください.

ただし、データには、提供したおもちゃの例よりも多くの情報が含まれていることを願っています。つまり、個人ごとに複数の観測があります (以下のように)。個人/ユニットごとに 1 つの観測だけでは、モデルを特定することは事実上不可能です。」

次に、彼は次の例を提供しました。

library(rstan)

set.seed(1)
x1 <- rnorm(100)  
x2 <- rnorm(100)
z <- 1 + 2*x1 + 3*x2 + 0.1 * rt(100, 7)
pr <- 1/(1+exp(-z))       
y <- rbinom(100,10,pr)  

df <- list(N=100, y=y, x1=x1, x2=x2, nu = 7)

# Stan code
model1 <- '
data {                          
   int<lower=0> N;          
   int<lower=0,upper=10> y[N];  
   vector[N] x1;         
   vector[N] x2;
   real nu;
}
parameters {
   real beta_0;     
   real beta_1;        
   real beta_2; 
   real<lower=0> sigma_e;
   vector[N] e;
}
model {
   e ~ student_t(nu, 0, sigma_e);
   sigma_e ~ cauchy(0, 1);
   y ~ binomial_logit(10, beta_0 + beta_1 * x1 + beta_2 * x2 + e);
}
'
# Run the model
fit <- stan(model_code = model1, data = df, iter = 4000, chains = 2)
print(fit)

ボブ・カーペンターもこの質問について簡単にコメントしました。

「[...] はい、階層設定でも同じことができますが、正常に近づくにつれてスケールが無限大になると、自由度のモデリングが難しくなる可能性があるため、注意が必要です。」

Bernd の質問に答えて、Lucy ~ bernoulli_logit(10...はモデル コードに次のように記述した理由を明らかにしました。

「私が提供したサンプル コードでは、10 がサンプル サイズです。おもちゃのデータには、個体 / ユニットごとに複数の観測値 (つまり、ユニットごとに 10 個の観測値) が含まれていることに気付いたかもしれません。

Stan のマニュアルには、関数の引数とサンプリング ステートメントに関する広範な情報も記載されています。」

于 2014-10-21T20:57:09.150 に答える
1

Stan 2.4リファレンスマニュアルの26ページ：

y ~ bernoulli(Phi( beta_0 + beta_1 * x1 + beta_2 * x2 ))

一般的な解決策はy ~ bernoulli(link_function(eta))、link_functionたとえばPhiです。bernoulli_logitたまたま、この機能をラップし、数値的に安定している特別な関数があります。

この理由が明確でない場合は、一般化された線形モデルを読むことをお勧めします。ウィキペディアのページはそれほど悪いレビューではありません。

于 2014-10-21T04:00:29.230 に答える