I'm trying to fit my data to Inverse-Gaussian distribution. The model runs with no error, but the results are not satisfactory. I get R.hat values bigger than 1 and n.eff in single digits. What did I do wrong/how can I improve my results?
It's hard to say what is wrong if you don't give a reproducible example - i.e. one with data.that readers can try out for themselves.
Having said that, your problem closely resembles the example from chapter 8 of Bayesian Modeling Using WinBUGS by Ioannis Ntzoufras. I have reproduced this example below and it runs fine. This leads me to conclude that the problem lies with your data (which you have not shared). Specifically, you have 18 covariates, which is a lot, If they are highly correlated then you will get poor performance from Gibbs sampling, which updates each beta[i] parameter individually.
The response variable P[i] is insurance premium and the covariates are risk classes.I'll consider reducing the covariates but please let me know if you can get fine results by using the 18 covariates.Thanks!
You can centre the covariates by defining X like this
X <- Data[,1:17]
x.bar <-apply(X,2,mean)
X <-sweep(X,2, x.bar, FUN="-")
X <-cbind(1, X)
xbar <-c(1, xbar)
You will still get correct estimates for beta[2] up to beta[18]. The intercept parameter beta[1] is different, but if you want to you can recover the original intercept parameter as
alpha <- inprod(beta, x.bar)
With this reparameterization, your model runs well with a thinning interval of (10, i.e. 10000 iterations to get 1000 samples)
Martyn
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The reparameterization has managed to lower the values of Rhat close to 1, but some n.eff values are still low compared to the number of simulations saved. Here's the updated model:
Hi,
I'm trying to fit my data to Inverse-Gaussian distribution. The model runs with no error, but the results are not satisfactory. I get R.hat values bigger than 1 and n.eff in single digits. What did I do wrong/how can I improve my results?
Below are the codes:
~~~~~~~~~~
modelText <- "
data{
C <- 10000
for(i in 1:N){
zeros[i] <- 0
}
}
model{
for(i in 1:N){
zeros[i] ~ dpois(zeros.mean[i])
zeros.mean[i] <- -l[i] + C
#log-likelihood
l[i] <- 0.5(log(lambda) - log(23.14) - 3log(P[i])) - 0.5 lambda pow( (P[i]- mu[i])/mu[i],2)/P[i]
#log-link + linear predictor
log(mu[i]) <- inprod(X[i,],beta[])
}
#priors
for(j in 1:18){
beta[j] ~ dnorm(0.0, 0.001)
}
lambda ~ dgamma(0.01,0.01)
s2 <- 1/lambda
Deviance <- -2sum(l[1:N])
}"
~~~~~~~~~~~~
Thanks in advance.
It's hard to say what is wrong if you don't give a reproducible example - i.e. one with data.that readers can try out for themselves.
Having said that, your problem closely resembles the example from chapter 8 of Bayesian Modeling Using WinBUGS by Ioannis Ntzoufras. I have reproduced this example below and it runs fine. This leads me to conclude that the problem lies with your data (which you have not shared). Specifically, you have 18 covariates, which is a lot, If they are highly correlated then you will get poor performance from Gibbs sampling, which updates each
beta[i]parameter individually.Model
Data
Thanks for the quick response.The data is as attached.Below is the full code.
~~~~~~
Data <- read.csv("Data.csv")
Define data matrix
X <- cbind(1,Data[,c(1:17)])
attach(X)
P <- Data$P
N <- length(Data$P)
IGaussian.data <- list("P","N","X")
Inverse-Gaussian using the zero trick
modelText <- "
data{
C <- 10000
for(i in 1:N){
zeros[i] <- 0
}
}
model{
for(i in 1:N){
zeros[i] ~ dpois(zeros.mean[i])
zeros.mean[i] <- -l[i] + C
#log-likelihood
l[i] <- 0.5(log(lambda) - log(23.14) - 3log(P[i])) - 0.5 lambda pow( (P[i]-mu[i])/mu[i],2)/P[i]
#log-link + linear predictor
log(mu[i]) <- inprod(X[i,],beta[])
}
#priors
for(j in 1:18){
beta[j] ~ dnorm(0.0, 0.001)
}
lambda ~ dgamma(0.01,0.01)
s2 <- 1/lambda
Deviance <- -2sum(l[1:N])
}"
writeLines(modelText,"IGaussian.txt")
Define parameters
IGaussian.params <- c(paste("beta[",i=1:18,"]",sep=""),"Deviance")
Define starting values
IGaussian.inits <- function(){list("beta"=rep(-0.1,18))}
~~~~~~~
The response variable P[i] is insurance premium and the covariates are risk classes.I'll consider reducing the covariates but please let me know if you can get fine results by using the 18 covariates.Thanks!
Last edit: Fadzli Fuzi 2014-10-10
You can centre the covariates by defining X like this
You will still get correct estimates for
beta[2]up tobeta[18]. The intercept parameterbeta[1]is different, but if you want to you can recover the original intercept parameter asWith this reparameterization, your model runs well with a thinning interval of (10, i.e. 10000 iterations to get 1000 samples)
Martyn
The reparameterization has managed to lower the values of Rhat close to 1, but some n.eff values are still low compared to the number of simulations saved. Here's the updated model:
~~~~~~~~~~~~~~
modelText <- "
data{
C <- 10000
for(i in 1:N){
zeros[i] <- 0
}
}
model{
for(i in 1:N){
zeros[i] ~ dpois(zeros.mean[i])
zeros.mean[i] <- -l[i] + C
#log-likelihood
l[i] <- 0.5(log(lambda) - log(23.14) - 3log(P[i])) - 0.5 lambda pow( (P[i]-mu[i])/mu[i],2)/P[i]
#log-link + linear predictor
log(mu[i]) <- inprod(xbar,beta)
}
#priors
for(j in 1:18){
beta[j] ~ dnorm(0.0, 0.001)
}
lambda ~ dgamma(0.01,0.01)
s2 <- 1/lambda
Deviance <- -2sum(l[1:N])
}"
~~~~~~~~~
Fadzli