update_RJ: update_RJ

In BAREB: A Bayesian Repulsive Biclustering Model for Periodontal Data

update_RJ

Description

Update the number of site level clusters for each patient level cluster via RJMCMC

Usage

update_RJ(
  w,
  K,
  Gamma,
  Beta,
  E,
  Z,
  X,
  R,
  mu,
  mu_star,
  Y,
  delta,
  c,
  sigma_square,
  C,
  S,
  theta,
  tau,
  q,
  m,
  T0,
  hyper_delta = 1
)

Arguments

w	The weights for site level clusters. Should be a matrix with S rows.
K	The number of site level cluster for each patient level cluster. Should be a vector of length S
Gamma	The linear coefficients for site level covariates. Should be a 3-dimensional array.
Beta	The linear coefficients for patient level covariates. Should be a matrix with S rows
E	A vector that records the current clustering membership.
Z	The design matrix for site level clusters.
X	The design matrix for patient level clusters.
R	The current site level clustering membership
mu	The current CAL mean matrix
mu_star	The latent value matrix for missingness model
Y	The observed CAL value matrix
delta	The missing indicator matrix
c	The linear coefficients for missingness model
sigma_square	The current noise variance
C	The DPP related kernel matrices. Should be an array of 3 dimensions
S	The number of patient level clusters.
theta	The DPP hyper-parameter for site level
tau	A fixed DPP hyper-parameter, which we suggest high value, say 10^5
q	The number of site level covariates
m	The number of sites
T0	The number of teeth
hyper_delta	The hyper-parameter with default value being 1

Value

A list with following updated parameters:

K	The numbers of site level clusters
w	The weights for site level clusters
Gamma	Linear coefficients for site level covariates
R	The site level clusteirng membership
C	The DPP related kernel matrice

Examples

library(BAREB)
data("obs")
X <- obs$X
Y <- obs$Y
Z <- obs$Z
delta <- obs$delta
data("truth")

set.seed(1)

n<-80
m<-168
T0<-28
q<-3
p<-3
S<-3
theta1 <- theta2 <- 5
tau <- 100000
D<-matrix(0, nrow = T0, ncol = m)
for(i in 1:T0){
  indi<-1:6
  indi<-indi+6*(i-1)
  D[i,indi]<-rep(1/6,6)
}
nu_gamma<-0.05
nu_beta <- 0.05
data("Init")
Beta0 <- Init$Beta
Gamma0 <- Init$Gamma

Niter<-10
record <- NULL
record$E<-matrix(NA,nrow = Niter, ncol = n)
record$R<-array(NA, dim = c(Niter, S, m))
record$Gamma <-array(NA,dim = c(Niter, 10, q, S))
record$Beta <- array(NA,dim = c(Niter, S, p))
record$K <- matrix(NA,nrow = Niter, ncol = S)
record$sigma_square <-rep(NA,Niter)
record$theta1<-rep(0,Niter)
record$theta2<-rep(0,Niter)
record$c<-matrix(0,nrow =  Niter,ncol = T0)
record$mu<-array(NA,dim = c(Niter,n,m))
record$w_beta <- array(NA, dim = c(Niter, S))
record$w <- array(NA, dim = c(Niter, S, 10))
set.seed(1)
E<-sample.int(S,n,TRUE)
Beta <- matrix(NA,nrow = S, ncol = p)
Beta[1,] <- Beta[2,] <- Beta[3,]  <- Beta0
Beta<- Beta + matrix(rnorm(S*p, 0, 1) , nrow = S, ncol = p)
Gamma <- array(NA,dim = c(10,q,S))
Gamma[1,,1] <- Gamma[2,,1] <- Gamma[3,,1] <- Gamma0
Gamma[,,2] <- Gamma[,,3]   <- Gamma[,,1]
Gamma <- Gamma + array(rnorm(10*q*S, 0, 5), dim = c(10, q, S))
K <- rep(3,S)
R <- matrix(NA,nrow = S, ncol = m)
R[1,]<- R[2,]<- R[3,]  <- sample.int(3,m,TRUE)
mu<-updatemu(R,Z,X,Gamma,K,Beta,E,m,n,p,q)
mu_star<-updatemustar(mu,rep(0.01,2),n,T0,D)
z_star<-updateZstar(mu_star,delta,n,T0)
sigma_square <- 10
C<-array(NA,dim = c(10,10,S))
w<-matrix(NA,nrow = S, ncol = 10)
w_beta<-rep(1/S,S)
for(i in 1:S){
 C[1:K[i],1:K[i],i]<-updateC(Gamma[1:K[i],,i],theta2,tau)
  w[i, 1:K[i]]<-rep(1/K[i],K[i])
}
c<-c(0,0.01)
start <- Sys.time()
for(iter in 1:Niter){
  c<-updatec(z_star, mu,D, 100,sigma_square, n, T0)
  mu_star<-updatemustar(mu,c,n,T0,D)
  z_star<-updateZstar(mu_star,delta,n, T0)
  w <- update_w(K, R, S)
  R <- updateR(w, Gamma, Beta,
               Y, Z, delta, mu, mu_star, c[2], S,
               sigma_square, K, E, X,
               m, n, q, p, T0)
  for(i in unique(E)){
    ind<-sort(unique(R[i,]))
    KK<-length(ind)
    Gamma_temp<-Gamma[ind,,i]
    Gamma[,,i]<-NA
    Gamma[1:KK,,i]<-Gamma_temp
    w_temp<-w[i,ind]
    w_temp<-w_temp/sum(w_temp)
    w[i,]<-NA
    w[i,1:KK]<-w_temp
    for(k in 1:KK){
      R[i,which(R[i,]==ind[k])]<-k
    }
    K[i]<-KK
  }
  mu<-updatemu(R,Z,X,Gamma,K,Beta,E,m,n,p,q)
  mu_star<-updatemustar(mu,c,n,T0,D)
  step <- array(rnorm(max(K) * S *q, 0, nu_gamma),dim=c( max(K), q,S))
  run<- array(runif(max(K) * S *q, 0, 1),dim=c( max(K), q,S))
  A<-updateGamma(X,Y, Z, delta, Beta, Gamma, E, R, S, K , mu, mu_star, sigma_square, rep(c,T0),
                 step,  run,  n,  m, T0,  p,  q,  D,theta2,  tau)
  Gamma<-A$Gamma
  mu<-A$mu
  mu_star<-A$mustar
  for(i in 1:S){
    if(K[i]==1){
      Gammai = t(as.matrix(Gamma[1:K[i],,i]))
      C[1:K[i],1:K[i],i]<-updateC(Gammai,theta2,tau)
    }
    else{
      C[1:K[i],1:K[i],i]<-updateC(Gamma[1:K[i],,i],theta2,tau)
    }
  }
  A<-update_RJ(w, K, Gamma,Beta, E,
               Z, X, R, mu, mu_star, Y, delta, c,sigma_square, C,
               S, theta2, tau, q, m, T0)
  K<-A$K
  w<-A$w
  Gamma<-A$Gamma
  R<-A$R
  C<-A$C
  mu<-updatemu(R,Z,X,Gamma,K,Beta,E,m,n,p,q)
  mu_star<-updatemustar(mu,rep(c,T0),n,T0,D)
  C_beta<-updateC(Beta,theta1,tau)
  step<-matrix(rnorm(S*p,0,nu_beta),nrow = S)
  runif<-matrix(runif(S*p,0,1),nrow = S)
  A<- updateBeta( X,Y, Z, delta,
                  Beta, Gamma, E,R,S,K,mu_star, mu,sigma_square,rep(c,T0), C_beta, step,runif,
                  n, m, T0,  p,  q,  D,
                  theta1, tau)
  Beta<-A$Beta
  mu<-A$mu
  mu_star<-A$mustar
  record$Beta[iter,,]<-Beta
  A<-updateE( Beta,Gamma, w_beta, X,Y,Z,delta,E, R, S,K, mu, mu_star,sigma_square,rep(c,T0),
              n,  m,  T0,  p,  q,  D)
  E<-A$E
  mu<-A$mu
  mu_star<-A$mustar
  K<-A$Ds
  w_beta<- update_w_beta(S, w_beta, E)
  theta1<-update_theta_beta(theta1,tau,Beta)
  theta2<-update_theta_gamma(theta2,tau,Gamma,S,K)
  sigma_square <- update_sigma_squre(Y,mu)
}

BAREB documentation built on April 12, 2025, 1:53 a.m.