[R] Bayesian question (problem using adapt)

[cmdicaso at ncsu.edu]

In the following code I have created the posterior density for a Bayesian
survival model with four parameters. However, when I try to use the adapt
function to perform integration in four dimensions (on my old version of R
I get an error message saying that I have applied a non-function, although
the function does work when I type kernel2(param0, theta0), or on the
newer version of R the computer crashes.

I have attached the following R code:

hazard1 <- function(x, param) {
 if(min(x, param)<=0) {return("function undefined")}
 return(param[1]*x^(param[2] -1))
 }

 gm <- function( theta, param) {
 dgamma(param[1], shape = .1, scale = .1) *
 dgamma(param[2], shape = .1, scale = .1) *
 dgamma(theta[1], shape = .1, scale = .1) *
 dgamma(theta[2], shape = .1, scale = .1)
 }


 chazard1 <-function(x, param) {
 if(min(x, param) <=0) {return("function undefined")}
 param[1]*x^param[2] / param[2]
 }

 hazard2 <-function(x, theta, param) {
 if(min(x, theta, param) <= 0) {return("function undefined")}
 ratio = 1/(theta[1]*exp(-chazard1(x, param)) +
theta[2]*(1-exp(-chazard1(x, param))))
 return(ratio*hazard1(x, param))
 }

 chazard2 <- function(x, theta, param) {
 if(min(x, theta, param) <=0) {return("function undefined")}
 val=c()
 for (i in 1:length(x)) {
 val=c(val, integrate(hazard2, lower=0, upper=x, theta=theta,
param=param)$value)}
 return(val)
 }


 param0 = c(1.11, .833)
 theta0 = c(1.11, .833)


 yn1 <-read.table("yang1.txt", col.names = c("trtmt", "deathday", "reason"))
 yn2 <-read.table("yang2.txt", col.names = c("trtmt", "deathday", "reason"))



 for (i in 1:length(yn1$reason)) {
 if (yn1$reason[i] == -1) (yn1$reason[i] = 0) else (yn1$reason[i] = 1) }

 for (i in 1:length(yn2$reason)) {
 if (yn2$reason[i] == -1) (yn2$reason[i] = 0) else (yn2$reason[i] = 1) }

 x1 = yn1$deathday / (1*10^3);
 d1 = yn1$reason;
 x2 = yn2$deathday / (1* 10^3);
 d2 = yn2$reason;



 loglik <- function( theta, param) { sum(d1*log(hazard1(x1, param)) -
chazard1(x1,

param) + d2*log(hazard2(x2, theta, param)) - chazard2(x2, theta, param))}
/ 10000



 kernel1 <- function(theta, param) { exp(loglik(theta, param)) + log(10000)}

 kernel2 <- function(theta, param) {( exp(loglik(theta, param)) +
log(10000)) * gm(theta, param)  }

 kernel2(theta0, param0)



 adapt(4, lower = c(0, 0, 0, 0), upper = c(1, 1, 1, 1), functn =
kernel2(theta0, param0))

So I am trying to integrate the posterior density with four different
parameters. But I cannot do anything to get the adapt function to work.