R/be_kmr.R

In apeterson91/bekmr:

#'Fits a Kernel Machine Regression logistic regression model
#'
#' @param formula similar to \code{\link[stats]{glm}}
#' @param failures
#' @param Z matrix of pairwise probabilities
#' @param iter_max maximum number of iterations
#' @param burnin number of iterations to use for burn-in
#'
#' @export
be_kmr <- function(formula,data,
                   Z,
                   iter_max = 2E3,
                   burnin = 1E3,
                   shape.sigma2 = 2,
                   scale.sigma2 = 1,
                   phi = 0.1
                   )
{


  call <- match.call(expand.dots = TRUE)
  mf <- match.call(expand.dots = FALSE)
  mf$formula <- formula
  m <- match(c("formula"),table = names(mf), nomatch=0L)
  mf <- mf[c(1L,m)]
  mf$data <- data
  mf$drop.unused.levels <- TRUE
  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf,parent.frame())
  mt <- attr(mf,"terms")
  if(is.empty.model(mt))
    stop("No intercept or predictors specified.",.call = FALSE)
  y <- model.response(mf,type = "any")
  X <- model.matrix(mt,mf)
  if(attr(mt,"intercept"))
    X <- X[,2:ncol(X),drop=F]

	## Paste code that turns formula into X and successes/failures vectors
	init.beta <- coef(glm(formula,data=data,family=binomial(link='logit')))[2:(ncol(X)+1)]
	init.h.z <- rep(0,ncol(Z))
	init.sigma.2 <-mean(apply(Z,1,var))
	init.phi <- phi
	# init.phi <- runif(1,min(as.numeric(fields::rdist(Z))),max(as.numeric(fields::rdist(Z))))

  ## proposal var
	delta.h <- 1 / sqrt(ncol(Z))
	delta.beta <- 1 / sqrt(ncol(X))


	## sample containers
	beta_samples <- matrix(0,nrow=iter_max,ncol=ncol(X))
	colnames(beta_samples) <- colnames(X)
	h_samples <- matrix(0,ncol=nrow(Z),nrow=iter_max)
	colnames(h_samples) <- paste0("h_",1:ncol(h_samples))
	sigma_sq_samples <- matrix(0,ncol=1,nrow=iter_max)
	phi_samples <- matrix(0,ncol=1,nrow=iter_max)


	n <- rowSums(y)
	y <- y[,1]
	p <- ncol(X)
	Dist.mat <- fields::rdist(Z)
	Dist.mat.sqrd <- Dist.mat*Dist.mat

	## prior information

	iV <- diag(ncol(X)) / 9
	m.beta <- rep(0,p)

	lower.phi <- min(as.numeric(fields::rdist(Z)))
	upper.phi <- max(as.numeric(fields::rdist(Z)))

	comp.gau.icorrfct <- function(Dist.mat.sqrd,phi){
		R.phi <- exp(-Dist.mat.sqrd/phi)
		iR.phi <- solve(R.phi)
		rm(R.phi)
		return(iR.phi)
	}

	accept.h <- 0
	accept.beta <- 0
	for(k in 1:iter_max){
		if(k==1){
			init.iR.phi <- comp.gau.icorrfct(Dist.mat.sqrd,init.phi)
			new.h.z <- sample.h(y,init.h.z,X,init.beta,n,delta.h,init.iR.phi)
			new.beta <- sample.beta(y,new.h.z,X,init.beta,p,delta.beta,m.beta,iV,n)
			new.sigma.2 <- sample.sigma.2(new.h.z,n,init.iR.phi,shape.sigma2,scale.sigma2)
			new.phi <- sample.phi(new.h.z,init.phi,init.iR.phi,Dist.mat.sqrd,lower.phi,upper.phi,n)
			h.z <- init.h.z
			beta <- init.beta
			sigma.2 <- init.iR.phi
			phi <- init.phi
		}

		if(k!=1){
			iR.phi <- comp.gau.icorrfct(Dist.mat.sqrd,phi)
			new.h.z <- sample.h(y,h.z,X,beta,n,delta.h,iR.phi)
			new.beta <- sample.beta(y,new.h.z,X,beta,p,delta.beta,m.beta,iV,n)
			new.sigma.2 <- sample.sigma.2(new.h.z,n,iR.phi,shape.sigma2,scale.sigma2)
			new.phi <- init.phi # new.phi <- sample.phi(new.h.z,phi,iR.phi,Dist.mat.sqrd,lower.phi,upper.phi,n)
		}

		ifelse(sum(new.h.z==h.z)!=length(n),accept.h <- accept.h+1,accept.h <- accept.h)
		ifelse(sum(new.beta==beta)!=p,accept.beta <- accept.beta+1,accept.beta <- accept.beta)

		h.z <- new.h.z
		beta <- new.beta
		sigma.2 <- new.sigma.2
		phi <- new.phi

		beta_samples[k,] <- beta
		h_samples[k,] <- h.z
		sigma_sq_samples[k] <- sigma.2
		phi_samples[k] <- new.phi

		### the treshold for acceptance can be lowered to match what suggested
		### in various papers; also the factors by which the sd of the proposals
		### are rescaled can be changed
		if((k <= burnin) & ((k%%100)==0)){
			if((accept.h/100) >= 0.4){
				delta.h <- delta.h*1.5
			}
			if((accept.h/100) < 0.1){
				delta.h <- delta.h/1.5
			}
			if((accept.beta/100) >= 0.4){
				delta.beta <- delta.beta*1.5
			}
			if((accept.beta/100) < 0.1){
				delta.beta <- delta.beta/1.5
			}
		  accept.h <- 0
		  accept.beta <- 0
		}

		if(k%%100==0)
		  print(paste("Iteration ",k,sep=""))
	}

	return(list(beta_samples = beta_samples[(burnin+1):iter_max,,drop=F],
				phi_samples = phi_samples[(burnin+1):iter_max,,drop=F],
				sigma_sq_samples = sigma_sq_samples[(burnin+1):iter_max,,drop=F],
				h_samples = h_samples[(burnin+1):iter_max,,drop=F],
				accept.h = accept.h,
				accept.beta = accept.beta))

}