Description Usage Arguments Value Examples
Update the number of site level clusters for each patient level cluster via RJMCMC
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
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 |
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 |
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)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.