# [R] problem that arises after using the new version of "BRugs"

moumita chatterjee mcmoumita8 at gmail.com
Fri Jan 18 11:55:48 CET 2013

```Respected Sir,
With reference to my mail to you and the reply mail
by you dated 9th and 16th January, 2013, I am sending the reproducible code
in the attached document named " MODIFIED ANS ". I am also attaching the
txt file named "hazModel", which is required to save in my documents folder
to run the program. The file also contains the error message along with the
code.
I would be very much obliged if
you kindly give me a solution for this question.
Thanking you.
Moumita Chatterjee
Research Scholar
University Of Calcutta.
-------------- next part --------------
x<-c(13,5,8,0,11)
delta<-c(1,1,0,1,0)

ZOSull<-
function (x, range.x, intKnots, drv = 0)
{
if (drv > 2)
stop("splines not smooth enough for more than 2 derivatives")
library(splines)
if (missing(range.x))
range.x <- c(1.05 * min(x) - 0.05 * max(x), 1.05 * max(x) -
0.05 * min(x))
if (missing(intKnots)) {
numIntKnots <- min(length(unique(x)), 35)
intKnots <- quantile(unique(x), seq(0, 1, length = (numIntKnots +
2))[-c(1, (numIntKnots + 2))])
}
numIntKnots <- length(intKnots)
allKnots <- c(rep(range.x, 4), intKnots, rep(range.x,
4))
K <- length(intKnots)
L <- 3 * (K + 8)
xtilde <- (rep(allKnots, each = 3)[-c(1, (L - 1), L)] + rep(allKnots,
each = 3)[-c(1, 2, L)])/2
wts <- rep(diff(allKnots), each = 3) * rep(c(1, 4, 1)/6,
K + 7)
Bdd <- spline.des(allKnots, xtilde, derivs = rep(2, length(xtilde)),
outer.ok = TRUE)\$design
Omega <- t(Bdd * wts) %*% Bdd
eigOmega <- eigen(Omega)
indsZ <- 1:(numIntKnots + 2)
UZ <- eigOmega\$vectors[, indsZ]
LZ <- t(t(UZ)/sqrt(eigOmega\$values[indsZ]))
indsX <- (numIntKnots + 3):(numIntKnots + 4)
UX <- eigOmega\$vectors[, indsX]
L <- cbind(UX, LZ)
stabCheck <- t(crossprod(L, t(crossprod(L, Omega))))
if (sum(stabCheck^2) > 1.0001 * (numIntKnots + 2))
print("WARNING: NUMERICAL INSTABILITY ARISING\\\n              FROM SPECTRAL DECOMPOSITION")
B <- spline.des(allKnots, x, derivs = rep(drv, length(x)),
outer.ok = TRUE)\$design
Z <- B %*% LZ
attr(Z, "range.x") <- range.x
attr(Z, "intKnots") <- intKnots
return(Z)
}

BRugsMCMC<-
function (data, inits, parametersToSave, nBurnin, nIter, nThin,
modelFile)
{

if ((nBurnin < 100) | (nIter < 100))
stop("currently only working for chains longer than 100")
if ((100 * round(nBurnin/100) != nBurnin) | (100 * round(nIter/100) !=
nIter))
warning("chain lengths not multiples of 100; truncation may occur.")

modelCheck(modelFile)

modelData(bugsData(data))

modelCompile(numChains = 1)

modelInits(bugsInits(inits))

burninIncSize <- round(nBurnin/(nThin * 100))
for (i in 1:100) {
modelUpdate(burninIncSize, thin = nThin)

}
samplesSet(parametersToSave)
iterIncSize <- round(nIter/(nThin * 100))
for (i in 1:100) {
modelUpdate(iterIncSize, thin = nThin)

}

sims <- samplesStats("*")
return(list(Stats = sims))
}

require(BRugs)

numKnots <- 25

range.x <- c(0,202)
ox <- order(x)
x <- x[ox]
delta <- delta[ox]
n <- length(x)
xTil <- as.numeric(names(table(x)))
cnts <- as.numeric(table(x))
nTil <- length(xTil)
deltaTil <- rep(NA, nTil)
for (iTil in 1:nTil) deltaTil[iTil] <- sum(delta[(1:n)[x ==
xTil[iTil]]])
sd.xTil <- sd(xTil)
xTil <- xTil/sd.xTil
range.x <- range.x/sd.xTil
d1xTil <- xTil[2:nTil] - xTil[1:(nTil - 1)]
d2xTil <- xTil[3:nTil] - xTil[1:(nTil - 2)]
rcsrCnts <- rev(cumsum(rev(cnts)))
trapLef <- 2 * cnts * xTil + (sum(cnts) - cnts) *
(xTil + xTil)
trapMid <- cnts[2:(nTil - 1)] * d1xTil[-length(d1xTil)] +
rcsrCnts[1:(nTil - 2)] * d2xTil
trapRig <- cnts[nTil] * d1xTil[length(d1xTil)]
oTil <- log(c(trapLef, trapMid, trapRig)/2)
X <- cbind(rep(1, nTil), xTil)
numIntKnots <- 25
intKnots <- quantile(xTil, seq(0, 1, length = numIntKnots +
2)[-c(1, numIntKnots + 2)])

Z <- ZOSull(xTil, intKnots = intKnots, range.x = range.x)
numKnots <- ncol(Z)
allData <- list(n = nTil, numKnots = numKnots, x = xTil,
delta = deltaTil, offs = oTil, Z = Z)
parInits <- list(list(beta0 = 0.1, beta1 = 0.1, u = rep(0.1,
numKnots), tauU = 1))

#################################################
# the file hazModel.txt must be in the
# My Documents folder for the next step to run, I am attaching this with my mail.
#################################################

BRugsMCMC(data = allData, inits = parInits, parametersToSave = c("beta0",
"beta1", "u", "tauU"), nBurnin =150, nIter =150,
nThin = 15, modelFile = "hazModel.txt")

betaMCMC <- rbind(samplesSample("beta0"), samplesSample("beta1"))

# An error is coming in this line. THe error window shows that BUGSHE~1.EXE has stopped working.#

# the error message is #
model has probably not yet been updated
model has probably not yet been updated
Error in handleRes(res) : NA
running command '"C:/Users/MOU/Documents/R/win-library/2.15/BRugs/exec/BugsHelper.exe" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S/trash" "file7f819516f3f.bug" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S/cmds.txt" "2"' had status 255

-------------- next part --------------
model

{

for(i in 1:n)

{

log(mu[i]) <- beta0 + beta1*x[i] + inprod(u[],Z[i,]) + offs[i]

delta[i] ~ dpois(mu[i])

}

for (k in 1:numKnots)

{

u[k] ~  dnorm(0,tauU)

}

beta0 ~ dnorm(0,1.0E-8)

beta1 ~ dnorm(0,1.0E-8)

tauU   ~ dgamma(0.01,0.01)

sigU <- 1/sqrt(tauU)

}
```