[R] Metropolis code help

From: francogrex <francogrex_at_mail.com>
Date: Fri, 01 Jun 2007 00:18:34 -0700 (PDT)

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.

______________________________________________
R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Fri 01 Jun 2007 - 07:27:46 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 01 Jun 2007 - 08:31:45 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.