[R] Metropolis code help

[francogrex]

Dears, I have the below code for metropolis of the GLM logit (logistic
regression) using a flat prior. Can someone help me modify the prior so that
the model becomes hierarchical by using a flat prior for mu and sigma, the
derived density for beta ~ N(mu, sigma^2)? Actually I took my code from a
teacher that posted on the internet and modified it to the GLM logit but I
can't adapt it to the hierarchical version. 
Here is the original code of the teacher with both flat prior on betas and a
hierarchical version:
www.stats.uwo.ca/faculty/murdoch/458/metropolis.r 


Below is My code with a flat prior on beta only (I'd like also to have the
hierarchical version!)
X<- cbind(1,DF$nsaid,DF$diuretic,DF$diuretic*DF$nsaid)
y<- DF$Var3

Metropolis <- function(logtarget, start, R = 1000, sd
= 1) {
    parmcount <- length(start)
    sims <- matrix(NA, nrow=R, ncol = parmcount)
    colnames(sims) <- names(start)

    sims[1,] <- start
    oldlogalpha <- logtarget(start)
    accepts <- 0

    for (i in 2:R) {
	jump <- rnorm(parmcount, mean=0, sd=sd)
	y <- sims[i-1,] + jump

    	newlogalpha <- logtarget(y)
	if (log(runif(1)) < newlogalpha - oldlogalpha) {
	    sims[i,] <- y
	    oldlogalpha <- newlogalpha
	    accepts <- accepts + 1
	} else {
	    sims[i,] <- sims[i-1,]
	}
    }
    cat('Accepted ',100*accepts/(R-1),'%\n')
    sims
}

# Use the binomial likelihood 
logitll=function(beta,y,X)
{
X<- cbind(1,DF$nsaid,DF$diuretic,DF$diuretic*DF$nsaid)
y<- DF$Var3
lF1=plogis(X%*%as.vector(beta),log.p=TRUE)
lF2=plogis(-X%*%as.vector(beta),log.p=TRUE)
sum(y*lF1+(1-y)*lF2)
}

# Use a uniform prior for p
logprior <- function(beta,y,X) 0

# The log posterior is the sum.  It's the target of our MCMC run
logposterior <- function(beta,y,X) logprior(beta,y,X)+logitll(beta,y,X)

start <- c(0,0,0,0)
sims <- Metropolis(logposterior, start, 10000, sd=1)

se_sims <- apply(sims, 2, sd)

sims <- Metropolis(logposterior, start, 10000,sd=se_sims)

cbind(rbind(mean(sims[1001:10000,1]),mean(sims[1001:10000,2]),mean(sims[1001:10000,3]),mean(sims[1001:10000,4])),
rbind(sd(sims[1001:10000,1]),sd(sims[1001:10000,2]),sd(sims[1001:10000,3]),sd(sims[1001:10000,4])))

Thanks in advance.

-- 
View this message in context: http://www.nabble.com/Metropolis-code-help-tf3850742.html#a10907938
Sent from the R help mailing list archive at Nabble.com.