Nothing
#' update_w_beta
#'
#' This function updates the weights for each patient level cluster
#'
#' @param S The number of patient level clusters
#' @param E A vector that records the current clustering membership.
#' @param hyper_delta The hyper-parameter with default value being 1
#' @return The updated weights for each patient level cluster
#' @seealso \link{update_RJ} for a complete example for all functions in this package.
#' @examples
#' #Suppose we know the number of patient level cluster is 4
#' #Suppose the current clustering membership indicates 3 patients in cluster 1,
#' #2 patients in cluster 2, 3 patinets in cluster 3, 1 patient in cluster 4
#' #Use the default value, 1, for hyper-parameter
#' update_w_beta(S=4,E=c(1,1,1,2,2,3,3,3,4))
#'
#'
#' #To change the hyper-parameter to, for example 2
#' update_w_beta(S=4,E=c(1,1,1,2,2,3,3,3,4),hyper_delta = 2)
#' @export
update_w_beta <- function(S, E, hyper_delta = 1){
w_star = rep(NA, S)
for (l in 1:S){
Dir_par = hyper_delta + sum((E == l))
w_star[l] = rgamma(1, shape = Dir_par, rate = 1)
}
w_star = w_star/sum(w_star)
return(w_star)
}
#' update_w
#'
#' This function updates the weights for site level clusters
#'
#' @param K The number of site level cluster for each patient level cluster. Should be a vector of length \code{S}
#' @param R The current site level clustering membership. Should be a matrix with \code{S} rows.
#' @param S The number of patient level clusters.
#' @param hyper_delta The hyper-parameter with default value being 1
#' @return The updated weights for site level clusters. Should be a matrix with \code{S} rows.
#' @seealso \link{update_RJ} for a complete example for all functions in this package.
#' @examples
#' #Suppose we know the number of patient level cluster is 2,
#' #one has 2 site level clusters and one has 3.
#' #Use the default value, 1, for hyper-parameter
#' update_w(K=c(2,3),R=matrix(c(1,2,2,1,1,2,3,2),nrow=2,byrow=TRUE), S=2)
#'
#' #To change the hyper-parameter to, for example 2
#' update_w(K=c(2,3),R=matrix(c(1,2,2,1,1,2,3,2),nrow=2,byrow=TRUE), S=2, hyper_delta = 2)
#' @details
#' It returns a matrix with 10 columns.
#' For example, first patient cluster has 2 site level clusters. The first row's first 2 elements give the weights for site level clusters in patient cluster 1. Last 8 elements are NA's
#' @export
update_w <- function(K, R, S, hyper_delta = 1){
#S is the number of patient level clusters
#K should be a vector of length S
ans <- matrix(NA, nrow = S, ncol = 10)
for(i in 1:S){
w_star = rep(NA, K[i])
for (l in 1:K[i]){
Dir_par = hyper_delta + sum((R[i,] == l))
w_star[l] = rgamma(1, shape = Dir_par, rate = 1)
}
w_star = w_star/sum(w_star)
ans[i, 1:K[i]]<- w_star
}
return(ans)
}
#' update_RJ
#'
#' Update the number of site level clusters for each patient level cluster via RJMCMC
#'
#' @param w The weights for site level clusters. Should be a matrix with S rows.
#' @param K The number of site level cluster for each patient level cluster. Should be a vector of length S
#' @param Gamma The linear coefficients for site level covariates. Should be a 3-dimensional array.
#' @param Beta The linear coefficients for patient level covariates. Should be a matrix with S rows
#' @param E A vector that records the current clustering membership.
#' @param Z The design matrix for site level clusters.
#' @param X The design matrix for patient level clusters.
#' @param R The current site level clustering membership
#' @param mu The current CAL mean matrix
#' @param mu_star The latent value matrix for missingness model
#' @param Y The observed CAL value matrix
#' @param delta The missing indicator matrix
#' @param c The linear coefficients for missingness model
#' @param sigma_square The current noise variance
#' @param C The DPP related kernel matrices. Should be an array of 3 dimensions
#' @param S The number of patient level clusters.
#' @param theta The DPP hyper-parameter for site level
#' @param tau A fixed DPP hyper-parameter, which we suggest high value, say 10^5
#' @param q The number of site level covariates
#' @param m The number of sites
#' @param T0 The number of teeth
#' @param hyper_delta The hyper-parameter with default value being 1
#' @return A list with following updated parameters:
#' @return \item{K}{The numbers of site level clusters}
#' @return \item{w}{The weights for site level clusters}
#' @return \item{Gamma}{Linear coefficients for site level covariates}
#' @return \item{R}{The site level clusteirng membership}
#' @return \item{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)
#' }
#'
#' @export
update_RJ <- function(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){
for(i in 1:S){
wi<- w[i, 1:K[i]]
Gammai<- Gamma[1:K[i],,i]
Betai<- Beta[i,]
indi<- which(E == i)
Ri<-R[i,]
Ki<-K[i]
Ci<- C[1:Ki,1:Ki, i]
if(Ki==1){
Gammai = t(as.matrix(Gammai))
Ci = as.matrix(Ci)
wi = as.vector(wi)
}
if(length(indi) == 0){
A<-RJi_empty(wi, Ki, Gammai, Ri, Ci,theta, tau, m, q, hyper_delta)
}
else{
Yi<-Y[indi,]
deltai<-delta[indi,]
mui<-mu[indi,]
mu_stari<-mu_star[indi,]
Xi<-X[indi,]
if(length(indi)==1){
Yi <- t(as.matrix(Yi))
deltai <- t(as.matrix(deltai))
mui <- t(as.matrix(mui))
mu_stari <- t(as.matrix(mu_stari))
Xi <- t(as.matrix(Xi))
}
A<-RJi(wi, Ki, as.matrix(Gammai), Betai, Xi, Yi, Z, Ri, deltai, mui[1,], mu_stari, 0.01, sigma_square, Ci,
theta, tau, m, length(indi),q,T0, hyper_delta=1)
}
K[i]<-A$K
w[i, 1:K[i]]<- A$w
Gamma[,,i]<-NA
Gamma[1:K[i],,i]<-A$Gamma
R[i,]<-A$R
C[1:K[i],1:K[i], i]<-A$C
}
return(list(K = K, w = w, Gamma = Gamma, R = R, C = C))
}
#' update_sigma_squre
#'
#' Update the noise variance
#'
#' @param Y The CAL observation matrix, with missing values
#' @param mu Current estimated mean matrix for CAL
#' @param a The hyper-parameter with default value being 1
#' @param b The hyper-parameter with default value being 1
#' @return Updated noise variance
#' @seealso \link{update_RJ} for a complete example for all functions in this package.
#' @export
update_sigma_squre<-function(Y, mu, a =1,b=1){
temp <- na.omit(as.vector(Y-mu))
s <- sum(temp^2)
sigma2<-1/rgamma(1,a+length(temp)/2,b+s/2)
return (sigma2)
}
#' const
#'
#' The function to compute the normalizing constant for Determinantal Point Process (DPP)
#'
#' @param theta The parameter that controls how repulsive the DPP is
#' @param tau The parameter we fixed at a large number
#' @return The computed normalizing constant
#' @export
const <- function(theta, tau){
a = 1/(4*tau^2)
b = 1/(theta^2)
c = sqrt(a^2 + 2*a*b)
sum = a + b + c
lambda = rep(0, 100)
for (i in 1:100){
eig1 = sqrt(2*a/sum)*(b/sum)^(i-1)
if (eig1 > 0.000001) {lambda[i] = eig1}
else{
break
}
}
return(prod(lambda+1))
}
#' likeli_theta
#'
#' The function computes the likelihood of the hyperparameter
#'
#' @param theta The parameter that controls how repulsive the DPP is
#' @param tau The parameter we fixed at a large number
#' @param Beta The matrix that records the samples
#' @return The likelihood
#' @export
likeli_theta<-function(theta,tau, Beta){
#function to minimize in order to update theta
if(is.matrix(Beta)){
C<-updateC(Beta,theta,tau)
Likeli<-log(det(C))-log(const(theta,tau))
}
else{
C<-kernelC(Beta,Beta,theta,tau)
Likeli<-log(C)-log(const(theta,tau))
}
return(Likeli)
}
#' update_theta_beta
#'
#' Update the DPP hyper-parameter for patient level
#'
#' @param theta The DPP hyper-parameter for patient level
#' @param tau A fixed DPP hyper-parameter, which we suggest high value, say 10^5
#' @param Beta The linear coefficients for patient level covariates. Should be a matrix with S rows
#' @param sig The hyper-parameter with default value being 10
#' @return updated DPP hyper-parameter for patient level
#' @seealso \link{update_RJ} for a complete example for all functions in this package.
#' @export
update_theta_beta<-function(theta, tau, Beta,sig = 10){
theta.new<- theta + (rnorm(1,0,0.1))
ratio <- (likeli_theta(theta.new,tau,Beta) + dnorm(theta.new,0,sig,log=TRUE) -
likeli_theta(theta,tau,Beta) - dnorm(theta,0,sig,log=TRUE))
if(ratio>log(runif(1))){
return(theta.new)
}
else{
return(theta)
}
}
#' update_theta_gamma
#'
#' Update the DPP hyper-parameter for site level
#'
#' @param theta The DPP hyper-parameter for patient level
#' @param tau A fixed DPP hyper-parameter, which we suggest high value, say 10^5
#' @param Gamma The linear coefficients for site level covariates. Should be a 3-dimensional array.
#' @param S The number of patient level clusters
#' @param Ds The number for site level clusters for each patient level cluster
#' @param sig The hyper-parameter with default value being 10
#' @return updated DPP hyper-parameter for site level
#' @seealso \link{update_RJ} for a complete example for all functions in this package.
#' @export
update_theta_gamma<-function(theta,tau,Gamma, S, Ds, sig = 10 ){
theta.new<- theta + (rnorm(1,0,0.1))
ratio <- 0
for(i in 1:S){
Gammatemp<-Gamma[1:Ds[i],,i]
ratio <- ratio + (likeli_theta(theta.new,tau,Gammatemp) + dnorm(theta.new,0,sig,log=TRUE) -
likeli_theta(theta,tau,Gammatemp) - dnorm(theta,0,sig,log=TRUE))
}
if(ratio>log(runif(1))){
return(theta.new)
}
else{
return(theta)
}
}
#tested
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.