##### November 1st 2014, Mathias Laplace Approximation Simulation ##### Bayesian logistic regression # Function to calculate the likelihood function loglik <- function(y,Xsub,betasub){ return(sum(y*Xsub%*%betasub-log(1+exp(Xsub%*%betasub)))) } # Monte Carlo to get Bayesian Integrand montecarlo <- function(y,Xsub,q,s){ estloglik <- NULL for(i in 1:s){ betasub <- matrix(rnorm(q,0,1),ncol=1) estloglik <- c(estloglik,loglik(y,Xsub,betasub)) } return(log(mean(exp(estloglik+n/2)))) } # Function to calculate determinant of Hessian Hessian <- function(Xsub,betasub){ res <- matrix(0,ncol(Xsub),ncol(Xsub)) for(i in 1:nrow(Xsub)){ res <- res+t(matrix(Xsub[i,],nrow=1))%*%matrix(Xsub[i,],nrow=1)*as.numeric(exp(matrix(Xsub[i,],nrow=1)%*%betasub)/((1+exp(matrix(Xsub[i,],nrow=1)%*%betasub))^2)) } return(res) } # Run simulation for various n, p, and q=1 runsimOne <- function(n,p,q,s,iter){ set.seed(iter) beta <- matrix(c(rep(2,q),rep(0,p-q)),nrow=p) X <- matrix(rnorm(n*p,0,1),n,p) prob <- exp(X%*%beta)/(1+exp(X%*%beta)) y <- rbinom(n,1,prob) Bayes <- Bayes2 <- laplace <- NULL for(j in 1:p){ Xsub <- matrix(X[,j],ncol=1) Bayes <- c(Bayes,montecarlo(y,Xsub,q,s=s)) Bayes2 <- c(Bayes2,montecarlo(y,Xsub,q,s=s)) estbeta <- matrix(summary(glm(y~-1+Xsub,family="binomial"))$coef[,1],ncol=1) laplace <- c(laplace,log(exp(loglik(y,Xsub,estbeta)+n/2)*exp(sum(-estbeta^2/2))*1/sqrt(det(Hessian(Xsub,estbeta))))) } return(list(Bayes=Bayes,Bayes2=Bayes2,laplace=laplace)) } # Run simulation for various n, p, and q=2 runsimTwo <- function(n,p,q,s,iter){ set.seed(iter) beta <- matrix(c(rep(2,q),rep(0,p-q)),nrow=p) X <- matrix(rnorm(n*p,0,1),n,p) prob <- exp(X%*%beta)/(1+exp(X%*%beta)) y <- rbinom(n,1,prob) Bayes <- Bayes2 <- laplace <- NULL for(j in 1:(p-1)){ for(k in (j+1):p){ Xsub <- X[,c(j,k)] Bayes <- c(Bayes,montecarlo(y,Xsub,q,s=s)) Bayes2 <- c(Bayes2,montecarlo(y,Xsub,q,s=s)) estbeta <- matrix(summary(glm(y~-1+Xsub,family="binomial"))$coef[,1],ncol=1) laplace <- c(laplace,log(exp(loglik(y,Xsub,estbeta)+n/2)*exp(sum(-estbeta^2/2))*1/sqrt(det(Hessian(Xsub,estbeta))))) } } return(list(Bayes=Bayes,Bayes2=Bayes2,laplace=laplace)) } # Run simulation for various n, p, and q=3 runsimThree <- function(n,p,q,s,iter){ set.seed(iter) beta <- matrix(c(rep(2,q),rep(0,p-q)),nrow=p) X <- matrix(rnorm(n*p,0,1),n,p) prob <- exp(X%*%beta)/(1+exp(X%*%beta)) y <- rbinom(n,1,prob) Bayes <- Bayes2 <- laplace <- NULL for(j in 1:(p-2)){ for(k in (j+1):(p-1)){ for(z in (k+1):p){ Xsub <- X[,c(j,k,z)] Bayes <- c(Bayes,montecarlo(y,Xsub,q,s=s)) Bayes2 <- c(Bayes2,montecarlo(y,Xsub,q,s=s)) estbeta <- matrix(summary(glm(y~-1+Xsub,family="binomial"))$coef[,1],ncol=1) laplace <- c(laplace,log(exp(loglik(y,Xsub,estbeta)+n/2)*exp(sum(-estbeta^2/2))*1/sqrt(det(Hessian(Xsub,estbeta))))) } } } return(list(Bayes=Bayes,Bayes2=Bayes2,laplace=laplace)) }