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.
