############# November 13, 2009 ############# Jean-Michel Marin and Christian P. Robert ############# On resolving the Savage--Dickey paradox ################################################################################ probitll=function(beta,y,X) { # LOG-LIKELIHOOD OF THE PROBIT MODEL if (is.matrix(beta)==F) beta=as.matrix(t(beta)) niter=dim(beta)[1] n=length(y) k=dim(beta)[2] lF1=pnorm(X%*%t(beta),log=TRUE) lF2=pnorm(-X%*%t(beta),log=TRUE) lala=matrix(rep(y,niter),n) lulu=lala*lF1+(1-lala)*lF2 colSums(lulu) } ################################################################################ probitlpost=function(beta,y,X) { # LOG-POSTERIOR OF THE PROBIT MODEL if (is.matrix(beta)==F) beta=as.matrix(t(beta)) niter=dim(beta)[1] n=length(y) k=dim(beta)[2] lF1=pnorm(X%*%t(beta),log=TRUE) lF2=pnorm(-X%*%t(beta),log=TRUE) lala=matrix(rep(y,niter),n) lulu=lala*lF1+(1-lala)*lF2 colSums(lulu)+dmvnorm(beta,mean=rep(0,k),sigma=n*ginv(t(X)%*%X),log=TRUE) } ################################################################################ rntneg=function(n,mu,sigma2) { # PRODUCE n SAMPLES FROM A SPECIFIED RIGHT-TRUNCATED TO 0 NORMAL DISTRIBUTION qnorm(runif(n)*pnorm(0,mu,sqrt(sigma2)),mu,sqrt(sigma2)) } ################################################################################ rntpos=function(n,mu,sigma2) { # PRODUCE n SAMPLES FROM THE SPECIFIED LEFT-TRUNCATED TO 0 NORMAL DISTRIBUTION -rntneg(n,-mu,sigma2) } ################################################################################ gibbsprobit=function(Niter,y,X,Sigmapost) { # GIBBS SAMPLING FOR THE PROBIT MODEL # CENTERED GAUSSIAN PRIOR ON THE REGRESSION COEFFICIENTS WITH COVARIANCE MATRIX # Sigmaprior IN THAT CASE Sigmapost=ginv(t(X)%*%X+ginv(Sigmaprior)) p=dim(X)[2] n=length(y) beta=matrix(0,Niter,p) Mu=matrix(0,Niter,p) z=rep(0,n) mod=summary(glm(y~-1+X,family=binomial(link="probit"))) beta[1,]=as.vector(mod$coefficient[,1]) moan=X%*%beta[1,] z[y==1]=rntpos(sum(y==1),moan[y==1],1) z[y==0]=rntneg(sum(y==0),moan[y==0],1) Mu[1,]=Sigmapost%*%t(X)%*%z for (i in 2:Niter) { beta[i,]=rmvnorm(1,mean=Mu[i-1,],sigma=Sigmapost) moan=X%*%beta[i,] z[y==1]=rntpos(sum(y==1),moan[y==1],1) z[y==0]=rntneg(sum(y==0),moan[y==0],1) Mu[i,]=Sigmapost%*%t(X)%*%z } list(beta=beta,mu=Mu) } ################################################################################ gibbsprobitsd1=function(Niter,y,Sigmapost) { # SPECIAL GIBBS SAMPLER FOR CALCULATING THE VERDINELLI-WASSERMAN (1995) SAVAGE-DICKEY # APPROXIMATION OF THE BAYES FACTOR WHICH CORERSPONDS OF TESTING THE NULLITY OF # THE COEFFICIENT ASSOCIATED TO REGRESSOR 3 p=dim(X0)[2] n=length(y) beta=matrix(0,Niter,p) z=rep(0,n) mod=summary(glm(y~-1+X0,family=binomial(link="probit"))) beta[1,]=as.vector(mod$coefficient[,1]) moan=X0%*%beta[1,] z[y==1]=rntpos(sum(y==1),moan[y==1],1) z[y==0]=rntneg(sum(y==0),moan[y==0],1) Mu=Sigmapost%*%t(X1)%*%z Mucond=Mu[1:2]-Mu[3]/Sigmapost[3,3]*matrix(Sigmapost[1:2,3],2,1) Sigmacond=Sigmapost[1:2,1:2]-matrix(Sigmapost[1:2,3],2,1)%*%matrix(Sigmapost[3,1:2],1,2)/ Sigmapost[3,3] for (i in 2:Niter) { beta[i,]=rmvnorm(1,mean=Mucond,sigma=Sigmacond) moan=X0%*%beta[i,] z[y==1]=rntpos(sum(y==1),moan[y==1],1) z[y==0]=rntneg(sum(y==0),moan[y==0],1) Mu=Sigmapost%*%t(X1)%*%z Mucond=Mu[1:2]-Mu[3]/Sigmapost[3,3]*matrix(Sigmapost[1:2,3],2,1) } beta } ################################################################################ library(MASS) library(mvtnorm) y=dget("y.R") X0=dget("X0.R") X1=dget("X1.R") n=length(y) Niter=20000 NMC=100 model1=summary(glm(y~-1+X0,family=binomial(link="probit"))) model2=summary(glm(y~-1+X1,family=binomial(link="probit"))) ############# MARIN-ROBERT (2009) DICKEY-SAVAGE APPROXIMATION bfsd1=rep(0,NMC) Sigmaprior1=n*ginv(t(X1)%*%X1) Sigmapriortilde=matrix(0,3,3) Sigmapriortilde[1:2,1:2]=n*ginv(t(X0)%*%X0) Sigmapriortilde[3,3]=Sigmaprior1[3,3] Sigmapost1=ginv(t(X1)%*%X1+ginv(Sigmaprior1)) Sigmaposttilde=ginv(t(X1)%*%X1+ginv(Sigmapriortilde)) for (i in 1:NMC) { gibbs1=gibbsprobit(Niter,y,X1,Sigmapost1) gibbs2=gibbsprobit(Niter,y,X1,Sigmaposttilde) mcond2=gibbs2$mu[,3]+matrix(Sigmaposttilde[3,1:2],1,2)%*%ginv(Sigmaposttilde[1:2,1:2])%*% t(gibbs2$beta[,1:2]-gibbs2$mu[,1:2]) vcond2=Sigmaposttilde[3,3]-matrix(Sigmaposttilde[3,1:2],1,2)%*%ginv(Sigmaposttilde[1:2,1:2])%*% matrix(Sigmaposttilde[1:2,3],2,1) bfsd1[i]=mean(dnorm(0,mcond2,sqrt(vcond2)))* mean(exp(dmvnorm(gibbs1$beta,rep(0,3),Sigmapriortilde,log=TRUE)- dmvnorm(gibbs1$beta,rep(0,3),Sigmaprior1,log=TRUE))) print(i) } bfsd1=bfsd1/dnorm(0,sd=sqrt(Sigmaprior1[3,3])) ############# VERDINELLI-WASSERMAN (1995) DICKEY-SAVAGE APPROXIMATION bfsd2=rep(0,NMC) Sigmaprior0=n*ginv(t(X0)%*%X0) Sigmaprior1=n*ginv(t(X1)%*%X1) Sigmapost1=ginv(t(X1)%*%X1+ginv(Sigmaprior1)) for (i in 1:NMC) { gibbs1=gibbsprobit(Niter,y,X1,Sigmapost1) gibbs2=gibbsprobitsd1(Niter,y,Sigmapost1) mcond1=gibbs1$mu[,3]+matrix(Sigmapost1[3,1:2],1,2)%*%ginv(Sigmapost1[1:2,1:2])%*% t(gibbs1$beta[,1:2]-gibbs1$mu[,1:2]) vcond1=Sigmapost1[3,3]-matrix(Sigmapost1[3,1:2],1,2)%*%ginv(Sigmapost1[1:2,1:2])%*% matrix(Sigmapost1[1:2,3],2,1) bfsd2[i]=mean(dnorm(0,mcond1,sqrt(vcond1)))* mean(exp(dmvnorm(gibbs2,rep(0,2),Sigmaprior0,log=TRUE)- dmvnorm(cbind(gibbs2,rep(0,Niter)),rep(0,3),Sigmaprior1,log=TRUE))) print(i) } ############# CHIB APPROXIMATION Sigmaprior0=n*ginv(t(X0)%*%X0) Sigmapost0=ginv(t(X0)%*%X0+ginv(Sigmaprior0)) Sigmaprior1=n*ginv(t(X1)%*%X1) Sigmapost1=ginv(t(X1)%*%X1+ginv(Sigmaprior1)) bfchi=rep(0,NMC) for (i in 1:NMC) { gibbs1=gibbsprobit(Niter,y,X0,Sigmapost0) gibbs2=gibbsprobit(Niter,y,X1,Sigmapost1) bfchi[i]=mean(exp(dmvnorm(t(t(gibbs2$mu)-model2$coeff[,1]),mean=rep(0,3),sigma=Sigmapost1,log=TRUE)-probitlpost(model2$coeff[,1],y,X1)))/ mean(exp(dmvnorm(t(t(gibbs1$mu)-model1$coeff[,1]),mean=rep(0,2),sigma=Sigmapost0,log=TRUE)-probitlpost(model1$coeff[,1],y,X0))) print(i) } ############# BRIDGE SAMPLING APPROXIMATION Niter=100 Sigmaprior0=n*ginv(t(X0)%*%X0) Sigmapost0=ginv(t(X0)%*%X0+ginv(Sigmaprior0)) Sigmaprior1=n*ginv(t(X1)%*%X1) Sigmapost1=ginv(t(X1)%*%X1+ginv(Sigmaprior1)) bfbs=rep(0,NMC) cova=model2$cov.unscaled expecta=model2$coeff[,1] meanw=function(beta) { if (is.matrix(beta)==F) beta=as.matrix(t(beta)) niter=dim(beta)[1] lala=matrix(rep(expecta[1:2],niter),2) t(expecta[3]+t(cova[1:2,3])%*%ginv(cova[1:2,1:2])%*%(t(beta)-lala)) } covw=cova[3,3]-t(cova[1:2,3])%*%ginv(cova[1:2,1:2])%*%cova[1:2,3] for (i in 1:NMC) { gibbs1=gibbsprobit(Niter,y,X0,Sigmapost0) gibbs2=gibbsprobit(Niter,y,X1,Sigmapost1) pseudo=rnorm(Niter,meanw(gibbs1$beta),sqrt(covw)) gibbs1p=cbind(gibbs1$beta,pseudo) bfbs[i]=mean(exp(probitlpost(gibbs2$beta[,1:2],y,X0)+dnorm(gibbs2$beta[,3],meanw(gibbs2$beta[,1:2]),sqrt(covw),log=TRUE))/ (dmvnorm(gibbs2$beta,expecta,cova)+dnorm(gibbs2$beta[,3],expecta[3],cova[3,3])))/ mean(exp(probitlpost(probit1p,y,X1))/(dmvnorm(probit1p,expecta,cova)+dnorm(pseudo,expecta[3],cova[3,3]))) print(i) } ############# IMPORTANCE SAMPLING APPROXIMATION bfis=rep(0,NMC) for (i in 1:NMC) { is1=rmvnorm(Niter,mean=model1$coeff[,1],sigma=model1$cov.unscaled) is2=rmvnorm(Niter,mean=model2$coeff[,1],sigma=model2$cov.unscaled) bfis[i]=mean(exp(probitlpost(is1,y,X0)-dmvnorm(is1,mean=model1$coeff[,1],sigma=model1$cov.unscaled,log=TRUE)))/ mean(exp(probitlpost(is2,y,X1)-dmvnorm(is2,mean=model2$coeff[,1],sigma=model2$cov.unscaled,log=TRUE))) print(i) } ################################################################################