# [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])))

```--
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