R/sampling_fcts_KMR.R

In apeterson91/bekmr:

### Authored by Veronica Berrocal and Adam Peterson
### Functions to sample parameters according
### to a MCMC algorithm
### The model we consider is: for i=1,...,n
###     logit(y_i) = h(z_i) + X_i*beta
### where z_i = (z_{i1},z_{i2},..,z_{i(i-1)},z_{i(i+1)},..,z_{in}).
###
### The priors that we are using are:
###
###   h(z) ~ GP( 0, sigma^2*R(phi))
###
###  with R(phi) correlation function and phi
###  covariance parameter (range parameter).
###
### sigma^2 ~ InverseGamma(shape.sigma2,scale.sigma2)
### phi ~ Uniform(lower.phi,upper.phi)
###
### beta~Normal(m.beta,V)
###
### Notation:
###  y: nx1 vector of 0-1 for the n schools
###  X: nxp design matrix of covariates with first column being the intercept
###  beta: px1 vector of regression coefficients
###  z: nx(n-1) matrix with each row the vector of co-clustering probabilities for a school
###  h.z: nx1 vector with the values of the Gaussian process at "coordinates" z_1, z_2,..z_n
###  sigma.2: marginal variance of the Gaussian process prior on h
###  phi: range parameter of the Gaussian covariance function used to model the covariance
###       of the Gaussian process
###  m.beta: px1 prior mean of the multivariate normal distribution for the vector of
###          regression coefficients beta
###  V: pxp prior covariance matrix of the multivariate normal distribution for the vector of
###          regression coefficients beta
###  iV: inverse of the prior covariance matrix V for the beta vector

#' This is to compute the distances among the z vectors for each school
#' z.mat is a n by (n-1) matrix where each row is the (n-1) vector of co-clustering probabilities
#' Specifically, the i-th row of z.mat is the vector of co-clustering probabilities for school i
# n <- length(y)
# Dist.mat <- rdist(z.mat)
# Dist.mat.sqrd <- Dist.mat*Dist.mat

#' This function computes the inverse of the matrix obtained using a Gaussian correlation function
comp.gau.icorrfct <- function(Dist.mat.sqrd,phi){
    R.phi <- exp(-Dist.mat.sqrd/(phi^2))
    iR.phi <- solve(R.phi)
    return(iR.phi)
}

#' This function computes the log likelihood for the binary data
#' for the n schools, when the probability of obesity at school i
#' is linked through a logistic regression to covariates X_i (the i-th row
#' of the design matrix X refers to the p covariates - including the intercept -
#' relative to school i), with h(z_i) being the kernel machine regression value at school i
binary.logreg.loglik <- function(y,h.z,X,beta,n){

    num.loglik <- sum(h.z+as.numeric(X%*%beta)[which(y!=0)])
    den.loglik <- sum(log(1+exp(h.z+as.numeric(X%*%beta))))
    bin.loglik <- num.loglik-den.loglik
    return(bin.loglik)
}


#' This function computes the log likelihood for the binomial data using a logit link function
#'
#' @param y vector of length n where y[i] = #obese students for school i
#' @param h.z vector  of length n where h[i] is the "cluster" effect for school i
#' @param X
#' @param beta
binomial.logreg.loglik <- function(y,h.z,X,beta,n){

	eta <- h.z + as.numeric(X%*%beta) ## calculate linear predictor eta
	mu <- binomial(link='logit')$linkinv(eta) ## take sigmoid of linear predictor eta
	bin.loglik <-  sum(dbinom(x = y, size = n, prob = mu, log = TRUE)) ## calculate loglikelihood
	return(bin.loglik)

}

#' This function samples a new vector h.z via Metropolis Hastings
#' The proposal distribution is a Multivariate Normal centered
#' at the current value of the h.z with a diagonal covariance matrix
#' with proposal sdev delta.h
#' The value of delta.h can be changed depending on the acceptance rate of the proposed values
#' during the MCMC algorithm and during burnin
sample.h <- function(y,h.z,X,beta,n,delta.h,iR.phi){

    # sample a proposed nx1 vector, h*, for h
    h.star.z <- MASS::mvrnorm(1,h.z,diag((delta.h^2),length(n) ))

    ### numerator of the Metropolis Hasting ratio on the log scale:
    ### this is the sum of two parts - the log of the binary likelihood
    ### and the log of the multivariate normal prior
    ### both evaluated when h is equal to the proposed value
    num.metrop.ratio.pt1 <- binomial.logreg.loglik(y,h.star.z,X,beta,n)
    num.metrop.ratio.pt2 <- -0.5*as.numeric(t(matrix(h.star.z,nrow=length(n),ncol=1))%*%iR.phi%*%matrix(h.star.z,nrow=length(n),ncol=1))
    num.metrop.ratio <- num.metrop.ratio.pt1+num.metrop.ratio.pt2

    ### denominator of the Metropolis Hasting ratio on the log scale:
    ### this is the sum of two parts - the log of the binary likelihood
    ### and the log of the multivariate normal prior
    ### both evaluated when h is the current value
    den.metrop.ratio.pt1 <- binomial.logreg.loglik(y,h.z,X,beta,n)
    den.metrop.ratio.pt2 <- -0.5*as.numeric(t(matrix(h.z,nrow=length(n),ncol=1))%*%iR.phi%*%matrix(h.z,nrow=length(n),ncol=1))
    den.metrop.ratio <- den.metrop.ratio.pt1+den.metrop.ratio.pt2

    ratio <- min(1,exp(num.metrop.ratio-den.metrop.ratio))
    u <- runif(1,0,1)
    ifelse(u <= ratio,new.h.z <- h.star.z,new.h.z <- h.z)
    return(new.h.z)
}


#' This function samples a new vector beta via Metropolis Hastings
#' The proposal distribution is a Multivariate Normal centered
#' at the current value of the beta with a diagonal covariance matrix
#' with proposal sdev delta.beta
#' The value of delta.beta can be changed depending on the acceptance rate of the proposed values.
#' during the MCMC algorithm and during burnin
sample.beta <- function(y,h.z,X,beta,p,delta.beta,m.beta,iV,n){

    # sample a proposed px1 vector, beta*, for beta
    beta.star <- MASS::mvrnorm(1,beta,diag((delta.beta^2),p))

    ### numerator of the Metropolis Hasting ratio on the log scale:
    ### this is the sum of two parts - the log of the binary likelihood
    ### and the log of the multivariate normal prior
    ### both evaluated when beta is equal to the proposed value
    num.metrop.ratio.pt1 <- binomial.logreg.loglik(y,h.z,X,beta.star,n)
    num.metrop.ratio.pt2 <- -0.5*as.numeric(t(matrix((beta.star-m.beta),nrow=p,ncol=1))%*%iV%*%matrix((beta.star-m.beta),nrow=p,ncol=1))
    num.metrop.ratio <- num.metrop.ratio.pt1+num.metrop.ratio.pt2

    ### denominator of the Metropolis Hasting ratio on the log scale:
    ### this is the sum of two parts - the log of the binary likelihood
    ### and the log of the multivariate normal prior
    ### both evaluated when beta is the current value
    den.metrop.ratio.pt1 <- binomial.logreg.loglik(y,h.z,X,beta,n)
    den.metrop.ratio.pt2 <- -0.5*as.numeric(t(matrix((beta-m.beta),nrow=p,ncol=1))%*%iV%*%matrix((beta-m.beta),nrow=p,ncol=1))
    den.metrop.ratio <- den.metrop.ratio.pt1+den.metrop.ratio.pt2

    ratio <- min(1,exp(num.metrop.ratio-den.metrop.ratio))
    u <- runif(1,0,1)
    ifelse(u <= ratio,new.beta <- beta.star,new.beta <- beta)
    return(new.beta)
}


#' This function samples a new value for sigma.2
sample.sigma.2 <- function(h.z,n,iR.phi,shape.sigma2,scale.sigma2){

    post.shape.sigma2 <- shape.sigma2+(length(n)/2)
    post.scale.sigma2 <- (0.5*as.numeric(t(matrix(h.z,nrow=length(n),ncol=1))%*%iR.phi%*%(matrix(h.z,nrow=length(n),ncol=1))))+scale.sigma2

    new.sigma2 <- 1/rgamma(1,post.shape.sigma2,post.scale.sigma2)
    return(new.sigma2)
}


#' This function samples a new value of the range parameter phi via Metropolis Hastings
#' The proposal distribution is a Uniform distribution on the interval (lower.phi,upper.phi)
#' that does not depend on the current value of phi
sample.phi <- function(h.z,phi,iR.phi,Dist.mat.sqrd,lower.phi,upper.phi,n){

    # sample a proposed value, phi*, for phi
    phi.star <- runif(1,lower.phi,upper.phi)
    # this is the inverse of the Gaussian correlation matrix corresponding to the proposed
    # value for phi
    iR.phi.star <- comp.gau.icorrfct(Dist.mat.sqrd,phi.star)

    ### numerator of the Metropolis Hasting ratio on the log scale:
    ### this is the log of multivariate normal prior for h evaluated when
    ### the precision matrix is corresponding to the proposed value of phi, phi.star
    num.metrop.ratio <- (0.5*determinant(iR.phi.star)$modulus)-0.5*as.numeric(t(matrix(h.z,nrow = length(n),ncol=1))%*%iR.phi.star%*%matrix(h.z,nrow = length(n),ncol=1))

    ### denominator of the Metropolis Hasting ratio on the log scale:
    ### this is the log of multivariate normal prior for h evaluated when
    ### the precision matrix is corresponding to the current value of phi
    den.metrop.ratio <- (0.5*determinant(iR.phi)$modulus)-0.5*as.numeric(t(matrix(h.z,nrow=length(n),ncol=1))%*%iR.phi%*%matrix(h.z,nrow = length(n),ncol=1))

    ratio <- min(1,exp(num.metrop.ratio-den.metrop.ratio))
    u <- runif(1,0,1)
    ifelse(u <= ratio,new.phi <- phi.star,new.phi <- phi)
    return(new.phi)
}