Как написать R-код для mcmc логистической модели изменений

Question

Я хочу написать код mcmc для логической модели точки изменения, моя модель выглядит следующим образом:

logit{p(x)}=beta0+beta1*doseA+(beta2-beta1)*(doseA-w)*I(doseA>w) +beta3*doseB+(beta4-beta3)*(doseB-v)*I(doseB>v)

, в которой оцениваются следующие параметры: beta0-beta4,w,v;

x,doseA,doseB предоставлены ли данные

У меня есть код Metroplis Hasting для более простой модели:

logit{p(x)}=beta0+beta1*doseA+beta2*doseB

, но я не знаю, как обращаться с параметрами точки изменения в моей модели, интересно,они должны быть отобраны вместе с другими параметрами или разделены.Кроме того, модель до и после точки изменения имеет разное количество параметров, поэтому они используют разные априоры, но мне интересно, стоит ли разделять мои данные для этого.

Я пытался использовать rstan для mcmc, но этобыл очень медленнымЯ надеюсь, что смогу написать код mcmc сам или использовать несколько пакетов, таких как "mcmc" или "HI"

# dose-toxicity probability model
#  input: doseA and doseB are vectors with same size, 
#        doseA[i] and doseB[i] comprise a dose regimen
#  output: p is a vector, p[i] is the prob. of doseA[i] and doseB[i] 
probfun = function(doseA,doseB,param){
  p = param[1]+param[2]*doseA+param[3]*doseB
  p = exp(p)
  p = 1-1/(1+p)
  return(p)
}

# log-likelihood of toxicity probability
#used by function logpostT
loglike = function(doseA,doseB,param,x,n){
  mu = param[1]+param[2]*doseA+param[3]*doseB
  lk = sum(x*mu-n*log(1+exp(mu)))       
  return(lk)
}

# Cauchy prior distribution
prior_c = function(b0,mub0,scaleb0){
  dcauchy(b0,mub0,scaleb0)
}

# posterior core for toxicity
#used by function mhT , to compute parameter posterior 
# input:doseA,doseB,x,n are input data, vectors with same size,
#       muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2 are hyperparameters for prior 
logpostT = function(doseA,doseB,param,x,n,muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2){
  h = loglike(doseA,doseB,param,x,n)+log(prior_c(param[1],muT0,scaleT0))+
    log(prior_c(param[2],shapeT1,rateT1))+log(prior_c(param[3],shapeT2,rateT2))
  return(h)
}
# candidate(proposal) density for toxicity
#used by function mhT , to provide new parameters for a new iteration
candiT = function(param,jsdb0,jsdb1,jsdb2){  
  a = rnorm(1,param[1],jsdb0)
  b = rnorm(1,param[2],jsdb1)
  c = rnorm(1,param[3],jsdb2)
  return(c(a,b,c))
}

# acceptance ratio in MCMC for toxicity
#used by function mhT , to compute r in Metropolis-Hasting algorithm, likelihood_new/likelihood_old*prior_new/prior_old
ratioT = function(doseA,doseB,x,n,b,bs,muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2){
  r = exp(sum(x*((bs[1]-b[1])+doseA*(bs[2]-b[2])+doseB*(bs[3]-b[3])))-       
            sum(n*(log(1+exp(bs[1]+doseA*bs[2]+doseB*bs[3]))-log(1+exp(b[1]+doseA*b[2]+doseB*b[3])))))*
    prior_c(bs[1],muT0,scaleT0)/prior_c(b[1],muT0,scaleT0)*
    prior_c(bs[2],shapeT1,rateT1)/prior_c(b[2],shapeT1,rateT1)*prior_c(bs[3],shapeT2,rateT2)/prior_c(b[3],shapeT2,rateT2)
  if(is.na(r)==T) r=Inf
  return(r)
}

# Metropolis-Hasting algorithm for toxicity
mhT = function(doseA,doseB,x,n,muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2,jsdb0,jsdb1,jsdb2,size,burn,thin){

  start = c(-5,5,5)

  b_pre = start
  n_var = length(b_pre)
  ac = matrix(0,size,n_var)
  draws = matrix(NA,size,n_var)
  post = c()       

  for (i in 1:size) {
    b_cand = candiT(b_pre,jsdb0,jsdb1,jsdb2)
    b_old = b_pre

    #sequential updated every parameter in each iteration
    for (j in 1:n_var){
      #give the b_star based on the b_old and b_cand
      if (j==1) {
        b_star = c(b_cand[j],b_old[(j+1):n_var])
      } else if (j==n_var) {
        b_star = c(b_old[1:(j-1)],b_cand[j])
      } else {
        b_star = c(b_old[1:(j-1)],b_cand[j],b_old[(j+1):n_var])
      }
      #calculate ratio and determine whether to accept or not
      r = ratioT(doseA,doseB,x,n,b_old,b_star,muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2)
      if (runif(1)<r) {
        b_new = b_star
        ac[i,j]=1
      }else {
        b_new = b_old
        ac[i,j]=0
      }
      b_old = b_new 
    }

    draws[i,] = b_new
    post[i] = logpostT(doseA,doseB,draws[i,],x,n,muT0,scaleT0,shapeT1,rateT1,shapeT2,rateT2)
    b_pre = b_new
  }

  return(list(draws=draws[seq(floor(size*burn)+1,size,thin),],post=post[seq(floor(size*burn)+1,size,thin)],ac=apply(ac,2,sum)))
}