биномиальная регрессия с групповыми эффектами - PullRequest
Новая
       23

биномиальная регрессия с групповыми эффектами

0 голосов
/ 27 мая 2020

Я пытаюсь построить модель биномиальной регрессии с помощью rstan.

Цель состоит в том, чтобы получить размер эффекта bt между двумя условиями X в группе t в целом и размер эффекта btg для подгрупп g . 

library(rstan)

df <- data.frame(hits=c(36,1261,36,1261,49,1248,17,7670,25,759,29,755),trials=c(118,53850,184,53784,209,53759,118,53850,184,53784,209,53759)
                ,X=rep(c(1,0),6),g=rep(rep(1:3, each=2),2),t=rep(1:2,each=6),tg=rep(1:6,each=2) )

stanIn <- list(Nt=length(unique(df$t)), #number of groups t
            Nc=length(df$t), #number of rows
            Ng=length(unique(df$g)), #number of subgroups g
            Ntg=length(unique(df$tg)), #number of t and g combinations
            N=df$trials, 
            n=df$hits,
            X=df$X, #condition 1 or 0
            t=df$t, #index of groups
            g=df$g, #index of subgroups
            tg=df$tg) #index of combinations between t and g

model <- stan(data=stanIn, file="minimal.stan", chains = 4)

с minimal.stan , как показано ниже. 

data {
  int<lower=1> Nt; 
  int<lower=1> Nc; 
  int<lower=1> Ng; 
  int<lower=1> Ntg; 
  int<lower=1> N[Nc];
  int<lower=0> n[Nc]; 
  int<lower=0,upper=1> X[Nc]; 
  int<lower=1> t[Nc]; 
  int<lower=0> g[Nc]; 
  int<lower=1> tg[Nc]; 
}

parameters {
  real at[Nt]; // group intercepts 
  real bt[Nt]; // group slopes  
  real btg[Ntg]; // subgroup slopes 
  real atg[Ntg]; // subgroup intercept
}

transformed parameters {
  vector[Nc] theta; // binomial probabilities

  for (i in 1:Nc) { // linear model
    theta[i] = inv_logit( (atg[tg[i]]+at[t[i]] ) + (bt[t[i]]+btg[tg[i]]) * X[i]);
    //theta[i] = inv_logit(at[t[i]] + bt[t[i]] * X[i]); //group effect
    //theta[i] = inv_logit(atg[tg[i]] + btg[tg[i]] * X[i]); //subgroup effects
  }
}

model {
  at ~ normal(0.0, 20.0);
  bt ~ normal(0.0, 20.0);  
  atg ~ normal(0.0, 20.0);
  btg ~ normal(0.0, 20.0);
  n ~ binomial(N, theta);
}

Я могу смоделировать общий групповой эффект с помощью первой прокомментированной строки (в преобразованных параметрах ) и эффекты подгруппы со второй строкой комментария. Идея заключалась в том, чтобы объединить как групповой эффект, так и отклонение от него для отдельной группы (первая линия).

Однако это дает действительно странные результаты для bt и btg (A), и я ожидал чего-то более похожего на (B) (я не могу воссоздать поведение видно в A с минимальным примером, это происходит только в полном наборе данных.)

enter image description here

Если это не очевидно из типа вопроса, я совершенно новичок в моделировании statisti c и подозреваю, что у меня есть концептуальная ошибка. Так что я был бы признателен за любой намек по этой проблеме или источник, чтобы прочитать об этих вещах (кажется обычным делом, но я ничего не нашел).

1 Ответ

0 голосов
/ 28 мая 2020

Я вижу несколько проблем. Наиболее прямо модель не идентифицируется, поскольку параметры atg и btg включают at и bt. Я изменил их на ag и bg, поскольку они являются параметрами вашей подгруппы, и проиндексировал их как таковые, как показано ниже: 

library(rstan)

df <- data.frame(hits   = c(36, 1261, 36, 1261, 49, 1248, 
                            17, 7670, 25, 759, 29, 755),
                 trials = c(118, 53850, 184, 53784, 209, 53759, 
                            118, 53850, 184, 53784, 209, 53759),
                 X  = rep(c(1,0), 6),
                 g  = rep(rep(1:3, each=2), 2),
                 t  = rep(1:2, each=6),
                 tg = rep(1:6, each=2) )

stanIn <- list(Nt = length(unique(df$t)), #number of groups t
               Nc = length(df$t),         #number of rows
               Ng = length(unique(df$g)), #number of subgroups g
               N  = df$trials, 
               n  = df$hits,
               X  = df$X,                 #condition 1 or 0
               t  = df$t,                 #index of groups
               g  = df$g)                 #index of subgroups

model <- stan(data = stanIn, file = "minimal.stan", 
              cores = 4, chains = 4,
              control = list(max_treedepth = 14))

с этими модифицированными образцами модели Stan без проблем: 

data {
  int<lower=1> Nt; 
  int<lower=1> Nc; 
  int<lower=1> Ng; 
  int<lower=1> N[Nc];
  int<lower=0> n[Nc]; 
  vector<lower=0,upper=1>[Nc] X; 
  int<lower=1> t[Nc]; 
  int<lower=0> g[Nc]; 
}

parameters {
  vector<offset=0, multiplier=20>[Nt] at; // group intercepts 
  vector<offset=0, multiplier=20>[Nt] bt; // group slopes  
  vector<offset=0, multiplier=20>[Ng] ag; // subgroup intercepts
  vector<offset=0, multiplier=20>[Ng] bg; // subgroup slopes
}

model {
  at ~ normal(0.0, 20.0);
  bt ~ normal(0.0, 20.0);  
  ag ~ normal(0.0, 20.0);
  bg ~ normal(0.0, 20.0);
  n  ~ binomial_logit(N, ag[g] + at[t] + (bt[t] + bg[g]) .* X);
}

generated quantities {
  vector[Nc] theta = inv_logit(ag[g] + at[t] + (bt[t] + bg[g]) .* X);
}

Следует отметить, что мне пришлось использовать высокий max_treedepth, чтобы соответствовать модели, но мне трудно комментировать это, не понимая данных. Я также переместил theta в generated quantities на случай, если вам нужны эти вычисления, но binomial_logit обрабатывает это напрямую.

Я также установил нецентрированные параметры, используя offset и multiplier, поэтому что у стандартного сэмплера лучшее пространство параметров для выборки. Наконец, я перекодировал циклы как векторы.

...