update_RJ: update_RJ

Description Usage Arguments Value Examples

View source: R/updates.R

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
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 March 26, 2020, 7:36 p.m.

Related to update_RJ in BAREB...